Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:40:18 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
8b150b4af9 Andr*0002       subroutine lwrio (nymd,nhms,bi,bj,myid,istrip,npcs,
                0003      .                  low_level,mid_level,
082e38725b Andr*0004      .                  im,jm,lm,
5cede9ba3f Andr*0005      .                  pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,
                0006      .                  co2,cfc11,cfc12,cfc22,methane,n2o,emissivity,
1662f365b2 Andr*0007      .                  tgz,radlwg,st4,dst4,
                0008      .                  dtradlw,dlwdtg,dtradlwc,lwgclr,
032fb71841 Andr*0009      .                  nlwcld,cldlw,clwmo,nlwlz,lwlz,
7d59649d1b Andr*0010      .                  lpnt,imstturb,qliqave,fccave,landtype)
1662f365b2 Andr*0011 
                0012       implicit none
                0013 
                0014 c Input Variables
                0015 c ---------------
8b150b4af9 Andr*0016       integer nymd,nhms,istrip,npcs,bi,bj,myid
d2aaec7e02 Andr*0017       integer mid_level,low_level
1662f365b2 Andr*0018       integer im,jm,lm        
a456aa407c Andr*0019       _RL pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1)
                0020       _RL dpres(im,jm,lm),pkht(im,jm,lm+1),pkz(im,jm,lm)
                0021       _RL tz(im,jm,lm),qz(im,jm,lm),oz(im,jm,lm)
                0022       _RL co2,cfc11,cfc12,cfc22,methane(lm),n2o(lm)    
                0023       _RL emissivity(im,jm,10)
                0024       _RL tgz(im,jm),radlwg(im,jm),st4(im,jm),dst4(im,jm)     
                0025       _RL dtradlw(im,jm,lm),dlwdtg (im,jm,lm) 
                0026       _RL dtradlwc(im,jm,lm),lwgclr(im,jm)     
5cede9ba3f Andr*0027       integer nlwcld,nlwlz
a456aa407c Andr*0028       _RL cldlw(im,jm,lm),clwmo(im,jm,lm),lwlz(im,jm,lm)
5cede9ba3f Andr*0029       logical lpnt
                0030       integer imstturb
a456aa407c Andr*0031       _RL qliqave(im,jm,lm),fccave(im,jm,lm)
5cede9ba3f Andr*0032       integer landtype(im,jm)
1662f365b2 Andr*0033 
                0034 c Local Variables
                0035 c ---------------
                0036       integer i,j,l,n,nn
                0037 
a456aa407c Andr*0038       _RL cldtot (im,jm,lm)
                0039       _RL cldmxo (im,jm,lm)
                0040 
                0041       _RL pl(istrip,lm)
                0042       _RL pk(istrip,lm)
                0043       _RL pke(istrip,lm)
                0044       _RL ple(istrip,lm+1)
                0045 
                0046       _RL ADELPL(ISTRIP,lm)
                0047       _RL dtrad(istrip,lm),dtradc(istrip,lm)
                0048       _RL OZL(ISTRIP,lm),TZL(ISTRIP,lm)
                0049       _RL SHZL(ISTRIP,lm),CLRO(ISTRIP,lm)
                0050       _RL CLMO(ISTRIP,lm)
                0051       _RL flx(ISTRIP,lm+1),flxclr(ISTRIP,lm+1)
                0052       _RL cldlz(istrip,lm)
                0053       _RL dfdts(istrip,lm+1),dtdtg(istrip,lm)
                0054 
                0055       _RL emiss(istrip,10)
                0056       _RL taual(istrip,lm,10)
                0057       _RL ssaal(istrip,lm,10)
                0058       _RL asyal(istrip,lm,10)
                0059       _RL cwc(istrip,lm,3)
                0060       _RL reff(istrip,lm,3)
                0061       _RL tauc(istrip,lm,3)
                0062 
                0063       _RL SGMT4(ISTRIP),TSURF(ISTRIP),dsgmt4(ISTRIP)
5cede9ba3f Andr*0064       integer lwi(istrip)
                0065 
032fb71841 Andr*0066       _RL tmpstrip(istrip,lm)
                0067       _RL tmpimjm(im,jm,lm)
f1eddbb0f1 Andr*0068       _RL tempor1(im,jm),tempor2(im,jm)
032fb71841 Andr*0069 
a456aa407c Andr*0070       _RL getcon,secday,convrt
8b150b4af9 Andr*0071 #ifdef ALLOW_DIAGNOSTICS
                0072       logical  diagnostics_is_on
                0073       external diagnostics_is_on
                0074       _RL tmpdiag(im,jm)
                0075 #endif
5cede9ba3f Andr*0076 
                0077       logical high,  trace, cldwater
032fb71841 Andr*0078 c     data high /.true./
                0079 c     data trace /.true./
                0080       data high /.false./
                0081       data trace /.false./
1662f365b2 Andr*0082       data cldwater /.false./
                0083 
                0084 C **********************************************************************
                0085 C ****                     INITIALIZATION                           ****
                0086 C **********************************************************************
                0087  
                0088       SECDAY = GETCON('SDAY')
                0089       CONVRT = GETCON('GRAVITY') / ( 100.0 * GETCON('CP') )
                0090 
                0091 c Adjust cloud fractions and cloud liquid water due to moist turbulence
                0092 c ---------------------------------------------------------------------
                0093       if(imstturb.ne.0) then
                0094         do L =1,lm
                0095         do j =1,jm
                0096         do i =1,im
e7070f537c Cons*0097            cldtot(i,j,L)=min(1.0 _d 0,max(cldlw(i,j,L),fccave(i,j,L)/
                0098      $          imstturb))
                0099            cldmxo(i,j,L) =  min( 1.0 _d 0,      clwmo(i,j,L) )
1662f365b2 Andr*0100            lwlz(i,j,L) =  lwlz(i,j,L) + qliqave(i,j,L)/imstturb
                0101         enddo
                0102         enddo
                0103         enddo
                0104       else
                0105         do L =1,lm
                0106         do j =1,jm
                0107         do i =1,im
32362b8fa8 Cons*0108          cldtot(i,j,L) =  min( 1.0 _d 0,cldlw(i,j,L) )
                0109          cldmxo(i,j,L) =  min( 1.0 _d 0,clwmo(i,j,L) )
1662f365b2 Andr*0110         enddo
                0111         enddo
                0112         enddo
                0113       endif
                0114 
                0115 C***********************************************************************
                0116 C ****                     LOOP OVER STRIPS                         ****
                0117 C **********************************************************************
                0118 
                0119       DO 1000 NN=1,NPCS
                0120 
                0121 C **********************************************************************
                0122 C ****                   VARIABLE INITIALIZATION                    ****
                0123 C **********************************************************************
                0124 
                0125        CALL STRIP (OZ,OZL      ,im*jm,ISTRIP,lm   ,NN)
                0126        CALL STRIP (PLZE,ple    ,im*jm,ISTRIP,lm+1 ,NN)
                0127        CALL STRIP (PLZ ,pl     ,im*jm,ISTRIP,lm   ,NN)
                0128        CALL STRIP (PKZ ,pk     ,im*jm,ISTRIP,lm   ,NN)
                0129        CALL STRIP (PKHT,pke    ,im*jm,ISTRIP,lm   ,NN)
                0130        CALL STRIP (TZ,TZL      ,im*jm,ISTRIP,lm   ,NN)
                0131        CALL STRIP (qz,SHZL     ,im*jm,ISTRIP,lm   ,NN)
                0132        CALL STRIP (cldtot,CLRO ,im*jm,ISTRIP,lm   ,NN)
                0133        CALL STRIP (cldmxo,CLMO ,im*jm,ISTRIP,lm   ,NN)
                0134        CALL STRIP (lwlz,cldlz  ,im*jm,ISTRIP,lm   ,NN)
                0135        CALL STRIP (tgz,tsurf   ,im*jm,ISTRIP,1    ,NN)
                0136 
                0137        CALL STRIP (emissivity,emiss,im*jm,ISTRIP,10,NN)
                0138 
                0139        call stripitint (landtype,lwi,im*jm,im*jm,istrip,1,nn)
                0140 
032fb71841 Andr*0141        DO L = 1,lm
                0142        DO I = 1,ISTRIP
                0143         ADELPL(I,L) = convrt / ( ple(I,L+1)-ple(I,L) )
                0144        ENDDO
1662f365b2 Andr*0145        ENDDO
                0146  
                0147 C Compute Clouds
                0148 C --------------
                0149        IF(NLWCLD.NE.0)THEN
                0150        DO L = 1,lm
                0151        DO I = 1,ISTRIP
32362b8fa8 Cons*0152         CLRO(I,L) = min( 1.0 _d 0,clro(i,L) )
                0153         CLMO(I,L) = min( 1.0 _d 0,clmo(i,L) )
1662f365b2 Andr*0154        ENDDO
                0155        ENDDO
                0156        ENDIF
                0157  
124bb68fab Andr*0158 C Convert to Temperature from Fizhi Theta
                0159 C ---------------------------------------
1662f365b2 Andr*0160       DO L = 1,lm
                0161       DO I = 1,ISTRIP
                0162       TZL(I,L) = TZL(I,L)*pk(I,L)
                0163       ENDDO
                0164       ENDDO
                0165       DO I = 1,ISTRIP
                0166 C To reduce longwave heating in lowest model layer
                0167       TZL(I,lm) = ( 2*tzl(i,lm)+tsurf(i) )/3.0  
                0168       ENDDO
                0169 
                0170 C Compute Optical Thicknesses
                0171 C ---------------------------
                0172       call opthk ( tzl,pl,ple,cldlz,clro,clmo,lwi,istrip,1,lm,tauc )
                0173  
                0174       do n = 1,3
                0175       do L = 1,lm
                0176       do i = 1,istrip
                0177       tauc(i,L,n) = tauc(i,L,n)*0.75
                0178       enddo
                0179       enddo
                0180       enddo
                0181 
                0182 C Set the optical thickness, single scattering albedo and assymetry factor
                0183 C     for aerosols equal to zero to omit the contribution of aerosols
                0184 c-------------------------------------------------------------------------
                0185       do n = 1,10
                0186       do L = 1,lm
                0187       do i = 1,istrip
                0188       taual(i,L,n) = 0.
                0189       ssaal(i,L,n) = 0.
                0190       asyal(i,L,n) = 0.
                0191       enddo
                0192       enddo
                0193       enddo
                0194 
                0195 C Set cwc and reff to zero - when cldwater is .false.
                0196 c----------------------------------------------------
                0197       do n = 1,3
                0198       do L = 1,lm
                0199       do i = 1,istrip
                0200        cwc(i,L,n) = 0.
                0201       reff(i,L,n) = 0.
                0202       enddo
                0203       enddo
                0204       enddo
                0205 
                0206 C **********************************************************************
                0207 C ****                CALL M-D Chou RADIATION                       ****
                0208 C **********************************************************************
                0209  
                0210       call irrad ( istrip,1,1,lm,ple,tzl,shzl,ozl,tsurf,co2,
                0211      .             n2o,methane,cfc11,cfc12,cfc22,emiss,
                0212      .             cldwater,cwc,tauc,reff,clro,mid_level,low_level,
                0213      .             taual,ssaal,asyal,
                0214      .             high,trace,flx,flxclr,dfdts,sgmt4 )
                0215  
                0216 C **********************************************************************
                0217 C ****                   HEATING RATES                              ****
                0218 C **********************************************************************
                0219 
                0220        do L = 1,lm
                0221        do i = 1,istrip
                0222          dtrad(i,L) = (   flx(i,L)-   flx(i,L+1))*adelpl(i,L)
154b1b9078 Andr*0223          tmpstrip(i,L) = flx(i,L)
1662f365b2 Andr*0224          dtdtg(i,L) = ( dfdts(i,L)- dfdts(i,L+1))*adelpl(i,L)
                0225         dtradc(i,L) = (flxclr(i,L)-flxclr(i,L+1))*adelpl(i,L)
                0226        enddo
                0227        enddo
                0228 
                0229 C **********************************************************************
                0230 C ****              NET UPWARD LONGWAVE FLUX  (W/M**2)              ****
                0231 C **********************************************************************
                0232 
                0233        DO I = 1,ISTRIP
                0234         flx   (i,1)    = -flx   (i,1)
                0235         flxclr(i,1)    = -flxclr(i,1)
                0236         flx   (i,lm+1) = -flx   (i,lm+1)
                0237         flxclr(i,lm+1) = -flxclr(i,lm+1)
                0238               sgmt4(i) = - sgmt4(i)
                0239              dsgmt4(i) = - dfdts(i,lm+1)
                0240        ENDDO
                0241 
                0242        CALL PASTE ( flx   (1,lm+1),RADLWG,ISTRIP,im*jm,1,NN )
                0243        CALL PASTE ( flxclr(1,lm+1),LWGCLR,ISTRIP,im*jm,1,NN )
                0244 
                0245        CALL PASTE (  sgmt4, st4,ISTRIP,im*jm,1,NN )
                0246        CALL PASTE ( dsgmt4,dst4,ISTRIP,im*jm,1,NN )
                0247 
                0248 C **********************************************************************
                0249 C ****            PASTE AND BUMP SOME DIAGNOSTICS                   ****
                0250 C **********************************************************************
                0251 
f1eddbb0f1 Andr*0252       CALL PASTE(flx(1,1),tempor1,ISTRIP,im*jm,1,NN)
                0253       CALL PASTE(flxclr(1,1),tempor2,ISTRIP,im*jm,1,NN)
1662f365b2 Andr*0254 
                0255 C **********************************************************************
                0256 C ****                 TENDENCY UPDATES                             ****
                0257 C **********************************************************************
                0258 
                0259       DO L = 1,lm
                0260       DO I = 1,ISTRIP
032fb71841 Andr*0261       DTRAD (I,L) = ple(i,lm+1) * DTRAD (I,L)/pk(I,L)
                0262       DTRADC(I,L) = ple(i,lm+1) * DTRADC(I,L)/pk(I,L)
                0263        dtdtg(I,L) = ple(i,lm+1) * dtdtg (I,L)/pk(I,L)
1662f365b2 Andr*0264       ENDDO
                0265       ENDDO
032fb71841 Andr*0266       CALL PASTE ( tmpstrip ,tmpimjm ,ISTRIP,im*jm,lm,NN )
1662f365b2 Andr*0267       CALL PASTE ( DTRAD ,DTRADLW ,ISTRIP,im*jm,lm,NN )
                0268       CALL PASTE ( DTRADC,DTRADLWC,ISTRIP,im*jm,lm,NN )
                0269       CALL PASTE ( dtdtg ,dlwdtg  ,ISTRIP,im*jm,lm,NN )
                0270 
                0271  1000 CONTINUE
                0272 
                0273 C **********************************************************************
                0274 C ****                     BUMP DIAGNOSTICS                         ****
                0275 C **********************************************************************
                0276 
8b150b4af9 Andr*0277 #ifdef ALLOW_DIAGNOSTICS
                0278       if(diagnostics_is_on('TGRLW   ',myid) ) then
                0279        call diagnostics_fill(tgz,'TGRLW   ',0,1,3,bi,bj,myid)
1662f365b2 Andr*0280       endif
                0281 
                0282       do L = 1,lm
                0283 
8b150b4af9 Andr*0284        if(diagnostics_is_on('TLW     ',myid) ) then
                0285         do j = 1,jm
                0286         do i = 1,im
                0287          tmpdiag(i,j) = tz(i,j,L)*pkz(i,j,L)
                0288         enddo
                0289         enddo
                0290         call diagnostics_fill(tmpdiag,'TLW     ',L,1,3,bi,bj,myid)
                0291        endif
1662f365b2 Andr*0292 
8b150b4af9 Andr*0293        if(diagnostics_is_on('SHRAD   ',myid) ) then
                0294         do j = 1,jm
                0295         do i = 1,im
                0296          tmpdiag(i,j) = qz(i,j,L)*1000.
                0297         enddo
                0298         enddo
                0299         call diagnostics_fill(tmpdiag,'SHRAD   ',L,1,3,bi,bj,myid)
                0300        endif
                0301 
                0302        if(diagnostics_is_on('OZLW    ',myid) ) then
                0303         do j = 1,jm
                0304         do i = 1,im
                0305          tmpdiag(i,j) = oz(i,j,L)
                0306         enddo
                0307         enddo
                0308         call diagnostics_fill(tmpdiag,'OZLW    ',L,1,3,bi,bj,myid)
                0309        endif
032fb71841 Andr*0310 
f1eddbb0f1 Andr*0311       enddo
8b150b4af9 Andr*0312 
                0313       if(diagnostics_is_on('OLR     ',myid) ) then
                0314        call diagnostics_fill(tempor1,'OLR     ',0,1,3,bi,bj,myid)
f1eddbb0f1 Andr*0315       endif
                0316 
8b150b4af9 Andr*0317       if(diagnostics_is_on('OLRCLR  ',myid) ) then
                0318        call diagnostics_fill(tempor2,'OLRCLR  ',0,1,3,bi,bj,myid)
