C $Header: /u/gcmpack/MITgcm/verification/fizhi-cs-aqualev20/code/fizhi_swrad.F,v 1.1 2006/04/03 20:55:14 molod Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
subroutine SWRIO (nymd,nhms,bi,bj,ndswr,myid,istrip,npcs,
. low_level,mid_level,im,jm,lm,
. pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,co2,
. albvisdr,albvisdf,albnirdr,albnirdf,
. dtradsw,dtswclr,radswg,swgclr,
. fdifpar,fdirpar,osr,osrclr,
. ptop,nswcld,cldsw,cswmo,nswlz,swlz,
. lpnt,imstturb,qliqave,fccave,landtype,xlats,xlons)
implicit none
c Input Variables
c ---------------
integer nymd,nhms,bi,bj,ndswr,myid,istrip,npcs
integer mid_level,low_level
integer im,jm,lm
_RL ptop
_RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm)
_RL pkht(im,jm,lm+1),pkz(im,jm,lm)
_RL tz(im,jm,lm),qz(im,jm,lm)
_RL oz(im,jm,lm)
_RL co2
_RL albvisdr(im,jm),albvisdf(im,jm),albnirdr(im,jm)
_RL albnirdf(im,jm)
_RL radswg(im,jm),swgclr(im,jm),fdifpar(im,jm),fdirpar(im,jm)
_RL osr(im,jm),osrclr(im,jm),dtradsw(im,jm,lm)
_RL dtswclr(im,jm,lm)
integer nswcld,nswlz
_RL cldsw(im,jm,lm),cswmo(im,jm,lm),swlz(im,jm,lm)
logical lpnt
integer imstturb
_RL qliqave(im,jm,lm),fccave(im,jm,lm)
integer landtype(im,jm)
_RL xlats(im,jm),xlons(im,jm)
c Local Variables
c ---------------
integer i,j,L,nn,nsecf
integer ntmstp,nymd2,nhms2
_RL getcon,grav,cp,undef
_RL ra,alf,reffw,reffi,tminv
parameter ( reffw = 10.0 )
parameter ( reffi = 65.0 )
_RL tdry(im,jm,lm)
_RL TEMP1(im,jm)
_RL TEMP2(im,jm)
_RL zenith (im,jm)
_RL cldtot (im,jm,lm)
_RL cldmxo (im,jm,lm)
_RL totcld (im,jm)
_RL cldlow (im,jm)
_RL cldmid (im,jm)
_RL cldhi (im,jm)
_RL taulow (im,jm)
_RL taumid (im,jm)
_RL tauhi (im,jm)
_RL tautype(im,jm,lm,3)
_RL tau(im,jm,lm)
_RL albedo(im,jm)
_RL PK(ISTRIP,lm)
_RL qzl(ISTRIP,lm),CLRO(ISTRIP,lm)
_RL TZL(ISTRIP,lm)
_RL OZL(ISTRIP,lm)
_RL PLE(ISTRIP,lm+1)
_RL COSZ(ISTRIP)
_RL dpstrip(ISTRIP,lm)
_RL albuvdr(ISTRIP),albuvdf(ISTRIP)
_RL albirdr(ISTRIP),albirdf(ISTRIP)
_RL difpar (ISTRIP),dirpar (ISTRIP)
_RL fdirir(istrip),fdifir(istrip)
_RL fdiruv(istrip),fdifuv(istrip)
_RL flux(istrip,lm+1)
_RL fluxclr(istrip,lm+1)
_RL dtsw(istrip,lm)
_RL dtswc(istrip,lm)
_RL taul (istrip,lm)
_RL reff (istrip,lm,2)
_RL tauc (istrip,lm,2)
_RL taua (istrip,lm)
_RL tstrip (istrip)
#ifdef ALLOW_DIAGNOSTICS
logical diagnostics_is_on
external
_RL tmpdiag(im,jm),tmpdiag2(im,jm)
#endif
logical first
data first /.true./
C **********************************************************************
C **** INITIALIZATION ****
C **********************************************************************
grav = getcon('GRAVITY')
cp = getcon('CP')
undef = getcon('UNDEF')
NTMSTP = nsecf(NDSWR)
TMINV = 1./float(ntmstp)
C Compute Temperature from Theta
C ------------------------------
do L=1,lm
do j=1,jm
do i=1,im
tdry(i,j,L) = tz(i,j,L)*pkz(i,j,L)
enddo
enddo
enddo
if (first .and. myid.eq.1 ) then
print *
print *,'Low-Level Clouds are Grouped between levels: ',
. lm,' and ',low_level
print *,'Mid-Level Clouds are Grouped between levels: ',
. low_level-1,' and ',mid_level
print *
first = .false.
endif
C **********************************************************************
C **** CALCULATE COSINE OF THE ZENITH ANGLE ****
C **********************************************************************
CALL ASTRO ( 20040321, NHMS, XLATS,XLONS, im*jm, TEMP1,RA )
NYMD2 = NYMD
NHMS2 = NHMS
CALL TICK ( NYMD2, NHMS2, NTMSTP )
CALL ASTRO ( 20040321, NHMS2, XLATS,XLONS, im*jm, TEMP2,RA )
do j = 1,jm
do i = 1,im
zenith(I,j) = TEMP1(I,j) + TEMP2(I,j)
IF( zenith(I,j) .GT. 0.0 ) THEN
zenith(I,j) = (2./3.)*( zenith(I,j)-TEMP2(I,j)*TEMP1(I,j)
. / zenith(I,j) )
ENDIF
ENDDO
ENDDO
C **********************************************************************
c **** Compute Two-Dimension Total Cloud Fraction (0-1) ****
C **********************************************************************
c Initialize Clear Sky Probability for Low, Mid, and High Clouds
c --------------------------------------------------------------
do j =1,jm
do i =1,im
cldlow(i,j) = 0.0
cldmid(i,j) = 0.0
cldhi (i,j) = 0.0
enddo
enddo
c Adjust cloud fractions and cloud liquid water due to moist turbulence
c ---------------------------------------------------------------------
if(imstturb.ne.0) then
do L =1,lm
do j =1,jm
do i =1,im
cldtot(i,j,L)=min(1.0 _d 0,max(cldsw(i,j,L),fccave(i,j,L)/
$ imstturb))
cldmxo(i,j,L)=min(1.0 _d 0,cswmo(i,j,L))
swlz(i,j,L)=swlz(i,j,L)+qliqave(i,j,L)/imstturb
enddo
enddo
enddo
else
do L =1,lm
do j =1,jm
do i =1,im
cldtot(i,j,L) = min( 1.0 _d 0,cldsw(i,j,L) )
cldmxo(i,j,L) = min( 1.0 _d 0,cswmo(i,j,L) )
enddo
enddo
enddo
endif
c Compute 3-D Cloud Fractions
c ---------------------------
if( nswcld.ne.0 ) then
do L = 1,lm
do j = 1,jm
do i = 1,im
c Compute Low-Mid-High Maximum Overlap Cloud Fractions
c ----------------------------------------------------
if( L.lt.mid_level ) then
cldhi (i,j) = max( cldtot(i,j,L),cldhi (i,j) )
elseif( L.ge.mid_level .and.
. L.lt.low_level ) then
cldmid(i,j) = max( cldtot(i,j,L),cldmid(i,j) )
elseif( L.ge.low_level ) then
cldlow(i,j) = max( cldtot(i,j,L),cldlow(i,j) )
endif
enddo
enddo
enddo
endif
c Totcld => Product of Clear Sky Probabilities for Low, Mid, and High Clouds
c --------------------------------------------------------------------------
do j = 1,jm
do i = 1,im
totcld(i,j) = 1.0 - (1.-cldhi (i,j))
. * (1.-cldmid(i,j))
. * (1.-cldlow(i,j))
enddo
enddo
#ifdef ALLOW_DIAGNOSTICS
c Compute Cloud Diagnostics
c -------------------------
if(diagnostics_is_on('CLDFRC ',myid) ) then
call DIAGNOSTICS_FILL(totcld,'CLDFRC ',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('CLDRAS ',myid) ) then
do L=1,lm
do j=1,jm
do i=1,im
tmpdiag(i,j) = cswmo(i,j,L)
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'CLDRAS ',L,1,3,bi,bj,myid)
enddo
endif
if(diagnostics_is_on('CLDTOT ',myid) ) then
do L=1,lm
do j=1,jm
do i=1,im
tmpdiag(i,j) = cldtot(i,j,L)
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'CLDTOT ',L,1,3,bi,bj,myid)
enddo
endif
if(diagnostics_is_on('CLDLOW ',myid) ) then
call DIAGNOSTICS_FILL(cldlow,'CLDLOW ',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('CLDMID ',myid) ) then
call DIAGNOSTICS_FILL(cldmid,'CLDMID ',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('CLDHI ',myid) ) then
call DIAGNOSTICS_FILL(cldhi,'CLDHI ',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('LZRAD ',myid) ) then
do L=1,lm
do j=1,jm
do i=1,im
tmpdiag(i,j) = swlz(i,j,L) * 1.0e6
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'LZRAD ',L,1,3,bi,bj,myid)
enddo
endif
c Albedo Diagnostics
c ------------------
if(diagnostics_is_on('ALBVISDR',myid) ) then
call DIAGNOSTICS_FILL(albvisdr,'ALBVISDR',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('ALBVISDF',myid) ) then
call DIAGNOSTICS_FILL(albvisdf,'ALBVISDF',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('ALBNIRDR',myid) ) then
call DIAGNOSTICS_FILL(albnirdr,'ALBNIRDR',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('ALBNIRDF',myid) ) then
call DIAGNOSTICS_FILL(albnirdf,'ALBNIRDF',0,1,3,bi,bj,myid)
endif
#endif
C Compute Optical Thicknesses and Diagnostics
C -------------------------------------------
call OPTHK(tdry,plz,plze,swlz,cldtot,cldmxo,landtype,im,jm,lm,
. tautype)
do L = 1,lm
do j = 1,jm
do i = 1,im
tau(i,j,L) = tautype(i,j,L,1)+tautype(i,j,L,2)+tautype(i,j,L,3)
enddo
enddo
enddo
#ifdef ALLOW_DIAGNOSTICS
if(diagnostics_is_on('TAUAVE ',myid) ) then
do L=1,lm
do j=1,jm
do i=1,im
tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'TAUAVE ',L,1,3,bi,bj,myid)
enddo
endif
if(diagnostics_is_on('TAUCLD ',myid) .and.
. diagnostics_is_on('TAUCLDC ',myid) ) then
do L=1,lm
do j=1,jm
do i=1,im
if( cldtot(i,j,L).ne.0.0 ) then
tmpdiag(i,j) = tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L))
tmpdiag2(i,j) = 1.
else
tmpdiag(i,j) = 0.
tmpdiag2(i,j) = 0.
endif
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'TAUCLD ',L,1,3,bi,bj,myid)
call DIAGNOSTICS_FILL(tmpdiag,'TAUCLDC ',L,1,3,bi,bj,myid)
enddo
endif
c Compute Low, Mid, and High Cloud Optical Depth Diagnostics
c ----------------------------------------------------------
if(diagnostics_is_on('TAULOW ',myid) .and.
. diagnostics_is_on('TAULOWC ',myid) ) then
do j = 1,jm
do i = 1,im
taulow(i,j) = 0.0
if( cldlow(i,j).ne.0.0 ) then
do L = low_level,lm
taulow(i,j) = taulow(i,j) + tau(i,j,L)
enddo
tmpdiag2(i,j) = 1.
else
tmpdiag(i,j) = 0.
endif
enddo
enddo
call DIAGNOSTICS_FILL(taulow,'TAULOW ',0,1,3,bi,bj,myid)
call DIAGNOSTICS_FILL(tmpdiag2,'TAULOWC ',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('TAUMID ',myid) .and.
. diagnostics_is_on('TAUMIDC ',myid) ) then
do j = 1,jm
do i = 1,im
taumid(i,j) = 0.0
if( cldmid(i,j).ne.0.0 ) then
do L = mid_level,low_level+1
taumid(i,j) = taumid(i,j) + tau(i,j,L)
enddo
tmpdiag2(i,j) = 1.
else
tmpdiag(i,j) = 0.
endif
enddo
enddo
call DIAGNOSTICS_FILL(taumid,'TAUMID ',0,1,3,bi,bj,myid)
call DIAGNOSTICS_FILL(tmpdiag2,'TAUMIDC ',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('TAUHI ',myid) .and.
. diagnostics_is_on('TAUHIC ',myid) ) then
do j = 1,jm
do i = 1,im
tauhi(i,j) = 0.0
if( cldhi(i,j).ne.0.0 ) then
do L = 1,mid_level+1
tauhi(i,j) = tauhi(i,j) + tau(i,j,L)
enddo
tmpdiag2(i,j) = 1.
else
tmpdiag(i,j) = 0.
endif
enddo
enddo
call DIAGNOSTICS_FILL(tauhi,'TAUHI ',0,1,3,bi,bj,myid)
call DIAGNOSTICS_FILL(tmpdiag2,'TAUHIC ',0,1,3,bi,bj,myid)
endif
#endif
C***********************************************************************
C **** LOOP OVER GLOBAL REGIONS ****
C **********************************************************************
do 1000 nn = 1,npcs
C **********************************************************************
C **** VARIABLE INITIALIZATION ****
C **********************************************************************
CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN )
CALL STRIP ( plze, ple ,im*jm,ISTRIP,lm+1,NN)
CALL STRIP ( pkz , pk ,im*jm,ISTRIP,lm ,NN)
CALL STRIP ( dpres,dpstrip,im*jm,ISTRIP,lm ,NN)
CALL STRIP ( tdry, tzl ,im*jm,ISTRIP,lm ,NN)
CALL STRIP ( qz , qzl ,im*jm,ISTRIP,lm ,NN)
CALL STRIP ( oz , ozl ,im*jm,ISTRIP,lm ,NN)
CALL STRIP ( tau , taul ,im*jm,ISTRIP,lm ,NN)
CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN )
CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN )
CALL STRIP ( albnirdr,albirdr,im*jm,ISTRIP,1,NN )
CALL STRIP ( albnirdf,albirdf,im*jm,ISTRIP,1,NN )
call STRIP ( cldtot,clro,im*jm,istrip,lm,nn )
c Partition Tau between Water and Ice Particles
c ---------------------------------------------
do L= 1,lm
do i= 1,istrip
alf = min( max((tzl(i,l)-253.15)/20.,0. _d 0) ,1. _d 0)
taua(i,L) = 0.
if( alf.ne.0.0 .and. alf.ne.1.0 ) then
tauc(i,L,1) = taul(i,L)/(1.+alf/(1-alf) * (reffi/reffw*0.80) )
tauc(i,L,2) = taul(i,L)/(1.+(1-alf)/alf * (reffw/reffi*1.25) )
endif
if( alf.eq.0.0 ) then
tauc(i,L,1) = taul(i,L)
tauc(i,L,2) = 0.0
endif
if( alf.eq.1.0 ) then
tauc(i,L,1) = 0.0
tauc(i,L,2) = taul(i,L)
endif
reff(i,L,1) = reffi
reff(i,L,2) = reffw
enddo
enddo
call SORAD ( istrip,1,1,lm,ple,tzl,qzl,ozl,co2,
. tauc,reff,clro,mid_level,low_level,taua,
. albirdr,albirdf,albuvdr,albuvdf,cosz,
. flux,fluxclr,fdirir,fdifir,dirpar,difpar,
. fdiruv,fdifuv )
C **********************************************************************
C **** Compute Mass-Weighted Theta Tendencies from Fluxes ****
C **********************************************************************
do l=1,lm
do i=1,istrip
alf = grav*(ple(i,Lm+1)-ptop)/(cp*dpstrip(i,L)*100)
dtsw (i,L) = alf*( flux (i,L)-flux (i,L+1) )/pk(i,L)
dtswc(i,L) = alf*( fluxclr(i,L)-fluxclr(i,L+1) )/pk(i,L)
enddo
enddo
call PASTE ( dtsw , dtradsw ,istrip,im*jm,lm,nn )
call PASTE ( dtswc, dtswclr ,istrip,im*jm,lm,nn )
call PASTE ( flux (1,1),osr ,istrip,im*jm,1,nn )
call PASTE ( fluxclr(1,1),osrclr,istrip,im*jm,1,nn )
call PASTE ( flux (1,lm+1),radswg,istrip,im*jm,1,nn )
call PASTE ( fluxclr(1,lm+1),swgclr,istrip,im*jm,1,nn )
call PASTE ( dirpar,fdirpar,istrip,im*jm,1,nn )
call PASTE ( difpar,fdifpar,istrip,im*jm,1,nn )
c Calculate Mean Albedo
c ---------------------
do i=1,istrip
if( cosz(i).gt.0.0 ) then
tstrip(i) = 1.0 - flux(i,lm+1)/
. ( fdirir(i)+fdifir(i)+dirpar(i)+difpar(i) + fdiruv(i)+fdifuv(i) )
if( tstrip(i).lt.0.0 ) tstrip(i) = undef
else
tstrip(i) = undef
endif
enddo
call PASTE ( tstrip,albedo,istrip,im*jm,1,nn )
1000 CONTINUE
#ifdef ALLOW_DIAGNOSTICS
c Mean Albedo Diagnostic
c ----------------------
if(diagnostics_is_on('ALBEDO ',myid) .and.
. diagnostics_is_on('ALBEDOC ',myid) ) then
do j=1,jm
do i=1,im
if( albedo(i,j).ne.undef ) then
tmpdiag(i,j) = albedo(i,j)
tmpdiag2(i,j) = 1.
else
tmpdiag(i,j) = 0.
tmpdiag2(i,j) = 0.
endif
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'ALBEDO ',0,1,3,bi,bj,myid)
call DIAGNOSTICS_FILL(tmpdiag2,'ALBEDOC ',0,1,3,bi,bj,myid)
endif
#endif
C **********************************************************************
C **** ZERO-OUT OR BUMP COUNTERS ****
C **********************************************************************
nswlz = 0
nswcld = 0
imstturb = 0
do L = 1,lm
do j = 1,jm
do i = 1,im
fccave(i,j,L) = 0.
qliqave(i,j,L) = 0.
enddo
enddo
enddo
return
end
subroutine OPTHK ( tl,pl,ple,lz,cf,cfm,lwi,im,jm,lm,tau )
C***********************************************************************
C
C PURPOSE:
C ========
C Compute cloud optical thickness using temperature and
C cloud fractions.
C
C INPUT:
C ======
C tl ......... Temperature at Model Layers (K)
C pl ......... Model Layer Pressures (mb)
C ple ........ Model Edge Pressures (mb)
C lz ......... Diagnosed Convective and Large-Scale Cloud Water (g/g)
C cf ......... Total Cloud Fraction (Random + Maximum Overlap)
C cfm ........ Maximum Overlap Cloud Fraction
C lwi ........ Surface Land Type
C im ......... Longitudinal dimension
C jm ......... Latitudinal dimension
C lm ......... Vertical dimension
C
C OUTPUT:
C =======
C tau ........ Cloud Optical Thickness (non-dimensional)
C tau(im,jm,lm,1): Suspended Ice
C tau(im,jm,lm,2): Suspended Water
C tau(im,jm,lm,3): Raindrops
C
C***********************************************************************
implicit none
integer im,jm,lm,i,j,L
_RL tl(im,jm,lm)
_RL pl(im,jm,lm)
_RL ple(im,jm,lm+1)
_RL lz(im,jm,lm)
_RL cf(im,jm,lm)
_RL cfm(im,jm,lm)
_RL tau(im,jm,lm,3)
integer lwi(im,jm)
_RL dp, alf, fracls, fraccu
_RL tauice, tauh2o, tauras
c Compute Cloud Optical Depths
c ----------------------------
do L=1,lm
do j=1,jm
do i=1,im
alf = min( max(( tl(i,j,L)-233.15)/20.,0. _d 0), 1. _d 0)
dp = ple(i,j,L+1)-ple(i,j,L)
tau(i,j,L,1) = 0.0
tau(i,j,L,2) = 0.0
tau(i,j,L,3) = 0.0
if( cf(i,j,L).ne.0.0 ) then
c Determine fraction of large-scale and cumulus clouds
c ----------------------------------------------------
fracls = ( cf(i,j,L)-cfm(i,j,L) )/cf(i,j,L)
fraccu = 1.0-fracls
c Define tau for large-scale ice and water clouds
c Note: tauice is scaled between (0.02 & .2) for: 0 < lz < 2 mg/kg over Land
c Note: tauh2o is scaled between (0.20 & 20) for: 0 < lz < 5 mg/kg over Land
c Note: tauh2o is scaled between (0.20 & 12) for: 0 < lz < 50 mg/kg over Ocean
c -------------------------------------------------------------------------------
c Large-Scale Ice
c ---------------
tauice = max( 0.0002 _d 0, 0.002*min(500*lz(i,j,L)*1000,
$ 1.0 _d 0) )
tau(i,j,L,1) = fracls*(1-alf)*tauice*dp
c Large-Scale Water
c -----------------
C Over Land
if( lwi(i,j).le.10 ) then
tauh2o = max( 0.0020 _d 0, 0.200*min(200*lz(i,j,L)*1000,
$ 1.0 _d 0) )
tau(i,j,L,3) = fracls*alf*tauh2o*dp
else
C Non-Precipitation Clouds Over Ocean
if( lz(i,j,L).eq.0.0 ) then
tauh2o = .12
tau(i,j,L,2) = fracls*alf*tauh2o*dp
else
C Over Ocean
tauh2o = max( 0.0020 _d 0, 0.120*min( 20*lz(i,j,L)*1000,
$ 1.0 _d 0) )
tau(i,j,L,3) = fracls*alf*tauh2o*dp
endif
endif
c Sub-Grid Convective
c -------------------
if( tl(i,j,L).gt.210.0 ) then
tauras = .16
tau(i,j,L,3) = tau(i,j,L,3) + fraccu*tauras*dp
else
tauras = tauice
tau(i,j,L,1) = tau(i,j,L,1) + fraccu*tauras*dp
endif
c Define tau for large-scale and cumulus clouds
c ---------------------------------------------
ccc tau(i,j,L) = ( fracls*( alf*tauh2o + (1-alf)*tauice )
ccc . + fraccu*tauras )*dp
endif
enddo
enddo
enddo
return
end
subroutine SORAD(m,n,ndim,np,pl,ta,wa,oa,co2,
* taucld,reff,fcld,ict,icb,
* taual,rsirbm,rsirdf,rsuvbm,rsuvdf,cosz,
* flx,flc,fdirir,fdifir,fdirpar,fdifpar,
* fdiruv,fdifuv)
c********************************************************************
c
c This routine computes solar fluxes due to the absoption by water
c vapor, ozone, co2, o2, clouds, and aerosols and due to the
c scattering by clouds, aerosols, and gases.
c
c This is a vectorized code. It computes the fluxes simultaneous for
c (m x n) soundings, which is a subset of the (m x ndim) soundings.
c In a global climate model, m and ndim correspond to the numbers of
c grid boxes in the zonal and meridional directions, respectively.
c
c Ice and liquid cloud particles are allowed to co-exist in any of the
c np layers. Two sets of cloud parameters are required as inputs, one
c for ice paticles and the other for liquid particles. These parameters
c are optical thickness (taucld) and effective particle size (reff).
c
c If no information is available for reff, a default value of
c 10 micron for liquid water and 75 micron for ice can be used.
c
c Clouds are grouped into high, middle, and low clouds separated by the
c level indices ict and icb. For detail, see the subroutine cldscale.
c
c----- Input parameters:
c units size
c number of soundings in zonal direction (m) n/d 1
c number of soundings in meridional direction (n) n/d 1
c maximum number of soundings in n/d 1
c meridional direction (ndim)
c number of atmospheric layers (np) n/d 1
c level pressure (pl) mb m*ndim*(np+1)
c layer temperature (ta) k m*ndim*np
c layer specific humidity (wa) gm/gm m*ndim*np
c layer ozone concentration (oa) gm/gm m*ndim*np
c co2 mixing ratio by volumn (co2) parts/part 1
c cloud optical thickness (taucld) n/d m*ndim*np*2
c index 1 for ice particles
c index 2 for liquid drops
c effective cloud-particle size (reff) micrometer m*ndim*np*2
c index 1 for ice particles
c index 2 for liquid drops
c cloud amount (fcld) fraction m*ndim*np
c level index separating high and middle n/d 1
c clouds (ict)
c level index separating middle and low clouds n/d 1
c clouds (icb)
c aerosol optical thickness (taual) n/d m*ndim*np
c solar ir surface albedo for beam fraction m*ndim
c radiation (rsirbm)
c solar ir surface albedo for diffuse fraction m*ndim
c radiation (rsirdf)
c uv + par surface albedo for beam fraction m*ndim
c radiation (rsuvbm)
c uv + par surface albedo for diffuse fraction m*ndim
c radiation (rsuvdf)
c cosine of solar zenith angle (cosz) n/d m*ndim
c
c----- Output parameters
c
c all-sky flux (downward minus upward) (flx) fraction m*ndim*(np+1)
c clear-sky flux (downward minus upward) (flc) fraction m*ndim*(np+1)
c all-sky direct downward ir (0.7-10 micron)
c flux at the surface (fdirir) fraction m*ndim
c all-sky diffuse downward ir flux at
c the surface (fdifir) fraction m*ndim
c all-sky direct downward par (0.4-0.7 micron)
c flux at the surface (fdirpar) fraction m*ndim
c all-sky diffuse downward par flux at
c the surface (fdifpar) fraction m*ndim
*
c----- Notes:
c
c (1) The unit of flux is fraction of the incoming solar radiation
c at the top of the atmosphere. Therefore, fluxes should
c be equal to flux multiplied by the extra-terrestrial solar
c flux and the cosine of solar zenith angle.
c (2) Clouds and aerosols can be included in any layers by specifying
c fcld(i,j,k), taucld(i,j,k,*) and taual(i,j,k), k=1,np.
c For an atmosphere without clouds and aerosols,
c set fcld(i,j,k)=taucld(i,j,k,*)=taual(i,j,k)=0.0.
c (3) Aerosol single scattering albedos and asymmetry
c factors are specified in the subroutines solir and soluv.
c (4) pl(i,j,1) is the pressure at the top of the model, and
c pl(i,j,np+1) is the surface pressure.
c (5) the pressure levels ict and icb correspond approximately
c to 400 and 700 mb.
c
c**************************************************************************
implicit none
c-----Explicit Inline Directives
#ifdef CRAY
#ifdef f77
cfpp$ expand (expmn)
#endif
#endif
_RL expmn
c-----input parameters
integer m,n,ndim,np,ict,icb
_RL pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np)
_RL taucld(m,ndim,np,2),reff(m,ndim,np,2)
_RL fcld(m,ndim,np),taual(m,ndim,np)
_RL rsirbm(m,ndim),rsirdf(m,ndim),
* rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2
c-----output parameters
_RL flx(m,ndim,np+1),flc(m,ndim,np+1)
_RL fdirir(m,ndim),fdifir(m,ndim)
_RL fdirpar(m,ndim),fdifpar(m,ndim)
_RL fdiruv(m,ndim),fdifuv(m,ndim)
c-----temporary array
integer i,j,k
_RL cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np)
_RL dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np)
_RL swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1)
_RL sdf(m,n),sclr(m,n),csm(m,n),x
c-----------------------------------------------------------------
do j= 1, n
do i= 1, m
swh(i,j,1)=0.
so2(i,j,1)=0.
c-----csm is the effective secant of the solar zenith angle
c see equation (12) of Lacis and Hansen (1974, JAS)
csm(i,j)=35./sqrt(1224.*cosz(i,j)*cosz(i,j)+1.)
enddo
enddo
do k= 1, np
do j= 1, n
do i= 1, m
c-----compute layer thickness and pressure-scaling function.
c indices for the surface level and surface layer
c are np+1 and np, respectively.
dp(i,j,k)=pl(i,j,k+1)-pl(i,j,k)
scal(i,j,k)=dp(i,j,k)*(.5*(pl(i,j,k)+pl(i,j,k+1))/300.)**.8
c-----compute scaled water vapor amount, unit is g/cm**2
wh(i,j,k)=1.02*wa(i,j,k)*scal(i,j,k)*
* (1.+0.00135*(ta(i,j,k)-240.))
swh(i,j,k+1)=swh(i,j,k)+wh(i,j,k)
c-----compute ozone amount, unit is (cm-atm)stp.
oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7
enddo
enddo
enddo
c-----scale cloud optical thickness in each layer from taucld (with
c cloud amount fcld) to tauclb and tauclf (with cloud amount cc).
c tauclb is the scaled optical thickness for beam radiation and
c tauclf is for diffuse radiation.
call CLDSCALE(m,n,ndim,np,cosz,fcld,taucld,ict,icb,
* cc,tauclb,tauclf)
c-----initialize fluxes for all-sky (flx), clear-sky (flc), and
c flux reduction (df)
do k=1, np+1
do j=1, n
do i=1, m
flx(i,j,k)=0.
flc(i,j,k)=0.
df(i,j,k)=0.
enddo
enddo
enddo
c-----compute solar ir fluxes
call SOLIR (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff,ict,icb
* ,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir)
c-----compute and update uv and par fluxes
call SOLUV (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff,ict,icb
* ,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc
* ,fdirpar,fdifpar,fdiruv,fdifuv)
c-----compute scaled amount of o2 (so2), unit is (cm-atm)stp.
do k= 1, np
do j= 1, n
do i= 1, m
so2(i,j,k+1)=so2(i,j,k)+165.22*scal(i,j,k)
enddo
enddo
enddo
c-----compute flux reduction due to oxygen following
c chou (J. climate, 1990). The fraction 0.0287 is the
c extraterrestrial solar flux in the o2 bands.
do k= 2, np+1
do j= 1, n
do i= 1, m
x=so2(i,j,k)*csm(i,j)
df(i,j,k)=df(i,j,k)+0.0287*(1.-expmn(-0.00027*sqrt(x)))
enddo
enddo
enddo
c-----compute scaled amounts for co2 (so2). unit is (cm-atm)stp.
do k= 1, np
do j= 1, n
do i= 1, m
so2(i,j,k+1)=so2(i,j,k)+co2*789.*scal(i,j,k)
enddo
enddo
enddo
c-----compute and update flux reduction due to co2 following
c chou (J. Climate, 1990)
call FLXCO2(m,n,np,so2,swh,csm,df)
c-----adjust for the effect of o2 cnd co2 on clear-sky fluxes.
do k= 2, np+1
do j= 1, n
do i= 1, m
flc(i,j,k)=flc(i,j,k)-df(i,j,k)
enddo
enddo
enddo
c-----adjust for the all-sky fluxes due to o2 and co2. It is
c assumed that o2 and co2 have no effects on solar radiation
c below clouds. clouds are assumed randomly overlapped.
do j=1,n
do i=1,m
sdf(i,j)=0.0
sclr(i,j)</