C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_lwrad.F,v 1.26 2010/03/16 00:19:33 jmc Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
subroutine LWRIO (nymd,nhms,bi,bj,myid,istrip,npcs,
. low_level,mid_level,
. im,jm,lm,
. pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,
. co2,cfc11,cfc12,cfc22,methane,n2o,emissivity,
. tgz,radlwg,st4,dst4,
. dtradlw,dlwdtg,dtradlwc,lwgclr,
. nlwcld,cldlw,clwmo,nlwlz,lwlz,
. lpnt,imstturb,qliqave,fccave,landtype)
implicit none
c Input Variables
c ---------------
integer nymd,nhms,istrip,npcs,bi,bj,myid
integer mid_level,low_level
integer im,jm,lm
_RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1)
_RL dpres(im,jm,lm),pkht(im,jm,lm+1),pkz(im,jm,lm)
_RL tz(im,jm,lm),qz(im,jm,lm),oz(im,jm,lm)
_RL co2,cfc11,cfc12,cfc22,methane(lm),n2o(lm)
_RL emissivity(im,jm,10)
_RL tgz(im,jm),radlwg(im,jm),st4(im,jm),dst4(im,jm)
_RL dtradlw(im,jm,lm),dlwdtg (im,jm,lm)
_RL dtradlwc(im,jm,lm),lwgclr(im,jm)
integer nlwcld,nlwlz
_RL cldlw(im,jm,lm),clwmo(im,jm,lm),lwlz(im,jm,lm)
logical lpnt
integer imstturb
_RL qliqave(im,jm,lm),fccave(im,jm,lm)
integer landtype(im,jm)
c Local Variables
c ---------------
integer i,j,l,n,nn
_RL cldtot (im,jm,lm)
_RL cldmxo (im,jm,lm)
_RL pl(istrip,lm)
_RL pk(istrip,lm)
_RL pke(istrip,lm)
_RL ple(istrip,lm+1)
_RL ADELPL(ISTRIP,lm)
_RL dtrad(istrip,lm),dtradc(istrip,lm)
_RL OZL(ISTRIP,lm),TZL(ISTRIP,lm)
_RL SHZL(ISTRIP,lm),CLRO(ISTRIP,lm)
_RL CLMO(ISTRIP,lm)
_RL flx(ISTRIP,lm+1),flxclr(ISTRIP,lm+1)
_RL cldlz(istrip,lm)
_RL dfdts(istrip,lm+1),dtdtg(istrip,lm)
_RL emiss(istrip,10)
_RL taual(istrip,lm,10)
_RL ssaal(istrip,lm,10)
_RL asyal(istrip,lm,10)
_RL cwc(istrip,lm,3)
_RL reff(istrip,lm,3)
_RL tauc(istrip,lm,3)
_RL SGMT4(ISTRIP),TSURF(ISTRIP),dsgmt4(ISTRIP)
integer lwi(istrip)
_RL tmpstrip(istrip,lm)
_RL tmpimjm(im,jm,lm)
_RL tempor1(im,jm),tempor2(im,jm)
_RL getcon,secday,convrt
#ifdef ALLOW_DIAGNOSTICS
logical diagnostics_is_on
external
_RL tmpdiag(im,jm)
#endif
logical high, trace, cldwater
c data high /.true./
c data trace /.true./
data high /.false./
data trace /.false./
data cldwater /.false./
C **********************************************************************
C **** INITIALIZATION ****
C **********************************************************************
SECDAY = GETCON('SDAY')
CONVRT = GETCON('GRAVITY') / ( 100.0 * GETCON('CP') )
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(cldlw(i,j,L),fccave(i,j,L)/
$ imstturb))
cldmxo(i,j,L) = min( 1.0 _d 0, clwmo(i,j,L) )
lwlz(i,j,L) = lwlz(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,cldlw(i,j,L) )
cldmxo(i,j,L) = min( 1.0 _d 0,clwmo(i,j,L) )
enddo
enddo
enddo
endif
C***********************************************************************
C **** LOOP OVER STRIPS ****
C **********************************************************************
DO 1000 NN=1,NPCS
C **********************************************************************
C **** VARIABLE INITIALIZATION ****
C **********************************************************************
CALL STRIP (OZ,OZL ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (PLZE,ple ,im*jm,ISTRIP,lm+1 ,NN)
CALL STRIP (PLZ ,pl ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (PKZ ,pk ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (PKHT,pke ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (TZ,TZL ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (qz,SHZL ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (cldtot,CLRO ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (cldmxo,CLMO ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (lwlz,cldlz ,im*jm,ISTRIP,lm ,NN)
CALL STRIP (tgz,tsurf ,im*jm,ISTRIP,1 ,NN)
CALL STRIP (emissivity,emiss,im*jm,ISTRIP,10,NN)
call STRIPITINT (landtype,lwi,im*jm,im*jm,istrip,1,nn)
DO L = 1,lm
DO I = 1,ISTRIP
ADELPL(I,L) = convrt / ( ple(I,L+1)-ple(I,L) )
ENDDO
ENDDO
C Compute Clouds
C --------------
IF(NLWCLD.NE.0)THEN
DO L = 1,lm
DO I = 1,ISTRIP
CLRO(I,L) = min( 1.0 _d 0,clro(i,L) )
CLMO(I,L) = min( 1.0 _d 0,clmo(i,L) )
ENDDO
ENDDO
ENDIF
C Convert to Temperature from Fizhi Theta
C ---------------------------------------
DO L = 1,lm
DO I = 1,ISTRIP
TZL(I,L) = TZL(I,L)*pk(I,L)
ENDDO
ENDDO
DO I = 1,ISTRIP
C To reduce longwave heating in lowest model layer
TZL(I,lm) = ( 2*tzl(i,lm)+tsurf(i) )/3.0
ENDDO
C Compute Optical Thicknesses
C ---------------------------
call OPTHK ( tzl,pl,ple,cldlz,clro,clmo,lwi,istrip,1,lm,tauc )
do n = 1,3
do L = 1,lm
do i = 1,istrip
tauc(i,L,n) = tauc(i,L,n)*0.75
enddo
enddo
enddo
C Set the optical thickness, single scattering albedo and assymetry factor
C for aerosols equal to zero to omit the contribution of aerosols
c-------------------------------------------------------------------------
do n = 1,10
do L = 1,lm
do i = 1,istrip
taual(i,L,n) = 0.
ssaal(i,L,n) = 0.
asyal(i,L,n) = 0.
enddo
enddo
enddo
C Set cwc and reff to zero - when cldwater is .false.
c----------------------------------------------------
do n = 1,3
do L = 1,lm
do i = 1,istrip
cwc(i,L,n) = 0.
reff(i,L,n) = 0.
enddo
enddo
enddo
C **********************************************************************
C **** CALL M-D Chou RADIATION ****
C **********************************************************************
call IRRAD ( istrip,1,1,lm,ple,tzl,shzl,ozl,tsurf,co2,
. n2o,methane,cfc11,cfc12,cfc22,emiss,
. cldwater,cwc,tauc,reff,clro,mid_level,low_level,
. taual,ssaal,asyal,
. high,trace,flx,flxclr,dfdts,sgmt4 )
C **********************************************************************
C **** HEATING RATES ****
C **********************************************************************
do L = 1,lm
do i = 1,istrip
dtrad(i,L) = ( flx(i,L)- flx(i,L+1))*adelpl(i,L)
tmpstrip(i,L) = flx(i,L)
dtdtg(i,L) = ( dfdts(i,L)- dfdts(i,L+1))*adelpl(i,L)
dtradc(i,L) = (flxclr(i,L)-flxclr(i,L+1))*adelpl(i,L)
enddo
enddo
C **********************************************************************
C **** NET UPWARD LONGWAVE FLUX (W/M**2) ****
C **********************************************************************
DO I = 1,ISTRIP
flx (i,1) = -flx (i,1)
flxclr(i,1) = -flxclr(i,1)
flx (i,lm+1) = -flx (i,lm+1)
flxclr(i,lm+1) = -flxclr(i,lm+1)
sgmt4(i) = - sgmt4(i)
dsgmt4(i) = - dfdts(i,lm+1)
ENDDO
CALL PASTE ( flx (1,lm+1),RADLWG,ISTRIP,im*jm,1,NN )
CALL PASTE ( flxclr(1,lm+1),LWGCLR,ISTRIP,im*jm,1,NN )
CALL PASTE ( sgmt4, st4,ISTRIP,im*jm,1,NN )
CALL PASTE ( dsgmt4,dst4,ISTRIP,im*jm,1,NN )
C **********************************************************************
C **** PASTE AND BUMP SOME DIAGNOSTICS ****
C **********************************************************************
CALL PASTE(flx(1,1),tempor1,ISTRIP,im*jm,1,NN)
CALL PASTE(flxclr(1,1),tempor2,ISTRIP,im*jm,1,NN)
C **********************************************************************
C **** TENDENCY UPDATES ****
C **********************************************************************
DO L = 1,lm
DO I = 1,ISTRIP
DTRAD (I,L) = ple(i,lm+1) * DTRAD (I,L)/pk(I,L)
DTRADC(I,L) = ple(i,lm+1) * DTRADC(I,L)/pk(I,L)
dtdtg(I,L) = ple(i,lm+1) * dtdtg (I,L)/pk(I,L)
ENDDO
ENDDO
CALL PASTE ( tmpstrip ,tmpimjm ,ISTRIP,im*jm,lm,NN )
CALL PASTE ( DTRAD ,DTRADLW ,ISTRIP,im*jm,lm,NN )
CALL PASTE ( DTRADC,DTRADLWC,ISTRIP,im*jm,lm,NN )
CALL PASTE ( dtdtg ,dlwdtg ,ISTRIP,im*jm,lm,NN )
1000 CONTINUE
C **********************************************************************
C **** BUMP DIAGNOSTICS ****
C **********************************************************************
#ifdef ALLOW_DIAGNOSTICS
if(diagnostics_is_on('TGRLW ',myid) ) then
call DIAGNOSTICS_FILL(tgz,'TGRLW ',0,1,3,bi,bj,myid)
endif
do L = 1,lm
if(diagnostics_is_on('TLW ',myid) ) then
do j = 1,jm
do i = 1,im
tmpdiag(i,j) = tz(i,j,L)*pkz(i,j,L)
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'TLW ',L,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('SHRAD ',myid) ) then
do j = 1,jm
do i = 1,im
tmpdiag(i,j) = qz(i,j,L)*1000.
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'SHRAD ',L,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('OZLW ',myid) ) then
do j = 1,jm
do i = 1,im
tmpdiag(i,j) = oz(i,j,L)
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'OZLW ',L,1,3,bi,bj,myid)
endif
enddo
if(diagnostics_is_on('OLR ',myid) ) then
call DIAGNOSTICS_FILL(tempor1,'OLR ',0,1,3,bi,bj,myid)
endif
if(diagnostics_is_on('OLRCLR ',myid) ) then
call DIAGNOSTICS_FILL(tempor2,'OLRCLR ',0,1,3,bi,bj,myid)
endif
#endif
C **********************************************************************
C **** Increment Diagnostics Counters and Zero-Out Cloud Info ****
C **********************************************************************
nlwlz = 0
nlwcld = 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
c********************** November 26, 1997 **************************
subroutine IRRAD (m,n,ndim,np,pl,ta,wa,oa,ts,co2,
* n2o,ch4,cfc11,cfc12,cfc22,emiss,
* cldwater,cwc,taucl,reff,fcld,ict,icb,
* taual,ssaal,asyal,
* high,trace,flx,flc,dfdts,st4)
c***********************************************************************
c
c Version IR-5 (September, 1998)
c
c New features of this version are:
c (1) The effect of aerosol scattering on LW fluxes is included.
c (2) Absorption and scattering of rain drops are included.
c
c***********************************************************************
c
c Version IR-4 (October, 1997)
c
c New features of this version are:
c (1) The surface is treated as non-black. The surface
c emissivity, emiss, is an input parameter
c (2) The effect of cloud scattering on LW fluxes is included
c
c*********************************************************************
c
c This routine computes ir fluxes due to water vapor, co2, o3,
c trace gases (n2o, ch4, cfc11, cfc12, cfc22, co2-minor),
c clouds, and aerosols.
c
c This is a vectorized code. It computes fluxes simultaneously for
c (m x n) soundings, which is a subset of (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 Some detailed descriptions of the radiation routine are given in
c Chou and Suarez (1994).
c
c Ice and liquid cloud particles are allowed to co-exist in any of the
c np layers.
c
c If no information is available for the effective cloud particle size,
c reff, default values of 10 micron for liquid water and 75 micron
c for ice can be used.
c
c The maximum-random assumption is applied for cloud overlapping.
c clouds are grouped into high, middle, and low clouds separated by the
c level indices ict and icb. Within each of the three groups, clouds
c are assumed maximally overlapped, and the cloud cover of a group is
c the maximum cloud cover of all the layers in the group. clouds among
c the three groups are assumed randomly overlapped. The indices ict and
c icb correpond approximately to the 400 mb and 700 mb levels.
c
c Aerosols are allowed to be in any of the np layers. Aerosol optical
c properties can be specified as functions of height and spectral band.
c
c There are options for computing fluxes:
c
c If cldwater=.true., taucl is computed from cwc and reff as a
c function of height and spectral band.
c If cldwater=.false., taucl must be given as input to the radiation
c routine. It is independent of spectral band.
c
c If high = .true., transmission functions in the co2, o3, and the
c three water vapor bands with strong absorption are computed using
c table look-up. cooling rates are computed accurately from the
c surface up to 0.01 mb.
c If high = .false., transmission functions are computed using the
c k-distribution method with linear pressure scaling for all spectral
c bands and gases. cooling rates are not calculated accurately for
c pressures less than 20 mb. the computation is faster with
c high=.false. than with high=.true.
c If trace = .true., absorption due to n2o, ch4, cfcs, and the
c two minor co2 bands in the window region is included.
c If trace = .false., absorption in minor bands is neglected.
c
c The IR spectrum is divided into nine bands:
c
c band wavenumber (/cm) absorber
c
c 1 0 - 340 h2o
c 2 340 - 540 h2o
c 3 540 - 800 h2o,cont,co2,n2o
c 4 800 - 980 h2o,cont
c co2,f11,f12,f22
c 5 980 - 1100 h2o,cont,o3
c co2,f11
c 6 1100 - 1215 h2o,cont
c n2o,ch4,f12,f22
c 7 1215 - 1380 h2o
c n2o,ch4
c 8 1380 - 1900 h2o
c 9 1900 - 3000 h2o
c
c In addition, a narrow band in the 15 micrometer region is added to
c compute flux reduction due to n2o
c
c 10 540 - 620 h2o,cont,co2,n2o
c
c Band 3 (540-800/cm) is further divided into 3 sub-bands :
c
c subband wavenumber (/cm)
c
c 1 540 - 620
c 2 620 - 720
c 3 720 - 800
c
c---- Input parameters units size
c
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
c meridional direction (ndim>=n) n/d 1
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) g/g m*ndim*np
c layer ozone mixing ratio by mass (oa) g/g m*ndim*np
c surface temperature (ts) k m*ndim
c co2 mixing ratio by volumn (co2) pppv 1
c n2o mixing ratio by volumn (n2o) pppv np
c ch4 mixing ratio by volumn (ch4) pppv np
c cfc11 mixing ratio by volumn (cfc11) pppv 1
c cfc12 mixing ratio by volumn (cfc12) pppv 1
c cfc22 mixing ratio by volumn (cfc22) pppv 1
c surface emissivity (emiss) fraction m*ndim*10
c input option for cloud optical thickness n/d 1
c cldwater="true" if cwc is provided
c cldwater="false" if taucl is provided
c cloud water mixing ratio (cwc) gm/gm m*ndim*np*3
c index 1 for ice particles
c index 2 for liquid drops
c index 3 for rain drops
c cloud optical thickness (taucl) n/d m*ndim*np*3
c index 1 for ice particles
c index 2 for liquid drops
c index 3 for rain drops
c effective cloud-particle size (reff) micrometer m*ndim*np*3
c index 1 for ice particles
c index 2 for liquid drops
c index 3 for rain 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 n/d 1
c clouds (icb)
c aerosol optical thickness (taual) n/d m*ndim*np*10
c aerosol single-scattering albedo (ssaal) n/d m*ndim*np*10
c aerosol asymmetry factor (asyal) n/d m*ndim*np*10
c high (see explanation above) 1
c trace (see explanation above) 1
c
c Data used in table look-up for transmittance calculations:
c
c c1 , c2, c3: for co2 (band 3)
c o1 , o2, o3: for o3 (band 5)
c h11,h12,h13: for h2o (band 1)
c h21,h22,h23: for h2o (band 2)
c h81,h82,h83: for h2o (band 8)
c
c---- output parameters
c
c net downward flux, all-sky (flx) w/m**2 m*ndim*(np+1)
c net downward flux, clear-sky (flc) w/m**2 m*ndim*(np+1)
c sensitivity of net downward flux
c to surface temperature (dfdts) w/m**2/k m*ndim*(np+1)
c emission by the surface (st4) w/m**2 m*ndim
c
c Notes:
c
c (1) Water vapor continuum absorption is included in 540-1380 /cm.
c (2) Scattering is parameterized for clouds and aerosols.
c (3) Diffuse cloud and aerosol transmissions are computed
c from exp(-1.66*tau).
c (4) If there are no clouds, flx=flc.
c (5) plevel(1) is the pressure at the top of the model atmosphere,
c and plevel(np+1) is the surface pressure.
c (6) Downward flux is positive and upward flux is negative.
c (7) dfdts is negative because upward flux is defined as negative.
c (8) For questions and coding errors, contact with Ming-Dah Chou,
c Code 913, NASA/Goddard Space Flight Center, Greenbelt, MD 20771.
c Phone: 301-286-4012, Fax: 301-286-1759,
c e-mail: chou@climate.gsfc.nasa.gov
c
c-----parameters defining the size of the pre-computed tables for transmittance
c calculations using table look-up.
c
c "nx" is the number of intervals in pressure
c "no" is the number of intervals in o3 amount
c "nc" is the number of intervals in co2 amount
c "nh" is the number of intervals in h2o amount
c "nt" is the number of copies to be made from the pre-computed
c transmittance tables to reduce "memory-bank conflict"
c in parallel machines and, hence, enhancing the speed of
c computations using table look-up.
c If such advantage does not exist, "nt" can be set to 1.
c***************************************************************************
cfpp$ expand (h2oexps)
cfpp$ expand (conexps)
cfpp$ expand (co2exps)
cfpp$ expand (n2oexps)
cfpp$ expand (ch4exps)
cfpp$ expand (comexps)
cfpp$ expand (cfcexps)
cfpp$ expand (b10exps)
cfpp$ expand (tablup)
cfpp$ expand (h2okdis)
cfpp$ expand (co2kdis)
cfpp$ expand (n2okdis)
cfpp$ expand (ch4kdis)
cfpp$ expand (comkdis)
cfpp$ expand (cfckdis)
cfpp$ expand (b10kdis)
implicit none
integer nx,no,nc,nh,nt
integer i,j,k,ip,iw,it,ib,ik,iq,isb,k1,k2
parameter (nx=26,no=21,nc=24,nh=31,nt=7)
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),
* ts(m,ndim)
_RL co2,n2o(np),ch4(np),cfc11,cfc12,cfc22,emiss(m,ndim,10)
_RL cwc(m,ndim,np,3),taucl(m,ndim,np,3),reff(m,ndim,np,3),
* fcld(m,ndim,np)
_RL taual(m,ndim,np,10),ssaal(m,ndim,np,10),asyal(m,ndim,np,10)
logical cldwater,high,trace
c---- output parameters ------
_RL flx(m,ndim,np+1),flc(m,ndim,np+1),dfdts(m,ndim,np+1),
* st4(m,ndim)
c---- static data -----
_RL cb(5,10)
_RL xkw(9),aw(9),bw(9),pm(9),fkw(6,9),gkw(6,3),xke(9)
_RL aib(3,10),awb(4,10),aiw(4,10),aww(4,10),aig(4,10),awg(4,10)
integer ne(9),mw(9)
c---- temporary arrays -----
_RL pa(m,n,np),dt(m,n,np)
_RL sh2o(m,n,np+1),swpre(m,n,np+1),swtem(m,n,np+1)
_RL sco3(m,n,np+1),scopre(m,n,np+1),scotem(m,n,np+1)
_RL dh2o(m,n,np),dcont(m,n,np),dco2(m,n,np),do3(m,n,np)
_RL dn2o(m,n,np),dch4(m,n,np)
_RL df11(m,n,np),df12(m,n,np),df22(m,n,np)
_RL th2o(m,n,6),tcon(m,n,3),tco2(m,n,6,2)
_RL tn2o(m,n,4),tch4(m,n,4),tcom(m,n,2)
_RL tf11(m,n),tf12(m,n),tf22(m,n)
_RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
_RL n2oexp(m,n,np,4),ch4exp(m,n,np,4),comexp(m,n,np,2)
_RL f11exp(m,n,np), f12exp(m,n,np), f22exp(m,n,np)
_RL clr(m,n,0:np+1),fclr(m,n)
_RL blayer(m,n,0:np+1),dlayer(m,n,np+1),dbs(m,n)
_RL clrlw(m,n),clrmd(m,n),clrhi(m,n)
_RL cwp(m,n,np,3)
_RL trant(m,n),tranal(m,n),transfc(m,n,np+1),trantcr(m,n,np+1)
_RL flxu(m,n,np+1),flxd(m,n,np+1),flcu(m,n,np+1),flcd(m,n,np+1)
_RL rflx(m,n,np+1),rflc(m,n,np+1)
logical oznbnd,co2bnd,h2otbl,conbnd,n2obnd
logical ch4bnd,combnd,f11bnd,f12bnd,f22bnd,b10bnd
_RL c1 (nx,nc,nt),c2 (nx,nc,nt),c3 (nx,nc,nt)
_RL o1 (nx,no,nt),o2 (nx,no,nt),o3 (nx,no,nt)
_RL h11(nx,nh,nt),h12(nx,nh,nt),h13(nx,nh,nt)
_RL h21(nx,nh,nt),h22(nx,nh,nt),h23(nx,nh,nt)
_RL h81(nx,nh,nt),h82(nx,nh,nt),h83(nx,nh,nt)
_RL dp,xx,p1,dwe,dpe,a1,b1,fk1,a2,b2,fk2
_RL w1,w2,w3,g1,g2,g3,ww,gg,ff,taux,reff1,reff2
c-----the following coefficients (equivalent to table 2 of
c chou and suarez, 1995) are for computing spectrally
c integrated planck fluxes using eq. (22)
data cb/
1 -2.6844e-1,-8.8994e-2, 1.5676e-3,-2.9349e-6, 2.2233e-9,
2 3.7315e+1,-7.4758e-1, 4.6151e-3,-6.3260e-6, 3.5647e-9,
3 3.7187e+1,-3.9085e-1,-6.1072e-4, 1.4534e-5,-1.6863e-8,
4 -4.1928e+1, 1.0027e+0,-8.5789e-3, 2.9199e-5,-2.5654e-8,
5 -4.9163e+1, 9.8457e-1,-7.0968e-3, 2.0478e-5,-1.5514e-8,
6 -4.7107e+1, 8.9130e-1,-5.9735e-3, 1.5596e-5,-9.5911e-9,
7 -5.4041e+1, 9.5332e-1,-5.7784e-3, 1.2555e-5,-2.9377e-9,
8 -6.9233e+0,-1.5878e-1, 3.9160e-3,-2.4496e-5, 4.9301e-8,
9 1.1483e+2,-2.2376e+0, 1.6394e-2,-5.3672e-5, 6.6456e-8,
* 1.9668e+1,-3.1161e-1, 1.2886e-3, 3.6835e-7,-1.6212e-9/
c-----xkw are the absorption coefficients for the first
c k-distribution function due to water vapor line absorption
c (tables 4 and 7). units are cm**2/g
data xkw / 29.55 , 4.167e-1, 1.328e-2, 5.250e-4,
* 5.25e-4, 9.369e-3, 4.719e-2, 1.320e-0, 5.250e-4/
c-----xke are the absorption coefficients for the first
c k-distribution function due to water vapor continuum absorption
c (table 6). units are cm**2/g
data xke / 0.00, 0.00, 27.40, 15.8,
* 9.40, 7.75, 0.0, 0.0, 0.0/
c-----mw are the ratios between neighboring absorption coefficients
c for water vapor line absorption (tables 4 and 7).
data mw /6,6,8,6,6,8,9,6,16/
c-----aw and bw (table 3) are the coefficients for temperature scaling
c in eq. (25).
data aw/ 0.0021, 0.0140, 0.0167, 0.0302,
* 0.0307, 0.0195, 0.0152, 0.0008, 0.0096/
data bw/ -1.01e-5, 5.57e-5, 8.54e-5, 2.96e-4,
* 2.86e-4, 1.108e-4, 7.608e-5, -3.52e-6, 1.64e-5/
data pm/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.77, 0.5, 1.0, 1.0/
c-----fkw is the planck-weighted k-distribution function due to h2o
c line absorption given in table 4 of Chou and Suarez (1995).
c the k-distribution function for the third band, fkw(*,3), is not used
data fkw / 0.2747,0.2717,0.2752,0.1177,0.0352,0.0255,
2 0.1521,0.3974,0.1778,0.1826,0.0374,0.0527,
3 6*1.00,
4 0.4654,0.2991,0.1343,0.0646,0.0226,0.0140,
5 0.5543,0.2723,0.1131,0.0443,0.0160,0.0000,
6 0.5955,0.2693,0.0953,0.0335,0.0064,0.0000,
7 0.1958,0.3469,0.3147,0.1013,0.0365,0.0048,
8 0.0740,0.1636,0.4174,0.1783,0.1101,0.0566,
9 0.1437,0.2197,0.3185,0.2351,0.0647,0.0183/
c-----gkw is the planck-weighted k-distribution function due to h2o
c line absorption in the 3 subbands (800-720,620-720,540-620 /cm)
c of band 3 given in table 7. Note that the order of the sub-bands
c is reversed.
data gkw/ 0.1782,0.0593,0.0215,0.0068,0.0022,0.0000,
2 0.0923,0.1675,0.0923,0.0187,0.0178,0.0000,
3 0.0000,0.1083,0.1581,0.0455,0.0274,0.0041/
c-----ne is the number of terms used in each band to compute water vapor
c continuum transmittance (table 6).
data ne /0,0,3,1,1,1,0,0,0/
c
c-----coefficients for computing the extinction coefficient
c for cloud ice particles
c polynomial fit: y=a1+a2/x**a3; x is in m**2/gm
c
data aib / -0.44171, 0.62951, 0.06465,
2 -0.13727, 0.61291, 0.28962,
3 -0.01878, 1.67680, 0.79080,
4 -0.01896, 1.06510, 0.69493,
5 -0.04788, 0.88178, 0.54492,
6 -0.02265, 1.57390, 0.76161,
7 -0.01038, 2.15640, 0.89045,
8 -0.00450, 2.51370, 0.95989,
9 -0.00044, 3.15050, 1.03750,
* -0.02956, 1.44680, 0.71283/
c
c-----coefficients for computing the extinction coefficient
c for cloud liquid drops
c polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
c
data awb / 0.08641, 0.01769, -1.5572e-3, 3.4896e-5,
2 0.22027, 0.00997, -1.8719e-3, 5.3112e-5,
3 0.39252, -0.02817, 7.2931e-4, -3.8151e-6,
4 0.39975, -0.03426, 1.2884e-3, -1.7930e-5,
5 0.34021, -0.02805, 1.0654e-3, -1.5443e-5,
6 0.15587, 0.00371, -7.7705e-4, 2.0547e-5,
7 0.05518, 0.04544, -4.2067e-3, 1.0184e-4,
8 0.12724, 0.04751, -5.2037e-3, 1.3711e-4,
9 0.30390, 0.01656, -3.5271e-3, 1.0828e-4,
* 0.63617, -0.06287, 2.2350e-3, -2.3177e-5/
c
c-----coefficients for computing the single-scattering albedo
c for cloud ice particles
c polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
c
data aiw/ 0.17201, 1.2229e-2, -1.4837e-4, 5.8020e-7,
2 0.81470, -2.7293e-3, 9.7816e-8, 5.7650e-8,
3 0.54859, -4.8273e-4, 5.4353e-6, -1.5679e-8,
4 0.39218, 4.1717e-3, - 4.8869e-5, 1.9144e-7,
5 0.71773, -3.3640e-3, 1.9713e-5, -3.3189e-8,
6 0.77345, -5.5228e-3, 4.8379e-5, -1.5151e-7,
7 0.74975, -5.6604e-3, 5.6475e-5, -1.9664e-7,
8 0.69011, -4.5348e-3, 4.9322e-5, -1.8255e-7,
9 0.83963, -6.7253e-3, 6.1900e-5, -2.0862e-7,
* 0.64860, -2.8692e-3, 2.7656e-5, -8.9680e-8/
c
c-----coefficients for computing the single-scattering albedo
c for cloud liquid drops
c polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
c
data aww/ -7.8566e-2, 8.0875e-2, -4.3403e-3, 8.1341e-5,
2 -1.3384e-2, 9.3134e-2, -6.0491e-3, 1.3059e-4,
3 3.7096e-2, 7.3211e-2, -4.4211e-3, 9.2448e-5,
4 -3.7600e-3, 9.3344e-2, -5.6561e-3, 1.1387e-4,
5 0.40212, 7.8083e-2, -5.9583e-3, 1.2883e-4,
6 0.57928, 5.9094e-2, -5.4425e-3, 1.2725e-4,
7 0.68974, 4.2334e-2, -4.9469e-3, 1.2863e-4,
8 0.80122, 9.4578e-3, -2.8508e-3, 9.0078e-5,
9 1.02340, -2.6204e-2, 4.2552e-4, 3.2160e-6,
* 0.05092, 7.5409e-2, -4.7305e-3, 1.0121e-4/
c
c-----coefficients for computing the asymmetry factor for cloud ice particles
c polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
c
data aig / 0.57867, 1.0135e-2, -1.1142e-4, 4.1537e-7,
2 0.72259, 3.1149e-3, -1.9927e-5, 5.6024e-8,
3 0.76109, 4.5449e-3, -4.6199e-5, 1.6446e-7,
4 0.86934, 2.7474e-3, -3.1301e-5, 1.1959e-7,
5 0.89103, 1.8513e-3, -1.6551e-5, 5.5193e-8,
6 0.86325, 2.1408e-3, -1.6846e-5, 4.9473e-8,
7 0.85064, 2.5028e-3, -2.0812e-5, 6.3427e-8,
8 0.86945, 2.4615e-3, -2.3882e-5, 8.2431e-8,
9 0.80122, 3.1906e-3, -2.4856e-5, 7.2411e-8,
* 0.73290, 4.8034e-3, -4.4425e-5, 1.4839e-7/
c
c-----coefficients for computing the asymmetry factor for cloud liquid drops
c polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
c
data awg / -0.51930, 0.20290, -1.1747e-2, 2.3868e-4,
2 -0.22151, 0.19708, -1.2462e-2, 2.6646e-4,
3 0.14157, 0.14705, -9.5802e-3, 2.0819e-4,
4 0.41590, 0.10482, -6.9118e-3, 1.5115e-4,
5 0.55338, 7.7016e-2, -5.2218e-3, 1.1587e-4,
6 0.61384, 6.4402e-2, -4.6241e-3, 1.0746e-4,
7 0.67891, 4.8698e-2, -3.7021e-3, 9.1966e-5,
8 0.78169, 2.0803e-2, -1.4749e-3, 3.9362e-5,
9 0.93218, -3.3425e-2, 2.9632e-3, -6.9362e-5,
* 0.01649, 0.16561, -1.0723e-2, 2.3220e-4/
c
c-----include tables used in the table look-up for co2 (band 3),
c o3 (band 5), and h2o (bands 1, 2, and 7) transmission functions.
logical first
data first /.true./
#include "h2o-tran3.h"
#include "co2-tran3.h"
#include "o3-tran3.h"
c save c1,c2,c3,o1,o2,o3
c save h11,h12,h13,h21,h22,h23,h81,h82,h83
if (first) then
c-----tables co2 and h2o are only used with 'high' option
if (high) then
do iw=1,nh
do ip=1,nx
h11(ip,iw,1)=1.0-h11(ip,iw,1)
h21(ip,iw,1)=1.0-h21(ip,iw,1)
h81(ip,iw,1)=1.0-h81(ip,iw,1)
enddo
enddo
do iw=1,nc
do ip=1,nx
c1(ip,iw,1)=1.0-c1(ip,iw,1)
enddo
enddo
c-----copy tables to enhance the speed of co2 (band 3), o3 (band 5),
c and h2o (bands 1, 2, and 8 only) transmission calculations
c using table look-up.
c-----tables are replicated to avoid memory bank conflicts
do it=2,nt
do iw=1,nc
do ip=1,nx
c1 (ip,iw,it)= c1(ip,iw,1)
c2 (ip,iw,it)= c2(ip,iw,1)
c3 (ip,iw,it)= c3(ip,iw,1)
enddo
enddo
do iw=1,nh
do ip=1,nx
h11(ip,iw,it)=h11(ip,iw,1)
h12(ip,iw,it)=h12(ip,iw,1)
h13(ip,iw,it)=h13(ip,iw,1)
h21(ip,iw,it)=h21(ip,iw,1)
h22(ip,iw,it)=h22(ip,iw,1)
h23(ip,iw,it)=h23(ip,iw,1)
h81(ip,iw,it)=h81(ip,iw,1)
h82(ip,iw,it)=h82(ip,iw,1)
h83(ip,iw,it)=h83(ip,iw,1)
enddo
enddo
enddo
endif
c-----always use table look-up for ozone transmittance
do iw=1,no
do ip=1,nx
o1(ip,iw,1)=1.0-o1(ip,iw,1)
enddo
enddo
do it=2,nt
do iw=1,no
do ip=1,nx
o1 (ip,iw,it)= o1(ip,iw,1)
o2 (ip,iw,it)= o2(ip,iw,1)
o3 (ip,iw,it)= o3(ip,iw,1)
enddo
enddo
enddo
first=.false.
endif
c-----set the pressure at the top of the model atmosphere
c to 1.0e-4 if it is zero
do j=1,n
do i=1,m
if (pl(i,j,1) .eq. 0.0) then
pl(i,j,1)=1.0e-4
endif
enddo
enddo
c-----compute layer pressure (pa) and layer temperature minus 250K (dt)
do k=1,np
do j=1,n
do i=1,m
pa(i,j,k)=0.5*(pl(i,j,k)+pl(i,j,k+1))
dt(i,j,k)=ta(i,j,k)-250.0
enddo
enddo
enddo
c-----compute layer absorber amount
c dh2o : water vapor amount (g/cm**2)
c dcont: scaled water vapor amount for continuum absorption
c (g/cm**2)
c dco2 : co2 amount (cm-atm)stp
c do3 : o3 amount (cm-atm)stp
c dn2o : n2o amount (cm-atm)stp
c dch4 : ch4 amount (cm-atm)stp
c df11 : cfc11 amount (cm-atm)stp
c df12 : cfc12 amount (cm-atm)stp
c df22 : cfc22 amount (cm-atm)stp
c the factor 1.02 is equal to 1000/980
c factors 789 and 476 are for unit conversion
c the factor 0.001618 is equal to 1.02/(.622*1013.25)
c the factor 6.081 is equal to 1800/296
do k=1,np
do j=1,n
do i=1,m
dp = pl(i,j,k+1)-pl(i,j,k)
dh2o (i,j,k) = 1.02*wa(i,j,k)*dp+1.e-10
do3 (i,j,k) = 476.*oa(i,j,k)*dp+1.e-10
dco2 (i,j,k) = 789.*co2*dp+1.e-10
dch4 (i,j,k) = 789.*ch4(k)*dp+1.e-10
dn2o (i,j,k) = 789.*n2o(k)*dp+1.e-10
df11 (i,j,k) = 789.*cfc11*dp+1.e-10
df12 (i,j,k) = 789.*cfc12*dp+1.e-10
df22 (i,j,k) = 789.*cfc22*dp+1.e-10
c-----compute scaled water vapor amount for h2o continuum absorption
c following eq. (43).
xx=pa(i,j,k)*0.001618*wa(i,j,k)*wa(i,j,k)*dp
dcont(i,j,k) = xx*exp(1800./ta(i,j,k)-6.081)+1.e-10
enddo
enddo
enddo
c-----compute column-integrated h2o amoumt, h2o-weighted pressure
c and temperature. it follows eqs. (37) and (38).
if (high) then
call COLUMN(m,n,np,pa,dt,dh2o,sh2o,swpre,swtem)
endif
c-----compute layer cloud water amount (gm/m**2)
c index is 1 for ice, 2 for waterdrops and 3 for raindrops.
c Rain optical thickness is 0.00307 /(gm/m**2).
c It is for a specific drop size distribution provided by Q. Fu.
if (cldwater) then
do k=1,np
do j=1,n
do i=1,m
dp =pl(i,j,k+1)-pl(i,j,k)
cwp(i,j,k,1)=1.02*10000.0*cwc(i,j,k,1)*dp
cwp(i,j,k,2)=1.02*10000.0*cwc(i,j,k,2)*dp
cwp(i,j,k,3)=1.02*10000.0*cwc(i,j,k,3)*dp
taucl(i,j,k,3)=0.00307*cwp(i,j,k,3)
enddo
enddo
enddo
endif
c-----the surface (np+1) is treated as a layer filled with black clouds.
c clr is the equivalent clear fraction of a layer
c transfc is the transmttance between the surface and a pressure level.
c trantcr is the clear-sky transmttance between the surface and a
c pressure level.
do j=1,n
do i=1,m
clr(i,j,0) = 1.0
clr(i,j,np+1) = 0.0
st4(i,j) = 0.0
transfc(i,j,np+1)=1.0
trantcr(i,j,np+1)=1.0
enddo
enddo
c-----initialize fluxes
do k=1,np+1
do j=1,n
do i=1,m
flx(i,j,k) = 0.0
flc(i,j,k) = 0.0
dfdts(i,j,k)= 0.0
rflx(i,j,k) = 0.0
rflc(i,j,k) = 0.0
enddo
enddo
enddo
c-----integration over spectral bands
do 1000 ib=1,10
c-----if h2otbl, compute h2o (line) transmittance using table look-up.
c if conbnd, compute h2o (continuum) transmittance in bands 3-6.
c if co2bnd, compute co2 transmittance in band 3.
c if oznbnd, compute o3 transmittance in band 5.
c if n2obnd, compute n2o transmittance in bands 6 and 7.
c if ch4bnd, compute ch4 transmittance in bands 6 and 7.
c if combnd, compute co2-minor transmittance in bands 4 and 5.
c if f11bnd, compute cfc11 transmittance in bands 4 and 5.
c if f12bnd, compute cfc12 transmittance in bands 4 and 6.
c if f22bnd, compute cfc22 transmittance in bands 4 and 6.
c if b10bnd, compute flux reduction due to n2o in band 10.
h2otbl=high.and.(ib.eq.1.or.ib.eq.2.or.ib.eq.8)
conbnd=ib.ge.3.and.ib.le.6
co2bnd=ib.eq.3
oznbnd=ib.eq.5
n2obnd=ib.eq.6.or.ib.eq.7
ch4bnd=ib.eq.6.or.ib.eq.7
combnd=ib.eq.4.or.ib.eq.5
f11bnd=ib.eq.4.or.ib.eq.5
f12bnd=ib.eq.4.or.ib.eq.6
f22bnd=ib.eq.4.or.ib.eq.6
b10bnd=ib.eq.10
c-----blayer is the spectrally integrated planck flux of the mean layer
c temperature derived from eq. (22)
c the fitting for the planck flux is valid in the range 160-345 K.
do k=1,np
do j=1,n
do i=1,m
blayer(i,j,k)=ta(i,j,k)*(ta(i,j,k)*(ta(i,j,k)
* *(ta(i,j,k)*cb(5,ib)+cb(4,ib))+cb(3,ib))
* +cb(2,ib))+cb(1,ib)
enddo
enddo
enddo
c-----the earth surface, with an index "np+1", is treated as a layer
do j=1,n
do i=1,m
blayer(i,j,0) = 0.0
blayer(i,j,np+1)= ( ts(i,j)*(ts(i,j)*(ts(i,j)
* *(ts(i,j)*cb(5,ib)+cb(4,ib))+cb(3,ib))
* +cb(2,ib))+cb(1,ib) )*emiss(i,j,ib)
c-----dbs is the derivative of the surface emission with respect to
c surface temperature (eq. 59).
dbs(i,j)=(ts(i,j)*(ts(i,j)*(ts(i,j)*4.*cb(5,ib)
* +3.*cb(4,ib))+2.*cb(3,ib))+cb(2,ib))*emiss(i,j,ib)
enddo
enddo
do k=1,np+1
do j=1,n
do i=1,m
dlayer(i,j,k)=blayer(i,j,k-1)-blayer(i,j,k)
enddo
enddo
enddo
c-----compute column-integrated absorber amoumt, absorber-weighted
c pressure and temperature for co2 (band 3) and o3 (band 5).
c it follows eqs. (37) and (38).
c-----this is in the band loop to save storage
if (high .and. co2bnd) then
call COLUMN(m,n,np,pa,dt,dco2,sco3,scopre,scotem)
endif
if (oznbnd) then
call COLUMN(m,n,np,pa,dt,do3,sco3,scopre,scotem)
endif
c-----compute cloud optical thickness
if (cldwater) then
do k= 1, np
do j= 1, n
do i= 1, m
taucl(i,j,k,1)=cwp(i,j,k,1)*(aib(1,ib)+aib(2,ib)/
* reff(i,j,k,1)**aib(3,ib))
taucl(i,j,k,2)=cwp(i,j,k,2)*(awb(1,ib)+(awb(2,ib)+(awb(3,ib)
* +awb(4,ib)*reff(i,j,k,2))*reff(i,j,k,2))*reff(i,j,k,2))
enddo
enddo
enddo
endif
c-----compute cloud single-scattering albedo and asymmetry factor for
c a mixture of ice particles, liquid drops, and rain drops
c single-scattering albedo and asymmetry factor of rain are set
c to 0.54 and 0.95, respectively.
do k= 1, np
do j= 1, n
do i= 1, m
clr(i,j,k) = 1.0
taux=taucl(i,j,k,1)+taucl(i,j,k,2)+taucl(i,j,k,3)
if (taux.gt.0.02 .and. fcld(i,j,k).gt.0.01) then
reff1=min(reff(i,j,k,1),130. _d 0)
reff2=min(reff(i,j,k,2),20.0 _d 0)
w1=taucl(i,j,k,1)*(aiw(1,ib)+(aiw(2,ib)+(aiw(3,ib)
* +aiw(4,ib)*reff1)*reff1)*reff1)
w2=taucl(i,j,k,2)*(aww(1,ib)+(aww(2,ib)+(aww(3,ib)
* +aww(4,ib)*reff2)*reff2)*reff2)
w3=taucl(i,j,k,3)*0.54
ww=(w1+w2+w3)/taux
g1=w1*(aig(1,ib)+(aig(2,ib)+(aig(3,ib)
* +aig(4,ib)*reff1)*reff1)*reff1)
g2=w2*(awg(1,ib)+(awg(2,ib)+(awg(3,ib)
* +awg(4,ib)*reff2)*reff2)*reff2)
g3=w3*0.95
gg=(g1+g2+g3)/(w1+w2+w3)
c-----parameterization of LW scattering
c compute forward-scattered fraction and scale optical thickness
ff=0.5+(0.3739+(0.0076+0.1185*gg)*gg)*gg
taux=taux*(1.-ww*ff)
c-----compute equivalent cloud-free fraction, clr, for each layer
c the cloud diffuse transmittance is approximated by using
c a diffusivity factor of 1.66.
clr(i,j,k)=1.0-(fcld(i,j,k)*(1.0-exp(-1.66*taux)))
endif
enddo
enddo
enddo
c-----compute the exponential terms (eq. 32) at each layer due to
c water vapor line absorption when k-distribution is used
if (.not.h2otbl .and. .not.b10bnd) then
call H2OEXPS(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)
endif
c-----compute the exponential terms (eq. 46) at each layer due to
c water vapor continuum absorption
if (conbnd) then
call CONEXPS(ib,m,n,np,dcont,xke,conexp)
endif
c-----compute the exponential terms (eq. 32) at each layer due to
c co2 absorption
if (.not.high .and. co2bnd) then
call CO2EXPS(m,n,np,dco2,pa,dt,co2exp)
endif
c***** for trace gases *****
if (trace) then
c-----compute the exponential terms at each layer due to n2o absorption
if (n2obnd) then
call N2OEXPS(ib,m,n,np,dn2o,pa,dt,n2oexp)
endif
c-----compute the exponential terms at each layer due to ch4 absorption
if (ch4bnd) then
call CH4EXPS(ib,m,n,np,dch4,pa,dt,ch4exp)
endif
c-----compute the exponential terms due to co2 minor absorption
if (combnd) then
call COMEXPS(ib,m,n,np,dco2,dt,comexp)
endif
c-----compute the exponential terms due to cfc11 absorption
if (f11bnd) then
a1 = 1.26610e-3
b1 = 3.55940e-6
fk1 = 1.89736e+1
a2 = 8.19370e-4
b2 = 4.67810e-6
fk2 = 1.01487e+1
call CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df11,dt,f11exp)
endif
c-----compute the exponential terms due to cfc12 absorption
if (f12bnd) then
a1 = 8.77370e-4
b1 =-5.88440e-6
fk1 = 1.58104e+1
a2 = 8.62000e-4
b2 =-4.22500e-6
fk2 = 3.70107e+1
call CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df12,dt,f12exp)
endif
c-----compute the exponential terms due to cfc22 absorption
if (f22bnd) then
a1 = 9.65130e-4
b1 = 1.31280e-5
fk1 = 6.18536e+0
a2 =-3.00010e-5
b2 = 5.25010e-7
fk2 = 3.27912e+1
call CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df22,dt,f22exp)
endif
c-----compute the exponential terms at each layer in band 10 due to
c h2o line and continuum, co2, and n2o absorption
if (b10bnd) then
call B10EXPS(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt
* ,h2oexp,conexp,co2exp,n2oexp)
endif
endif
c-----compute transmittances for regions between levels k1 and k2
c and update the fluxes at the two levels.
c-----initialize fluxes
do k=1,np+1
do j=1,n
do i=1,m
flxu(i,j,k) = 0.0
flxd(i,j,k) = 0.0
flcu(i,j,k) = 0.0
flcd(i,j,k) = 0.0
enddo
enddo
enddo
do 2000 k1=1,np
c-----initialize fclr, th2o, tcon, tco2, and tranal
c fclr is the clear fraction of the line-of-sight
c clrlw, clrmd, and clrhi are the clear-sky fractions of the
c low, middle, and high troposphere, respectively
c tranal is the aerosol transmission function
do j=1,n
do i=1,m
fclr(i,j) = 1.0
clrlw(i,j) = 1.0
clrmd(i,j) = 1.0
clrhi(i,j) = 1.0
tranal(i,j)= 1.0
enddo
enddo
c-----for h2o line transmission
if (.not. h2otbl) then
do ik=1,6
do j=1,n
do i=1,m
th2o(i,j,ik)=1.0
enddo
enddo
enddo
endif
c-----for h2o continuum transmission
if (conbnd) then
do iq=1,3
do j=1,n
do i=1,m
tcon(i,j,iq)=1.0
enddo
enddo
enddo
endif
c-----for co2 transmission using k-distribution method.
c band 3 is divided into 3 sub-bands, but sub-bands 3a and 3c
c are combined in computing the co2 transmittance.
if (.not.high .and. co2bnd) then
do isb=1,2
do ik=1,6
do j=1,n
do i=1,m
tco2(i,j,ik,isb)=1.0
enddo
enddo
enddo
enddo
endif
c***** for trace gases *****
if (trace) then
c-----for n2o transmission using k-distribution method.
if (n2obnd) then
do ik=1,4
do j=1,n
do i=1,m
tn2o(i,j,ik)=1.0
enddo
enddo
enddo
endif
c-----for ch4 transmission using k-distribution method.
if (ch4bnd) then
do ik=1,4
do j=1,n
do i=1,m
tch4(i,j,ik)=1.0
enddo
enddo
enddo
endif
c-----for co2-minor transmission using k-distribution method.
if (combnd) then
do ik=1,2
do j=1,n
do i=1,m
tcom(i,j,ik)=1.0
enddo
enddo
enddo
endif
c-----for cfc-11 transmission using k-distribution method.
if (f11bnd) then
do j=1,n
do i=1,m
tf11(i,j)=1.0
enddo
enddo
endif
c-----for cfc-12 transmission using k-distribution method.
if (f12bnd) then
do j=1,n
do i=1,m
tf12(i,j)=1.0
enddo
enddo
endif
c-----for cfc-22 transmission when using k-distribution method.
if (f22bnd) then
do j=1,n
do i=1,m
tf22(i,j)=1.0
enddo
enddo
endif
c-----for the transmission in band 10 using k-distribution method.
if (b10bnd) then
do ik=1,6
do j=1,n
do i=1,m
th2o(i,j,ik)=1.0
tco2(i,j,ik,1)=1.0
enddo
enddo
enddo
do iq=1,3
do j=1,n
do i=1,m
tcon(i,j,iq)=1.0
enddo
enddo
enddo
do ik=1,4
do j=1,n
do i=1,m
tn2o(i,j,ik)=1.0
enddo
enddo
enddo
endif
endif
c-----loop over the bottom level of the region (k2)
do 3000 k2=k1+1,np+1
c-----initialize total transmission function (trant)
do j=1,n
do i=1,m
trant(i,j)=1.0
enddo
enddo
if (h2otbl) then
w1=-8.0
p1=-2.0
dwe=0.3
dpe=0.2
c-----compute water vapor transmittance using table look-up
if (ib.eq.1) then
call TABLUP(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
* w1,p1,dwe,dpe,h11,h12,h13,trant)
endif
if (ib.eq.2) then
call TABLUP(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
* w1,p1,dwe,dpe,h21,h22,h23,trant)
endif
if (ib.eq.8) then
call TABLUP(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
* w1,p1,dwe,dpe,h81,h82,h83,trant)
endif
else
c-----compute water vapor transmittance using k-distribution
if (.not.b10bnd) then
call H2OKDIS(ib,m,n,np,k2-1,fkw,gkw,ne,h2oexp,conexp,
* th2o,tcon,trant)
endif
endif
if (co2bnd) then
if (high) then
c-----compute co2 transmittance using table look-up method
w1=-4.0
p1=-2.0
dwe=0.3
dpe=0.2
call TABLUP(k1,k2,m,n,np,nx,nc,nt,sco3,scopre,scotem,
* w1,p1,dwe,dpe,c1,c2,c3,trant)
else
c-----compute co2 transmittance using k-distribution method
call CO2KDIS(m,n,np,k2-1,co2exp,tco2,trant)
endif
endif
c-----All use table look-up to compute o3 transmittance.
if (oznbnd) then
w1=-6.0
p1=-2.0
dwe=0.3
dpe=0.2
call TABLUP(k1,k2,m,n,np,nx,no,nt,sco3,scopre,scotem,
* w1,p1,dwe,dpe,o1,o2,o3,trant)
endif
c***** for trace gases *****
if (trace) then
c-----compute n2o transmittance using k-distribution method
if (n2obnd) then
call N2OKDIS(ib,m,n,np,k2-1,n2oexp,tn2o,trant)
endif
c-----compute ch4 transmittance using k-distribution method
if (ch4bnd) then
call CH4KDIS(ib,m,n,np,k2-1,ch4exp,tch4,trant)
endif
c-----compute co2-minor transmittance using k-distribution method
if (combnd) then
call COMKDIS(ib,m,n,np,k2-1,comexp,tcom,trant)
endif
c-----compute cfc11 transmittance using k-distribution method
if (f11bnd) then
call CFCKDIS(m,n,np,k2-1,f11exp,tf11,trant)
endif
c-----compute cfc12 transmittance using k-distribution method
if (f12bnd) then
call CFCKDIS(m,n,np,k2-1,f12exp,tf12,trant)
endif
c-----compute cfc22 transmittance using k-distribution method
if (f22bnd) then
call CFCKDIS(m,n,np,k2-1,f22exp,tf22,trant)
endif
c-----compute transmittance in band 10 using k-distribution method
c here, trant is the change in transmittance due to n2o absorption
if (b10bnd) then
call B10KDIS(m,n,np,k2-1,h2oexp,conexp,co2exp,n2oexp
* ,th2o,tcon,tco2,tn2o,trant)
endif
endif
c
c-----include aerosol effect
c
do j=1,n
do i=1,m
ff=0.5+(0.3739+(0.0076+0.1185*asyal(i,j,k2-1,ib))
* *asyal(i,j,k2-1,ib))*asyal(i,j,k2-1,ib)
taux=taual(i,j,k2-1,ib)*(1.-ssaal(i,j,k2-1,ib)*ff)
tranal(i,j)=tranal(i,j)*exp(-1.66*taux)
trant (i,j)=trant(i,j) *tranal(i,j)
enddo
enddo
c-----Identify the clear-sky fractions clrhi, clrmd and clrlw, in the
c high, middle and low cloud groups.
c fclr is the clear-sky fraction between levels k1 and k2 assuming
c the three cloud groups are randomly overlapped.
do j=1,n
do i=1,m
if( k2 .le. ict ) then
clrhi(i,j)=min(clr(i,j,k2-1),clrhi(i,j))
elseif( k2 .gt. ict .and. k2 .le. icb ) then
clrmd(i,j)=min(clr(i,j,k2-1),clrmd(i,j))
elseif( k2 .gt. icb ) then
clrlw(i,j)=min(clr(i,j,k2-1),clrlw(i,j))
endif
fclr(i,j)=clrlw(i,j)*clrmd(i,j)*clrhi(i,j)
enddo
enddo
c-----compute upward and downward fluxes (band 1-9).
c add "boundary" terms to the net downward flux.
c these are the first terms on the right-hand-side of
c eqs. (56a) and (56b). downward fluxes are positive.
if (.not. b10bnd) then
if (k2 .eq. k1+1) then
do j=1,n
do i=1,m
c-----compute upward and downward fluxes for clear-sky situation
flcu(i,j,k1)=flcu(i,j,k1)-blayer(i,j,k1)
flcd(i,j,k2)=flcd(i,j,k2)+blayer(i,j,k1)
c-----compute upward and downward fluxes for all-sky situation
flxu(i,j,k1)=flxu(i,j,k1)-blayer(i,j,k1)
flxd(i,j,k2)=flxd(i,j,k2)+blayer(i,j,k1)
enddo
enddo
endif
c-----add flux components involving the four layers above and below
c the levels k1 and k2. it follows eqs. (56a) and (56b).
do j=1,n
do i=1,m
xx=trant(i,j)*dlayer(i,j,k2)
flcu(i,j,k1) =flcu(i,j,k1)+xx
flxu(i,j,k1) =flxu(i,j,k1)+xx*fclr(i,j)
xx=trant(i,j)*dlayer(i,j,k1)
flcd(i,j,k2) =flcd(i,j,k2)+xx
flxd(i,j,k2) =flxd(i,j,k2)+xx*fclr(i,j)
enddo
enddo
else
c-----flux reduction due to n2o in band 10
if (trace) then
do j=1,n
do i=1,m
rflx(i,j,k1) = rflx(i,j,k1)+trant(i,j)*fclr(i,j)
* *dlayer(i,j,k2)
rflx(i,j,k2) = rflx(i,j,k2)+trant(i,j)*fclr(i,j)
* *dlayer(i,j,k1)
rflc(i,j,k1) = rflc(i,j,k1)+trant(i,j)
* *dlayer(i,j,k2)
rflc(i,j,k2) = rflc(i,j,k2)+trant(i,j)
* *dlayer(i,j,k1)
enddo
enddo
endif
endif
3000 continue
c-----transmission between level k1 and the surface
do j=1,n
do i=1,m
trantcr(i,j,k1) =trant(i,j)
transfc(i,j,k1) =trant(i,j)*fclr(i,j)
enddo
enddo
c-----compute the partial derivative of fluxes with respect to
c surface temperature (eq. 59)
if (trace .or. (.not. b10bnd) ) then
do j=1,n
do i=1,m
dfdts(i,j,k1) =dfdts(i,j,k1)-dbs(i,j)*transfc(i,j,k1)
enddo
enddo
endif
2000 continue
if (.not. b10bnd) then
c-----add contribution from the surface to the flux terms at the surface
do j=1,n
do i=1,m
flcu(i,j,np+1)=flcu(i,j,np+1)-blayer(i,j,np+1)
flxu(i,j,np+1)=flxu(i,j,np+1)-blayer(i,j,np+1)
st4(i,j)=st4(i,j)-blayer(i,j,np+1)
dfdts(i,j,np+1)=dfdts(i,j,np+1)-dbs(i,j)
enddo
enddo
c-----add the flux component reflected by the surface
c note: upward flux is negative
do k=1, np+1
do j=1,n
do i=1,m
flcu(i,j,k)=flcu(i,j,k)-
* flcd(i,j,np+1)*trantcr(i,j,k)*(1.-emiss(i,j,ib))
flxu(i,j,k)=flxu(i,j,k)-
* flxd(i,j,np+1)*transfc(i,j,k)*(1.-emiss(i,j,ib))
enddo
enddo
enddo
endif
c-----summation of fluxes over spectral bands
do k=1,np+1
do j=1,n
do i=1,m
flc(i,j,k)=flc(i,j,k)+flcd(i,j,k)+flcu(i,j,k)
flx(i,j,k)=flx(i,j,k)+flxd(i,j,k)+flxu(i,j,k)
enddo
enddo
enddo
1000 continue
c-----adjust fluxes due to n2o absorption in band 10
do k=1,np+1
do j=1,n
do i=1,m
flc(i,j,k)=flc(i,j,k)+rflc(i,j,k)
flx(i,j,k)=flx(i,j,k)+rflx(i,j,k)
enddo
enddo
enddo
return
end
c***********************************************************************
subroutine COLUMN (m,n,np,pa,dt,sabs0,sabs,spre,stem)
c***********************************************************************
c-----compute column-integrated (from top of the model atmosphere)
c absorber amount (sabs), absorber-weighted pressure (spre) and
c temperature (stem).
c computations of spre and stem follows eqs. (37) and (38).
c
c--- input parameters
c number of soundings in zonal direction (m)
c number of soundings in meridional direction (n)
c number of atmospheric layers (np)
c layer pressure (pa)
c layer temperature minus 250K (dt)
c layer absorber amount (sabs0)
c
c--- output parameters
c column-integrated absorber amount (sabs)
c column absorber-weighted pressure (spre)
c column absorber-weighted temperature (stem)
c
c--- units of pa and dt are mb and k, respectively.
c units of sabs are g/cm**2 for water vapor and (cm-atm)stp
c for co2 and o3
c***********************************************************************
implicit none
integer m,n,np,i,j,k
c---- input parameters -----
_RL pa(m,n,np),dt(m,n,np),sabs0(m,n,np)
c---- output parameters -----
_RL sabs(m,n,np+1),spre(m,n,np+1),stem(m,n,np+1)
c*********************************************************************
do j=1,n
do i=1,m
sabs(i,j,1)=0.0
spre(i,j,1)=0.0
stem(i,j,1)=0.0
enddo
enddo
do k=1,np
do j=1,n
do i=1,m
sabs(i,j,k+1)=sabs(i,j,k)+sabs0(i,j,k)
spre(i,j,k+1)=spre(i,j,k)+pa(i,j,k)*sabs0(i,j,k)
stem(i,j,k+1)=stem(i,j,k)+dt(i,j,k)*sabs0(i,j,k)
enddo
enddo
enddo
return
end
c**********************************************************************
subroutine H2OEXPS(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)
c**********************************************************************
c compute exponentials for water vapor line absorption
c in individual layers.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c layer water vapor amount for line absorption (dh2o)
c layer pressure (pa)
c layer temperature minus 250K (dt)
c absorption coefficients for the first k-distribution
c function due to h2o line absorption (xkw)
c coefficients for the temperature and pressure scaling (aw,bw,pm)
c ratios between neighboring absorption coefficients for
c h2o line absorption (mw)
c
c---- output parameters
c 6 exponentials for each layer (h2oexp)
c**********************************************************************
implicit none
integer ib,m,n,np,i,j,k,ik
c---- input parameters ------
_RL dh2o(m,n,np),pa(m,n,np),dt(m,n,np)
c---- output parameters -----
_RL h2oexp(m,n,np,6)
c---- static data -----
integer mw(9)
_RL xkw(9),aw(9),bw(9),pm(9)
c---- temporary arrays -----
_RL xh,xh1
c**********************************************************************
c note that the 3 sub-bands in band 3 use the same set of xkw, aw,
c and bw, therefore, h2oexp for these sub-bands are identical.
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
c-----xh is the scaled water vapor amount for line absorption
c computed from (27)
xh = dh2o(i,j,k)*(pa(i,j,k)/500.)**pm(ib)
1 * ( 1.+(aw(ib)+bw(ib)* dt(i,j,k))*dt(i,j,k) )
c-----h2oexp is the water vapor transmittance of the layer k
c due to line absorption
h2oexp(i,j,k,1) = exp(-xh*xkw(ib))
enddo
enddo
enddo
do ik=2,6
if (mw(ib).eq.6) then
do k=1,np
do j=1,n
do i=1,m
xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
h2oexp(i,j,k,ik) = xh*xh*xh
enddo
enddo
enddo
elseif (mw(ib).eq.8) then
do k=1,np
do j=1,n
do i=1,m
xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
xh = xh*xh
h2oexp(i,j,k,ik) = xh*xh
enddo
enddo
enddo
elseif (mw(ib).eq.9) then
do k=1,np
do j=1,n
do i=1,m
xh=h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
xh1 = xh*xh
h2oexp(i,j,k,ik) = xh*xh1
enddo
enddo
enddo
else
do k=1,np
do j=1,n
do i=1,m
xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
xh = xh*xh
xh = xh*xh
h2oexp(i,j,k,ik) = xh*xh
enddo
enddo
enddo
endif
enddo
return
end
c**********************************************************************
subroutine CONEXPS(ib,m,n,np,dcont,xke,conexp)
c**********************************************************************
c compute exponentials for continuum absorption in individual layers.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c layer scaled water vapor amount for continuum absorption (dcont)
c absorption coefficients for the first k-distribution function
c due to water vapor continuum absorption (xke)
c
c---- output parameters
c 1 or 3 exponentials for each layer (conexp)
c**********************************************************************
implicit none
integer ib,m,n,np,i,j,k,iq
c---- input parameters ------
_RL dcont(m,n,np)
c---- updated parameters -----
_RL conexp(m,n,np,3)
c---- static data -----
_RL xke(9)
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
conexp(i,j,k,1) = exp(-dcont(i,j,k)*xke(ib))
enddo
enddo
enddo
if (ib .eq. 3) then
c-----the absorption coefficients for sub-bands 3b (iq=2) and 3a (iq=3)
c are, respectively, two and three times the absorption coefficient
c for sub-band 3c (iq=1) (table 6).
do iq=2,3
do k=1,np
do j=1,n
do i=1,m
conexp(i,j,k,iq) = conexp(i,j,k,iq-1) *conexp(i,j,k,iq-1)
enddo
enddo
enddo
enddo
endif
return
end
c**********************************************************************
subroutine CO2EXPS(m,n,np,dco2,pa,dt,co2exp)
c**********************************************************************
c compute co2 exponentials for individual layers.
c
c---- input parameters
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c layer co2 amount (dco2)
c layer pressure (pa)
c layer temperature minus 250K (dt)
c
c---- output parameters
c 6 exponentials for each layer (co2exp)
c**********************************************************************
implicit none
integer m,n,np,i,j,k
c---- input parameters -----
_RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)
c---- output parameters -----
_RL co2exp(m,n,np,6,2)
c---- temporary arrays -----
_RL xc
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
c-----compute the scaled co2 amount from eq. (27) for band-wings
c (sub-bands 3a and 3c).
xc = dco2(i,j,k)*(pa(i,j,k)/300.0)**0.5
1 *(1.+(0.0182+1.07e-4*dt(i,j,k))*dt(i,j,k))
c-----six exponentials by powers of 8 (table 7).
co2exp(i,j,k,1,1)=exp(-xc*2.656e-5)
xc=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)
xc=xc*xc
co2exp(i,j,k,2,1)=xc*xc
xc=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)
xc=xc*xc
co2exp(i,j,k,3,1)=xc*xc
xc=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)
xc=xc*xc
co2exp(i,j,k,4,1)=xc*xc
xc=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)
xc=xc*xc
co2exp(i,j,k,5,1)=xc*xc
xc=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)
xc=xc*xc
co2exp(i,j,k,6,1)=xc*xc
c-----compute the scaled co2 amount from eq. (27) for band-center
c region (sub-band 3b).
xc = dco2(i,j,k)*(pa(i,j,k)/30.0)**0.85
1 *(1.+(0.0042+2.00e-5*dt(i,j,k))*dt(i,j,k))
co2exp(i,j,k,1,2)=exp(-xc*2.656e-3)
xc=co2exp(i,j,k,1,2)*co2exp(i,j,k,1,2)
xc=xc*xc
co2exp(i,j,k,2,2)=xc*xc
xc=co2exp(i,j,k,2,2)*co2exp(i,j,k,2,2)
xc=xc*xc
co2exp(i,j,k,3,2)=xc*xc
xc=co2exp(i,j,k,3,2)*co2exp(i,j,k,3,2)
xc=xc*xc
co2exp(i,j,k,4,2)=xc*xc
xc=co2exp(i,j,k,4,2)*co2exp(i,j,k,4,2)
xc=xc*xc
co2exp(i,j,k,5,2)=xc*xc
xc=co2exp(i,j,k,5,2)*co2exp(i,j,k,5,2)
xc=xc*xc
co2exp(i,j,k,6,2)=xc*xc
enddo
enddo
enddo
return
end
c**********************************************************************
subroutine N2OEXPS(ib,m,n,np,dn2o,pa,dt,n2oexp)
c**********************************************************************
c compute n2o exponentials for individual layers.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c layer n2o amount (dn2o)
c layer pressure (pa)
c layer temperature minus 250K (dt)
c
c---- output parameters
c 2 or 4 exponentials for each layer (n2oexp)
c**********************************************************************
implicit none
integer ib,m,n,np,i,j,k
c---- input parameters -----
_RL dn2o(m,n,np),pa(m,n,np),dt(m,n,np)
c---- output parameters -----
_RL n2oexp(m,n,np,4)
c---- temporary arrays -----
_RL xc,xc1,xc2
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
c-----four exponential by powers of 21 for band 6
if (ib.eq.6) then
xc=dn2o(i,j,k)*(1.+(1.9297e-3+4.3750e-6*dt(i,j,k))*dt(i,j,k))
n2oexp(i,j,k,1)=exp(-xc*6.31582e-2)
xc=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
xc1=xc*xc
xc2=xc1*xc1
n2oexp(i,j,k,2)=xc*xc1*xc2
c-----four exponential by powers of 8 for band 7
else
xc=dn2o(i,j,k)*(pa(i,j,k)/500.0)**0.48
* *(1.+(1.3804e-3+7.4838e-6*dt(i,j,k))*dt(i,j,k))
n2oexp(i,j,k,1)=exp(-xc*5.35779e-2)
xc=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
xc=xc*xc
n2oexp(i,j,k,2)=xc*xc
xc=n2oexp(i,j,k,2)*n2oexp(i,j,k,2)
xc=xc*xc
n2oexp(i,j,k,3)=xc*xc
xc=n2oexp(i,j,k,3)*n2oexp(i,j,k,3)
xc=xc*xc
n2oexp(i,j,k,4)=xc*xc
endif
enddo
enddo
enddo
return
end
c**********************************************************************
subroutine CH4EXPS(ib,m,n,np,dch4,pa,dt,ch4exp)
c**********************************************************************
c compute ch4 exponentials for individual layers.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c layer ch4 amount (dch4)
c layer pressure (pa)
c layer temperature minus 250K (dt)
c
c---- output parameters
c 1 or 4 exponentials for each layer (ch4exp)
c**********************************************************************
implicit none
integer ib,m,n,np,i,j,k
c---- input parameters -----
_RL dch4(m,n,np),pa(m,n,np),dt(m,n,np)
c---- output parameters -----
_RL ch4exp(m,n,np,4)
c---- temporary arrays -----
_RL xc
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
c-----four exponentials for band 6
if (ib.eq.6) then
xc=dch4(i,j,k)*(1.+(1.7007e-2+1.5826e-4*dt(i,j,k))*dt(i,j,k))
ch4exp(i,j,k,1)=exp(-xc*5.80708e-3)
c-----four exponentials by powers of 12 for band 7
else
xc=dch4(i,j,k)*(pa(i,j,k)/500.0)**0.65
* *(1.+(5.9590e-4-2.2931e-6*dt(i,j,k))*dt(i,j,k))
ch4exp(i,j,k,1)=exp(-xc*6.29247e-2)
xc=ch4exp(i,j,k,1)*ch4exp(i,j,k,1)*ch4exp(i,j,k,1)
xc=xc*xc
ch4exp(i,j,k,2)=xc*xc
xc=ch4exp(i,j,k,2)*ch4exp(i,j,k,2)*ch4exp(i,j,k,2)
xc=xc*xc
ch4exp(i,j,k,3)=xc*xc
xc=ch4exp(i,j,k,3)*ch4exp(i,j,k,3)*ch4exp(i,j,k,3)
xc=xc*xc
ch4exp(i,j,k,4)=xc*xc
endif
enddo
enddo
enddo
return
end
c**********************************************************************
subroutine COMEXPS(ib,m,n,np,dcom,dt,comexp)
c**********************************************************************
c compute co2-minor exponentials for individual layers.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c layer co2 amount (dcom)
c layer temperature minus 250K (dt)
c
c---- output parameters
c 2 exponentials for each layer (comexp)
c**********************************************************************
implicit none
integer ib,m,n,np,i,j,k
c---- input parameters -----
_RL dcom(m,n,np),dt(m,n,np)
c---- output parameters -----
_RL comexp(m,n,np,2)
c---- temporary arrays -----
_RL xc,xc1,xc2
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
c-----two exponentials by powers of 60 for band 4
if (ib.eq.4) then
xc=dcom(i,j,k)*(1.+(3.5775e-2+4.0447e-4*dt(i,j,k))*dt(i,j,k))
comexp(i,j,k,1)=exp(-xc*2.18947e-5)
xc=comexp(i,j,k,1)*comexp(i,j,k,1)*comexp(i,j,k,1)
xc=xc*xc
xc1=xc*xc
xc=xc1*xc1
xc=xc*xc
comexp(i,j,k,2)=xc*xc1
c-----two exponentials by powers of 44 for band 5
else
xc=dcom(i,j,k)*(1.+(3.4268e-2+3.7401e-4*dt(i,j,k))*dt(i,j,k))
comexp(i,j,k,1)=exp(-xc*4.74570e-5)
xc=comexp(i,j,k,1)*comexp(i,j,k,1)
xc1=xc*xc
xc2=xc1*xc1
xc=xc2*xc2
xc=xc*xc
comexp(i,j,k,2)=xc1*xc2*xc
endif
enddo
enddo
enddo
return
end
c**********************************************************************
subroutine CFCEXPS(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,dcfc,dt,cfcexp)
c**********************************************************************
c compute cfc(-11, -12, -22) exponentials for individual layers.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c parameters for computing the scaled cfc amounts
c for temperature scaling (a1,b1,a2,b2)
c the absorption coefficients for the
c first k-distribution function due to cfcs (fk1,fk2)
c layer cfc amounts (dcfc)
c layer temperature minus 250K (dt)
c
c---- output parameters
c 1 exponential for each layer (cfcexp)
c**********************************************************************
implicit none
integer ib,m,n,np,i,j,k
c---- input parameters -----
_RL dcfc(m,n,np),dt(m,n,np)
c---- output parameters -----
_RL cfcexp(m,n,np)
c---- static data -----
_RL a1,b1,fk1,a2,b2,fk2
c---- temporary arrays -----
_RL xf
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
c-----compute the scaled cfc amount (xf) and exponential (cfcexp)
if (ib.eq.4) then
xf=dcfc(i,j,k)*(1.+(a1+b1*dt(i,j,k))*dt(i,j,k))
cfcexp(i,j,k)=exp(-xf*fk1)
else
xf=dcfc(i,j,k)*(1.+(a2+b2*dt(i,j,k))*dt(i,j,k))
cfcexp(i,j,k)=exp(-xf*fk2)
endif
enddo
enddo
enddo
return
end
c**********************************************************************
subroutine B10EXPS(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt
* ,h2oexp,conexp,co2exp,n2oexp)
c**********************************************************************
c compute band3a exponentials for individual layers.
c
c---- input parameters
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of layers (np)
c layer h2o amount for line absorption (dh2o)
c layer h2o amount for continuum absorption (dcont)
c layer co2 amount (dco2)
c layer n2o amount (dn2o)
c layer pressure (pa)
c layer temperature minus 250K (dt)
c
c---- output parameters
c
c exponentials for each layer (h2oexp,conexp,co2exp,n2oexp)
c**********************************************************************
implicit none
integer m,n,np,i,j,k
c---- input parameters -----
_RL dh2o(m,n,np),dcont(m,n,np),dn2o(m,n,np)
_RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)
c---- output parameters -----
_RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
* ,n2oexp(m,n,np,4)
c---- temporary arrays -----
_RL xx,xx1,xx2,xx3
c**********************************************************************
do k=1,np
do j=1,n
do i=1,m
c-----compute the scaled h2o-line amount for the sub-band 3a
c table 3, Chou et al. (JAS, 1993)
xx=dh2o(i,j,k)*(pa(i,j,k)/500.0)
1 *(1.+(0.0149+6.20e-5*dt(i,j,k))*dt(i,j,k))
c-----six exponentials by powers of 8
c the constant 0.10624 is equal to 1.66*0.064
h2oexp(i,j,k,1)=exp(-xx*0.10624)
xx=h2oexp(i,j,k,1)*h2oexp(i,j,k,1)
xx=xx*xx
h2oexp(i,j,k,2)=xx*xx
xx=h2oexp(i,j,k,2)*h2oexp(i,j,k,2)
xx=xx*xx
h2oexp(i,j,k,3)=xx*xx
xx=h2oexp(i,j,k,3)*h2oexp(i,j,k,3)
xx=xx*xx
h2oexp(i,j,k,4)=xx*xx
xx=h2oexp(i,j,k,4)*h2oexp(i,j,k,4)
xx=xx*xx
h2oexp(i,j,k,5)=xx*xx
xx=h2oexp(i,j,k,5)*h2oexp(i,j,k,5)
xx=xx*xx
h2oexp(i,j,k,6)=xx*xx
c-----compute the scaled co2 amount for the sub-band 3a
c table 1, Chou et al. (JAS, 1993)
xx=dco2(i,j,k)*(pa(i,j,k)/300.0)**0.5
1 *(1.+(0.0179+1.02e-4*dt(i,j,k))*dt(i,j,k))
c-----six exponentials by powers of 8
c the constant 2.656e-5 is equal to 1.66*1.60e-5
co2exp(i,j,k,1,1)=exp(-xx*2.656e-5)
xx=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)
xx=xx*xx
co2exp(i,j,k,2,1)=xx*xx
xx=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)
xx=xx*xx
co2exp(i,j,k,3,1)=xx*xx
xx=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)
xx=xx*xx
co2exp(i,j,k,4,1)=xx*xx
xx=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)
xx=xx*xx
co2exp(i,j,k,5,1)=xx*xx
xx=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)
xx=xx*xx
co2exp(i,j,k,6,1)=xx*xx
c-----one exponential of h2o continuum for sub-band 3a
c tabl 5 of Chou et. al. (JAS, 1993)
c the constant 1.04995e+2 is equal to 1.66*63.25
conexp(i,j,k,1)=exp(-dcont(i,j,k)*1.04995e+2)
c-----compute the scaled n2o amount for sub-band 3a
xx=dn2o(i,j,k)*(1.+(1.4476e-3+3.6656e-6*dt(i,j,k))*dt(i,j,k))
c-----two exponential2 by powers of 58
n2oexp(i,j,k,1)=exp(-xx*0.25238)
xx=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
xx1=xx*xx
xx1=xx1*xx1
xx2=xx1*xx1
xx3=xx2*xx2
n2oexp(i,j,k,2)=xx*xx1*xx2*xx3
enddo
enddo
enddo
return
end
c**********************************************************************
subroutine TABLUP(k1,k2,m,n,np,nx,nh,nt,sabs,spre,stem,w1,p1,
* dwe,dpe,coef1,coef2,coef3,tran)
c**********************************************************************
c compute water vapor, co2 and o3 transmittances between levels
c k1 and k2 for m x n soundings, using table look-up.
c
c calculations follow Eq. (40) of Chou and Suarez (1994)
c
c---- input ---------------------
c indices for pressure levels (k1 and k2)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of atmospheric layers (np)
c number of pressure intervals in the table (nx)
c number of absorber amount intervals in the table (nh)
c number of tables copied (nt)
c column-integrated absorber amount (sabs)
c column absorber amount-weighted pressure (spre)
c column absorber amount-weighted temperature (stem)
c first value of absorber amount (log10) in the table (w1)
c first value of pressure (log10) in the table (p1)
c size of the interval of absorber amount (log10) in the table (dwe)
c size of the interval of pressure (log10) in the table (dpe)
c pre-computed coefficients (coef1, coef2, and coef3)
c
c---- updated ---------------------
c transmittance (tran)
c
c Note:
c (1) units of sabs are g/cm**2 for water vapor and
c (cm-atm)stp for co2 and o3.
c (2) units of spre and stem are, respectively, mb and K.
c (3) there are nt identical copies of the tables (coef1, coef2, and
c coef3). the prupose of using the multiple copies of tables is
c to increase the speed in parallel (vectorized) computations.
C if such advantage does not exist, nt can be set to 1.
c
c**********************************************************************
implicit none
integer k1,k2,m,n,np,nx,nh,nt,i,j
c---- input parameters -----
_RL w1,p1,dwe,dpe
_RL sabs(m,n,np+1),spre(m,n,np+1),stem(m,n,np+1)
_RL coef1(nx,nh,nt),coef2(nx,nh,nt),coef3(nx,nh,nt)
c---- update parameter -----
_RL tran(m,n)
c---- temporary variables -----
_RL x1,x2,x3,we,pe,fw,fp,pa,pb,pc,ax,ba,bb,t1,ca,cb,t2
integer iw,ip,nn
c**********************************************************************
do j=1,n
do i=1,m
nn=mod(i,nt)+1
x1=sabs(i,j,k2)-sabs(i,j,k1)
x2=(spre(i,j,k2)-spre(i,j,k1))/x1
x3=(stem(i,j,k2)-stem(i,j,k1))/x1
we=(log10(x1)-w1)/dwe
pe=(log10(x2)-p1)/dpe
we=max(we,w1-2.*dwe)
pe=max(pe,p1)
iw=int(we+1.5)
ip=int(pe+1.5)
iw=min(iw,nh-1)
iw=max(iw, 2)
ip=min(ip,nx-1)
ip=max(ip, 1)
fw=we-float(iw-1)
fp=pe-float(ip-1)
c-----linear interpolation in pressure
pa = coef1(ip,iw-1,nn)*(1.-fp)+coef1(ip+1,iw-1,nn)*fp
pb = coef1(ip,iw, nn)*(1.-fp)+coef1(ip+1,iw, nn)*fp
pc = coef1(ip,iw+1,nn)*(1.-fp)+coef1(ip+1,iw+1,nn)*fp
c-----quadratic interpolation in absorber amount for coef1
ax = (-pa*(1.-fw)+pc*(1.+fw)) *fw*0.5 + pb*(1.-fw*fw)
c-----linear interpolation in absorber amount for coef2 and coef3
ba = coef2(ip,iw, nn)*(1.-fp)+coef2(ip+1,iw, nn)*fp
bb = coef2(ip,iw+1,nn)*(1.-fp)+coef2(ip+1,iw+1,nn)*fp
t1 = ba*(1.-fw) + bb*fw
ca = coef3(ip,iw, nn)*(1.-fp)+coef3(ip+1,iw, nn)*fp
cb = coef3(ip,iw+1,nn)*(1.-fp)+coef3(ip+1,iw+1,nn)*fp
t2 = ca*(1.-fw) + cb*fw
c-----update the total transmittance between levels k1 and k2
tran(i,j)= (ax + (t1+t2*x3) * x3)*tran(i,j)
enddo
enddo
return
end
c**********************************************************************
subroutine H2OKDIS(ib,m,n,np,k,fkw,gkw,ne,h2oexp,conexp,
* th2o,tcon,tran)
c**********************************************************************
c compute water vapor transmittance between levels k1 and k2 for
c m x n soundings, using the k-distribution method.
c
c computations follow eqs. (34), (46), (50) and (52).
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of levels (np)
c current level (k)
c planck-weighted k-distribution function due to
c h2o line absorption (fkw)
c planck-weighted k-distribution function due to
c h2o continuum absorption (gkw)
c number of terms used in each band to compute water vapor
c continuum transmittance (ne)
c exponentials for line absorption (h2oexp)
c exponentials for continuum absorption (conexp)
c
c---- updated parameters
c transmittance between levels k1 and k2 due to
c water vapor line absorption (th2o)
c transmittance between levels k1 and k2 due to
c water vapor continuum absorption (tcon)
c total transmittance (tran)
c
c**********************************************************************
implicit none
integer ib,m,n,np,k,i,j
c---- input parameters ------
_RL conexp(m,n,np,3),h2oexp(m,n,np,6)
integer ne(9)
_RL fkw(6,9),gkw(6,3)
c---- updated parameters -----
_RL th2o(m,n,6),tcon(m,n,3),tran(m,n)
c---- temporary arrays -----
_RL trnth2o
c-----tco2 are the six exp factors between levels k1 and k2
c tran is the updated total transmittance between levels k1 and k2
c-----th2o is the 6 exp factors between levels k1 and k2 due to
c h2o line absorption.
c-----tcon is the 3 exp factors between levels k1 and k2 due to
c h2o continuum absorption.
c-----trnth2o is the total transmittance between levels k1 and k2 due
c to both line and continuum absorption computed from eq. (52).
do j=1,n
do i=1,m
th2o(i,j,1) = th2o(i,j,1)*h2oexp(i,j,k,1)
th2o(i,j,2) = th2o(i,j,2)*h2oexp(i,j,k,2)
th2o(i,j,3) = th2o(i,j,3)*h2oexp(i,j,k,3)
th2o(i,j,4) = th2o(i,j,4)*h2oexp(i,j,k,4)
th2o(i,j,5) = th2o(i,j,5)*h2oexp(i,j,k,5)
th2o(i,j,6) = th2o(i,j,6)*h2oexp(i,j,k,6)
enddo
enddo
if (ne(ib).eq.0) then
do j=1,n
do i=1,m
trnth2o =(fkw(1,ib)*th2o(i,j,1)
* + fkw(2,ib)*th2o(i,j,2)
* + fkw(3,ib)*th2o(i,j,3)
* + fkw(4,ib)*th2o(i,j,4)
* + fkw(5,ib)*th2o(i,j,5)
* + fkw(6,ib)*th2o(i,j,6))
tran(i,j)=tran(i,j)*trnth2o
enddo
enddo
elseif (ne(ib).eq.1) then
do j=1,n
do i=1,m
tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)
trnth2o =(fkw(1,ib)*th2o(i,j,1)
* + fkw(2,ib)*th2o(i,j,2)
* + fkw(3,ib)*th2o(i,j,3)
* + fkw(4,ib)*th2o(i,j,4)
* + fkw(5,ib)*th2o(i,j,5)
* + fkw(6,ib)*th2o(i,j,6))*tcon(i,j,1)
tran(i,j)=tran(i,j)*trnth2o
enddo
enddo
else
do j=1,n
do i=1,m
tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)
tcon(i,j,2)= tcon(i,j,2)*conexp(i,j,k,2)
tcon(i,j,3)= tcon(i,j,3)*conexp(i,j,k,3)
trnth2o = ( gkw(1,1)*th2o(i,j,1)
* + gkw(2,1)*th2o(i,j,2)
* + gkw(3,1)*th2o(i,j,3)
* + gkw(4,1)*th2o(i,j,4)
* + gkw(5,1)*th2o(i,j,5)
* + gkw(6,1)*th2o(i,j,6) ) * tcon(i,j,1)
* + ( gkw(1,2)*th2o(i,j,1)
* + gkw(2,2)*th2o(i,j,2)
* + gkw(3,2)*th2o(i,j,3)
* + gkw(4,2)*th2o(i,j,4)
* + gkw(5,2)*th2o(i,j,5)
* + gkw(6,2)*th2o(i,j,6) ) * tcon(i,j,2)
* + ( gkw(1,3)*th2o(i,j,1)
* + gkw(2,3)*th2o(i,j,2)
* + gkw(3,3)*th2o(i,j,3)
* + gkw(4,3)*th2o(i,j,4)
* + gkw(5,3)*th2o(i,j,5)
* + gkw(6,3)*th2o(i,j,6) ) * tcon(i,j,3)
tran(i,j)=tran(i,j)*trnth2o
enddo
enddo
endif
return
end
c**********************************************************************
subroutine CO2KDIS(m,n,np,k,co2exp,tco2,tran)
c**********************************************************************
c compute co2 transmittances between levels k1 and k2 for
c m x n soundings, using the k-distribution method with linear
c pressure scaling. computations follow eq. (34).
c
c---- input parameters
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of levels (np)
c current level (k)
c exponentials for co2 absorption (co2exp)
c
c---- updated parameters
c transmittance between levels k1 and k2 due to co2 absorption
c for the various values of the absorption coefficient (tco2)
c total transmittance (tran)
c
c**********************************************************************
implicit none
integer m,n,np,k,i,j
c---- input parameters -----
_RL co2exp(m,n,np,6,2)
c---- updated parameters -----
_RL tco2(m,n,6,2),tran(m,n)
c---- temporary arrays -----
_RL xc
c-----tco2 is the 6 exp factors between levels k1 and k2.
c xc is the total co2 transmittance given by eq. (53).
do j=1,n
do i=1,m
c-----band-wings
tco2(i,j,1,1)=tco2(i,j,1,1)*co2exp(i,j,k,1,1)
xc= 0.1395 *tco2(i,j,1,1)
tco2(i,j,2,1)=tco2(i,j,2,1)*co2exp(i,j,k,2,1)
xc=xc+0.1407 *tco2(i,j,2,1)
tco2(i,j,3,1)=tco2(i,j,3,1)*co2exp(i,j,k,3,1)
xc=xc+0.1549 *tco2(i,j,3,1)
tco2(i,j,4,1)=tco2(i,j,4,1)*co2exp(i,j,k,4,1)
xc=xc+0.1357 *tco2(i,j,4,1)
tco2(i,j,5,1)=tco2(i,j,5,1)*co2exp(i,j,k,5,1)
xc=xc+0.0182 *tco2(i,j,5,1)
tco2(i,j,6,1)=tco2(i,j,6,1)*co2exp(i,j,k,6,1)
xc=xc+0.0220 *tco2(i,j,6,1)
c-----band-center region
tco2(i,j,1,2)=tco2(i,j,1,2)*co2exp(i,j,k,1,2)
xc=xc+0.0766 *tco2(i,j,1,2)
tco2(i,j,2,2)=tco2(i,j,2,2)*co2exp(i,j,k,2,2)
xc=xc+0.1372 *tco2(i,j,2,2)
tco2(i,j,3,2)=tco2(i,j,3,2)*co2exp(i,j,k,3,2)
xc=xc+0.1189 *tco2(i,j,3,2)
tco2(i,j,4,2)=tco2(i,j,4,2)*co2exp(i,j,k,4,2)
xc=xc+0.0335 *tco2(i,j,4,2)
tco2(i,j,5,2)=tco2(i,j,5,2)*co2exp(i,j,k,5,2)
xc=xc+0.0169 *tco2(i,j,5,2)
tco2(i,j,6,2)=tco2(i,j,6,2)*co2exp(i,j,k,6,2)
xc=xc+0.0059 *tco2(i,j,6,2)
tran(i,j)=tran(i,j)*xc
enddo
enddo
return
end
c**********************************************************************
subroutine N2OKDIS(ib,m,n,np,k,n2oexp,tn2o,tran)
c**********************************************************************
c compute n2o transmittances between levels k1 and k2 for
c m x n soundings, using the k-distribution method with linear
c pressure scaling.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of levels (np)
c current level (k)
c exponentials for n2o absorption (n2oexp)
c
c---- updated parameters
c transmittance between levels k1 and k2 due to n2o absorption
c for the various values of the absorption coefficient (tn2o)
c total transmittance (tran)
c
c**********************************************************************
implicit none
integer ib,m,n,np,k,i,j
c---- input parameters -----
_RL n2oexp(m,n,np,4)
c---- updated parameters -----
_RL tn2o(m,n,4),tran(m,n)
c---- temporary arrays -----
_RL xc
c-----tn2o is the 2 exp factors between levels k1 and k2.
c xc is the total n2o transmittance
do j=1,n
do i=1,m
c-----band 6
if (ib.eq.6) then
tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
xc= 0.940414*tn2o(i,j,1)
tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
xc=xc+0.059586*tn2o(i,j,2)
c-----band 7
else
tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
xc= 0.561961*tn2o(i,j,1)
tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
xc=xc+0.138707*tn2o(i,j,2)
tn2o(i,j,3)=tn2o(i,j,3)*n2oexp(i,j,k,3)
xc=xc+0.240670*tn2o(i,j,3)
tn2o(i,j,4)=tn2o(i,j,4)*n2oexp(i,j,k,4)
xc=xc+0.058662*tn2o(i,j,4)
endif
tran(i,j)=tran(i,j)*xc
enddo
enddo
return
end
c**********************************************************************
subroutine CH4KDIS(ib,m,n,np,k,ch4exp,tch4,tran)
c**********************************************************************
c compute ch4 transmittances between levels k1 and k2 for
c m x n soundings, using the k-distribution method with
c linear pressure scaling.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of levels (np)
c current level (k)
c exponentials for ch4 absorption (ch4exp)
c
c---- updated parameters
c transmittance between levels k1 and k2 due to ch4 absorption
c for the various values of the absorption coefficient (tch4)
c total transmittance (tran)
c
c**********************************************************************
implicit none
integer ib,m,n,np,k,i,j
c---- input parameters -----
_RL ch4exp(m,n,np,4)
c---- updated parameters -----
_RL tch4(m,n,4),tran(m,n)
c---- temporary arrays -----
_RL xc
c-----tch4 is the 2 exp factors between levels k1 and k2.
c xc is the total ch4 transmittance
do j=1,n
do i=1,m
c-----band 6
if (ib.eq.6) then
tch4(i,j,1)=tch4(i,j,1)*ch4exp(i,j,k,1)
xc= tch4(i,j,1)
c-----band 7
else
tch4(i,j,1)=tch4(i,j,1)*ch4exp(i,j,k,1)
xc= 0.610650*tch4(i,j,1)
tch4(i,j,2)=tch4(i,j,2)*ch4exp(i,j,k,2)
xc=xc+0.280212*tch4(i,j,2)
tch4(i,j,3)=tch4(i,j,3)*ch4exp(i,j,k,3)
xc=xc+0.107349*tch4(i,j,3)
tch4(i,j,4)=tch4(i,j,4)*ch4exp(i,j,k,4)
xc=xc+0.001789*tch4(i,j,4)
endif
tran(i,j)=tran(i,j)*xc
enddo
enddo
return
end
c**********************************************************************
subroutine COMKDIS(ib,m,n,np,k,comexp,tcom,tran)
c**********************************************************************
c compute co2-minor transmittances between levels k1 and k2
c for m x n soundings, using the k-distribution method
c with linear pressure scaling.
c
c---- input parameters
c spectral band (ib)
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of levels (np)
c current level (k)
c exponentials for co2-minor absorption (comexp)
c
c---- updated parameters
c transmittance between levels k1 and k2 due to co2-minor absorption
c for the various values of the absorption coefficient (tcom)
c total transmittance (tran)
c
c**********************************************************************
implicit none
integer ib,m,n,np,k,i,j
c---- input parameters -----
_RL comexp(m,n,np,2)
c---- updated parameters -----
_RL tcom(m,n,2),tran(m,n)
c---- temporary arrays -----
_RL xc
c-----tcom is the 2 exp factors between levels k1 and k2.
c xc is the total co2-minor transmittance
do j=1,n
do i=1,m
c-----band 4
if (ib.eq.4) then
tcom(i,j,1)=tcom(i,j,1)*comexp(i,j,k,1)
xc= 0.972025*tcom(i,j,1)
tcom(i,j,2)=tcom(i,j,2)*comexp(i,j,k,2)
xc=xc+0.027975*tcom(i,j,2)
c-----band 5
else
tcom(i,j,1)=tcom(i,j,1)*comexp(i,j,k,1)
xc= 0.961324*tcom(i,j,1)
tcom(i,j,2)=tcom(i,j,2)*comexp(i,j,k,2)
xc=xc+0.038676*tcom(i,j,2)
endif
tran(i,j)=tran(i,j)*xc
enddo
enddo
return
end
c**********************************************************************
subroutine CFCKDIS(m,n,np,k,cfcexp,tcfc,tran)
c**********************************************************************
c compute cfc-(11,12,22) transmittances between levels k1 and k2
c for m x n soundings, using the k-distribution method with
c linear pressure scaling.
c
c---- input parameters
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of levels (np)
c current level (k)
c exponentials for cfc absorption (cfcexp)
c
c---- updated parameters
c transmittance between levels k1 and k2 due to cfc absorption
c for the various values of the absorption coefficient (tcfc)
c total transmittance (tran)
c
c**********************************************************************
implicit none
integer m,n,np,k,i,j
c---- input parameters -----
_RL cfcexp(m,n,np)
c---- updated parameters -----
_RL tcfc(m,n),tran(m,n)
c-----tcfc is the exp factors between levels k1 and k2.
do j=1,n
do i=1,m
tcfc(i,j)=tcfc(i,j)*cfcexp(i,j,k)
tran(i,j)=tran(i,j)*tcfc(i,j)
enddo
enddo
return
end
c**********************************************************************
subroutine B10KDIS(m,n,np,k,h2oexp,conexp,co2exp,n2oexp
* ,th2o,tcon,tco2,tn2o,tran)
c**********************************************************************
c
c compute h2o (line and continuum),co2,n2o transmittances between
c levels k1 and k2 for m x n soundings, using the k-distribution
c method with linear pressure scaling.
c
c---- input parameters
c number of grid intervals in zonal direction (m)
c number of grid intervals in meridional direction (n)
c number of levels (np)
c current level (k)
c exponentials for h2o line absorption (h2oexp)
c exponentials for h2o continuum absorption (conexp)
c exponentials for co2 absorption (co2exp)
c exponentials for n2o absorption (n2oexp)
c
c---- updated parameters
c transmittance between levels k1 and k2 due to h2o line absorption
c for the various values of the absorption coefficient (th2o)
c transmittance between levels k1 and k2 due to h2o continuum
c absorption for the various values of the absorption
c coefficient (tcon)
c transmittance between levels k1 and k2 due to co2 absorption
c for the various values of the absorption coefficient (tco2)
c transmittance between levels k1 and k2 due to n2o absorption
c for the various values of the absorption coefficient (tn2o)
c total transmittance (tran)
c
c**********************************************************************
implicit none
integer m,n,np,k,i,j
c---- input parameters -----
_RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
* ,n2oexp(m,n,np,4)
c---- updated parameters -----
_RL th2o(m,n,6),tcon(m,n,3),tco2(m,n,6,2),tn2o(m,n,4)
* ,tran(m,n)
c---- temporary arrays -----
_RL xx
c-----initialize tran
do j=1,n
do i=1,m
tran(i,j)=1.0
enddo
enddo
c-----for h2o line
do j=1,n
do i=1,m
th2o(i,j,1)=th2o(i,j,1)*h2oexp(i,j,k,1)
xx= 0.3153*th2o(i,j,1)
th2o(i,j,2)=th2o(i,j,2)*h2oexp(i,j,k,2)
xx=xx+0.4604*th2o(i,j,2)
th2o(i,j,3)=th2o(i,j,3)*h2oexp(i,j,k,3)
xx=xx+0.1326*th2o(i,j,3)
th2o(i,j,4)=th2o(i,j,4)*h2oexp(i,j,k,4)
xx=xx+0.0798*th2o(i,j,4)
th2o(i,j,5)=th2o(i,j,5)*h2oexp(i,j,k,5)
xx=xx+0.0119*th2o(i,j,5)
tran(i,j)=tran(i,j)*xx
enddo
enddo
c-----for h2o continuum
do j=1,n
do i=1,m
tcon(i,j,1)=tcon(i,j,1)*conexp(i,j,k,1)
tran(i,j)=tran(i,j)*tcon(i,j,1)
enddo
enddo
c-----for co2
do j=1,n
do i=1,m
tco2(i,j,1,1)=tco2(i,j,1,1)*co2exp(i,j,k,1,1)
xx= 0.2673*tco2(i,j,1,1)
tco2(i,j,2,1)=tco2(i,j,2,1)*co2exp(i,j,k,2,1)
xx=xx+ 0.2201*tco2(i,j,2,1)
tco2(i,j,3,1)=tco2(i,j,3,1)*co2exp(i,j,k,3,1)
xx=xx+ 0.2106*tco2(i,j,3,1)
tco2(i,j,4,1)=tco2(i,j,4,1)*co2exp(i,j,k,4,1)
xx=xx+ 0.2409*tco2(i,j,4,1)
tco2(i,j,5,1)=tco2(i,j,5,1)*co2exp(i,j,k,5,1)
xx=xx+ 0.0196*tco2(i,j,5,1)
tco2(i,j,6,1)=tco2(i,j,6,1)*co2exp(i,j,k,6,1)
xx=xx+ 0.0415*tco2(i,j,6,1)
tran(i,j)=tran(i,j)*xx
enddo
enddo
c-----for n2o
do j=1,n
do i=1,m
tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
xx= 0.970831*tn2o(i,j,1)
tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
xx=xx+0.029169*tn2o(i,j,2)
tran(i,j)=tran(i,j)*(xx-1.0)
enddo
enddo
return
end