032fb71841 Andr*0319       endif
8b150b4af9 Andr*0320 #endif
1662f365b2 Andr*0321 
                0322 C **********************************************************************
                0323 C ****    Increment Diagnostics Counters and Zero-Out Cloud Info    ****
                0324 C **********************************************************************
                0325 
                0326       nlwlz    = 0
                0327       nlwcld   = 0
                0328       imstturb = 0
                0329 
                0330       do L = 1,lm
                0331       do j = 1,jm
                0332       do i = 1,im
                0333          fccave(i,j,L) = 0.
                0334         qliqave(i,j,L) = 0.
                0335       enddo
                0336       enddo
                0337       enddo
                0338 
                0339       return
                0340       end
                0341 c**********************   November 26, 1997   **************************
                0342 
                0343       subroutine irrad (m,n,ndim,np,pl,ta,wa,oa,ts,co2,
                0344      *                  n2o,ch4,cfc11,cfc12,cfc22,emiss,
                0345      *                  cldwater,cwc,taucl,reff,fcld,ict,icb,
                0346      *                  taual,ssaal,asyal,
                0347      *                  high,trace,flx,flc,dfdts,st4)
                0348 
                0349 c***********************************************************************
                0350 c
                0351 c                        Version IR-5 (September, 1998)
                0352 c
                0353 c  New features of this version are:
                0354 c   (1) The effect of aerosol scattering on LW fluxes is included.
                0355 c   (2) Absorption and scattering of rain drops are included.
                0356 c
                0357 c***********************************************************************
                0358 c
                0359 c                        Version IR-4 (October, 1997)
                0360 c
                0361 c  New features of this version are:
                0362 c   (1) The surface is treated as non-black.  The surface
                0363 c         emissivity, emiss, is an input parameter
                0364 c   (2) The effect of cloud scattering on LW fluxes is included
                0365 c
                0366 c*********************************************************************
                0367 c
                0368 c This routine computes ir fluxes due to water vapor, co2, o3,
                0369 c   trace gases (n2o, ch4, cfc11, cfc12, cfc22, co2-minor),
                0370 c   clouds, and aerosols.
                0371 c  
                0372 c This is a vectorized code.  It computes fluxes simultaneously for
                0373 c   (m x n) soundings, which is a subset of (m x ndim) soundings.
                0374 c   In a global climate model, m and ndim correspond to the numbers of
                0375 c   grid boxes in the zonal and meridional directions, respectively.
                0376 c
                0377 c Some detailed descriptions of the radiation routine are given in
                0378 c   Chou and Suarez (1994).
                0379 c
                0380 c Ice and liquid cloud particles are allowed to co-exist in any of the
                0381 c  np layers. 
                0382 c
                0383 c If no information is available for the effective cloud particle size,
                0384 c  reff, default values of 10 micron for liquid water and 75 micron
                0385 c  for ice can be used.
                0386 c
                0387 c The maximum-random assumption is applied for cloud overlapping.
                0388 c  clouds are grouped into high, middle, and low clouds separated by the
                0389 c  level indices ict and icb.  Within each of the three groups, clouds
                0390 c  are assumed maximally overlapped, and the cloud cover of a group is
                0391 c  the maximum cloud cover of all the layers in the group.  clouds among
                0392 c  the three groups are assumed randomly overlapped. The indices ict and
                0393 c  icb correpond approximately to the 400 mb and 700 mb levels.
                0394 c
                0395 c Aerosols are allowed to be in any of the np layers. Aerosol optical
                0396 c  properties can be specified as functions of height and spectral band.
                0397 c
                0398 c There are options for computing fluxes:
                0399 c
                0400 c   If cldwater=.true., taucl is computed from cwc and reff as a
                0401 c   function of height and spectral band. 
                0402 c   If cldwater=.false., taucl must be given as input to the radiation
                0403 c   routine. It is independent of spectral band.
                0404 c
                0405 c   If high = .true., transmission functions in the co2, o3, and the
                0406 c   three water vapor bands with strong absorption are computed using
                0407 c   table look-up.  cooling rates are computed accurately from the
                0408 c   surface up to 0.01 mb.
                0409 c   If high = .false., transmission functions are computed using the
                0410 c   k-distribution method with linear pressure scaling for all spectral
                0411 c   bands and gases.  cooling rates are not calculated accurately for
                0412 c   pressures less than 20 mb. the computation is faster with
                0413 c   high=.false. than with high=.true.
                0414 
                0415 c   If trace = .true., absorption due to n2o, ch4, cfcs, and the 
                0416 c   two minor co2 bands in the window region is included.
                0417 c   If trace = .false., absorption in minor bands is neglected.
                0418 c
                0419 c The IR spectrum is divided into nine bands:
                0420 c   
                0421 c   band     wavenumber (/cm)   absorber
                0422 c
                0423 c    1           0 - 340           h2o
                0424 c    2         340 - 540           h2o
                0425 c    3         540 - 800       h2o,cont,co2,n2o
                0426 c    4         800 - 980       h2o,cont
                0427 c                              co2,f11,f12,f22
                0428 c    5         980 - 1100      h2o,cont,o3
                0429 c                              co2,f11
                0430 c    6        1100 - 1215      h2o,cont
                0431 c                              n2o,ch4,f12,f22
                0432 c    7        1215 - 1380          h2o
                0433 c                              n2o,ch4
                0434 c    8        1380 - 1900          h2o
                0435 c    9        1900 - 3000          h2o
                0436 c
                0437 c In addition, a narrow band in the 15 micrometer region is added to
                0438 c    compute flux reduction due to n2o
                0439 c
                0440 c    10        540 - 620       h2o,cont,co2,n2o
                0441 c
                0442 c Band 3 (540-800/cm) is further divided into 3 sub-bands :
                0443 c
                0444 c   subband   wavenumber (/cm)
                0445 c
                0446 c    1          540 - 620
                0447 c    2          620 - 720
                0448 c    3          720 - 800
                0449 c
                0450 c---- Input parameters                               units    size
                0451 c
                0452 c   number of soundings in zonal direction (m)        n/d      1
                0453 c   number of soundings in meridional direction (n)   n/d      1
                0454 c   maximum number of soundings in
                0455 c                 meridional direction (ndim>=n)      n/d      1 
                0456 c   number of atmospheric layers (np)                 n/d      1
                0457 c   level pressure (pl)                               mb      m*ndim*(np+1)
                0458 c   layer temperature (ta)                            k       m*ndim*np
                0459 c   layer specific humidity (wa)                      g/g     m*ndim*np
                0460 c   layer ozone mixing ratio by mass (oa)             g/g     m*ndim*np
                0461 c   surface temperature (ts)                          k       m*ndim  
                0462 c   co2 mixing ratio by volumn (co2)                  pppv     1
                0463 c   n2o mixing ratio by volumn (n2o)                  pppv     np
                0464 c   ch4 mixing ratio by volumn (ch4)                  pppv     np
                0465 c   cfc11 mixing ratio by volumn (cfc11)              pppv     1
                0466 c   cfc12 mixing ratio by volumn (cfc12)              pppv     1
                0467 c   cfc22 mixing ratio by volumn (cfc22)              pppv     1
                0468 c   surface emissivity (emiss)                      fraction   m*ndim*10
                0469 c   input option for cloud optical thickness          n/d      1
                0470 c           cldwater="true" if cwc is provided
                0471 c           cldwater="false" if taucl is provided
                0472 c   cloud water mixing ratio (cwc)                   gm/gm   m*ndim*np*3
                0473 c           index 1 for ice particles
                0474 c           index 2 for liquid drops
                0475 c           index 3 for rain drops
                0476 c   cloud optical thickness (taucl)                   n/d    m*ndim*np*3
                0477 c           index 1 for ice particles
                0478 c           index 2 for liquid drops
                0479 c           index 3 for rain drops
                0480 c   effective cloud-particle size (reff)          micrometer m*ndim*np*3
                0481 c           index 1 for ice particles
                0482 c           index 2 for liquid drops
                0483 c           index 3 for rain drops
                0484 c   cloud amount (fcld)                             fraction  m*ndim*np
                0485 c   level index separating high and middle            n/d      1
                0486 c           clouds (ict)
                0487 c   level index separating middle and low             n/d      1
                0488 c           clouds (icb)
                0489 c   aerosol optical thickness (taual)                 n/d   m*ndim*np*10
                0490 c   aerosol single-scattering albedo (ssaal)          n/d   m*ndim*np*10
                0491 c   aerosol asymmetry factor (asyal)                  n/d   m*ndim*np*10
                0492 c   high (see explanation above)                               1
                0493 c   trace (see explanation above)                              1
                0494 c
                0495 c Data used in table look-up for transmittance calculations:
                0496 c
                0497 c   c1 , c2, c3: for co2 (band 3)
                0498 c   o1 , o2, o3: for  o3 (band 5)
                0499 c   h11,h12,h13: for h2o (band 1)
                0500 c   h21,h22,h23: for h2o (band 2)
                0501 c   h81,h82,h83: for h2o (band 8)
                0502 c
                0503 c---- output parameters
                0504 c
                0505 c   net downward flux, all-sky   (flx)             w/m**2  m*ndim*(np+1)
                0506 c   net downward flux, clear-sky (flc)             w/m**2  m*ndim*(np+1)
                0507 c   sensitivity of net downward flux  
                0508 c       to surface temperature (dfdts)            w/m**2/k m*ndim*(np+1)
                0509 c   emission by the surface (st4)                 w/m**2     m*ndim 
                0510 c 
                0511 c Notes: 
                0512 c
                0513 c   (1) Water vapor continuum absorption is included in 540-1380 /cm.
                0514 c   (2) Scattering is parameterized for clouds and aerosols.
                0515 c   (3) Diffuse cloud and aerosol transmissions are computed
                0516 c       from exp(-1.66*tau).
                0517 c   (4) If there are no clouds, flx=flc.
                0518 c   (5) plevel(1) is the pressure at the top of the model atmosphere,
                0519 c        and plevel(np+1) is the surface pressure.
                0520 c   (6) Downward flux is positive and upward flux is negative.
                0521 c   (7) dfdts is negative because upward flux is defined as negative.
                0522 c   (8) For questions and coding errors, contact with Ming-Dah Chou,
                0523 c       Code 913, NASA/Goddard Space Flight Center, Greenbelt, MD 20771.
                0524 c       Phone: 301-286-4012, Fax: 301-286-1759,
                0525 c       e-mail: chou@climate.gsfc.nasa.gov
                0526 c
                0527 c-----parameters defining the size of the pre-computed tables for transmittance
                0528 c     calculations using table look-up.
                0529 c
                0530 c     "nx" is the number of intervals in pressure
                0531 c     "no" is the number of intervals in o3 amount
                0532 c     "nc" is the number of intervals in co2 amount
                0533 c     "nh" is the number of intervals in h2o amount
                0534 c     "nt" is the number of copies to be made from the pre-computed
                0535 c          transmittance tables to reduce "memory-bank conflict"
                0536 c          in parallel machines and, hence, enhancing the speed of
                0537 c          computations using table look-up. 
                0538 c          If such advantage does not exist, "nt" can be set to 1.
                0539 c***************************************************************************
                0540 
                0541 cfpp$ expand (h2oexps)
                0542 cfpp$ expand (conexps)
                0543 cfpp$ expand (co2exps)
                0544 cfpp$ expand (n2oexps)
                0545 cfpp$ expand (ch4exps)
                0546 cfpp$ expand (comexps)
                0547 cfpp$ expand (cfcexps)
                0548 cfpp$ expand (b10exps)
                0549 cfpp$ expand (tablup)
                0550 cfpp$ expand (h2okdis)
                0551 cfpp$ expand (co2kdis)
                0552 cfpp$ expand (n2okdis)
                0553 cfpp$ expand (ch4kdis)
                0554 cfpp$ expand (comkdis)
                0555 cfpp$ expand (cfckdis)
                0556 cfpp$ expand (b10kdis)
                0557 
                0558       implicit none
                0559       integer nx,no,nc,nh,nt
                0560       integer i,j,k,ip,iw,it,ib,ik,iq,isb,k1,k2
                0561       parameter (nx=26,no=21,nc=24,nh=31,nt=7)
                0562 
                0563 c---- input parameters ------
                0564 
                0565       integer m,n,ndim,np,ict,icb
a456aa407c Andr*0566       _RL pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np),
1662f365b2 Andr*0567      *     ts(m,ndim)
a456aa407c Andr*0568       _RL co2,n2o(np),ch4(np),cfc11,cfc12,cfc22,emiss(m,ndim,10)
                0569       _RL cwc(m,ndim,np,3),taucl(m,ndim,np,3),reff(m,ndim,np,3),
1662f365b2 Andr*0570      *     fcld(m,ndim,np)
a456aa407c Andr*0571       _RL taual(m,ndim,np,10),ssaal(m,ndim,np,10),asyal(m,ndim,np,10)
1662f365b2 Andr*0572       logical cldwater,high,trace
                0573 
                0574 c---- output parameters ------
                0575 
a456aa407c Andr*0576       _RL flx(m,ndim,np+1),flc(m,ndim,np+1),dfdts(m,ndim,np+1),
1662f365b2 Andr*0577      *     st4(m,ndim)
                0578 
                0579 c---- static data -----
                0580 
a456aa407c Andr*0581       _RL cb(5,10)
                0582       _RL xkw(9),aw(9),bw(9),pm(9),fkw(6,9),gkw(6,3),xke(9)
                0583       _RL aib(3,10),awb(4,10),aiw(4,10),aww(4,10),aig(4,10),awg(4,10)
1662f365b2 Andr*0584       integer ne(9),mw(9)
                0585 
                0586 c---- temporary arrays -----
                0587 
a456aa407c Andr*0588       _RL pa(m,n,np),dt(m,n,np)
                0589       _RL sh2o(m,n,np+1),swpre(m,n,np+1),swtem(m,n,np+1)
                0590       _RL sco3(m,n,np+1),scopre(m,n,np+1),scotem(m,n,np+1)
                0591       _RL dh2o(m,n,np),dcont(m,n,np),dco2(m,n,np),do3(m,n,np)
                0592       _RL dn2o(m,n,np),dch4(m,n,np)
                0593       _RL df11(m,n,np),df12(m,n,np),df22(m,n,np)
                0594       _RL th2o(m,n,6),tcon(m,n,3),tco2(m,n,6,2)
                0595       _RL tn2o(m,n,4),tch4(m,n,4),tcom(m,n,2)
                0596       _RL tf11(m,n),tf12(m,n),tf22(m,n)
                0597       _RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
                0598       _RL n2oexp(m,n,np,4),ch4exp(m,n,np,4),comexp(m,n,np,2)
                0599       _RL f11exp(m,n,np),  f12exp(m,n,np),  f22exp(m,n,np)
                0600       _RL clr(m,n,0:np+1),fclr(m,n)
                0601       _RL blayer(m,n,0:np+1),dlayer(m,n,np+1),dbs(m,n)
                0602       _RL clrlw(m,n),clrmd(m,n),clrhi(m,n)
                0603       _RL cwp(m,n,np,3)
                0604       _RL trant(m,n),tranal(m,n),transfc(m,n,np+1),trantcr(m,n,np+1)
                0605       _RL flxu(m,n,np+1),flxd(m,n,np+1),flcu(m,n,np+1),flcd(m,n,np+1)
                0606       _RL rflx(m,n,np+1),rflc(m,n,np+1)
1662f365b2 Andr*0607 
                0608       logical oznbnd,co2bnd,h2otbl,conbnd,n2obnd
                0609       logical ch4bnd,combnd,f11bnd,f12bnd,f22bnd,b10bnd
                0610 
a456aa407c Andr*0611       _RL c1 (nx,nc,nt),c2 (nx,nc,nt),c3 (nx,nc,nt)
                0612       _RL o1 (nx,no,nt),o2 (nx,no,nt),o3 (nx,no,nt)
                0613       _RL h11(nx,nh,nt),h12(nx,nh,nt),h13(nx,nh,nt)
                0614       _RL h21(nx,nh,nt),h22(nx,nh,nt),h23(nx,nh,nt)
                0615       _RL h81(nx,nh,nt),h82(nx,nh,nt),h83(nx,nh,nt)
1662f365b2 Andr*0616 
a456aa407c Andr*0617       _RL dp,xx,p1,dwe,dpe,a1,b1,fk1,a2,b2,fk2
                0618       _RL w1,w2,w3,g1,g2,g3,ww,gg,ff,taux,reff1,reff2
