pkg/oasim/oasim_navaer.F

87dd4f7d5f Oliv*0001 #include "OASIM_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: OASIM_NAVAER
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE OASIM_NAVAER(
                0009      I  relhum,am,Vi,ws,
                0010      O  beta,eta,wa,asymp,
                0011      I  myThid)
                0012
                0013 C     !DESCRIPTION:
                0014 C     Computes aerosol parameters according to a simplified version
                0015 C     of the Navy marine aerosol model.
                0016
                0017 C     !USES:
                0018       IMPLICIT NONE
                0019
                0020 C     !INPUT PARAMETERS:
                0021       _RL relhum, am, Vi, ws
                0022       INTEGER myThid
                0023
                0024 C     !OUTPUT PARAMETERS:
                0025       _RL beta, eta, wa, asymp
                0026 CEOP
                0027
                0028 #ifdef ALLOW_OASIM
                0029
                0030 C     !LOCAL VARIABLES:
                0031       INTEGER n,i
                0032       _RL rnum, rden, frh, arg, rval, sumx, sumy, sumxy, sumx2
                0033       _RL rlrn, rldndr, gama, rlogc, alpha, cext
                0034       _RL a(3),ro(3),dndr(3),r(3),rlam
                0035       DATA ro /0.03 _d 0,0.24 _d 0,2.0 _d 0/
                0036       DATA r /0.1 _d 0,1.0 _d 0,10.0 _d 0/
                0037       DATA rlam /0.55 _d 0/
                0038 c
                0039 c  Relative humidity factor
                0040 c      IF (relhum .GE. 100.0 _d 0)relhum = 99.9 _d 0
                0041       relhum = MIN(99.9 _d 0,relhum)
                0042       rnum = 2.0 _d 0 - relhum/100.0 _d 0
                0043       rden = 6.0 _d 0*(1.0 _d 0-relhum/100.0 _d 0)
                0044       frh = (rnum/rden)**0.333 _d 0
                0045 c
                0046 c  Size distribution amplitude components
                0047       a(1) = 2000.0 _d 0*am*am
                0048       a(2) = 5.866 _d 0*(ws-2.2 _d 0)
                0049 c      IF (a(2) .LT. 0.5 _d 0)a(2) = 0.5 _d 0
                0050       a(2) = MAX(0.5 _d 0,a(2))
                0051       a(3) = 0.01527 _d 0*(ws-2.2 _d 0)*0.05 _d 0        !from Hughes 1987
                0052 c      IF (a(3) .LT. 1.4 _d -5)a(3) = 1.4 _d -5
                0053       a(3) = MAX(1.4 _d -5,a(3))
                0054 c
                0055 c  Compute size distribution at three selected radii according to
                0056 c  Navy method
                0057       DO n = 1,3
                0058        dndr(n) = 0.0
                0059        DO i = 1,3
                0060         rden = frh*ro(i)
                0061         arg = LOG(r(n)/rden)*LOG(r(n)/rden)
                0062         rval = a(i)*EXP(-arg)/frh
                0063         dndr(n) = dndr(n) + rval
                0064        ENDDO
                0065       ENDDO
                0066 c
                0067 c  Least squares approximation
                0068       sumx = 0.0 _d 0
                0069       sumy = 0.0 _d 0
                0070       sumxy = 0.0 _d 0
                0071       sumx2 = 0.0 _d 0
                0072       DO n = 1,3
                0073        rlrn = LOG10(r(n))
                0074        rldndr = LOG10(dndr(n))
                0075        sumx = sumx + rlrn
                0076        sumy = sumy + rldndr
                0077        sumxy = sumxy + rlrn*rldndr
                0078        sumx2 = sumx2 + rlrn*rlrn
                0079       ENDDO
                0080       gama = sumxy/sumx2
                0081       rlogc = sumy/3.0 _d 0 - gama*sumx/3.0 _d 0
                0082       alpha = -(gama+3.0 _d 0)
                0083       eta = -alpha
                0084 c
                0085 c  Compute beta
                0086       cext = 3.91 _d 0/Vi
                0087       beta = cext*rlam**alpha
                0088 c
                0089 c  Compute asymmetry parameter -- a function of alpha
                0090       IF (alpha .GT. 1.2 _d 0)THEN
                0091        asymp = 0.65 _d 0
                0092       ELSE IF (alpha .LT. 0.0 _d 0)THEN
                0093        asymp = 0.82 _d 0
                0094       ELSE
                0095        asymp = -0.14167 _d 0*alpha + 0.82 _d 0
                0096       ENDIF
                0097 c
                0098 c  Single scattering albedo at 550; function of RH
                0099       wa = (-0.0032 _d 0*am + 0.972 _d 0)*EXP(3.06 _d -4*relhum)
                0100
                0101 #endif /* ALLOW_OASIM */
                0102
                0103       RETURN
                0104       END