pkg/oasim/oasim_slingo.F

87dd4f7d5f Oliv*0001 #include "OASIM_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: OASIM_SLINGO
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE OASIM_SLINGO(
                0009      I                  rmu0,clwp,cre,
                0010      O                  Tcd,Tcs,
                0011      I                  myThid)
                0012
                0013 C     !DESCRIPTION:
                0014 C     Slingo's (1989) Delta-Eddington approximation for the two-stream
                0015 C     equations applied to clouds.
                0016
                0017 C     !USES:
                0018       IMPLICIT NONE
                0019 #include "SIZE.h"
                0020 #include "OASIM_SIZE.h"
                0021 #include "OASIM_INTERNAL.h"
                0022 #include "OASIM_SLINGO.h"
                0023
                0024 C     !INPUT PARAMETERS:
                0025 c     rmu0   :: Kasten's approx for cosine of solar zenith angle
                0026 c     clwp   :: liquid water path in cloud (g/m2)
                0027 c     cre    :: cloud droplet effective radius (um)
                0028       _RL rmu0, clwp, cre
                0029       INTEGER myThid
                0030
                0031 C     !OUTPUT PARAMETERS:
                0032 c     Tcd :: spectral transmittance for downwelling direct irradiance
                0033 c     Tcs :: spectral transmittance for downwelling diffuse irradiance
                0034       _RL Tcd(ncld),Tcs(ncld)
                0035 CEOP
                0036
                0037 #ifdef ALLOW_OASIM
                0038
                0039 C     !LOCAL VARIABLES:
                0040 C     Tcu :: spectral transmittance for upwelling irradiance
                0041       INTEGER izero, nc
                0042       _RL Tcu(ncld)
                0043       _RL b0, alpha1, alpha3, alpha2, alpha4, bmu0
                0044       _RL e, eps, f, g, gama1, gama2, oneomega, omega
                0045       _RL rdif, rden, re, rm, rnum, sqarg, tauc, tdb
                0046       _RL tdir, tdif, u2, val2, val1, val3
                0047 c
                0048 c  Compute re as funtion of cldtau and LWP according to eq. 1 in
                0049 c  Slingo.
                0050 c   tau is derived at this wavelength (0.6 um) in the ISCCP data set
                0051 c      re = clwp*bsl(9)/(cldtau - clwp*asl(9))
                0052 c      re = min(re,15.0)  !block high re -- produces excessive direct
                0053 c  Changes to the ISCCP-D2 data set make this relationship untenable
                0054 c  (excessive res are derived).  Instead choose a fixed re of 10 um
                0055 c  for ocean (Kiehl et al., 1998 -- J. Clim.)
                0056 c       re = 10.0
                0057 c  Paper by Han et al., 1994 (J.Clim.) show mean ocean cloud radius
                0058 c  = 11.8 um
                0059 c       re = 11.8
                0060 c  Mean of Kiehl and Han
                0061       re = (10.0 _d 0+11.8 _d 0)/2.0 _d 0
                0062 c
                0063 c  Compute spectral cloud characteristics
                0064 c   If MODIS re is available use it; otherwise use parameterized re above
                0065       IF (cre .GE. 0.0 _d 0)THEN   !use MODIS re
                0066        re = cre
                0067       ENDIF
                0068       izero = 0
                0069       DO nc = 1,ncld
                0070        tauc = clwp*(asl(nc)+bsl(nc)/re)
                0071        oneomega = csl(nc) + dsl(nc)*re
                0072        omega = 1.0 _d 0 - oneomega
                0073        g = esl(nc) + fsl(nc)*re
                0074        b0 = 3.0 _d 0/7.0 _d 0*(1.0 _d 0-g)
                0075        bmu0 = 0.5 _d 0 - 0.75 _d 0*rmu0*g/(1.0 _d 0+g)
                0076        f = g*g
                0077        U2 = U1*(1.0 _d 0-((1.0 _d 0-omega)/(7.0 _d 0*omega*b0)))
                0078        U2 = MAX(U2,0.0 _d 0)
                0079        alpha1 = U1*(1.0 _d 0-omega*(1.0 _d 0-b0))
                0080        alpha2 = U2*omega*b0
                0081        alpha3 = (1.0 _d 0-f)*omega*bmu0
                0082        alpha4 = (1.0 _d 0-f)*omega*(1.0 _d 0-bmu0)
                0083        sqarg = alpha1*alpha1 - alpha2*alpha2
                0084        sqarg = MAX(sqarg,1.0 _d -17)
                0085        eps = SQRT(sqarg)
                0086        rM = alpha2/(alpha1+eps)
                0087        E = EXP(-eps*tauc)
                0088        val1 = 1.0 _d 0 - omega*f
                0089        val2 = eps*eps*rmu0*rmu0
                0090        rnum = val1*alpha3 - rmu0*(alpha1*alpha3+alpha2*alpha4)
                0091        rden = val1*val1 - val2
                0092        gama1 = rnum/rden
                0093        rnum = -val1*alpha4 - rmu0*(alpha1*alpha4+alpha2*alpha3)
                0094        gama2 = rnum/rden
                0095        Tdb = EXP(-val1*tauc/rmu0)
                0096        val3 = 1.0 _d 0 - E*E*rM*rM
                0097        Rdif = rM*(1.0 _d 0-E*E)/val3
                0098        Tdif = E*(1.0 _d 0-rM*rM)/val3
                0099 c       Rdir = -gama2*Rdif - gama1*Tdb*Tdif + gama1
                0100        Tdir = -gama2*Tdif - gama1*Tdb*Rdif + gama2*Tdb
                0101 c       Tdir = max(Tdir,0.0 _d 0)
                0102        Tcd(nc) = Tdb
                0103        Tcs(nc) = Tdir
                0104        IF (Tcs(nc) .LT. 0.0 _d 0)izero=1
                0105 c       Tcu(nc) = Tdif
                0106       ENDDO
                0107       IF (izero .EQ. 1)THEN    !block negative diffuse irrads
                0108        DO nc = 1,ncld
                0109         Tcs(nc) = 0.0
                0110        ENDDO
                0111       ENDIF
                0112
                0113 #endif /* ALLOW_OASIM */
                0114
                0115       RETURN
                0116       END