1662f365b2 Andr*0619 
                0620 c-----the following coefficients (equivalent to table 2 of
                0621 c     chou and suarez, 1995) are for computing spectrally 
                0622 c     integrated planck fluxes using eq. (22)
                0623 
                0624        data cb/
                0625      1      -2.6844e-1,-8.8994e-2, 1.5676e-3,-2.9349e-6, 2.2233e-9,
                0626      2       3.7315e+1,-7.4758e-1, 4.6151e-3,-6.3260e-6, 3.5647e-9,
                0627      3       3.7187e+1,-3.9085e-1,-6.1072e-4, 1.4534e-5,-1.6863e-8,
                0628      4      -4.1928e+1, 1.0027e+0,-8.5789e-3, 2.9199e-5,-2.5654e-8,
                0629      5      -4.9163e+1, 9.8457e-1,-7.0968e-3, 2.0478e-5,-1.5514e-8,
                0630      6      -4.7107e+1, 8.9130e-1,-5.9735e-3, 1.5596e-5,-9.5911e-9,
                0631      7      -5.4041e+1, 9.5332e-1,-5.7784e-3, 1.2555e-5,-2.9377e-9,
                0632      8      -6.9233e+0,-1.5878e-1, 3.9160e-3,-2.4496e-5, 4.9301e-8,
                0633      9       1.1483e+2,-2.2376e+0, 1.6394e-2,-5.3672e-5, 6.6456e-8,
                0634      *       1.9668e+1,-3.1161e-1, 1.2886e-3, 3.6835e-7,-1.6212e-9/
                0635 
                0636 c-----xkw are the absorption coefficients for the first
                0637 c     k-distribution function due to water vapor line absorption
                0638 c     (tables 4 and 7).  units are cm**2/g    
                0639 
                0640       data xkw / 29.55  , 4.167e-1, 1.328e-2, 5.250e-4,
                0641      *           5.25e-4, 9.369e-3, 4.719e-2, 1.320e-0, 5.250e-4/
                0642 
                0643 c-----xke are the absorption coefficients for the first
                0644 c     k-distribution function due to water vapor continuum absorption
                0645 c     (table 6).  units are cm**2/g
                0646 
                0647       data xke /  0.00,   0.00,   27.40,   15.8,
                0648      *            9.40,   7.75,     0.0,    0.0,    0.0/
                0649 
                0650 c-----mw are the ratios between neighboring absorption coefficients
                0651 c     for water vapor line absorption (tables 4 and 7).
                0652 
                0653       data mw /6,6,8,6,6,8,9,6,16/
                0654 
                0655 c-----aw and bw (table 3) are the coefficients for temperature scaling
                0656 c     in eq. (25).
                0657 
                0658       data aw/ 0.0021, 0.0140, 0.0167, 0.0302,
                0659      *         0.0307, 0.0195, 0.0152, 0.0008, 0.0096/
                0660       data bw/ -1.01e-5, 5.57e-5, 8.54e-5, 2.96e-4,
                0661      *          2.86e-4, 1.108e-4, 7.608e-5, -3.52e-6, 1.64e-5/
                0662 
                0663       data pm/ 1.0, 1.0, 1.0, 1.0, 1.0, 0.77, 0.5, 1.0, 1.0/
                0664 
                0665 c-----fkw is the planck-weighted k-distribution function due to h2o
                0666 c     line absorption given in table 4 of Chou and Suarez (1995).
                0667 c     the k-distribution function for the third band, fkw(*,3), is not used
                0668 
                0669       data fkw / 0.2747,0.2717,0.2752,0.1177,0.0352,0.0255,
                0670      2           0.1521,0.3974,0.1778,0.1826,0.0374,0.0527,
                0671      3           6*1.00,
                0672      4           0.4654,0.2991,0.1343,0.0646,0.0226,0.0140,
                0673      5           0.5543,0.2723,0.1131,0.0443,0.0160,0.0000,
                0674      6           0.5955,0.2693,0.0953,0.0335,0.0064,0.0000,
                0675      7           0.1958,0.3469,0.3147,0.1013,0.0365,0.0048,
                0676      8           0.0740,0.1636,0.4174,0.1783,0.1101,0.0566,
                0677      9           0.1437,0.2197,0.3185,0.2351,0.0647,0.0183/
                0678 
                0679 c-----gkw is the planck-weighted k-distribution function due to h2o
                0680 c     line absorption in the 3 subbands (800-720,620-720,540-620 /cm)
                0681 c     of band 3 given in table 7.  Note that the order of the sub-bands
                0682 c     is reversed.
                0683 
                0684       data gkw/  0.1782,0.0593,0.0215,0.0068,0.0022,0.0000,
                0685      2           0.0923,0.1675,0.0923,0.0187,0.0178,0.0000,
                0686      3           0.0000,0.1083,0.1581,0.0455,0.0274,0.0041/
                0687 
                0688 c-----ne is the number of terms used in each band to compute water vapor
                0689 c     continuum transmittance (table 6).
                0690 
                0691       data ne /0,0,3,1,1,1,0,0,0/
                0692 c
                0693 c-----coefficients for computing the extinction coefficient
                0694 c     for cloud ice particles
                0695 c     polynomial fit: y=a1+a2/x**a3; x is in m**2/gm
                0696 c
                0697       data aib /  -0.44171,    0.62951,   0.06465,
                0698      2            -0.13727,    0.61291,   0.28962,
                0699      3            -0.01878,    1.67680,   0.79080,
                0700      4            -0.01896,    1.06510,   0.69493,
                0701      5            -0.04788,    0.88178,   0.54492,
                0702      6            -0.02265,    1.57390,   0.76161,
                0703      7            -0.01038,    2.15640,   0.89045, 
                0704      8            -0.00450,    2.51370,   0.95989,
                0705      9            -0.00044,    3.15050,   1.03750,
                0706      *            -0.02956,    1.44680,   0.71283/
                0707 c
                0708 c-----coefficients for computing the extinction coefficient
                0709 c     for cloud liquid drops
                0710 c     polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
                0711 c
                0712       data awb /   0.08641,    0.01769,    -1.5572e-3,   3.4896e-5,
                0713      2             0.22027,    0.00997,    -1.8719e-3,   5.3112e-5,
                0714      3             0.39252,   -0.02817,     7.2931e-4,  -3.8151e-6,
                0715      4             0.39975,   -0.03426,     1.2884e-3,  -1.7930e-5,
                0716      5             0.34021,   -0.02805,     1.0654e-3,  -1.5443e-5,
                0717      6             0.15587,    0.00371,    -7.7705e-4,   2.0547e-5,
                0718      7             0.05518,    0.04544,    -4.2067e-3,   1.0184e-4,
                0719      8             0.12724,    0.04751,    -5.2037e-3,   1.3711e-4,
                0720      9             0.30390,    0.01656,    -3.5271e-3,   1.0828e-4,
                0721      *             0.63617,   -0.06287,     2.2350e-3,  -2.3177e-5/
                0722 c
                0723 c-----coefficients for computing the single-scattering albedo
                0724 c     for cloud ice particles
                0725 c     polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
                0726 c
                0727       data aiw/    0.17201,    1.2229e-2,  -1.4837e-4,   5.8020e-7,
                0728      2             0.81470,   -2.7293e-3,   9.7816e-8,   5.7650e-8,
                0729      3             0.54859,   -4.8273e-4,   5.4353e-6,  -1.5679e-8,
                0730      4             0.39218,    4.1717e-3, - 4.8869e-5,   1.9144e-7,
                0731      5             0.71773,   -3.3640e-3,   1.9713e-5,  -3.3189e-8,
                0732      6             0.77345,   -5.5228e-3,   4.8379e-5,  -1.5151e-7,
                0733      7             0.74975,   -5.6604e-3,   5.6475e-5,  -1.9664e-7,
                0734      8             0.69011,   -4.5348e-3,   4.9322e-5,  -1.8255e-7,
                0735      9             0.83963,   -6.7253e-3,   6.1900e-5,  -2.0862e-7,
                0736      *             0.64860,   -2.8692e-3,   2.7656e-5,  -8.9680e-8/
                0737 c
                0738 c-----coefficients for computing the single-scattering albedo
                0739 c     for cloud liquid drops
                0740 c     polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
                0741 c
                0742       data aww/   -7.8566e-2,  8.0875e-2,  -4.3403e-3,   8.1341e-5,
                0743      2            -1.3384e-2,  9.3134e-2,  -6.0491e-3,   1.3059e-4,
                0744      3             3.7096e-2,  7.3211e-2,  -4.4211e-3,   9.2448e-5,
                0745      4            -3.7600e-3,  9.3344e-2,  -5.6561e-3,   1.1387e-4,
                0746      5             0.40212,    7.8083e-2,  -5.9583e-3,   1.2883e-4,
                0747      6             0.57928,    5.9094e-2,  -5.4425e-3,   1.2725e-4,
                0748      7             0.68974,    4.2334e-2,  -4.9469e-3,   1.2863e-4,
                0749      8             0.80122,    9.4578e-3,  -2.8508e-3,   9.0078e-5,
                0750      9             1.02340,   -2.6204e-2,   4.2552e-4,   3.2160e-6,
                0751      *             0.05092,    7.5409e-2,  -4.7305e-3,   1.0121e-4/ 
                0752 c
                0753 c-----coefficients for computing the asymmetry factor for cloud ice particles
                0754 c     polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
                0755 c
                0756       data aig /   0.57867,    1.0135e-2,  -1.1142e-4,   4.1537e-7,
                0757      2             0.72259,    3.1149e-3,  -1.9927e-5,   5.6024e-8,
                0758      3             0.76109,    4.5449e-3,  -4.6199e-5,   1.6446e-7,
                0759      4             0.86934,    2.7474e-3,  -3.1301e-5,   1.1959e-7,
                0760      5             0.89103,    1.8513e-3,  -1.6551e-5,   5.5193e-8,
                0761      6             0.86325,    2.1408e-3,  -1.6846e-5,   4.9473e-8,
                0762      7             0.85064,    2.5028e-3,  -2.0812e-5,   6.3427e-8,
                0763      8             0.86945,    2.4615e-3,  -2.3882e-5,   8.2431e-8,
                0764      9             0.80122,    3.1906e-3,  -2.4856e-5,   7.2411e-8,
                0765      *             0.73290,    4.8034e-3,  -4.4425e-5,   1.4839e-7/
                0766 c
                0767 c-----coefficients for computing the asymmetry factor for cloud liquid drops
                0768 c     polynomial fit: y=a1+a2*x+a3*x**2+a4*x**3; x is in m**2/gm
                0769 c
                0770       data awg /  -0.51930,    0.20290,    -1.1747e-2,   2.3868e-4,
                0771      2            -0.22151,    0.19708,    -1.2462e-2,   2.6646e-4,
                0772      3             0.14157,    0.14705,    -9.5802e-3,   2.0819e-4,
                0773      4             0.41590,    0.10482,    -6.9118e-3,   1.5115e-4,
                0774      5             0.55338,    7.7016e-2,  -5.2218e-3,   1.1587e-4,
                0775      6             0.61384,    6.4402e-2,  -4.6241e-3,   1.0746e-4,
                0776      7             0.67891,    4.8698e-2,  -3.7021e-3,   9.1966e-5,
                0777      8             0.78169,    2.0803e-2,  -1.4749e-3,   3.9362e-5,
                0778      9             0.93218,   -3.3425e-2,   2.9632e-3,  -6.9362e-5,
                0779      *             0.01649,    0.16561,    -1.0723e-2,   2.3220e-4/ 
                0780 c
                0781 c-----include tables used in the table look-up for co2 (band 3), 
                0782 c     o3 (band 5), and h2o (bands 1, 2, and 7) transmission functions.
                0783 
                0784       logical first
                0785       data first /.true./
                0786 
65b2c6ac12 Andr*0787 #include "h2o-tran3.h"
                0788 #include "co2-tran3.h"
                0789 #include "o3-tran3.h"
1662f365b2 Andr*0790 
032fb71841 Andr*0791 c     save c1,c2,c3,o1,o2,o3
                0792 c     save h11,h12,h13,h21,h22,h23,h81,h82,h83
1662f365b2 Andr*0793 
154b1b9078 Andr*0794       if (first) then
1662f365b2 Andr*0795 
                0796 c-----tables co2 and h2o are only used with 'high' option
                0797 
                0798        if (high) then
                0799 
                0800         do iw=1,nh
                0801          do ip=1,nx
                0802           h11(ip,iw,1)=1.0-h11(ip,iw,1)
                0803           h21(ip,iw,1)=1.0-h21(ip,iw,1)
                0804           h81(ip,iw,1)=1.0-h81(ip,iw,1)
                0805          enddo
                0806         enddo
                0807 
                0808         do iw=1,nc
                0809          do ip=1,nx
                0810           c1(ip,iw,1)=1.0-c1(ip,iw,1)
                0811          enddo
                0812         enddo
                0813 
                0814 c-----copy tables to enhance the speed of co2 (band 3), o3 (band 5),
                0815 c     and h2o (bands 1, 2, and 8 only) transmission calculations
                0816 c     using table look-up.
                0817 c-----tables are replicated to avoid memory bank conflicts
                0818 
                0819         do it=2,nt
                0820          do iw=1,nc
                0821           do ip=1,nx
                0822            c1 (ip,iw,it)= c1(ip,iw,1)
                0823            c2 (ip,iw,it)= c2(ip,iw,1)
                0824            c3 (ip,iw,it)= c3(ip,iw,1)
                0825           enddo
                0826          enddo
                0827          do iw=1,nh
                0828           do ip=1,nx
                0829            h11(ip,iw,it)=h11(ip,iw,1)
                0830            h12(ip,iw,it)=h12(ip,iw,1)
                0831            h13(ip,iw,it)=h13(ip,iw,1)
                0832            h21(ip,iw,it)=h21(ip,iw,1)
                0833            h22(ip,iw,it)=h22(ip,iw,1)
                0834            h23(ip,iw,it)=h23(ip,iw,1)
                0835            h81(ip,iw,it)=h81(ip,iw,1)
                0836            h82(ip,iw,it)=h82(ip,iw,1)
                0837            h83(ip,iw,it)=h83(ip,iw,1)
                0838           enddo
                0839          enddo
                0840         enddo
                0841 
                0842        endif
                0843 
                0844 c-----always use table look-up for ozone transmittance
                0845 
                0846         do iw=1,no
                0847          do ip=1,nx
                0848           o1(ip,iw,1)=1.0-o1(ip,iw,1)
                0849          enddo
                0850         enddo
                0851 
                0852         do it=2,nt
                0853          do iw=1,no
                0854           do ip=1,nx
                0855            o1 (ip,iw,it)= o1(ip,iw,1)
                0856            o2 (ip,iw,it)= o2(ip,iw,1)
                0857            o3 (ip,iw,it)= o3(ip,iw,1)
                0858           enddo
                0859          enddo
                0860         enddo
                0861 
154b1b9078 Andr*0862        first=.false.
1662f365b2 Andr*0863 
154b1b9078 Andr*0864       endif
1662f365b2 Andr*0865 
                0866 c-----set the pressure at the top of the model atmosphere
                0867 c     to 1.0e-4 if it is zero
                0868 
                0869       do j=1,n
                0870        do i=1,m
                0871           if (pl(i,j,1) .eq. 0.0) then
                0872               pl(i,j,1)=1.0e-4
                0873           endif
                0874        enddo
                0875       enddo
                0876 
                0877 c-----compute layer pressure (pa) and layer temperature minus 250K (dt)
                0878 
                0879       do k=1,np
                0880        do j=1,n
                0881         do i=1,m
                0882          pa(i,j,k)=0.5*(pl(i,j,k)+pl(i,j,k+1))
                0883          dt(i,j,k)=ta(i,j,k)-250.0
                0884         enddo
                0885        enddo
                0886       enddo
                0887 
                0888 c-----compute layer absorber amount
                0889 
                0890 c     dh2o : water vapor amount (g/cm**2)
                0891 c     dcont: scaled water vapor amount for continuum absorption
                0892 c            (g/cm**2)
                0893 c     dco2 : co2 amount (cm-atm)stp
                0894 c     do3  : o3 amount (cm-atm)stp
                0895 c     dn2o : n2o amount (cm-atm)stp
                0896 c     dch4 : ch4 amount (cm-atm)stp
                0897 c     df11 : cfc11 amount (cm-atm)stp
                0898 c     df12 : cfc12 amount (cm-atm)stp
                0899 c     df22 : cfc22 amount (cm-atm)stp
                0900 c     the factor 1.02 is equal to 1000/980
                0901 c     factors 789 and 476 are for unit conversion
                0902 c     the factor 0.001618 is equal to 1.02/(.622*1013.25) 
                0903 c     the factor 6.081 is equal to 1800/296
                0904 
                0905 
                0906       do k=1,np
                0907        do j=1,n
                0908         do i=1,m
                0909 
                0910          dp           = pl(i,j,k+1)-pl(i,j,k)
                0911          dh2o (i,j,k) = 1.02*wa(i,j,k)*dp+1.e-10
                0912          do3  (i,j,k) = 476.*oa(i,j,k)*dp+1.e-10
                0913          dco2 (i,j,k) = 789.*co2*dp+1.e-10
                0914          dch4 (i,j,k) = 789.*ch4(k)*dp+1.e-10
                0915          dn2o (i,j,k) = 789.*n2o(k)*dp+1.e-10
                0916          df11 (i,j,k) = 789.*cfc11*dp+1.e-10
                0917          df12 (i,j,k) = 789.*cfc12*dp+1.e-10
                0918          df22 (i,j,k) = 789.*cfc22*dp+1.e-10
                0919 
                0920 c-----compute scaled water vapor amount for h2o continuum absorption
                0921 c     following eq. (43).
                0922 
                0923          xx=pa(i,j,k)*0.001618*wa(i,j,k)*wa(i,j,k)*dp
                0924          dcont(i,j,k) = xx*exp(1800./ta(i,j,k)-6.081)+1.e-10
                0925 
                0926         enddo
                0927        enddo
                0928       enddo
                0929 
                0930 c-----compute column-integrated h2o amoumt, h2o-weighted pressure
                0931 c     and temperature.  it follows eqs. (37) and (38).
                0932 
                0933        if (high) then
                0934         call column(m,n,np,pa,dt,dh2o,sh2o,swpre,swtem)
                0935        endif
                0936 
                0937 c-----compute layer cloud water amount (gm/m**2)
                0938 c     index is 1 for ice, 2 for waterdrops and 3 for raindrops.
                0939 c     Rain optical thickness is 0.00307 /(gm/m**2).
                0940 c     It is for a specific drop size distribution provided by Q. Fu.
                0941 
                0942        if (cldwater) then
                0943         do k=1,np
                0944          do j=1,n
                0945           do i=1,m
                0946              dp          =pl(i,j,k+1)-pl(i,j,k)
                0947              cwp(i,j,k,1)=1.02*10000.0*cwc(i,j,k,1)*dp
                0948              cwp(i,j,k,2)=1.02*10000.0*cwc(i,j,k,2)*dp
                0949              cwp(i,j,k,3)=1.02*10000.0*cwc(i,j,k,3)*dp
                0950              taucl(i,j,k,3)=0.00307*cwp(i,j,k,3)
                0951           enddo
                0952          enddo
                0953         enddo
                0954        endif
                0955 
                0956 c-----the surface (np+1) is treated as a layer filled with black clouds.
                0957 c     clr is the equivalent clear fraction of a layer
                0958 c     transfc is the transmttance between the surface and a pressure level.
                0959 c     trantcr is the clear-sky transmttance between the surface and a
                0960 c     pressure level.
                0961  
                0962       do j=1,n
                0963        do i=1,m
                0964         clr(i,j,0)    = 1.0
                0965         clr(i,j,np+1) = 0.0
                0966         st4(i,j)      = 0.0
                0967         transfc(i,j,np+1)=1.0
                0968         trantcr(i,j,np+1)=1.0
                0969        enddo
                0970       enddo
                0971 
                0972 c-----initialize fluxes
                0973 
                0974       do k=1,np+1
                0975        do j=1,n
                0976         do i=1,m
                0977          flx(i,j,k)  = 0.0
                0978          flc(i,j,k)  = 0.0
                0979          dfdts(i,j,k)= 0.0
                0980          rflx(i,j,k) = 0.0
                0981          rflc(i,j,k) = 0.0
                0982         enddo
                0983        enddo
                0984       enddo
                0985 
                0986 c-----integration over spectral bands
                0987 
                0988       do 1000 ib=1,10
                0989 
                0990 c-----if h2otbl, compute h2o (line) transmittance using table look-up.
                0991 c     if conbnd, compute h2o (continuum) transmittance in bands 3-6.
                0992 c     if co2bnd, compute co2 transmittance in band 3.
                0993 c     if oznbnd, compute  o3 transmittance in band 5.
                0994 c     if n2obnd, compute n2o transmittance in bands 6 and 7.
                0995 c     if ch4bnd, compute ch4 transmittance in bands 6 and 7.
                0996 c     if combnd, compute co2-minor transmittance in bands 4 and 5.
                0997 c     if f11bnd, compute cfc11 transmittance in bands 4 and 5.
                0998 c     if f12bnd, compute cfc12 transmittance in bands 4 and 6.
                0999 c     if f22bnd, compute cfc22 transmittance in bands 4 and 6.
                1000 c     if b10bnd, compute flux reduction due to n2o in band 10.
                1001 
                1002        h2otbl=high.and.(ib.eq.1.or.ib.eq.2.or.ib.eq.8)
                1003        conbnd=ib.ge.3.and.ib.le.6
                1004        co2bnd=ib.eq.3
                1005        oznbnd=ib.eq.5
                1006        n2obnd=ib.eq.6.or.ib.eq.7
                1007        ch4bnd=ib.eq.6.or.ib.eq.7
                1008        combnd=ib.eq.4.or.ib.eq.5
                1009        f11bnd=ib.eq.4.or.ib.eq.5
                1010        f12bnd=ib.eq.4.or.ib.eq.6
                1011        f22bnd=ib.eq.4.or.ib.eq.6
                1012        b10bnd=ib.eq.10
                1013 
                1014 c-----blayer is the spectrally integrated planck flux of the mean layer
                1015 c     temperature derived from eq. (22)
                1016 c     the fitting for the planck flux is valid in the range 160-345 K.
                1017 
                1018        do k=1,np
                1019         do j=1,n
                1020          do i=1,m
                1021           blayer(i,j,k)=ta(i,j,k)*(ta(i,j,k)*(ta(i,j,k)
                1022      *                 *(ta(i,j,k)*cb(5,ib)+cb(4,ib))+cb(3,ib))
                1023      *                 +cb(2,ib))+cb(1,ib)
                1024          enddo
                1025         enddo
                1026        enddo
                1027 
3daafce20b Jean*1028 c-----the earth surface, with an index "np+1", is treated as a layer
1662f365b2 Andr*1029 
                1030        do j=1,n
                1031         do i=1,m
                1032          blayer(i,j,0)   = 0.0
                1033          blayer(i,j,np+1)= ( ts(i,j)*(ts(i,j)*(ts(i,j)
                1034      *                    *(ts(i,j)*cb(5,ib)+cb(4,ib))+cb(3,ib))
                1035      *                    +cb(2,ib))+cb(1,ib) )*emiss(i,j,ib)
                1036 
                1037 c-----dbs is the derivative of the surface emission with respect to
                1038 c     surface temperature (eq. 59).
                1039 
                1040          dbs(i,j)=(ts(i,j)*(ts(i,j)*(ts(i,j)*4.*cb(5,ib)
                1041      *           +3.*cb(4,ib))+2.*cb(3,ib))+cb(2,ib))*emiss(i,j,ib)
                1042 
                1043         enddo
                1044        enddo
                1045 
                1046        do k=1,np+1
                1047         do j=1,n
                1048          do i=1,m
                1049            dlayer(i,j,k)=blayer(i,j,k-1)-blayer(i,j,k)
                1050          enddo
                1051         enddo
                1052        enddo
                1053 
                1054 c-----compute column-integrated absorber amoumt, absorber-weighted
                1055 c     pressure and temperature for co2 (band 3) and o3 (band 5).
                1056 c     it follows eqs. (37) and (38).
                1057 
                1058 c-----this is in the band loop to save storage
                1059 
                1060       if (high .and. co2bnd) then
                1061         call column(m,n,np,pa,dt,dco2,sco3,scopre,scotem)
                1062       endif
                1063 
                1064       if (oznbnd) then
                1065         call column(m,n,np,pa,dt,do3,sco3,scopre,scotem)
                1066       endif
                1067 
                1068 c-----compute cloud optical thickness
                1069 
                1070       if (cldwater) then
                1071        do k= 1, np
                1072         do j= 1, n
                1073          do i= 1, m
                1074           taucl(i,j,k,1)=cwp(i,j,k,1)*(aib(1,ib)+aib(2,ib)/
                1075      *      reff(i,j,k,1)**aib(3,ib))
                1076           taucl(i,j,k,2)=cwp(i,j,k,2)*(awb(1,ib)+(awb(2,ib)+(awb(3,ib)
                1077      *      +awb(4,ib)*reff(i,j,k,2))*reff(i,j,k,2))*reff(i,j,k,2))
                1078          enddo
                1079         enddo
                1080        enddo
                1081       endif
                1082 
                1083 c-----compute cloud single-scattering albedo and asymmetry factor for
                1084 c     a mixture of ice particles, liquid drops, and rain drops
                1085 c     single-scattering albedo and asymmetry factor of rain are set
                1086 c     to 0.54 and 0.95, respectively.
                1087 
                1088        do k= 1, np
                1089         do j= 1, n
                1090          do i= 1, m
                1091 
                1092            clr(i,j,k)    = 1.0
                1093            taux=taucl(i,j,k,1)+taucl(i,j,k,2)+taucl(i,j,k,3)
                1094 
                1095           if (taux.gt.0.02 .and. fcld(i,j,k).gt.0.01) then
                1096 
32362b8fa8 Cons*1097            reff1=min(reff(i,j,k,1),130. _d 0)
                1098            reff2=min(reff(i,j,k,2),20.0 _d 0)
1662f365b2 Andr*1099 
                1100            w1=taucl(i,j,k,1)*(aiw(1,ib)+(aiw(2,ib)+(aiw(3,ib)
                1101      *       +aiw(4,ib)*reff1)*reff1)*reff1)
                1102            w2=taucl(i,j,k,2)*(aww(1,ib)+(aww(2,ib)+(aww(3,ib)
                1103      *       +aww(4,ib)*reff2)*reff2)*reff2)
                1104            w3=taucl(i,j,k,3)*0.54
                1105            ww=(w1+w2+w3)/taux
                1106 
                1107            g1=w1*(aig(1,ib)+(aig(2,ib)+(aig(3,ib)
                1108      *      +aig(4,ib)*reff1)*reff1)*reff1)
                1109            g2=w2*(awg(1,ib)+(awg(2,ib)+(awg(3,ib)
                1110      *      +awg(4,ib)*reff2)*reff2)*reff2)
                1111            g3=w3*0.95
                1112 
                1113            gg=(g1+g2+g3)/(w1+w2+w3)
                1114 
                1115 c-----parameterization of LW scattering
                1116 c     compute forward-scattered fraction and scale optical thickness
                1117 
                1118            ff=0.5+(0.3739+(0.0076+0.1185*gg)*gg)*gg
                1119            taux=taux*(1.-ww*ff)
                1120 
                1121 c-----compute equivalent cloud-free fraction, clr, for each layer
                1122 c     the cloud diffuse transmittance is approximated by using 
                1123 c     a diffusivity factor of 1.66.
                1124 
                1125            clr(i,j,k)=1.0-(fcld(i,j,k)*(1.0-exp(-1.66*taux)))
                1126 
                1127           endif
                1128 
                1129          enddo
                1130         enddo
                1131        enddo
                1132 
                1133 c-----compute the exponential terms (eq. 32) at each layer due to
                1134 c     water vapor line absorption when k-distribution is used
                1135 
                1136       if (.not.h2otbl .and. .not.b10bnd) then
                1137         call h2oexps(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)
                1138       endif
                1139 
                1140 c-----compute the exponential terms (eq. 46) at each layer due to
                1141 c     water vapor continuum absorption
                1142 
                1143       if (conbnd) then
                1144         call conexps(ib,m,n,np,dcont,xke,conexp)
                1145       endif
                1146 
                1147 c-----compute the exponential terms (eq. 32) at each layer due to
                1148 c     co2 absorption
                1149 
                1150       if (.not.high .and. co2bnd) then
                1151         call co2exps(m,n,np,dco2,pa,dt,co2exp)
                1152       endif
                1153 
                1154 c***** for trace gases *****
                1155 
                1156       if (trace) then
                1157 
                1158 c-----compute the exponential terms at each layer due to n2o absorption
                1159 
                1160        if (n2obnd) then
                1161         call n2oexps(ib,m,n,np,dn2o,pa,dt,n2oexp)
                1162        endif
                1163 
                1164 c-----compute the exponential terms at each layer due to ch4 absorption
                1165 
                1166        if (ch4bnd) then
                1167         call ch4exps(ib,m,n,np,dch4,pa,dt,ch4exp)
                1168        endif
                1169 
                1170 c-----compute the exponential terms due to co2 minor absorption
                1171 
                1172        if (combnd) then
                1173         call comexps(ib,m,n,np,dco2,dt,comexp)
                1174        endif
                1175 
                1176 c-----compute the exponential terms due to cfc11 absorption
                1177 
                1178        if (f11bnd) then
                1179             a1  = 1.26610e-3
                1180             b1  = 3.55940e-6
                1181             fk1 = 1.89736e+1
                1182             a2  = 8.19370e-4
                1183             b2  = 4.67810e-6
                1184             fk2 = 1.01487e+1
                1185         call cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df11,dt,f11exp)
                1186        endif
                1187 
                1188 c-----compute the exponential terms due to cfc12 absorption
                1189 
                1190        if (f12bnd) then
                1191             a1  = 8.77370e-4
                1192             b1  =-5.88440e-6
                1193             fk1 = 1.58104e+1
                1194             a2  = 8.62000e-4
                1195             b2  =-4.22500e-6
                1196             fk2 = 3.70107e+1
                1197         call cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df12,dt,f12exp)
                1198        endif
                1199 
                1200 c-----compute the exponential terms due to cfc22 absorption
                1201 
                1202        if (f22bnd) then
                1203             a1  = 9.65130e-4
                1204             b1  = 1.31280e-5
                1205             fk1 = 6.18536e+0
                1206             a2  =-3.00010e-5 
                1207             b2  = 5.25010e-7
                1208             fk2 = 3.27912e+1
                1209         call cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,df22,dt,f22exp)
                1210        endif
                1211 
                1212 c-----compute the exponential terms at each layer in band 10 due to
                1213 c     h2o line and continuum, co2, and n2o absorption
                1214 
                1215        if (b10bnd) then
                1216         call b10exps(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt
                1217      *              ,h2oexp,conexp,co2exp,n2oexp)
                1218        endif
                1219 
                1220       endif
                1221 
                1222 c-----compute transmittances for regions between levels k1 and k2
                1223 c     and update the fluxes at the two levels.
                1224 
                1225 
                1226 c-----initialize fluxes
                1227 
                1228       do k=1,np+1
                1229        do j=1,n
                1230         do i=1,m
                1231          flxu(i,j,k) = 0.0
                1232          flxd(i,j,k) = 0.0
                1233          flcu(i,j,k) = 0.0
                1234          flcd(i,j,k) = 0.0
                1235         enddo
                1236        enddo
                1237       enddo
                1238 
                1239       do 2000 k1=1,np
                1240 
                1241 c-----initialize fclr, th2o, tcon, tco2, and tranal
                1242 c     fclr is the clear fraction of the line-of-sight
                1243 c     clrlw, clrmd, and clrhi are the clear-sky fractions of the 
                1244 c      low, middle, and high troposphere, respectively
                1245 c     tranal is the aerosol transmission function
                1246 
                1247         do j=1,n
                1248          do i=1,m
                1249           fclr(i,j)  = 1.0
                1250           clrlw(i,j) = 1.0
                1251           clrmd(i,j) = 1.0
                1252           clrhi(i,j) = 1.0
                1253           tranal(i,j)= 1.0
                1254          enddo
                1255         enddo
                1256 
                1257 c-----for h2o line transmission
                1258 
                1259       if (.not. h2otbl) then
                1260         do ik=1,6
                1261          do j=1,n
                1262           do i=1,m
                1263            th2o(i,j,ik)=1.0
                1264           enddo
                1265          enddo
                1266         enddo
                1267       endif
                1268 
                1269 c-----for h2o continuum transmission
                1270 
                1271        if (conbnd) then
                1272          do iq=1,3
                1273           do j=1,n
                1274            do i=1,m
                1275             tcon(i,j,iq)=1.0
                1276            enddo
                1277           enddo
                1278          enddo
                1279        endif
                1280 
                1281 c-----for co2 transmission using k-distribution method.
                1282 c     band 3 is divided into 3 sub-bands, but sub-bands 3a and 3c
                1283 c     are combined in computing the co2 transmittance.
                1284 
                1285        if (.not.high .and. co2bnd) then
                1286          do isb=1,2
                1287           do ik=1,6
                1288            do j=1,n
                1289             do i=1,m
                1290              tco2(i,j,ik,isb)=1.0
                1291             enddo
                1292            enddo
                1293           enddo
                1294          enddo
                1295        endif
                1296 
                1297 c***** for trace gases *****
                1298 
                1299       if (trace) then
                1300 
                1301 c-----for n2o transmission using k-distribution method.
                1302 
                1303        if (n2obnd) then
                1304           do ik=1,4
                1305            do j=1,n
                1306             do i=1,m
                1307              tn2o(i,j,ik)=1.0
                1308             enddo
                1309            enddo
                1310           enddo
                1311        endif
                1312 
                1313 c-----for ch4 transmission using k-distribution method.
                1314 
                1315        if (ch4bnd) then
                1316           do ik=1,4
                1317            do j=1,n
                1318             do i=1,m
                1319              tch4(i,j,ik)=1.0
                1320             enddo
                1321            enddo
                1322           enddo
                1323        endif
                1324 
                1325 c-----for co2-minor transmission using k-distribution method.
                1326 
                1327        if (combnd) then
                1328           do ik=1,2
                1329            do j=1,n
                1330             do i=1,m
                1331              tcom(i,j,ik)=1.0
                1332             enddo
                1333            enddo
                1334           enddo
                1335        endif
                1336 
                1337 c-----for cfc-11 transmission using k-distribution method.
                1338 
                1339        if (f11bnd) then
                1340            do j=1,n
                1341             do i=1,m
                1342              tf11(i,j)=1.0
                1343             enddo
                1344            enddo
                1345        endif
                1346 
                1347 c-----for cfc-12 transmission using k-distribution method.
                1348 
                1349        if (f12bnd) then
                1350            do j=1,n
                1351             do i=1,m
                1352              tf12(i,j)=1.0
                1353             enddo
                1354            enddo
                1355        endif
                1356 
                1357 c-----for cfc-22 transmission when using k-distribution method.
                1358 
                1359        if (f22bnd) then
                1360            do j=1,n
                1361             do i=1,m
                1362              tf22(i,j)=1.0
                1363             enddo
                1364            enddo
                1365        endif
                1366 
                1367 c-----for the transmission in band 10 using k-distribution method.
                1368 
                1369        if (b10bnd) then
                1370           do ik=1,6
                1371            do j=1,n
                1372             do i=1,m
                1373               th2o(i,j,ik)=1.0
                1374               tco2(i,j,ik,1)=1.0
                1375             enddo
                1376            enddo
                1377           enddo
                1378           do iq=1,3
                1379            do j=1,n
                1380             do i=1,m
                1381               tcon(i,j,iq)=1.0
                1382             enddo
                1383            enddo
                1384           enddo
                1385           do ik=1,4
                1386            do j=1,n
                1387             do i=1,m
                1388               tn2o(i,j,ik)=1.0
                1389             enddo
                1390            enddo
                1391           enddo
                1392        endif
                1393 
                1394       endif
                1395 
                1396 c-----loop over the bottom level of the region (k2)
                1397 
                1398       do 3000 k2=k1+1,np+1
                1399 
                1400 c-----initialize total transmission function (trant)
                1401 
                1402           do j=1,n
                1403            do i=1,m
                1404             trant(i,j)=1.0
                1405            enddo
                1406           enddo
                1407 
                1408       if (h2otbl) then
                1409           w1=-8.0
                1410           p1=-2.0
                1411           dwe=0.3
                1412           dpe=0.2
                1413 
                1414 c-----compute water vapor transmittance using table look-up
                1415 
                1416           if (ib.eq.1) then
                1417            call tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
                1418      *                 w1,p1,dwe,dpe,h11,h12,h13,trant)
                1419           endif
                1420           if (ib.eq.2) then
                1421            call tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
                1422      *                 w1,p1,dwe,dpe,h21,h22,h23,trant)
                1423 
                1424           endif
                1425           if (ib.eq.8) then
                1426            call tablup(k1,k2,m,n,np,nx,nh,nt,sh2o,swpre,swtem,
                1427      *                 w1,p1,dwe,dpe,h81,h82,h83,trant)
                1428           endif
                1429 
                1430       else
                1431 
                1432 c-----compute water vapor transmittance using k-distribution
                1433 
                1434        if (.not.b10bnd) then
                1435         call h2okdis(ib,m,n,np,k2-1,fkw,gkw,ne,h2oexp,conexp,
                1436      *               th2o,tcon,trant)
                1437        endif
                1438 
                1439       endif
                1440 
                1441       if (co2bnd) then
                1442 
                1443         if (high) then
                1444 
                1445 c-----compute co2 transmittance using table look-up method
                1446 
                1447           w1=-4.0
                1448           p1=-2.0
                1449           dwe=0.3
                1450           dpe=0.2
                1451           call tablup(k1,k2,m,n,np,nx,nc,nt,sco3,scopre,scotem,
                1452      *                w1,p1,dwe,dpe,c1,c2,c3,trant)
046c3794c1 Andr*1453 
1662f365b2 Andr*1454        else
                1455 
                1456 c-----compute co2 transmittance using k-distribution method
                1457           call co2kdis(m,n,np,k2-1,co2exp,tco2,trant)
                1458 
                1459         endif
                1460 
                1461       endif
                1462 
                1463 c-----All use table look-up to compute o3 transmittance.
                1464 
                1465       if (oznbnd) then
                1466           w1=-6.0
                1467           p1=-2.0
                1468           dwe=0.3
                1469           dpe=0.2
                1470           call tablup(k1,k2,m,n,np,nx,no,nt,sco3,scopre,scotem,
                1471      *                w1,p1,dwe,dpe,o1,o2,o3,trant)
154b1b9078 Andr*1472 
1662f365b2 Andr*1473       endif
                1474 
                1475 c***** for trace gases *****
                1476 
                1477       if (trace) then
                1478 
                1479 c-----compute n2o transmittance using k-distribution method
                1480 
                1481        if (n2obnd) then
                1482           call n2okdis(ib,m,n,np,k2-1,n2oexp,tn2o,trant)
                1483        endif
                1484 
                1485 c-----compute ch4 transmittance using k-distribution method
                1486 
                1487        if (ch4bnd) then
                1488           call ch4kdis(ib,m,n,np,k2-1,ch4exp,tch4,trant)
                1489        endif
                1490 
                1491 c-----compute co2-minor transmittance using k-distribution method
                1492 
                1493        if (combnd) then
                1494           call comkdis(ib,m,n,np,k2-1,comexp,tcom,trant)
                1495        endif
                1496 
                1497 c-----compute cfc11 transmittance using k-distribution method
                1498 
                1499        if (f11bnd) then
                1500           call cfckdis(m,n,np,k2-1,f11exp,tf11,trant)
                1501        endif
                1502 
                1503 c-----compute cfc12 transmittance using k-distribution method
                1504 
                1505        if (f12bnd) then
                1506           call cfckdis(m,n,np,k2-1,f12exp,tf12,trant)
                1507        endif
                1508 
                1509 c-----compute cfc22 transmittance using k-distribution method
                1510 
                1511        if (f22bnd) then
                1512           call cfckdis(m,n,np,k2-1,f22exp,tf22,trant)
                1513        endif
                1514 
                1515 c-----compute transmittance in band 10 using k-distribution method
                1516 c     here, trant is the change in transmittance due to n2o absorption
                1517 
                1518        if (b10bnd) then
                1519           call b10kdis(m,n,np,k2-1,h2oexp,conexp,co2exp,n2oexp
                1520      *                ,th2o,tcon,tco2,tn2o,trant)
                1521        endif
                1522 
                1523       endif
                1524 c
                1525 c-----include aerosol effect
                1526 c
                1527       do j=1,n
                1528        do i=1,m
                1529          ff=0.5+(0.3739+(0.0076+0.1185*asyal(i,j,k2-1,ib))
                1530      *      *asyal(i,j,k2-1,ib))*asyal(i,j,k2-1,ib)
                1531          taux=taual(i,j,k2-1,ib)*(1.-ssaal(i,j,k2-1,ib)*ff)
                1532          tranal(i,j)=tranal(i,j)*exp(-1.66*taux)
                1533          trant (i,j)=trant(i,j) *tranal(i,j)
                1534        enddo
                1535       enddo
                1536 
                1537 c-----Identify the clear-sky fractions clrhi, clrmd and clrlw, in the
                1538 c     high, middle and low cloud groups.
                1539 c     fclr is the clear-sky fraction between levels k1 and k2 assuming
                1540 c     the three cloud groups are randomly overlapped.
                1541 
                1542       do j=1,n
                1543        do i=1,m
                1544         if( k2 .le. ict ) then
                1545             clrhi(i,j)=min(clr(i,j,k2-1),clrhi(i,j))
                1546         elseif( k2 .gt. ict .and. k2 .le. icb ) then
                1547             clrmd(i,j)=min(clr(i,j,k2-1),clrmd(i,j))
                1548         elseif( k2 .gt. icb ) then
                1549             clrlw(i,j)=min(clr(i,j,k2-1),clrlw(i,j))
                1550         endif
                1551             fclr(i,j)=clrlw(i,j)*clrmd(i,j)*clrhi(i,j)
                1552 
                1553        enddo
                1554       enddo
                1555 
                1556 c-----compute upward and downward fluxes (band 1-9).
                1557 c     add "boundary" terms to the net downward flux.
                1558 c     these are the first terms on the right-hand-side of
                1559 c     eqs. (56a) and (56b).  downward fluxes are positive.
                1560 
                1561       if (.not. b10bnd) then
                1562 
                1563        if (k2 .eq. k1+1) then
                1564 
                1565         do j=1,n
                1566          do i=1,m
                1567 
                1568 c-----compute upward and downward fluxes for clear-sky situation
                1569 
                1570           flcu(i,j,k1)=flcu(i,j,k1)-blayer(i,j,k1)
                1571           flcd(i,j,k2)=flcd(i,j,k2)+blayer(i,j,k1)
                1572 
                1573 c-----compute upward and downward fluxes for all-sky situation
                1574 
                1575           flxu(i,j,k1)=flxu(i,j,k1)-blayer(i,j,k1)
                1576           flxd(i,j,k2)=flxd(i,j,k2)+blayer(i,j,k1)
                1577 
                1578          enddo
                1579         enddo
                1580 
                1581        endif
                1582 
                1583 c-----add flux components involving the four layers above and below
                1584 c     the levels k1 and k2.  it follows eqs. (56a) and (56b).
                1585 
                1586        do j=1,n
                1587         do i=1,m
                1588          xx=trant(i,j)*dlayer(i,j,k2)
                1589          flcu(i,j,k1) =flcu(i,j,k1)+xx
                1590          flxu(i,j,k1) =flxu(i,j,k1)+xx*fclr(i,j)
                1591          xx=trant(i,j)*dlayer(i,j,k1)
                1592          flcd(i,j,k2) =flcd(i,j,k2)+xx
                1593          flxd(i,j,k2) =flxd(i,j,k2)+xx*fclr(i,j)
                1594         enddo
                1595        enddo
                1596 
                1597       else
                1598 
                1599 c-----flux reduction due to n2o in band 10
                1600 
                1601        if (trace) then
                1602 
                1603         do j=1,n
                1604          do i=1,m
                1605            rflx(i,j,k1) = rflx(i,j,k1)+trant(i,j)*fclr(i,j)
                1606      *                   *dlayer(i,j,k2)
                1607            rflx(i,j,k2) = rflx(i,j,k2)+trant(i,j)*fclr(i,j)
                1608      *                   *dlayer(i,j,k1)
                1609            rflc(i,j,k1) = rflc(i,j,k1)+trant(i,j)
                1610      *                   *dlayer(i,j,k2)
                1611            rflc(i,j,k2) = rflc(i,j,k2)+trant(i,j)
                1612      *                   *dlayer(i,j,k1)
                1613          enddo
                1614         enddo
                1615 
                1616        endif
                1617 
                1618       endif
                1619 
                1620  3000 continue
                1621 
                1622 
                1623 c-----transmission between level k1 and the surface
                1624 
                1625          do j=1,n
                1626           do i=1,m
                1627            trantcr(i,j,k1) =trant(i,j)
                1628            transfc(i,j,k1) =trant(i,j)*fclr(i,j)
                1629           enddo
                1630          enddo
                1631 
                1632 c-----compute the partial derivative of fluxes with respect to
                1633 c     surface temperature (eq. 59)
                1634 
                1635       if (trace .or. (.not. b10bnd) ) then
                1636 
                1637        do j=1,n
                1638         do i=1,m
                1639          dfdts(i,j,k1) =dfdts(i,j,k1)-dbs(i,j)*transfc(i,j,k1)
                1640         enddo
                1641        enddo
                1642 
                1643       endif
                1644 
                1645  2000 continue
                1646 
                1647       if (.not. b10bnd) then
                1648 
                1649 c-----add contribution from the surface to the flux terms at the surface
                1650 
                1651         do j=1,n
                1652          do i=1,m
                1653           flcu(i,j,np+1)=flcu(i,j,np+1)-blayer(i,j,np+1)
                1654           flxu(i,j,np+1)=flxu(i,j,np+1)-blayer(i,j,np+1)
                1655           st4(i,j)=st4(i,j)-blayer(i,j,np+1)
                1656           dfdts(i,j,np+1)=dfdts(i,j,np+1)-dbs(i,j)
                1657          enddo
                1658         enddo
                1659 
                1660 c-----add the flux component reflected by the surface
                1661 c     note: upward flux is negative
                1662 
                1663         do k=1, np+1
                1664          do j=1,n
                1665           do i=1,m
                1666            flcu(i,j,k)=flcu(i,j,k)-
                1667      *           flcd(i,j,np+1)*trantcr(i,j,k)*(1.-emiss(i,j,ib))
                1668            flxu(i,j,k)=flxu(i,j,k)-
                1669      *           flxd(i,j,np+1)*transfc(i,j,k)*(1.-emiss(i,j,ib))
                1670           enddo
                1671          enddo
                1672         enddo
                1673 
                1674       endif
                1675 
                1676 c-----summation of fluxes over spectral bands
                1677 
                1678       do k=1,np+1
                1679        do j=1,n
                1680         do i=1,m
                1681          flc(i,j,k)=flc(i,j,k)+flcd(i,j,k)+flcu(i,j,k)
                1682          flx(i,j,k)=flx(i,j,k)+flxd(i,j,k)+flxu(i,j,k)
                1683         enddo
                1684        enddo
                1685       enddo
                1686 
                1687  1000 continue
                1688 
                1689 c-----adjust fluxes due to n2o absorption in band 10
                1690 
                1691       do k=1,np+1
                1692        do j=1,n
                1693         do i=1,m
                1694          flc(i,j,k)=flc(i,j,k)+rflc(i,j,k)
                1695          flx(i,j,k)=flx(i,j,k)+rflx(i,j,k)
                1696         enddo
                1697        enddo
                1698       enddo
                1699 
                1700       return
                1701       end
                1702 c***********************************************************************
                1703       subroutine column (m,n,np,pa,dt,sabs0,sabs,spre,stem)
                1704 c***********************************************************************
                1705 c-----compute column-integrated (from top of the model atmosphere)
                1706 c     absorber amount (sabs), absorber-weighted pressure (spre) and
                1707 c     temperature (stem).
                1708 c     computations of spre and stem follows eqs. (37) and (38).
                1709 c
                1710 c--- input parameters
                1711 c   number of soundings in zonal direction (m)
                1712 c   number of soundings in meridional direction (n)
                1713 c   number of atmospheric layers (np)
                1714 c   layer pressure (pa)
                1715 c   layer temperature minus 250K (dt)
                1716 c   layer absorber amount (sabs0)
                1717 c
                1718 c--- output parameters
                1719 c   column-integrated absorber amount (sabs)
                1720 c   column absorber-weighted pressure (spre)
                1721 c   column absorber-weighted temperature (stem)
                1722 c
                1723 c--- units of pa and dt are mb and k, respectively.
                1724 c    units of sabs are g/cm**2 for water vapor and (cm-atm)stp
                1725 c    for co2 and o3
                1726 c***********************************************************************
                1727       implicit none
                1728       integer m,n,np,i,j,k
                1729 
                1730 c---- input parameters -----
                1731 
a456aa407c Andr*1732       _RL pa(m,n,np),dt(m,n,np),sabs0(m,n,np)
1662f365b2 Andr*1733 
                1734 c---- output parameters -----
                1735 
a456aa407c Andr*1736       _RL sabs(m,n,np+1),spre(m,n,np+1),stem(m,n,np+1)
1662f365b2 Andr*1737 
                1738 c*********************************************************************
                1739         do j=1,n
                1740          do i=1,m
                1741           sabs(i,j,1)=0.0
                1742           spre(i,j,1)=0.0
                1743           stem(i,j,1)=0.0
                1744          enddo
                1745         enddo
                1746 
                1747         do k=1,np
                1748          do j=1,n
                1749           do i=1,m
                1750            sabs(i,j,k+1)=sabs(i,j,k)+sabs0(i,j,k)
                1751            spre(i,j,k+1)=spre(i,j,k)+pa(i,j,k)*sabs0(i,j,k)
                1752            stem(i,j,k+1)=stem(i,j,k)+dt(i,j,k)*sabs0(i,j,k)
                1753           enddo
                1754          enddo
                1755         enddo
                1756 
                1757        return
                1758        end
                1759 c**********************************************************************
                1760       subroutine h2oexps(ib,m,n,np,dh2o,pa,dt,xkw,aw,bw,pm,mw,h2oexp)
                1761 c**********************************************************************
                1762 c   compute exponentials for water vapor line absorption
                1763 c   in individual layers.
                1764 c
                1765 c---- input parameters
                1766 c  spectral band (ib)
                1767 c  number of grid intervals in zonal direction (m)
                1768 c  number of grid intervals in meridional direction (n)
                1769 c  number of layers (np)
                1770 c  layer water vapor amount for line absorption (dh2o) 
                1771 c  layer pressure (pa)
                1772 c  layer temperature minus 250K (dt)
                1773 c  absorption coefficients for the first k-distribution
                1774 c     function due to h2o line absorption (xkw)
                1775 c  coefficients for the temperature and pressure scaling (aw,bw,pm)
                1776 c  ratios between neighboring absorption coefficients for
                1777 c     h2o line absorption (mw)
                1778 c
                1779 c---- output parameters
                1780 c  6 exponentials for each layer  (h2oexp)
                1781 c**********************************************************************
                1782       implicit none
                1783       integer ib,m,n,np,i,j,k,ik
                1784 
                1785 c---- input parameters ------
                1786 
a456aa407c Andr*1787       _RL dh2o(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*1788 
                1789 c---- output parameters -----
                1790 
a456aa407c Andr*1791       _RL h2oexp(m,n,np,6)
1662f365b2 Andr*1792 
                1793 c---- static data -----
                1794 
                1795       integer mw(9)
a456aa407c Andr*1796       _RL xkw(9),aw(9),bw(9),pm(9)
1662f365b2 Andr*1797 
                1798 c---- temporary arrays -----
                1799 
a456aa407c Andr*1800       _RL xh,xh1
1662f365b2 Andr*1801 
                1802 c**********************************************************************
                1803 c    note that the 3 sub-bands in band 3 use the same set of xkw, aw,
                1804 c    and bw,  therefore, h2oexp for these sub-bands are identical.
                1805 c**********************************************************************
                1806 
                1807         do k=1,np
                1808          do j=1,n
                1809           do i=1,m
                1810 
                1811 c-----xh is the scaled water vapor amount for line absorption
                1812 c     computed from (27)
                1813 
                1814            xh = dh2o(i,j,k)*(pa(i,j,k)/500.)**pm(ib)
                1815      1        * ( 1.+(aw(ib)+bw(ib)* dt(i,j,k))*dt(i,j,k) )
                1816 
                1817 c-----h2oexp is the water vapor transmittance of the layer k
                1818 c     due to line absorption
                1819 
                1820            h2oexp(i,j,k,1) = exp(-xh*xkw(ib))
                1821 
                1822           enddo
                1823          enddo
                1824         enddo
                1825 
                1826         do ik=2,6
                1827 
                1828          if (mw(ib).eq.6) then
                1829 
                1830           do k=1,np
                1831            do j=1,n
                1832             do i=1,m
                1833              xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
                1834              h2oexp(i,j,k,ik) = xh*xh*xh
                1835             enddo
                1836            enddo
                1837           enddo
                1838 
                1839         elseif (mw(ib).eq.8) then
                1840 
                1841           do k=1,np
                1842            do j=1,n
                1843             do i=1,m
                1844              xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
                1845              xh = xh*xh
                1846              h2oexp(i,j,k,ik) = xh*xh
                1847             enddo
                1848            enddo
                1849           enddo
                1850 
                1851         elseif (mw(ib).eq.9) then
                1852 
                1853           do k=1,np
                1854            do j=1,n
                1855             do i=1,m
                1856              xh=h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
                1857              xh1 = xh*xh
                1858              h2oexp(i,j,k,ik) = xh*xh1
                1859             enddo
                1860            enddo
                1861           enddo
                1862 
                1863         else
                1864 
                1865           do k=1,np
                1866            do j=1,n
                1867             do i=1,m
                1868              xh = h2oexp(i,j,k,ik-1)*h2oexp(i,j,k,ik-1)
                1869              xh = xh*xh
                1870              xh = xh*xh
                1871              h2oexp(i,j,k,ik) = xh*xh
                1872             enddo
                1873            enddo
                1874           enddo
                1875 
                1876         endif
                1877        enddo
                1878 
                1879       return
                1880       end
                1881 c**********************************************************************
                1882       subroutine conexps(ib,m,n,np,dcont,xke,conexp)
                1883 c**********************************************************************
                1884 c   compute exponentials for continuum absorption in individual layers.
                1885 c
                1886 c---- input parameters
                1887 c  spectral band (ib)
                1888 c  number of grid intervals in zonal direction (m)
                1889 c  number of grid intervals in meridional direction (n)
                1890 c  number of layers (np)
                1891 c  layer scaled water vapor amount for continuum absorption (dcont) 
                1892 c  absorption coefficients for the first k-distribution function
                1893 c     due to water vapor continuum absorption (xke)
                1894 c
                1895 c---- output parameters
                1896 c  1 or 3 exponentials for each layer (conexp)
                1897 c**********************************************************************
                1898       implicit none
                1899       integer ib,m,n,np,i,j,k,iq
                1900 
                1901 c---- input parameters ------
                1902 
a456aa407c Andr*1903       _RL dcont(m,n,np)
1662f365b2 Andr*1904 
                1905 c---- updated parameters -----
                1906 
a456aa407c Andr*1907       _RL conexp(m,n,np,3)
1662f365b2 Andr*1908 
                1909 c---- static data -----
                1910 
a456aa407c Andr*1911       _RL xke(9)
1662f365b2 Andr*1912 
                1913 c**********************************************************************
                1914 
                1915         do k=1,np
                1916          do j=1,n
                1917           do i=1,m
                1918            conexp(i,j,k,1) = exp(-dcont(i,j,k)*xke(ib))
                1919           enddo
                1920          enddo
                1921         enddo
                1922 
                1923        if (ib .eq. 3) then
                1924 
                1925 c-----the absorption coefficients for sub-bands 3b (iq=2) and 3a (iq=3)
                1926 c     are, respectively, two and three times the absorption coefficient
                1927 c     for sub-band 3c (iq=1) (table 6).
                1928 
                1929         do iq=2,3
                1930          do k=1,np
                1931           do j=1,n
                1932            do i=1,m
                1933             conexp(i,j,k,iq) = conexp(i,j,k,iq-1) *conexp(i,j,k,iq-1)
                1934            enddo
                1935           enddo
                1936          enddo
                1937         enddo
                1938 
                1939        endif
                1940 
                1941       return
                1942       end
                1943 c**********************************************************************
                1944       subroutine co2exps(m,n,np,dco2,pa,dt,co2exp)
                1945 c**********************************************************************
                1946 c   compute co2 exponentials for individual layers.
                1947 c
                1948 c---- input parameters
                1949 c  number of grid intervals in zonal direction (m)
                1950 c  number of grid intervals in meridional direction (n)
                1951 c  number of layers (np)
                1952 c  layer co2 amount (dco2)
                1953 c  layer pressure (pa)
                1954 c  layer temperature minus 250K (dt)
                1955 c
                1956 c---- output parameters
                1957 c  6 exponentials for each layer (co2exp)
                1958 c**********************************************************************
                1959       implicit none
                1960       integer m,n,np,i,j,k
                1961 
                1962 c---- input parameters -----
                1963 
a456aa407c Andr*1964       _RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*1965 
                1966 c---- output parameters -----
                1967 
a456aa407c Andr*1968       _RL co2exp(m,n,np,6,2)
1662f365b2 Andr*1969 
                1970 c---- temporary arrays -----
                1971 
a456aa407c Andr*1972       _RL xc
1662f365b2 Andr*1973 
                1974 c**********************************************************************
                1975 
                1976         do k=1,np
                1977          do j=1,n
                1978           do i=1,m
                1979 
                1980 c-----compute the scaled co2 amount from eq. (27) for band-wings
                1981 c     (sub-bands 3a and 3c).
                1982 
                1983            xc = dco2(i,j,k)*(pa(i,j,k)/300.0)**0.5
                1984      1             *(1.+(0.0182+1.07e-4*dt(i,j,k))*dt(i,j,k))
                1985 
                1986 c-----six exponentials by powers of 8 (table 7).
                1987 
                1988            co2exp(i,j,k,1,1)=exp(-xc*2.656e-5)
                1989 
                1990            xc=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)
                1991            xc=xc*xc
                1992            co2exp(i,j,k,2,1)=xc*xc
                1993 
                1994            xc=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)
                1995            xc=xc*xc
                1996            co2exp(i,j,k,3,1)=xc*xc
                1997 
                1998            xc=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)
                1999            xc=xc*xc
                2000            co2exp(i,j,k,4,1)=xc*xc
                2001 
                2002            xc=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)
                2003            xc=xc*xc
                2004            co2exp(i,j,k,5,1)=xc*xc
                2005 
                2006            xc=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)
                2007            xc=xc*xc
                2008            co2exp(i,j,k,6,1)=xc*xc
                2009 
                2010 c-----compute the scaled co2 amount from eq. (27) for band-center
                2011 c     region (sub-band 3b).
                2012 
                2013            xc = dco2(i,j,k)*(pa(i,j,k)/30.0)**0.85
                2014      1             *(1.+(0.0042+2.00e-5*dt(i,j,k))*dt(i,j,k))
                2015 
                2016            co2exp(i,j,k,1,2)=exp(-xc*2.656e-3)
                2017 
                2018            xc=co2exp(i,j,k,1,2)*co2exp(i,j,k,1,2)
                2019            xc=xc*xc
                2020            co2exp(i,j,k,2,2)=xc*xc
                2021 
                2022            xc=co2exp(i,j,k,2,2)*co2exp(i,j,k,2,2)
                2023            xc=xc*xc
                2024            co2exp(i,j,k,3,2)=xc*xc
                2025 
                2026            xc=co2exp(i,j,k,3,2)*co2exp(i,j,k,3,2)
                2027            xc=xc*xc
                2028            co2exp(i,j,k,4,2)=xc*xc
                2029 
                2030            xc=co2exp(i,j,k,4,2)*co2exp(i,j,k,4,2)
                2031            xc=xc*xc
                2032            co2exp(i,j,k,5,2)=xc*xc
                2033 
                2034            xc=co2exp(i,j,k,5,2)*co2exp(i,j,k,5,2)
                2035            xc=xc*xc
                2036            co2exp(i,j,k,6,2)=xc*xc
                2037 
                2038           enddo
                2039          enddo
                2040         enddo
                2041 
                2042       return
                2043       end
                2044 c**********************************************************************
                2045       subroutine n2oexps(ib,m,n,np,dn2o,pa,dt,n2oexp)
                2046 c**********************************************************************
                2047 c   compute n2o exponentials for individual layers.
                2048 c
                2049 c---- input parameters
                2050 c  spectral band (ib)
                2051 c  number of grid intervals in zonal direction (m)
                2052 c  number of grid intervals in meridional direction (n)
                2053 c  number of layers (np)
                2054 c  layer n2o amount (dn2o)
                2055 c  layer pressure (pa)
                2056 c  layer temperature minus 250K (dt)
                2057 c
                2058 c---- output parameters
                2059 c  2 or 4 exponentials for each layer (n2oexp)
                2060 c**********************************************************************
                2061       implicit none
                2062       integer ib,m,n,np,i,j,k
                2063 
                2064 c---- input parameters -----
                2065 
a456aa407c Andr*2066       _RL dn2o(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*2067 
                2068 c---- output parameters -----
                2069 
a456aa407c Andr*2070       _RL n2oexp(m,n,np,4)
1662f365b2 Andr*2071 
                2072 c---- temporary arrays -----
                2073 
a456aa407c Andr*2074       _RL xc,xc1,xc2
1662f365b2 Andr*2075 
                2076 c**********************************************************************
                2077 
                2078        do k=1,np
                2079         do j=1,n
                2080          do i=1,m
                2081 
                2082 c-----four exponential by powers of 21 for band 6
                2083 
                2084           if (ib.eq.6) then
                2085 
                2086            xc=dn2o(i,j,k)*(1.+(1.9297e-3+4.3750e-6*dt(i,j,k))*dt(i,j,k))
                2087            n2oexp(i,j,k,1)=exp(-xc*6.31582e-2)
                2088 
                2089            xc=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
                2090            xc1=xc*xc
                2091            xc2=xc1*xc1
                2092            n2oexp(i,j,k,2)=xc*xc1*xc2
                2093 
                2094 c-----four exponential by powers of 8 for band 7
                2095 
                2096           else
                2097 
                2098            xc=dn2o(i,j,k)*(pa(i,j,k)/500.0)**0.48
                2099      *        *(1.+(1.3804e-3+7.4838e-6*dt(i,j,k))*dt(i,j,k))
                2100            n2oexp(i,j,k,1)=exp(-xc*5.35779e-2)
                2101 
                2102            xc=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
                2103            xc=xc*xc
                2104            n2oexp(i,j,k,2)=xc*xc
                2105            xc=n2oexp(i,j,k,2)*n2oexp(i,j,k,2)
                2106            xc=xc*xc
                2107            n2oexp(i,j,k,3)=xc*xc
                2108            xc=n2oexp(i,j,k,3)*n2oexp(i,j,k,3)
                2109            xc=xc*xc
                2110            n2oexp(i,j,k,4)=xc*xc
                2111 
                2112           endif
                2113 
                2114          enddo
                2115         enddo
                2116        enddo
                2117 
                2118       return
                2119       end
                2120 c**********************************************************************
                2121       subroutine ch4exps(ib,m,n,np,dch4,pa,dt,ch4exp)
                2122 c**********************************************************************
                2123 c   compute ch4 exponentials for individual layers.
                2124 c
                2125 c---- input parameters
                2126 c  spectral band (ib)
                2127 c  number of grid intervals in zonal direction (m)
                2128 c  number of grid intervals in meridional direction (n)
                2129 c  number of layers (np)
                2130 c  layer ch4 amount (dch4)
                2131 c  layer pressure (pa)
                2132 c  layer temperature minus 250K (dt)
                2133 c
                2134 c---- output parameters
                2135 c  1 or 4 exponentials for each layer (ch4exp)
                2136 c**********************************************************************
                2137       implicit none
                2138       integer ib,m,n,np,i,j,k
                2139 
                2140 c---- input parameters -----
                2141 
a456aa407c Andr*2142       _RL dch4(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*2143 
                2144 c---- output parameters -----
                2145 
a456aa407c Andr*2146       _RL ch4exp(m,n,np,4)
1662f365b2 Andr*2147 
                2148 c---- temporary arrays -----
                2149 
a456aa407c Andr*2150       _RL xc
1662f365b2 Andr*2151 
                2152 c**********************************************************************
                2153 
                2154        do k=1,np
                2155         do j=1,n
                2156          do i=1,m
                2157 
                2158 c-----four exponentials for band 6
                2159 
                2160           if (ib.eq.6) then
                2161 
                2162            xc=dch4(i,j,k)*(1.+(1.7007e-2+1.5826e-4*dt(i,j,k))*dt(i,j,k))
                2163            ch4exp(i,j,k,1)=exp(-xc*5.80708e-3)
                2164 
                2165 c-----four exponentials by powers of 12 for band 7
                2166 
                2167           else
                2168 
                2169            xc=dch4(i,j,k)*(pa(i,j,k)/500.0)**0.65
                2170      *       *(1.+(5.9590e-4-2.2931e-6*dt(i,j,k))*dt(i,j,k))
                2171            ch4exp(i,j,k,1)=exp(-xc*6.29247e-2)
                2172 
                2173            xc=ch4exp(i,j,k,1)*ch4exp(i,j,k,1)*ch4exp(i,j,k,1)
                2174            xc=xc*xc
                2175            ch4exp(i,j,k,2)=xc*xc
                2176 
                2177            xc=ch4exp(i,j,k,2)*ch4exp(i,j,k,2)*ch4exp(i,j,k,2)
                2178            xc=xc*xc
                2179            ch4exp(i,j,k,3)=xc*xc
                2180 
                2181            xc=ch4exp(i,j,k,3)*ch4exp(i,j,k,3)*ch4exp(i,j,k,3)
                2182            xc=xc*xc
                2183            ch4exp(i,j,k,4)=xc*xc
                2184 
                2185           endif
                2186 
                2187          enddo
                2188         enddo
                2189        enddo
                2190 
                2191       return
                2192       end
                2193 c**********************************************************************
                2194       subroutine comexps(ib,m,n,np,dcom,dt,comexp)
                2195 c**********************************************************************
                2196 c   compute co2-minor exponentials for individual layers.
                2197 c
                2198 c---- input parameters
                2199 c  spectral band (ib)
                2200 c  number of grid intervals in zonal direction (m)
                2201 c  number of grid intervals in meridional direction (n)
                2202 c  number of layers (np)
                2203 c  layer co2 amount (dcom)
                2204 c  layer temperature minus 250K (dt)
                2205 c
                2206 c---- output parameters
                2207 c  2 exponentials for each layer (comexp)
                2208 c**********************************************************************
                2209       implicit none
                2210       integer ib,m,n,np,i,j,k
                2211 
                2212 c---- input parameters -----
                2213 
a456aa407c Andr*2214       _RL dcom(m,n,np),dt(m,n,np)
1662f365b2 Andr*2215 
                2216 c---- output parameters -----
                2217 
a456aa407c Andr*2218       _RL comexp(m,n,np,2)
1662f365b2 Andr*2219 
                2220 c---- temporary arrays -----
                2221 
a456aa407c Andr*2222       _RL xc,xc1,xc2
1662f365b2 Andr*2223 
                2224 c**********************************************************************
                2225 
                2226        do k=1,np
                2227         do j=1,n
                2228          do i=1,m
                2229 
                2230 c-----two exponentials by powers of 60 for band 4
                2231 
                2232           if (ib.eq.4) then
                2233 
                2234            xc=dcom(i,j,k)*(1.+(3.5775e-2+4.0447e-4*dt(i,j,k))*dt(i,j,k))
                2235            comexp(i,j,k,1)=exp(-xc*2.18947e-5)
                2236 
                2237            xc=comexp(i,j,k,1)*comexp(i,j,k,1)*comexp(i,j,k,1)
                2238            xc=xc*xc
                2239            xc1=xc*xc
                2240            xc=xc1*xc1
                2241            xc=xc*xc
                2242            comexp(i,j,k,2)=xc*xc1
                2243 
                2244 c-----two exponentials by powers of 44 for band 5
                2245 
                2246           else
                2247 
                2248            xc=dcom(i,j,k)*(1.+(3.4268e-2+3.7401e-4*dt(i,j,k))*dt(i,j,k))
                2249            comexp(i,j,k,1)=exp(-xc*4.74570e-5)
                2250 
                2251            xc=comexp(i,j,k,1)*comexp(i,j,k,1)
                2252            xc1=xc*xc
                2253            xc2=xc1*xc1
                2254            xc=xc2*xc2
                2255            xc=xc*xc
                2256            comexp(i,j,k,2)=xc1*xc2*xc
                2257 
                2258           endif
                2259 
                2260          enddo
                2261         enddo
                2262        enddo
                2263 
                2264       return
                2265       end
                2266 c**********************************************************************
                2267       subroutine cfcexps(ib,m,n,np,a1,b1,fk1,a2,b2,fk2,dcfc,dt,cfcexp)
                2268 c**********************************************************************
                2269 c   compute cfc(-11, -12, -22) exponentials for individual layers.
                2270 c
                2271 c---- input parameters
                2272 c  spectral band (ib)
                2273 c  number of grid intervals in zonal direction (m)
                2274 c  number of grid intervals in meridional direction (n)
                2275 c  number of layers (np)
                2276 c  parameters for computing the scaled cfc amounts
                2277 c             for temperature scaling (a1,b1,a2,b2)
                2278 c  the absorption coefficients for the
                2279 c     first k-distribution function due to cfcs (fk1,fk2)
                2280 c  layer cfc amounts (dcfc)
                2281 c  layer temperature minus 250K (dt)
                2282 c
                2283 c---- output parameters
                2284 c  1 exponential for each layer (cfcexp)
                2285 c**********************************************************************
                2286       implicit none
                2287       integer ib,m,n,np,i,j,k
                2288 
                2289 c---- input parameters -----
                2290 
a456aa407c Andr*2291       _RL dcfc(m,n,np),dt(m,n,np)
1662f365b2 Andr*2292 
                2293 c---- output parameters -----
                2294 
a456aa407c Andr*2295       _RL cfcexp(m,n,np)
1662f365b2 Andr*2296 
                2297 c---- static data -----
                2298 
a456aa407c Andr*2299       _RL a1,b1,fk1,a2,b2,fk2
1662f365b2 Andr*2300 
                2301 c---- temporary arrays -----
                2302 
a456aa407c Andr*2303       _RL xf
1662f365b2 Andr*2304 
                2305 c**********************************************************************
                2306 
                2307        do k=1,np
                2308         do j=1,n
                2309          do i=1,m
                2310 
                2311 c-----compute the scaled cfc amount (xf) and exponential (cfcexp)
                2312 
                2313           if (ib.eq.4) then
                2314            xf=dcfc(i,j,k)*(1.+(a1+b1*dt(i,j,k))*dt(i,j,k))
                2315            cfcexp(i,j,k)=exp(-xf*fk1)
                2316           else
                2317            xf=dcfc(i,j,k)*(1.+(a2+b2*dt(i,j,k))*dt(i,j,k))
                2318            cfcexp(i,j,k)=exp(-xf*fk2)
                2319           endif
                2320 
                2321          enddo
                2322         enddo
                2323        enddo
                2324 
                2325       return
                2326       end
                2327 c**********************************************************************
                2328       subroutine b10exps(m,n,np,dh2o,dcont,dco2,dn2o,pa,dt
                2329      *          ,h2oexp,conexp,co2exp,n2oexp)
                2330 c**********************************************************************
                2331 c   compute band3a exponentials for individual layers.
                2332 c
                2333 c---- input parameters
                2334 c  number of grid intervals in zonal direction (m)
                2335 c  number of grid intervals in meridional direction (n)
                2336 c  number of layers (np)
                2337 c  layer h2o amount for line absorption (dh2o)
                2338 c  layer h2o amount for continuum absorption (dcont)
                2339 c  layer co2 amount (dco2)
                2340 c  layer n2o amount (dn2o)
                2341 c  layer pressure (pa)
                2342 c  layer temperature minus 250K (dt)
                2343 c
                2344 c---- output parameters
                2345 c
                2346 c  exponentials for each layer (h2oexp,conexp,co2exp,n2oexp)
                2347 c**********************************************************************
                2348       implicit none
                2349       integer m,n,np,i,j,k
                2350 
                2351 c---- input parameters -----
                2352 
a456aa407c Andr*2353       _RL dh2o(m,n,np),dcont(m,n,np),dn2o(m,n,np)
                2354       _RL dco2(m,n,np),pa(m,n,np),dt(m,n,np)
1662f365b2 Andr*2355 
                2356 c---- output parameters -----
                2357 
a456aa407c Andr*2358       _RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
1662f365b2 Andr*2359      *    ,n2oexp(m,n,np,4)
                2360 
                2361 c---- temporary arrays -----
                2362 
a456aa407c Andr*2363       _RL xx,xx1,xx2,xx3
1662f365b2 Andr*2364 
                2365 c**********************************************************************
                2366 
                2367         do k=1,np
                2368          do j=1,n
                2369           do i=1,m
                2370 
                2371 c-----compute the scaled h2o-line amount for the sub-band 3a
                2372 c     table 3, Chou et al. (JAS, 1993)
                2373 
                2374            xx=dh2o(i,j,k)*(pa(i,j,k)/500.0)
                2375      1           *(1.+(0.0149+6.20e-5*dt(i,j,k))*dt(i,j,k))
                2376 
                2377 c-----six exponentials by powers of 8
                2378 c     the constant 0.10624 is equal to 1.66*0.064
                2379 
                2380            h2oexp(i,j,k,1)=exp(-xx*0.10624)
                2381 
                2382            xx=h2oexp(i,j,k,1)*h2oexp(i,j,k,1)
                2383            xx=xx*xx
                2384            h2oexp(i,j,k,2)=xx*xx
                2385 
                2386            xx=h2oexp(i,j,k,2)*h2oexp(i,j,k,2)
                2387            xx=xx*xx
                2388            h2oexp(i,j,k,3)=xx*xx
                2389 
                2390            xx=h2oexp(i,j,k,3)*h2oexp(i,j,k,3)
                2391            xx=xx*xx
                2392            h2oexp(i,j,k,4)=xx*xx
                2393 
                2394            xx=h2oexp(i,j,k,4)*h2oexp(i,j,k,4)
                2395            xx=xx*xx
                2396            h2oexp(i,j,k,5)=xx*xx
                2397 
                2398            xx=h2oexp(i,j,k,5)*h2oexp(i,j,k,5)
                2399            xx=xx*xx
                2400            h2oexp(i,j,k,6)=xx*xx
                2401 
                2402 c-----compute the scaled co2 amount for the sub-band 3a
                2403 c     table 1, Chou et al. (JAS, 1993)
                2404 
                2405            xx=dco2(i,j,k)*(pa(i,j,k)/300.0)**0.5
                2406      1           *(1.+(0.0179+1.02e-4*dt(i,j,k))*dt(i,j,k))
                2407 
                2408 c-----six exponentials by powers of 8
                2409 c     the constant 2.656e-5 is equal to 1.66*1.60e-5
                2410 
                2411            co2exp(i,j,k,1,1)=exp(-xx*2.656e-5)
                2412 
                2413            xx=co2exp(i,j,k,1,1)*co2exp(i,j,k,1,1)
                2414            xx=xx*xx
                2415            co2exp(i,j,k,2,1)=xx*xx
                2416 
                2417            xx=co2exp(i,j,k,2,1)*co2exp(i,j,k,2,1)
                2418            xx=xx*xx
                2419            co2exp(i,j,k,3,1)=xx*xx
                2420 
                2421            xx=co2exp(i,j,k,3,1)*co2exp(i,j,k,3,1)
                2422            xx=xx*xx
                2423            co2exp(i,j,k,4,1)=xx*xx
                2424 
                2425            xx=co2exp(i,j,k,4,1)*co2exp(i,j,k,4,1)
                2426            xx=xx*xx
                2427            co2exp(i,j,k,5,1)=xx*xx
                2428 
                2429            xx=co2exp(i,j,k,5,1)*co2exp(i,j,k,5,1)
                2430            xx=xx*xx
                2431            co2exp(i,j,k,6,1)=xx*xx
                2432 
                2433 c-----one exponential of h2o continuum for sub-band 3a
                2434 c     tabl 5 of Chou et. al. (JAS, 1993)
                2435 c     the constant 1.04995e+2 is equal to 1.66*63.25
                2436 
                2437            conexp(i,j,k,1)=exp(-dcont(i,j,k)*1.04995e+2)
                2438 
                2439 c-----compute the scaled n2o amount for sub-band 3a
                2440 
                2441            xx=dn2o(i,j,k)*(1.+(1.4476e-3+3.6656e-6*dt(i,j,k))*dt(i,j,k))
                2442 
                2443 c-----two exponential2 by powers of 58
                2444 
                2445            n2oexp(i,j,k,1)=exp(-xx*0.25238)
                2446 
                2447            xx=n2oexp(i,j,k,1)*n2oexp(i,j,k,1)
                2448            xx1=xx*xx
                2449            xx1=xx1*xx1
                2450            xx2=xx1*xx1
                2451            xx3=xx2*xx2
                2452            n2oexp(i,j,k,2)=xx*xx1*xx2*xx3
                2453 
                2454           enddo
                2455          enddo
                2456         enddo
                2457 
                2458       return
                2459       end
                2460 c**********************************************************************
                2461       subroutine tablup(k1,k2,m,n,np,nx,nh,nt,sabs,spre,stem,w1,p1,
                2462      *                  dwe,dpe,coef1,coef2,coef3,tran)
                2463 c**********************************************************************
                2464 c   compute water vapor, co2 and o3 transmittances between levels
                2465 c   k1 and k2 for m x n soundings, using table look-up.
                2466 c
                2467 c   calculations follow Eq. (40) of Chou and Suarez (1994)
                2468 c
                2469 c---- input ---------------------
                2470 c  indices for pressure levels (k1 and k2)
                2471 c  number of grid intervals in zonal direction (m)
                2472 c  number of grid intervals in meridional direction (n)
                2473 c  number of atmospheric layers (np)
                2474 c  number of pressure intervals in the table (nx)
                2475 c  number of absorber amount intervals in the table (nh)
                2476 c  number of tables copied (nt)
                2477 c  column-integrated absorber amount (sabs)
                2478 c  column absorber amount-weighted pressure (spre)
                2479 c  column absorber amount-weighted temperature (stem)
                2480 c  first value of absorber amount (log10) in the table (w1) 
                2481 c  first value of pressure (log10) in the table (p1) 
                2482 c  size of the interval of absorber amount (log10) in the table (dwe)
                2483 c  size of the interval of pressure (log10) in the table (dpe)
                2484 c  pre-computed coefficients (coef1, coef2, and coef3)
                2485 c
                2486 c---- updated ---------------------
                2487 c  transmittance (tran)
                2488 c
                2489 c  Note:
                2490 c   (1) units of sabs are g/cm**2 for water vapor and
                2491 c       (cm-atm)stp for co2 and o3.
                2492 c   (2) units of spre and stem are, respectively, mb and K.
                2493 c   (3) there are nt identical copies of the tables (coef1, coef2, and
                2494 c       coef3).  the prupose of using the multiple copies of tables is
                2495 c       to increase the speed in parallel (vectorized) computations.
                2496 C       if such advantage does not exist, nt can be set to 1.
                2497 c   
                2498 c**********************************************************************
                2499       implicit none
                2500       integer k1,k2,m,n,np,nx,nh,nt,i,j
                2501 
                2502 c---- input parameters -----
                2503 
a456aa407c Andr*2504       _RL w1,p1,dwe,dpe
                2505       _RL sabs(m,n,np+1),spre(m,n,np+1),stem(m,n,np+1)
                2506       _RL coef1(nx,nh,nt),coef2(nx,nh,nt),coef3(nx,nh,nt)
1662f365b2 Andr*2507 
                2508 c---- update parameter -----
                2509 
a456aa407c Andr*2510       _RL tran(m,n)
1662f365b2 Andr*2511 
                2512 c---- temporary variables -----
                2513 
a456aa407c Andr*2514       _RL x1,x2,x3,we,pe,fw,fp,pa,pb,pc,ax,ba,bb,t1,ca,cb,t2
1662f365b2 Andr*2515       integer iw,ip,nn
                2516 
                2517 c**********************************************************************
                2518 
                2519       do j=1,n
                2520        do i=1,m
                2521 
                2522         nn=mod(i,nt)+1
                2523 
                2524         x1=sabs(i,j,k2)-sabs(i,j,k1)
                2525         x2=(spre(i,j,k2)-spre(i,j,k1))/x1
                2526         x3=(stem(i,j,k2)-stem(i,j,k1))/x1
                2527 
                2528         we=(log10(x1)-w1)/dwe
                2529         pe=(log10(x2)-p1)/dpe
                2530 
                2531         we=max(we,w1-2.*dwe)
                2532         pe=max(pe,p1)
                2533 
                2534         iw=int(we+1.5)
                2535         ip=int(pe+1.5)
                2536 
                2537         iw=min(iw,nh-1)
                2538         iw=max(iw, 2)
                2539 
                2540         ip=min(ip,nx-1)
                2541         ip=max(ip, 1)
                2542 
                2543         fw=we-float(iw-1)
                2544         fp=pe-float(ip-1)
                2545 
                2546 c-----linear interpolation in pressure
                2547 
                2548         pa = coef1(ip,iw-1,nn)*(1.-fp)+coef1(ip+1,iw-1,nn)*fp
                2549         pb = coef1(ip,iw,  nn)*(1.-fp)+coef1(ip+1,iw,  nn)*fp
                2550         pc = coef1(ip,iw+1,nn)*(1.-fp)+coef1(ip+1,iw+1,nn)*fp
                2551 
                2552 c-----quadratic interpolation in absorber amount for coef1
                2553 
                2554         ax = (-pa*(1.-fw)+pc*(1.+fw)) *fw*0.5 + pb*(1.-fw*fw)
                2555 
                2556 c-----linear interpolation in absorber amount for coef2 and coef3
                2557 
                2558         ba = coef2(ip,iw,  nn)*(1.-fp)+coef2(ip+1,iw,  nn)*fp
                2559         bb = coef2(ip,iw+1,nn)*(1.-fp)+coef2(ip+1,iw+1,nn)*fp
                2560         t1 = ba*(1.-fw) + bb*fw
                2561 
                2562         ca = coef3(ip,iw,  nn)*(1.-fp)+coef3(ip+1,iw,  nn)*fp
                2563         cb = coef3(ip,iw+1,nn)*(1.-fp)+coef3(ip+1,iw+1,nn)*fp
                2564         t2 = ca*(1.-fw) + cb*fw
                2565 
                2566 c-----update the total transmittance between levels k1 and k2
                2567 
                2568         tran(i,j)= (ax + (t1+t2*x3) * x3)*tran(i,j)
                2569 
                2570        enddo
                2571       enddo
                2572 
                2573       return
                2574       end
                2575 c**********************************************************************
                2576       subroutine h2okdis(ib,m,n,np,k,fkw,gkw,ne,h2oexp,conexp,
                2577      *                   th2o,tcon,tran)
                2578 c**********************************************************************
                2579 c   compute water vapor transmittance between levels k1 and k2 for
                2580 c   m x n soundings, using the k-distribution method.
                2581 c
                2582 c   computations follow eqs. (34), (46), (50) and (52).
                2583 c
                2584 c---- input parameters
                2585 c  spectral band (ib)
                2586 c  number of grid intervals in zonal direction (m)
                2587 c  number of grid intervals in meridional direction (n)
                2588 c  number of levels (np)
                2589 c  current level (k)
                2590 c  planck-weighted k-distribution function due to
                2591 c    h2o line absorption (fkw)
                2592 c  planck-weighted k-distribution function due to
                2593 c    h2o continuum absorption (gkw)
                2594 c  number of terms used in each band to compute water vapor
                2595 c     continuum transmittance (ne)
                2596 c  exponentials for line absorption (h2oexp) 
                2597 c  exponentials for continuum absorption (conexp) 
                2598 c
                2599 c---- updated parameters
                2600 c  transmittance between levels k1 and k2 due to
                2601 c    water vapor line absorption (th2o)
                2602 c  transmittance between levels k1 and k2 due to
                2603 c    water vapor continuum absorption (tcon)
                2604 c  total transmittance (tran)
                2605 c
                2606 c**********************************************************************
                2607       implicit none
                2608       integer ib,m,n,np,k,i,j
                2609 
                2610 c---- input parameters ------
                2611 
a456aa407c Andr*2612       _RL conexp(m,n,np,3),h2oexp(m,n,np,6)
1662f365b2 Andr*2613       integer ne(9)
a456aa407c Andr*2614       _RL  fkw(6,9),gkw(6,3)
1662f365b2 Andr*2615 
                2616 c---- updated parameters -----
                2617 
a456aa407c Andr*2618       _RL th2o(m,n,6),tcon(m,n,3),tran(m,n)
1662f365b2 Andr*2619 
                2620 c---- temporary arrays -----
                2621 
a456aa407c Andr*2622       _RL trnth2o
1662f365b2 Andr*2623 
                2624 c-----tco2 are the six exp factors between levels k1 and k2 
                2625 c     tran is the updated total transmittance between levels k1 and k2
                2626 
                2627 c-----th2o is the 6 exp factors between levels k1 and k2 due to
                2628 c     h2o line absorption. 
                2629 
                2630 c-----tcon is the 3 exp factors between levels k1 and k2 due to
                2631 c     h2o continuum absorption.
                2632 
                2633 c-----trnth2o is the total transmittance between levels k1 and k2 due
                2634 c     to both line and continuum absorption computed from eq. (52).
                2635 
                2636          do j=1,n
                2637           do i=1,m
                2638            th2o(i,j,1) = th2o(i,j,1)*h2oexp(i,j,k,1)
                2639            th2o(i,j,2) = th2o(i,j,2)*h2oexp(i,j,k,2)
                2640            th2o(i,j,3) = th2o(i,j,3)*h2oexp(i,j,k,3)
                2641            th2o(i,j,4) = th2o(i,j,4)*h2oexp(i,j,k,4)
                2642            th2o(i,j,5) = th2o(i,j,5)*h2oexp(i,j,k,5)
                2643            th2o(i,j,6) = th2o(i,j,6)*h2oexp(i,j,k,6)
                2644           enddo
                2645          enddo
                2646 
                2647       if (ne(ib).eq.0) then
                2648 
                2649          do j=1,n
                2650           do i=1,m
                2651 
                2652            trnth2o      =(fkw(1,ib)*th2o(i,j,1)
                2653      *                  + fkw(2,ib)*th2o(i,j,2)
                2654      *                  + fkw(3,ib)*th2o(i,j,3)
                2655      *                  + fkw(4,ib)*th2o(i,j,4)
                2656      *                  + fkw(5,ib)*th2o(i,j,5)
                2657      *                  + fkw(6,ib)*th2o(i,j,6))
                2658 
                2659           tran(i,j)=tran(i,j)*trnth2o
                2660 
                2661           enddo
                2662          enddo
                2663 
                2664       elseif (ne(ib).eq.1) then
                2665 
                2666          do j=1,n
                2667           do i=1,m
                2668 
                2669            tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)
                2670 
                2671            trnth2o      =(fkw(1,ib)*th2o(i,j,1)
                2672      *                  + fkw(2,ib)*th2o(i,j,2)
                2673      *                  + fkw(3,ib)*th2o(i,j,3)
                2674      *                  + fkw(4,ib)*th2o(i,j,4)
                2675      *                  + fkw(5,ib)*th2o(i,j,5)
                2676      *                  + fkw(6,ib)*th2o(i,j,6))*tcon(i,j,1)
                2677 
                2678           tran(i,j)=tran(i,j)*trnth2o
                2679 
                2680           enddo
                2681          enddo
                2682 
                2683       else
                2684 
                2685          do j=1,n
                2686           do i=1,m
                2687 
                2688            tcon(i,j,1)= tcon(i,j,1)*conexp(i,j,k,1)
                2689            tcon(i,j,2)= tcon(i,j,2)*conexp(i,j,k,2)
                2690            tcon(i,j,3)= tcon(i,j,3)*conexp(i,j,k,3)
                2691 
                2692 
                2693            trnth2o      = (  gkw(1,1)*th2o(i,j,1)
                2694      *                     + gkw(2,1)*th2o(i,j,2)
                2695      *                     + gkw(3,1)*th2o(i,j,3)
                2696      *                     + gkw(4,1)*th2o(i,j,4)
                2697      *                     + gkw(5,1)*th2o(i,j,5)
                2698      *                     + gkw(6,1)*th2o(i,j,6) ) * tcon(i,j,1)
                2699      *                  + (  gkw(1,2)*th2o(i,j,1)
                2700      *                     + gkw(2,2)*th2o(i,j,2)
                2701      *                     + gkw(3,2)*th2o(i,j,3)
                2702      *                     + gkw(4,2)*th2o(i,j,4)
                2703      *                     + gkw(5,2)*th2o(i,j,5)
                2704      *                     + gkw(6,2)*th2o(i,j,6) ) * tcon(i,j,2)
                2705      *                  + (  gkw(1,3)*th2o(i,j,1)
                2706      *                     + gkw(2,3)*th2o(i,j,2)
                2707      *                     + gkw(3,3)*th2o(i,j,3)
                2708      *                     + gkw(4,3)*th2o(i,j,4)
                2709      *                     + gkw(5,3)*th2o(i,j,5)
                2710      *                     + gkw(6,3)*th2o(i,j,6) ) * tcon(i,j,3)
                2711 
                2712           tran(i,j)=tran(i,j)*trnth2o
                2713 
                2714           enddo
                2715          enddo
                2716 
                2717       endif
                2718 
                2719       return
                2720       end
                2721 c**********************************************************************
                2722       subroutine co2kdis(m,n,np,k,co2exp,tco2,tran)
                2723 c**********************************************************************
                2724 c   compute co2 transmittances between levels k1 and k2 for
                2725 c    m x n soundings, using the k-distribution method with linear
                2726 c    pressure scaling. computations follow eq. (34).
                2727 c
                2728 c---- input parameters
                2729 c   number of grid intervals in zonal direction (m)
                2730 c   number of grid intervals in meridional direction (n)
                2731 c   number of levels (np)
                2732 c   current level (k)
                2733 c   exponentials for co2 absorption (co2exp)
                2734 c
                2735 c---- updated parameters
                2736 c   transmittance between levels k1 and k2 due to co2 absorption
                2737 c     for the various values of the absorption coefficient (tco2)
                2738 c   total transmittance (tran)
                2739 c
                2740 c**********************************************************************
                2741       implicit none
                2742       integer m,n,np,k,i,j
                2743 
                2744 c---- input parameters -----
                2745 
a456aa407c Andr*2746       _RL co2exp(m,n,np,6,2)
1662f365b2 Andr*2747 
                2748 c---- updated parameters -----
                2749 
a456aa407c Andr*2750       _RL tco2(m,n,6,2),tran(m,n)
1662f365b2 Andr*2751 
                2752 c---- temporary arrays -----
                2753 
a456aa407c Andr*2754       _RL xc
1662f365b2 Andr*2755 
                2756 c-----tco2 is the 6 exp factors between levels k1 and k2. 
                2757 c     xc is the total co2 transmittance given by eq. (53).
                2758 
                2759          do j=1,n
                2760           do i=1,m
                2761 
                2762 c-----band-wings
                2763 
                2764            tco2(i,j,1,1)=tco2(i,j,1,1)*co2exp(i,j,k,1,1)
                2765            xc=   0.1395 *tco2(i,j,1,1)
                2766 
                2767            tco2(i,j,2,1)=tco2(i,j,2,1)*co2exp(i,j,k,2,1)
                2768            xc=xc+0.1407 *tco2(i,j,2,1)
                2769 
                2770            tco2(i,j,3,1)=tco2(i,j,3,1)*co2exp(i,j,k,3,1)
                2771            xc=xc+0.1549 *tco2(i,j,3,1)
                2772 
                2773            tco2(i,j,4,1)=tco2(i,j,4,1)*co2exp(i,j,k,4,1)
                2774            xc=xc+0.1357 *tco2(i,j,4,1)
                2775 
                2776            tco2(i,j,5,1)=tco2(i,j,5,1)*co2exp(i,j,k,5,1)
                2777            xc=xc+0.0182 *tco2(i,j,5,1)
                2778 
                2779            tco2(i,j,6,1)=tco2(i,j,6,1)*co2exp(i,j,k,6,1)
                2780            xc=xc+0.0220 *tco2(i,j,6,1)
                2781 
                2782 c-----band-center region
                2783 
                2784            tco2(i,j,1,2)=tco2(i,j,1,2)*co2exp(i,j,k,1,2)
                2785            xc=xc+0.0766 *tco2(i,j,1,2)
                2786 
                2787            tco2(i,j,2,2)=tco2(i,j,2,2)*co2exp(i,j,k,2,2)
                2788            xc=xc+0.1372 *tco2(i,j,2,2)
                2789 
                2790            tco2(i,j,3,2)=tco2(i,j,3,2)*co2exp(i,j,k,3,2)
                2791            xc=xc+0.1189 *tco2(i,j,3,2)
                2792 
                2793            tco2(i,j,4,2)=tco2(i,j,4,2)*co2exp(i,j,k,4,2)
                2794            xc=xc+0.0335 *tco2(i,j,4,2)
                2795 
                2796            tco2(i,j,5,2)=tco2(i,j,5,2)*co2exp(i,j,k,5,2)
                2797            xc=xc+0.0169 *tco2(i,j,5,2)
                2798 
                2799            tco2(i,j,6,2)=tco2(i,j,6,2)*co2exp(i,j,k,6,2)
                2800            xc=xc+0.0059 *tco2(i,j,6,2)
                2801 
                2802            tran(i,j)=tran(i,j)*xc
                2803 
                2804           enddo
                2805          enddo
                2806 
                2807       return
                2808       end
                2809 c**********************************************************************
                2810       subroutine n2okdis(ib,m,n,np,k,n2oexp,tn2o,tran)
                2811 c**********************************************************************
                2812 c   compute n2o transmittances between levels k1 and k2 for
                2813 c    m x n soundings, using the k-distribution method with linear
                2814 c    pressure scaling.
                2815 c
                2816 c---- input parameters
                2817 c   spectral band (ib)
                2818 c   number of grid intervals in zonal direction (m)
                2819 c   number of grid intervals in meridional direction (n)
                2820 c   number of levels (np)
                2821 c   current level (k)
                2822 c   exponentials for n2o absorption (n2oexp)
                2823 c
                2824 c---- updated parameters
                2825 c   transmittance between levels k1 and k2 due to n2o absorption
                2826 c     for the various values of the absorption coefficient (tn2o)
                2827 c   total transmittance (tran)
                2828 c
                2829 c**********************************************************************
                2830       implicit none
                2831       integer ib,m,n,np,k,i,j
                2832 
                2833 c---- input parameters -----
                2834 
a456aa407c Andr*2835       _RL n2oexp(m,n,np,4)
1662f365b2 Andr*2836 
                2837 c---- updated parameters -----
                2838 
a456aa407c Andr*2839       _RL tn2o(m,n,4),tran(m,n)
1662f365b2 Andr*2840 
                2841 c---- temporary arrays -----
                2842 
a456aa407c Andr*2843       _RL xc
1662f365b2 Andr*2844 
                2845 c-----tn2o is the 2 exp factors between levels k1 and k2. 
                2846 c     xc is the total n2o transmittance
                2847 
                2848         do j=1,n
                2849          do i=1,m
                2850 
                2851 c-----band 6
                2852 
                2853           if (ib.eq.6) then
                2854 
                2855            tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
                2856            xc=   0.940414*tn2o(i,j,1)
                2857 
                2858            tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
                2859            xc=xc+0.059586*tn2o(i,j,2)
                2860 
                2861 c-----band 7
                2862 
                2863           else
                2864 
                2865            tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
                2866            xc=   0.561961*tn2o(i,j,1)
                2867 
                2868            tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
                2869            xc=xc+0.138707*tn2o(i,j,2)
                2870 
                2871            tn2o(i,j,3)=tn2o(i,j,3)*n2oexp(i,j,k,3)
                2872            xc=xc+0.240670*tn2o(i,j,3)
                2873 
                2874            tn2o(i,j,4)=tn2o(i,j,4)*n2oexp(i,j,k,4)
                2875            xc=xc+0.058662*tn2o(i,j,4)
                2876 
                2877           endif
                2878 
                2879            tran(i,j)=tran(i,j)*xc
                2880 
                2881          enddo
                2882         enddo
                2883 
                2884       return
                2885       end
                2886 c**********************************************************************
                2887       subroutine ch4kdis(ib,m,n,np,k,ch4exp,tch4,tran)
                2888 c**********************************************************************
                2889 c   compute ch4 transmittances between levels k1 and k2 for
                2890 c    m x n soundings, using the k-distribution method with
                2891 c    linear pressure scaling.
                2892 c
                2893 c---- input parameters
                2894 c   spectral band (ib)
                2895 c   number of grid intervals in zonal direction (m)
                2896 c   number of grid intervals in meridional direction (n)
                2897 c   number of levels (np)
                2898 c   current level (k)
                2899 c   exponentials for ch4 absorption (ch4exp)
                2900 c
                2901 c---- updated parameters
                2902 c   transmittance between levels k1 and k2 due to ch4 absorption
                2903 c     for the various values of the absorption coefficient (tch4)
                2904 c   total transmittance (tran)
                2905 c
                2906 c**********************************************************************
                2907       implicit none
                2908       integer ib,m,n,np,k,i,j
                2909 
                2910 c---- input parameters -----
                2911 
a456aa407c Andr*2912       _RL ch4exp(m,n,np,4)
1662f365b2 Andr*2913 
                2914 c---- updated parameters -----
                2915 
a456aa407c Andr*2916       _RL tch4(m,n,4),tran(m,n)
1662f365b2 Andr*2917 
                2918 c---- temporary arrays -----
                2919 
a456aa407c Andr*2920       _RL xc
1662f365b2 Andr*2921 
                2922 c-----tch4 is the 2 exp factors between levels k1 and k2. 
                2923 c     xc is the total ch4 transmittance
                2924 
                2925         do j=1,n
                2926          do i=1,m
                2927 
                2928 c-----band 6
                2929 
                2930           if (ib.eq.6) then
                2931 
                2932            tch4(i,j,1)=tch4(i,j,1)*ch4exp(i,j,k,1)
                2933            xc= tch4(i,j,1)
                2934 
                2935 c-----band 7
                2936 
                2937           else
                2938 
                2939            tch4(i,j,1)=tch4(i,j,1)*ch4exp(i,j,k,1)
                2940            xc=   0.610650*tch4(i,j,1)
                2941 
                2942            tch4(i,j,2)=tch4(i,j,2)*ch4exp(i,j,k,2)
                2943            xc=xc+0.280212*tch4(i,j,2)
                2944 
                2945            tch4(i,j,3)=tch4(i,j,3)*ch4exp(i,j,k,3)
                2946            xc=xc+0.107349*tch4(i,j,3)
                2947 
                2948            tch4(i,j,4)=tch4(i,j,4)*ch4exp(i,j,k,4)
                2949            xc=xc+0.001789*tch4(i,j,4)
                2950 
                2951           endif
                2952 
                2953            tran(i,j)=tran(i,j)*xc
                2954 
                2955          enddo
                2956         enddo
                2957 
                2958       return
                2959       end
                2960 c**********************************************************************
                2961       subroutine comkdis(ib,m,n,np,k,comexp,tcom,tran)
                2962 c**********************************************************************
                2963 c  compute co2-minor transmittances between levels k1 and k2
                2964 c   for m x n soundings, using the k-distribution method
                2965 c   with linear pressure scaling.
                2966 c
                2967 c---- input parameters
                2968 c   spectral band (ib)
                2969 c   number of grid intervals in zonal direction (m)
                2970 c   number of grid intervals in meridional direction (n)
                2971 c   number of levels (np)
                2972 c   current level (k)
                2973 c   exponentials for co2-minor absorption (comexp)
                2974 c
                2975 c---- updated parameters
                2976 c   transmittance between levels k1 and k2 due to co2-minor absorption
                2977 c     for the various values of the absorption coefficient (tcom)
                2978 c   total transmittance (tran)
                2979 c
                2980 c**********************************************************************
                2981       implicit none
                2982       integer ib,m,n,np,k,i,j
                2983 
                2984 c---- input parameters -----
                2985 
a456aa407c Andr*2986       _RL comexp(m,n,np,2)
1662f365b2 Andr*2987 
                2988 c---- updated parameters -----
                2989 
a456aa407c Andr*2990       _RL tcom(m,n,2),tran(m,n)
1662f365b2 Andr*2991 
                2992 c---- temporary arrays -----
                2993 
a456aa407c Andr*2994       _RL xc
1662f365b2 Andr*2995 
                2996 c-----tcom is the 2 exp factors between levels k1 and k2. 
                2997 c     xc is the total co2-minor transmittance
                2998 
                2999          do j=1,n
                3000           do i=1,m
                3001 
                3002 c-----band 4
                3003 
                3004            if (ib.eq.4) then
                3005 
                3006             tcom(i,j,1)=tcom(i,j,1)*comexp(i,j,k,1)
                3007             xc=   0.972025*tcom(i,j,1)
                3008             tcom(i,j,2)=tcom(i,j,2)*comexp(i,j,k,2)
                3009             xc=xc+0.027975*tcom(i,j,2)
                3010 
                3011 c-----band 5
                3012 
                3013            else
                3014 
                3015             tcom(i,j,1)=tcom(i,j,1)*comexp(i,j,k,1)
                3016             xc=   0.961324*tcom(i,j,1)
                3017             tcom(i,j,2)=tcom(i,j,2)*comexp(i,j,k,2)
                3018             xc=xc+0.038676*tcom(i,j,2)
                3019 
                3020            endif
                3021 
                3022             tran(i,j)=tran(i,j)*xc
                3023 
                3024           enddo
                3025          enddo
                3026 
                3027       return
                3028       end
                3029 c**********************************************************************
                3030       subroutine cfckdis(m,n,np,k,cfcexp,tcfc,tran)
                3031 c**********************************************************************
                3032 c  compute cfc-(11,12,22) transmittances between levels k1 and k2
                3033 c   for m x n soundings, using the k-distribution method with
                3034 c   linear pressure scaling.
                3035 c
                3036 c---- input parameters
                3037 c   number of grid intervals in zonal direction (m)
                3038 c   number of grid intervals in meridional direction (n)
                3039 c   number of levels (np)
                3040 c   current level (k)
                3041 c   exponentials for cfc absorption (cfcexp)
                3042 c
                3043 c---- updated parameters
                3044 c   transmittance between levels k1 and k2 due to cfc absorption
                3045 c     for the various values of the absorption coefficient (tcfc)
                3046 c   total transmittance (tran)
                3047 c
                3048 c**********************************************************************
                3049       implicit none
                3050       integer m,n,np,k,i,j
                3051 
                3052 c---- input parameters -----
                3053 
a456aa407c Andr*3054       _RL cfcexp(m,n,np)
1662f365b2 Andr*3055 
                3056 c---- updated parameters -----
                3057 
a456aa407c Andr*3058       _RL tcfc(m,n),tran(m,n)
1662f365b2 Andr*3059 
                3060 c-----tcfc is the exp factors between levels k1 and k2. 
                3061 
                3062          do j=1,n
                3063           do i=1,m
                3064 
                3065             tcfc(i,j)=tcfc(i,j)*cfcexp(i,j,k)
                3066             tran(i,j)=tran(i,j)*tcfc(i,j)
                3067 
                3068           enddo
                3069          enddo
                3070 
                3071       return
                3072       end
                3073 c**********************************************************************
                3074       subroutine b10kdis(m,n,np,k,h2oexp,conexp,co2exp,n2oexp
                3075      *          ,th2o,tcon,tco2,tn2o,tran)
                3076 c**********************************************************************
                3077 c
                3078 c   compute h2o (line and continuum),co2,n2o transmittances between
                3079 c   levels k1 and k2 for m x n soundings, using the k-distribution
                3080 c   method with linear pressure scaling.
                3081 c
                3082 c---- input parameters
                3083 c   number of grid intervals in zonal direction (m)
                3084 c   number of grid intervals in meridional direction (n)
                3085 c   number of levels (np)
                3086 c   current level (k)
                3087 c   exponentials for h2o line absorption (h2oexp)
                3088 c   exponentials for h2o continuum absorption (conexp)
                3089 c   exponentials for co2 absorption (co2exp)
                3090 c   exponentials for n2o absorption (n2oexp)
                3091 c
                3092 c---- updated parameters
                3093 c   transmittance between levels k1 and k2 due to h2o line absorption
                3094 c     for the various values of the absorption coefficient (th2o)
                3095 c   transmittance between levels k1 and k2 due to h2o continuum
                3096 c     absorption for the various values of the absorption
                3097 c     coefficient (tcon)
                3098 c   transmittance between levels k1 and k2 due to co2 absorption
                3099 c     for the various values of the absorption coefficient (tco2)
                3100 c   transmittance between levels k1 and k2 due to n2o absorption
                3101 c     for the various values of the absorption coefficient (tn2o)
                3102 c   total transmittance (tran)
                3103 c
                3104 c**********************************************************************
                3105       implicit none
                3106       integer m,n,np,k,i,j
                3107 
                3108 c---- input parameters -----
                3109 
a456aa407c Andr*3110       _RL h2oexp(m,n,np,6),conexp(m,n,np,3),co2exp(m,n,np,6,2)
1662f365b2 Andr*3111      *    ,n2oexp(m,n,np,4)
                3112 
                3113 c---- updated parameters -----
                3114 
a456aa407c Andr*3115       _RL th2o(m,n,6),tcon(m,n,3),tco2(m,n,6,2),tn2o(m,n,4)
1662f365b2 Andr*3116      *    ,tran(m,n)
                3117 
                3118 c---- temporary arrays -----
                3119 
a456aa407c Andr*3120       _RL xx
1662f365b2 Andr*3121 
                3122 c-----initialize tran
                3123 
                3124        do j=1,n
                3125         do i=1,m
                3126            tran(i,j)=1.0
                3127         enddo
                3128        enddo
                3129 
                3130 c-----for h2o line
                3131 
                3132         do j=1,n
                3133          do i=1,m
                3134 
                3135            th2o(i,j,1)=th2o(i,j,1)*h2oexp(i,j,k,1)
                3136            xx=   0.3153*th2o(i,j,1)
                3137 
                3138            th2o(i,j,2)=th2o(i,j,2)*h2oexp(i,j,k,2)
                3139            xx=xx+0.4604*th2o(i,j,2)
                3140 
                3141            th2o(i,j,3)=th2o(i,j,3)*h2oexp(i,j,k,3)
                3142            xx=xx+0.1326*th2o(i,j,3)
                3143 
                3144            th2o(i,j,4)=th2o(i,j,4)*h2oexp(i,j,k,4)
                3145            xx=xx+0.0798*th2o(i,j,4)
                3146 
                3147            th2o(i,j,5)=th2o(i,j,5)*h2oexp(i,j,k,5)
                3148            xx=xx+0.0119*th2o(i,j,5)
                3149 
                3150            tran(i,j)=tran(i,j)*xx
                3151 
                3152          enddo
                3153         enddo
                3154 
                3155 c-----for h2o continuum
                3156 
                3157         do j=1,n
                3158          do i=1,m
                3159 
                3160            tcon(i,j,1)=tcon(i,j,1)*conexp(i,j,k,1)
                3161            tran(i,j)=tran(i,j)*tcon(i,j,1)
                3162 
                3163          enddo
                3164         enddo
                3165  
                3166 c-----for co2 
                3167 
                3168         do j=1,n
                3169          do i=1,m
                3170 
                3171            tco2(i,j,1,1)=tco2(i,j,1,1)*co2exp(i,j,k,1,1)
                3172            xx=    0.2673*tco2(i,j,1,1)
                3173 
                3174            tco2(i,j,2,1)=tco2(i,j,2,1)*co2exp(i,j,k,2,1)
                3175            xx=xx+ 0.2201*tco2(i,j,2,1)
                3176 
                3177            tco2(i,j,3,1)=tco2(i,j,3,1)*co2exp(i,j,k,3,1)
                3178            xx=xx+ 0.2106*tco2(i,j,3,1)
                3179 
                3180            tco2(i,j,4,1)=tco2(i,j,4,1)*co2exp(i,j,k,4,1)
                3181            xx=xx+ 0.2409*tco2(i,j,4,1)
                3182 
                3183            tco2(i,j,5,1)=tco2(i,j,5,1)*co2exp(i,j,k,5,1)
                3184            xx=xx+ 0.0196*tco2(i,j,5,1)
                3185 
                3186            tco2(i,j,6,1)=tco2(i,j,6,1)*co2exp(i,j,k,6,1)
                3187            xx=xx+ 0.0415*tco2(i,j,6,1)
                3188 
                3189            tran(i,j)=tran(i,j)*xx
                3190 
                3191          enddo
                3192         enddo
                3193 
                3194 c-----for n2o
                3195 
                3196         do j=1,n
                3197          do i=1,m
                3198 
                3199            tn2o(i,j,1)=tn2o(i,j,1)*n2oexp(i,j,k,1)
                3200            xx=   0.970831*tn2o(i,j,1)
                3201 
                3202            tn2o(i,j,2)=tn2o(i,j,2)*n2oexp(i,j,k,2)
                3203            xx=xx+0.029169*tn2o(i,j,2)
                3204 
                3205            tran(i,j)=tran(i,j)*(xx-1.0)
                3206 
                3207          enddo
                3208         enddo
                3209 
                3210       return
                3211       end