File: C:\NOAA\NEMS_11731\src\atmos\phys\radlw_main.f

1     !!!!!  ==========================================================  !!!!!
2     !!!!!              rrtm1 radiation package description             !!!!!
3     !!!!!  ==========================================================  !!!!!
4     !                                                                      !
5     !    the rrtm1 package includes these parts:                           !
6     !                                                                      !
7     !       'radlw_rrtm1_param.f'                                          !
8     !       'radlw_rrtm1_datatb.f'                                         !
9     !       'radlw_rrtm1_main.f'                                           !
10     !                                                                      !
11     !    the 'radlw_rrtm1_param.f' contains:                               !
12     !                                                                      !
13     !       'module_radlw_cntr_para'   -- control parameters set up        !
14     !       'module_radlw_parameters'  -- band parameters set up           !
15     !                                                                      !
16     !    the 'radlw_rrtm1_datatb.f' contains:                              !
17     !                                                                      !
18     !       'module_radlw_avplank'     -- plank flux data                  !
19     !       'module_radlw_cldprlw'     -- cloud property coefficients      !
20     !       'module_radlw_kgbnn'       -- absorption coeffients for 16     !
21     !                                     bands, where nn = 01-16          !
22     !                                                                      !
23     !    the 'radlw_rrtm1_main.f' contains:                                !
24     !                                                                      !
25     !       'module_radlw_main'        -- main lw radiation transfer       !
26     !                                                                      !
27     !    in the main module 'module_radlw_main' there are only two         !
28     !    externally callable subroutines:                                  !
29     !                                                                      !
30     !                                                                      !
31     !       'lwrad'     -- main rrtm1 lw radiation routine                 !
32     !          inputs:                                                     !
33     !           (pmid,pint,tmid,tint,qnm,o3mr,gasvmr,                      !
34     !            clouds,iauxil,aerosols,sfemis,sfgtmp,                     !
35     !            npts, nlay, nlp1, iflip, lprnt,                           !
36     !          outputs:                                                    !
37     !            hlwc,topflx,sfcflx,                                       !
38     !!         optional outputs:                                           !
39     !            HLW0,HLWB,FLXPRF)                                         !
40     !                                                                      !
41     !       'rlwinit'   -- initialization routine                          !
42     !          inputs:                                                     !
43     !           ( icwp, me, nlay, iovr, isubc )                            !
44     !          outputs:                                                    !
45     !           (none)                                                     !
46     !                                                                      !
47     !    all the lw radiation subprograms become contained subprograms     !
48     !    in module 'module_radlw_main' and many of them are not directly   !
49     !    accessable from places outside the module.                        !
50     !                                                                      !
51     !                                                                      !
52     !    derived data type constructs used:                                !
53     !                                                                      !
54     !     1. radiation flux at toa: (from module 'module_radlw_parameters')!
55     !          topflw_type   -  derived data type for toa rad fluxes       !
56     !            upfxc              total sky upward flux at toa           !
57     !            upfx0              clear sky upward flux at toa           !
58     !                                                                      !
59     !     2. radiation flux at sfc: (from module 'module_radlw_parameters')!
60     !          sfcflw_type   -  derived data type for sfc rad fluxes       !
61     !            upfxc              total sky upward flux at sfc           !
62     !            upfx0              clear sky upward flux at sfc           !
63     !            dnfxc              total sky downward flux at sfc         !
64     !            dnfx0              clear sky downward flux at sfc         !
65     !                                                                      !
66     !     3. radiation flux profiles(from module 'module_radlw_parameters')!
67     !          proflw_type    -  derived data type for rad vertical prof   !
68     !            upfxc              level upward flux for total sky        !
69     !            dnfxc              level downward flux for total sky      !
70     !            upfx0              level upward flux for clear sky        !
71     !            dnfx0              level downward flux for clear sky      !
72     !                                                                      !
73     !    external modules referenced:                                      !
74     !                                                                      !
75     !       'module machine'                                               !
76     !       'module physcons'                                              !
77     !                                                                      !
78     !    compilation sequence is:                                          !
79     !                                                                      !
80     !       'radlw_rrtm1_param.f'                                          !
81     !       'radlw_rrtm1_datatb.f'                                         !
82     !       'radlw_rrtm1_main.f'                                           !
83     !                                                                      !
84     !    and all should be put in front of routines that use lw modules    !
85     !                                                                      !
86     !                                                                      !
87     !                                                                      !
88     !    the original program declarations:                                !
89     !                                                                      !
90     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
91     !                                                                      !
92     ! Copyright 2002, 2003, Atmospheric & Environmental Research, Inc.(AER)!
93     ! This software may be used, copied, or redistributed as long as it is !
94     ! not sold and this copyright notice is reproduced on each copy made.  !
95     ! This model is provided as is without any express or implied warranties
96     !                      (http://www.rtweb.aer.com/)                     !
97     !                                                                      !
98     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
99     !                                                                      !
100     !                               rrtm                                   !
101     !                                                                      !
102     !                   rapid radiative transfer model                     !
103     !                                                                      !
104     !            atmospheric and environmental research, inc.              !
105     !                        840 memorial drive                            !
106     !                        cambridge, ma 02139                           !
107     !                                                                      !
108     !                           eli j. mlawer                              !
109     !                         steven j. taubman~                           !
110     !                         shepard a. clough                            !
111     !                                                                      !
112     !                         ~currently at gfdl                           !
113     !                                                                      !
114     !                       email:  mlawer@aer.com                         !
115     !                                                                      !
116     !        the authors wish to acknowledge the contributions of the      !
117     !        following people:  patrick d. brown, michael j. iacono,       !
118     !        ronald e. farren, luke chen, robert bergstrom.                !
119     !                                                                      !
120     ! ******************************************************************** !
121     !                                                                      !
122     !    ncep modifications history log:                                   !
123     !                                                                      !
124     !       nov 1999,  ken campana  -- received the original code from aer !
125     !                    (1998 ncar ccm version), updated to link up with  !
126     !                    ncep mrf model                                    !
127     !       jun 2000,  ken campana  -- added option to switch random and   !
128     !                    maximum/random cloud overlap                      !
129     !           2001,  s. moorthi   -- further updates for mrf model       !
130     !       may 2001,  yu-tai hou   -- updated on trace gases and cloud    !
131     !                    property based on rrtm_v3.0 codes.                !
132     !       dec 2001,  yu-tai hou   -- rewritten code into fortran 90 std  !
133     !                    set ncep radiation structure std that contains    !
134     !                    three plug-in compatable fortran program files:   !
135     !                    'radlw_param.f', 'radlw_datatb.f', 'radlw_main.f' !
136     !                    fixed bugs in subprograms taugb14, taugb2, etc.   !
137     !                    added out-of-bounds protections. (a detailed note !
138     !                    of up_to_date modifications/corrections by ncep   !
139     !                    was sent to aer in 2002)                          !
140     !       jun 2004,  yu-tai hou   -- add mike iacono's apr 2004          !
141     !                    modification of variable diffusivity angles.      !
142     !       apr 2005,  yu-tai hou   -- minor modifications on module       !
143     !                    structures include rain/snow effect (this version !
144     !                    of code was given back to aer in jun 2006)        !
145     !       mar 2007,  yu-tai hou   -- add aerosol effect for ncep models  !
146     !                    using the generallized aerosol optical property   !
147     !                    scheme for gfs model.                             !
148     !       apr 2007,  yu-tai hou   -- add spectral band heating as an     !
149     !                    optional output to support the 500 km gfs model's !
150     !                    upper stratospheric radiation calculations. and   !
151     !                    restructure optional outputs for easy access by   !
152     !                    different models.                                 !
153     !       mar 2009,  yu-tai hou   -- modified the program interface that !
154     !                    is consistant with newer radiation code structure.!
155     !                                                                      !
156     !                                                                      !
157     !!!!!  ==========================================================  !!!!!
158     !!!!!                       end descriptions                       !!!!!
159     !!!!!  ==========================================================  !!!!!
160     
161     
162     
163     !========================================!
164           module module_radlw_main           !
165     !........................................!
166     !
167           use machine,           only : kind_phys
168           use physcons,          only : con_g, con_cp, con_avgd, con_amd,   &
169          &                              con_amw, con_amo3
170     
171           use module_radlw_parameters
172           use module_radlw_cntr_para
173     !
174           implicit none
175     !
176           private
177     !
178     !  ...  version tag and last revision date
179     !     character(24), parameter :: VTAGLW='RRTM-LW v2.3g   Apr 2004'
180     !     character(24), parameter :: VTAGLW='RRTM-LW v2.3g   Mar 2007'
181           character(24), parameter :: VTAGLW='RRTM-LW v2.3g   Apr 2007'
182     
183     !  ---  constant values
184           real (kind=kind_phys) :: eps, oneminus, bpade, stpfac, wtnum      &
185          &,     co2fac, f_zero
186     
187           parameter (eps=1.0e-6,  oneminus=1.0-eps)
188           parameter (bpade=1.0/0.278)      ! pade approximation constant
189           parameter (stpfac=296./1013.)
190           parameter (wtnum=0.5)
191     !     parameter (avgdro=6.022e23)      ! avogadro constant  (1/mol)
192     !mji  parameter (secang=1.66)          ! diffusivity angle
193           parameter (co2fac=3.55e-4)       ! factor for cal. of co2mult
194           parameter (f_zero=0.0)
195     
196     !  ...  atomic weights for conversion from mass to volume mixing ratios
197           real (kind=kind_phys) :: amdw, amdo3
198     
199           parameter (amdw =con_amd/con_amw)
200           parameter (amdo3=con_amd/con_amo3)
201     
202     !  ...  band indices
203           integer :: nspa(NBANDS), nspb(NBANDS), ngb(NGPTLW)
204     
205           data nspa / 1, 1,10, 9, 9, 1, 9, 1,11, 1, 1, 9, 9, 1, 9, 9 /
206           data nspb / 1, 1, 5, 6, 5, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0 /
207           data ngb  /  8*1, 14*2, 16*3, 14*4, 16*5,  8*6, 12*7,  8*8,       &
208          &            12*9, 6*10, 8*11, 8*12, 4*13, 2*14, 2*15, 2*16 /
209     
210     !  ...  band wavenumber intervals
211     !     real (kind=kind_phys) :: wavenum1(NBANDS), wavenum2(NBANDS)
212     !     data wavenum1/                                                    &
213     !    &         10.,  250.,  500.,  630.,  700.,  820.,  980., 1080.,    &
214     !    &       1180., 1390., 1480., 1800., 2080., 2250., 2380., 2600. /
215     !     data wavenum2/                                                    &
216     !    &        250.,  500.,  630.,  700.,  820.,  980., 1080., 1180.,    &
217     !    &       1390., 1480., 1800., 2080., 2250., 2380., 2600., 3000. /
218           real (kind=kind_phys) :: delwave(NBANDS)
219           data delwave / 240., 250., 130.,  70., 120., 160., 100., 100.,    &
220          &               210.,  90., 320., 280., 170., 130., 220., 400. /
221     
222     !mji ... coefficients for variable diffusivity angle
223           real (kind=kind_phys), dimension(NBANDS) :: a0, a1, a2
224           data a0 /  1.66, 1.55, 1.58, 1.66, 1.54,1.454, 1.89, 1.33,        &
225          &          1.668, 1.66, 1.66, 1.66, 1.66, 1.66, 1.66, 1.66  /
226           data a1 /  0.00, 0.25, 0.22, 0.00, 0.13,0.446,-0.10, 0.40,        &
227          &         -0.006, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00  /
228           data a2 /  0.00,-12.0,-11.7, 0.00,-0.72,-0.243,0.19,-0.062,       &
229          &          0.414, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00  /
230     
231     !  ---  reference pressure and temperature
232           real (kind=kind_phys), dimension(59) :: pref, preflog, tref
233     
234     !  ...  these pressures are chosen such that the ln of the first one
235     !       has only a few non-zero digits (i.e. ln(pref(1)) = 6.96000) and
236     !       each subsequent ln(pref) differs from the previous one by 0.2.
237           data pref /                                                       &
238          &    1.05363e+03,8.62642e+02,7.06272e+02,5.78246e+02,4.73428e+02,  &
239          &    3.87610e+02,3.17348e+02,2.59823e+02,2.12725e+02,1.74164e+02,  &
240          &    1.42594e+02,1.16746e+02,9.55835e+01,7.82571e+01,6.40715e+01,  &
241          &    5.24573e+01,4.29484e+01,3.51632e+01,2.87892e+01,2.35706e+01,  &
242          &    1.92980e+01,1.57998e+01,1.29358e+01,1.05910e+01,8.67114e+00,  &
243          &    7.09933e+00,5.81244e+00,4.75882e+00,3.89619e+00,3.18993e+00,  &
244          &    2.61170e+00,2.13828e+00,1.75067e+00,1.43333e+00,1.17351e+00,  &
245          &    9.60789e-01,7.86628e-01,6.44036e-01,5.27292e-01,4.31710e-01,  &
246          &    3.53455e-01,2.89384e-01,2.36928e-01,1.93980e-01,1.58817e-01,  &
247          &    1.30029e-01,1.06458e-01,8.71608e-02,7.13612e-02,5.84256e-02,  &
248          &    4.78349e-02,3.91639e-02,3.20647e-02,2.62523e-02,2.14936e-02,  &
249          &    1.75975e-02,1.44076e-02,1.17959e-02,9.65769e-03 /
250           data preflog /                                                    &
251          &     6.9600e+00, 6.7600e+00, 6.5600e+00, 6.3600e+00, 6.1600e+00,  &
252          &     5.9600e+00, 5.7600e+00, 5.5600e+00, 5.3600e+00, 5.1600e+00,  &
253          &     4.9600e+00, 4.7600e+00, 4.5600e+00, 4.3600e+00, 4.1600e+00,  &
254          &     3.9600e+00, 3.7600e+00, 3.5600e+00, 3.3600e+00, 3.1600e+00,  &
255          &     2.9600e+00, 2.7600e+00, 2.5600e+00, 2.3600e+00, 2.1600e+00,  &
256          &     1.9600e+00, 1.7600e+00, 1.5600e+00, 1.3600e+00, 1.1600e+00,  &
257          &     9.6000e-01, 7.6000e-01, 5.6000e-01, 3.6000e-01, 1.6000e-01,  &
258          &    -4.0000e-02,-2.4000e-01,-4.4000e-01,-6.4000e-01,-8.4000e-01,  &
259          &    -1.0400e+00,-1.2400e+00,-1.4400e+00,-1.6400e+00,-1.8400e+00,  &
260          &    -2.0400e+00,-2.2400e+00,-2.4400e+00,-2.6400e+00,-2.8400e+00,  &
261          &    -3.0400e+00,-3.2400e+00,-3.4400e+00,-3.6400e+00,-3.8400e+00,  &
262          &    -4.0400e+00,-4.2400e+00,-4.4400e+00,-4.6400e+00 /
263     !  ...  these are the temperatures associated with the respective
264     !       pressures for the MLS standard atmosphere.
265           data tref /                                                       &
266          &     2.9420E+02, 2.8799E+02, 2.7894E+02, 2.6925E+02, 2.5983E+02,  &
267          &     2.5017E+02, 2.4077E+02, 2.3179E+02, 2.2306E+02, 2.1578E+02,  &
268          &     2.1570E+02, 2.1570E+02, 2.1570E+02, 2.1706E+02, 2.1858E+02,  &
269          &     2.2018E+02, 2.2174E+02, 2.2328E+02, 2.2479E+02, 2.2655E+02,  &
270          &     2.2834E+02, 2.3113E+02, 2.3401E+02, 2.3703E+02, 2.4022E+02,  &
271          &     2.4371E+02, 2.4726E+02, 2.5085E+02, 2.5457E+02, 2.5832E+02,  &
272          &     2.6216E+02, 2.6606E+02, 2.6999E+02, 2.7340E+02, 2.7536E+02,  &
273          &     2.7568E+02, 2.7372E+02, 2.7163E+02, 2.6955E+02, 2.6593E+02,  &
274          &     2.6211E+02, 2.5828E+02, 2.5360E+02, 2.4854E+02, 2.4348E+02,  &
275          &     2.3809E+02, 2.3206E+02, 2.2603E+02, 2.2000E+02, 2.1435E+02,  &
276          &     2.0887E+02, 2.0340E+02, 1.9792E+02, 1.9290E+02, 1.8809E+02,  &
277          &     1.8329E+02, 1.7849E+02, 1.7394E+02, 1.7212E+02 /
278     
279     !! ---  logical flags for optional output fields
280     
281           logical :: lhlwb  = .false.
282           logical :: lhlw0  = .false.
283           logical :: lflxprf= .false.
284     
285     !  ---  those data will be set up only once by "rlwinit"
286     
287     !  ...  fluxfac, heatfac are factors for fluxes (in w/m**2) and heating
288     !       rates (in k/day, or k/sec set by subroutine 'rlwinit')
289     !       semiss0 are default surface emissivity for each bands
290     
291           real (kind=kind_phys) :: fluxfac, heatfac, semiss0(NBANDS)
292     
293           real (kind=kind_phys), dimension(0:N5000) :: tau, tf, trans
294           real (kind=kind_phys), dimension(0:N200 ) :: corr1, corr2
295     
296     !  ...  iovrlw  is the clouds overlapping control flag
297     !            =0: random overlapping clouds
298     !            =1: maximum/random overlapping clouds
299     
300           integer :: iovrlw
301     
302           public lwrad, rlwinit
303     
304     
305     ! =================
306           contains
307     ! =================
308     
309     !-----------------------------------
310           subroutine lwrad                                                  &
311     !...................................
312     
313     !  ---  inputs:
314          &     ( pmid,pint,tmid,tint,qnm,o3mr,gasvmr,                       &
315          &       clouds,iauxil,aerosols,sfemis,sfgtmp,                      &
316          &       NPTS, NLAY, NLP1, iflip, lprnt,                            &
317     !  ---  outputs:
318          &       hlwc,topflx,sfcflx                                         &
319     !! ---  optional:
320          &,      HLW0,HLWB,FLXPRF                                           &
321          &     )
322     
323     !  ====================  defination of variables  ===================  !
324     !                                                                      !
325     !  input variables:                                                    !
326     !     pmid   (NPTS,NLAY)    - layer pressures (mb)                     !
327     !     pint   (NPTS,NLP1)    - interface pressures (mb)                 !
328     !     tmid   (NPTS,NLAY)    - layer temperature (k)                    !
329     !     tint   (NPTS,NLP1)    - interface temperatures (k)               !
330     !     qnm    (NPTS,NLAY)    - layer h2o mixing ratio (gm/gm)*see inside!
331     !     o3mr   (NPTS,NLAY)    - layer o3 mixing ratio (gm/gm) *see inside!
332     !     gasvmr (NPTS,NLAY,:)  - atmospheric gases amount:                !
333     !                       (check module_radiation_gases for definition)  !
334     !       gasvmr(:,:,1)   -      co2 volume mixing ratio                 !
335     !       gasvmr(:,:,2)   -      n2o volume mixing ratio                 !
336     !       gasvmr(:,:,3)   -      ch4 volume mixing ratio                 !
337     !       gasvmr(:,:,4)   -      o2  volume mixing ratio                 !
338     !       gasvmr(:,:,5)   -      co  volume mixing ratio                 !
339     !       gasvmr(:,:,6)   -      cfc11 volume mixing ratio               !
340     !       gasvmr(:,:,7)   -      cfc12 volume mixing ratio               !
341     !       gasvmr(:,:,8)   -      cfc22 volume mixing ratio               !
342     !       gasvmr(:,:,9)   -      ccl4  volume mixing ratio               !
343     !     clouds (NPTS,NLAY,:)  - layer cloud profiles:                    !
344     !                       (check module_radiation_clouds for definition) !
345     !                ---  for  iflagliq > 0  ---                           !
346     !       clouds(:,:,1)  -   layer total cloud fraction                  !
347     !       clouds(:,:,2)  -   layer cloud liq water path      (g/m**2)    !
348     !       clouds(:,:,3)  -   mean eff radius for liq cloud   (micron)    !
349     !       clouds(:,:,4)  -   layer cloud ice water path      (g/m**2)    !
350     !       clouds(:,:,5)  -   mean eff radius for ice cloud   (micron)    !
351     !       clouds(:,:,6)  -   layer rain drop water path      (g/m**2)    !
352     !       clouds(:,:,7)  -   mean eff radius for rain drop   (micron)    !
353     !       clouds(:,:,8)  -   layer snow flake water path     (g/m**2)    !
354     !   ** fu's scheme need to be normalized by snow density (g/m**3/1.0e6)!
355     !       clouds(:,:,9)  -   mean eff radius for snow flake  (micron)    !
356     !                ---  for  iflagliq = 0  ---                           !
357     !       clouds(:,:,1)  -   layer total cloud fraction                  !
358     !       clouds(:,:,2)  -   layer cloud optical depth                   !
359     !       clouds(:,:,3)  -   layer cloud single scattering albedo        !
360     !       clouds(:,:,4)  -   layer cloud asymmetry factor                !
361     !     iauxil(:)             - auxiliary special purpose array          !
362     !                       (*** not used in this version of code)         !
363     !     aerosols(NPTS,NLAY,NBANDS,:) - aerosol optical properties        !
364     !                       (check module_radiation_aerosols for definition!
365     !        (:,:,:,1)          - optical depth                            !
366     !        (:,:,:,2)          - single scattering albedo                 !
367     !        (:,:,:,3)          - asymmetry parameter                      !
368     !     sfemis (NPTS)         - surface emissivity                       !
369     !     sfgtmp (NPTS)         - surface ground temperature (k)           !
370     !     NPTS                  - total number of horizontal points        !
371     !     NLAY,NLP1             - total number of vertical layers, levels  !
372     !     iflip                 - control flag for in/out vertical index   !
373     !                             =0: index from toa to surface            !
374     !                             =1: index from surface to toa            !
375     !     lprnt                 - cntl flag for diagnostic print out       !
376     !                                                                      !
377     !  control parameters in module "module_radlw_cntr_para":              !
378     !     ilwrate               - heating rate unit selections             !
379     !                             =1: output in k/day                      !
380     !                             =2: output in k/second                   !
381     !     iaerlw                - control flag for aerosols                !
382     !                             =0: do not include aerosol effect        !
383     !                             >0: include aerosol effect               !
384     !     irgaslw               - control flag for rare gases              !
385     !                             (ch4,n2o,o2, etc.)                       !
386     !                             =0: do not include rare gases            !
387     !                             =1: include all rare gases               !
388     !     icfclw                - control flag for cfc gases               !
389     !                             =0: do not include cfc gases             !
390     !                             =1: include all cfc gases                !
391     !     iflagliq              - liq-cloud optical properties contrl flag !
392     !                             =0: input cld opt dep, ignor iflagice    !
393     !                             =1: input cwp,cip, (ccm2) ignor iflagice !
394     !                             =2: input cwp rew, (ccm3 method)         !
395     !                             =3: input cwp rew, hu and stamnes (1993) !
396     !     iflagice              - ice-cloud optical properties contrl flag !
397     !                       * * * if iflagliq .lt. 2, iflafice is ignored  !
398     !                             =0: input cip rei, (ccm3 method)         !
399     !                             =1: input cip rei, ebert and curry (1997)!
400     !                             =2: input cip rei, streamer (1996)       !
401     !                                                                      !
402     !  output variables:                                                   !
403     !     hlwc   (NPTS,NLAY)    - total sky heating rate (k/day or k/sec)  !
404     !     topflx (NPTS)         - radiation fluxes at top, component:      !
405     !                        (check module_radlw_paramters for definition) !
406     !        upfxc                 total sky upward flux at top (w/m2)     !
407     !        upfx0                 clear sky upward flux at top (w/m2)     !
408     !     sfcflx (NPTS)         - radiation fluxes at sfc, component:      !
409     !                        (check module_radlw_paramters for definition) !
410     !        upfxc                 total sky upward flux at sfc (w/m2)     !
411     !        upfx0                 clear sky upward flux at sfc (w/m2)     !
412     !        dnfxc                 total sky downward flux at sfc (w/m2)   !
413     !        dnfx0                 clear sky downward flux at sfc (w/m2)   !
414     !                                                                      !
415     !! optional output variables:                                          !
416     !     hlwb(NPTS,NLAY,NBANDS)- spectral band total sky heating rates    !
417     !     hlw0   (NPTS,NLAY)    - clear sky heating rate (k/day or k/sec)  !
418     !     flxprf (NPTS,NLP1)    - level radiative fluxes (w/m2), components!
419     !                        (check module_radlw_paramters for definition) !
420     !        upfxc                 total sky upward flux                   !
421     !        dnfxc                 total sky dnward flux                   !
422     !        upfx0                 clear sky upward flux                   !
423     !        dnfx0                 clear sky dnward flux                   !
424     !                                                                      !
425     !  module parameters, control and local variables:                     !
426     !     NBANDS                - number of longwave spectral bands        !
427     !     MAXGAS                - maximum number of absorbing gaseous      !
428     !     MAXXSEC               - maximum number of cross-sections         !
429     !     NGPTLW                - total number of g-point subintervals     !
430     !     NGnn   (nn=1-16)      - number of g-points in band nn            !
431     !     nspa,nspb(NBANDS)     - number of lower/upper ref atm's per band !
432     !     delwave(NBANDS)       - longwave band width (wavenumbers)        !
433     !     bpade                 - pade approximation constant (1/0.278)    !
434     !     iovrlw                - cloud overlapping control flag           !
435     !                             =0: random overlapping clouds            !
436     !                             =1: maximum/random overlapping clouds    !
437     !     pavel  (NLAY)         - layer pressures (mb)                     !
438     !     delp   (NLAY)         - layer pressure thickness (mb)            !
439     !     tavel  (NLAY)         - layer temperatures (k)                   !
440     !     tz     (0:NLAY)       - level (interface) temperatures (k)       !
441     !     semiss (NBANDS)       - surface emissivity for each band         !
442     !     wx     (NLAY,MAXXSEC) - cross-section molecules concentration    !
443     !     coldry (NLAY)         - dry air column amount                    !
444     !                                   (1.e-20*molecules/cm**2)           !
445     !     cldfrac(0:NLP1)       - layer cloud fraction                     !
446     !     taucloud(NBANDS,NLAY) - layer cloud optical depth for each band  !
447     !     taug   (NGPTLW,NLAY)  - gaseous optical depths                   !
448     !     pfrac  (NGPTLW,NLAY)  - planck fractions                         !
449     !     itr    (NGPTLW,NLAY)  - integer look-up table index              !
450     !     colamt (NLAY,MAXGAS)  - column amounts of absorbing gases        !
451     !                             1-MAXGAS are for watervapor, carbon      !
452     !                             dioxide, ozone, nitrous oxide, methane,  !
453     !                             oxigen, carbon monoxide, respectively    !
454     !                             (molecules/cm**2)                        !
455     !     pwvcm                 - column precipitable water vapor (cm)     !
456     !     secdiff(NBANDS)       - variable diffusivity angle defined as    !
457     !                             an exponential function of the column    !
458     !                             water amount in bands 2-3 and 5-9.       !
459     !                             this reduces the bias of several w/m2 in !
460     !                             downward surface flux in high water      !
461     !                             profiles caused by using the constant    !
462     !                             diffusivity angle of 1.66.         (mji) !
463     !     co2mult(NLAY)         - the factor used to multiply the ave co2  !
464     !                             abs coeff to get the added contribution  !
465     !                             to the optical depth relative to 355 ppm.!
466     !     facij  (NLAY)         - indicator of interpolation factors       !
467     !                             =0/1: indicate lower/higher temp & height!
468     !     selffac(NLAY)         - scale factor for self-continuum, equals  !
469     !                          (w.v. density)/(atm density at 296K,1013 mb)!
470     !     selffrac(NLAY)        - factor for temp interpolation of ref     !
471     !                             self-continuum data                      !
472     !     indself(NLAY)         - index of the lower two appropriate ref   !
473     !                             temp for the self-continuum interpolation!
474     !     laytrop,layswtch,laylow                                          !
475     !                           - layer at which switch is made from one   !
476     !                             combination of key species to another    !
477     !     totuflux(0:NLAY)      - upward longwave flux (w/m2)              !
478     !     totdflux(0:NLAY)      - downward longwave flux (w/m2)            !
479     !     totuclfl(0:NLAY)      - clear-sky upward longwave flux (w/m2)    !
480     !     totdclfl(0:NLAY)      - clear-sky downward longwave flux (w/m2)  !
481     !     fnet    (0:NLAY)      - net longwave flux (w/m2)                 !
482     !     fnetc   (0:NLAY)      - clear-sky net longwave flux (w/m2)       !
483     !                                                                      !
484     !                                                                      !
485     !  =====================    end of definitions    ===================  !
486     !
487           implicit none
488     
489     !  ---  inputs:
490           integer,  intent(in) :: NPTS, NLAY, NLP1, iflip, iauxil(:)
491     
492           logical,  intent(in) :: lprnt
493     
494           real (kind=kind_phys), dimension(:,:),  intent(in) :: pint, tint, &
495          &       pmid, tmid, qnm, o3mr
496     
497           real (kind=kind_phys), dimension(:,:,:),intent(in) :: gasvmr,     &
498          &       clouds
499     
500           real (kind=kind_phys), dimension(:),    intent(in) :: sfemis,     &
501          &       sfgtmp
502     
503           real (kind=kind_phys), dimension(:,:,:,:),intent(in) :: aerosols
504     
505     !  ---  outputs:
506           real (kind=kind_phys), dimension(:,:), intent(out) :: hlwc
507     
508           type (topflw_type),    dimension(:),   intent(out) :: topflx
509           type (sfcflw_type),    dimension(:),   intent(out) :: sfcflx
510     
511     !! ---  optional outputs:
512           real (kind=kind_phys),dimension(:,:,:),optional,intent(out):: hlwb
513           real (kind=kind_phys),dimension(:,:),optional,intent(out):: hlw0
514           type (proflw_type),   dimension(:,:),optional,intent(out):: flxprf
515     
516     !  ---  locals:
517           real (kind=kind_phys), dimension(0:NLP1) :: cldfrac
518     
519           real (kind=kind_phys), dimension(0:NLAY) :: totuflux, totdflux,   &
520          &       totuclfl, totdclfl, tz
521     
522           real (kind=kind_phys), dimension(NLAY)   :: htr, htrcl
523     
524           real (kind=kind_phys), dimension(NLAY)   :: pavel, tavel, delp,   &
525          &       taucl, cwp1, cip1, rew1, rei1, cda1, cda2, cda3, cda4,     &
526          &       coldry, co2mult, h2ovmr, o3vmr, fac00, fac01, fac10,       &
527          &       fac11, forfac, plog, selffac, selffrac, temcol
528     
529           real (kind=kind_phys) :: colamt(NLAY,MAXGAS), wx(NLAY,MAXXSEC),   &
530          &       taucloud(NBANDS,NLAY), pfrac(NGPTLW,NLAY), semiss(NBANDS), &
531          &       secdiff(NBANDS), tauaer(NBANDS,NLAY), htrb(NLAY,NBANDS)
532     
533           real (kind=kind_phys) :: fp, ft, ft1, tem0, tem1, tem2, pwvcm,    &
534          &       stemp
535     
536           integer, dimension(NLAY) :: jp, jt, jt1, indself
537           integer                  :: itr(NGPTLW,NLAY), laytrop, layswtch,  &
538          &                            laylow, jp1, j, k, k1, iplon
539     !
540     !===> ... begin here
541     !
542     
543           lhlwb  = present ( hlwb )
544           lhlw0  = present ( hlw0 )
545           lflxprf= present ( flxprf )
546     
547     !  ---  loop over horizontal NPTS profiles
548     
549           lab_do_iplon : do iplon = 1, NPTS
550     
551             if (sfemis(iplon) > eps .and. sfemis(iplon) <= 1.0) then  ! input surface emissivity
552               do j = 1, NBANDS
553                 semiss(j) = sfemis(iplon)
554               enddo
555             else                                                      ! use default values
556               do j = 1, NBANDS
557                 semiss(j) = semiss0(j)
558               enddo
559             endif
560     
561             stemp = sfgtmp(iplon)          ! surface ground temp
562     
563     !  ---  prepare atmospheric profile for use in rrtm
564     !       the vertical index of internal array is from surface to top
565     
566             if (iflip == 0) then        ! input from toa to sfc
567     
568               tem1 = 100.0 * con_g
569               tem2 = 1.0e-20 * 1.0e3 * con_avgd
570               tz(0) = tint(iplon,NLP1)
571     
572               do k = 1, NLAY
573                 k1 = NLP1 - k
574                 pavel(k)= pmid(iplon,k1)
575                 delp(k) = pint(iplon,k1+1) - pint(iplon,k1)
576                 tavel(k)= tmid(iplon,k1)
577                 tz(k)   = tint(iplon,k1)
578     
579     !  ---  set absorber amount
580     !test use
581     !           h2ovmr(k)= max(f_zero,qnm(iplon,k1)*amdw)                   ! input mass mixing ratio
582     !           h2ovmr(k)= max(f_zero,qnm(iplon,k1))                        ! input vol mixing ratio
583     !ncep model use
584                 h2ovmr(k)= max(f_zero,                                      &
585          &                     qnm(iplon,k1)*amdw/(1.0-qnm(iplon,k1)))      ! input specific humidity
586                 o3vmr (k)= max(f_zero,o3mr(iplon,k1)*amdo3)                 ! input mass mixing ratio
587     !test use   o3vmr (k)= max(f_zero,o3mr(iplon,k1))                       ! input vol mixing ratio
588     
589                 tem0 = (1.0 - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
590                 coldry(k) = tem2 * delp(k) / (tem1*tem0*(1.0 + h2ovmr(k)))
591                 temcol(k) = 1.0e-12 * coldry(k)
592     
593                 colamt(k,1) =                coldry(k)*h2ovmr(k)            ! h2o
594                 colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(iplon,k1,1))  ! co2
595                 colamt(k,3) =                coldry(k)*o3vmr(k)             ! o3
596               enddo
597     
598     !  ---  set aerosol optical properties
599     
600               if (iaerlw > 0) then
601                 do k = 1, NLAY
602                   k1 = NLP1 - k
603                   do j = 1, NBANDS
604                     tauaer(j,k) = aerosols(iplon,k1,j,1)                    &
605          &                      * (1.0 - aerosols(iplon,k1,j,2))
606                   enddo
607                 enddo
608               else
609                 tauaer(:,:) = f_zero
610               endif
611     
612               if (iflagliq > 0) then   ! use prognostic cloud method
613                 do k = 1, NLAY
614                   k1 = NLP1 - k
615                   cldfrac(k)= clouds(iplon,k1,1)
616                   cwp1 (k)  = clouds(iplon,k1,2)
617                   rew1 (k)  = clouds(iplon,k1,3)
618                   cip1 (k)  = clouds(iplon,k1,4)
619                   rei1 (k)  = clouds(iplon,k1,5)
620                   cda1 (k)  = clouds(iplon,k1,6)
621                   cda2 (k)  = clouds(iplon,k1,7)
622                   cda3 (k)  = clouds(iplon,k1,8)
623                   cda4 (k)  = clouds(iplon,k1,9)
624                 enddo
625               else                        ! use diagnostic cloud method
626                 do k = 1, NLAY
627                   k1 = NLP1 - k
628                   cldfrac(k)= clouds(iplon,k1,1)
629                   cda1(k)   = clouds(iplon,k1,2)
630                 enddo
631               endif                       ! end if_iflagliq
632     
633               cldfrac(0)    = 1.0         ! padding value only
634               cldfrac(NLP1) = f_zero      ! padding value only
635     
636             else                        ! input from sfc to toa
637     
638               tem1 = 100.0 * con_g
639               tem2 = 1.0e-20 * 1.0e3 * con_avgd
640               tz(0) = tint(iplon,1)
641     
642               do k = 1, NLAY
643                 pavel(k)= pmid(iplon,k)
644                 delp(k) = pint(iplon,k) - pint(iplon,k+1)
645                 tavel(k)= tmid(iplon,k)
646                 tz(k)   = tint(iplon,k+1)
647     
648     !  ---  set absorber amount
649     !test use
650     !           h2ovmr(k)= max(f_zero,qnm(iplon,k)*amdw)                    ! input mass mixing ratio
651     !           h2ovmr(k)= max(f_zero,qnm(iplon,k))                         ! input vol mixing ratio
652     !ncep model use
653                 h2ovmr(k)= max(f_zero,qnm(iplon,k)*amdw/(1.0-qnm(iplon,k))) ! input specific humidity
654                 o3vmr (k)= max(f_zero,o3mr(iplon,k)*amdo3)                  ! input mass mixing ratio
655     !test use   o3vmr (k)= max(f_zero,o3mr(iplon,k))                        ! input vol mixing ratio
656     
657                 tem0 = (1.0 - h2ovmr(k))*con_amd + h2ovmr(k)*con_amw
658                 coldry(k) = tem2 * delp(k) / (tem1*tem0*(1.0 + h2ovmr(k)))
659                 temcol(k) = 1.0e-12 * coldry(k)
660     
661                 colamt(k,1) =                coldry(k)*h2ovmr(k)           ! h2o
662                 colamt(k,2) = max(temcol(k), coldry(k)*gasvmr(iplon,k,1))  ! co2
663                 colamt(k,3) =                coldry(k)*o3vmr(k)            ! o3
664               enddo
665     
666     !  ---  set aerosol optical properties
667     
668               if (iaerlw > 0) then
669                 do k = 1, NLAY
670                   do j = 1, NBANDS
671                     tauaer(j,k) = aerosols(iplon,k,j,1)                     &
672          &                      * (1.0 - aerosols(iplon,k,j,2))
673                   enddo
674                 enddo
675               else
676                 tauaer(:,:) = f_zero
677               endif
678     
679               if (iflagliq > 0) then   ! use prognostic cloud method
680                 do k = 1, NLAY
681                   cldfrac(k)= clouds(iplon,k,1)
682                   cwp1 (k)  = clouds(iplon,k,2)
683                   rew1 (k)  = clouds(iplon,k,3)
684                   cip1 (k)  = clouds(iplon,k,4)
685                   rei1 (k)  = clouds(iplon,k,5)
686                   cda1 (k)  = clouds(iplon,k,6)
687                   cda2 (k)  = clouds(iplon,k,7)
688                   cda3 (k)  = clouds(iplon,k,8)
689                   cda4 (k)  = clouds(iplon,k,9)
690                 enddo
691               else                        ! use diagnostic cloud method
692                 do k = 1, NLAY
693                   cldfrac(k)= clouds(iplon,k,1)
694                   cda1(k)   = clouds(iplon,k,2)
695                 enddo
696               endif
697     
698               cldfrac(0)    = 1.0         ! padding value only
699               cldfrac(NLP1) = f_zero      ! padding value only
700     
701             endif                       ! if_iflip
702     
703     !  ---  set up col amount for rare gases, convert from volume mixing ratio to
704     !       molec/cm2 based on coldry (scaled to 1.0e-20) for use in rrtm
705     
706             if (iflip == 0) then        ! input from toa to sfc
707     
708               if (irgaslw == 1) then
709                 do k = 1, NLAY
710                   k1 = NLP1 - k
711                   colamt(k,4)=max(temcol(k), coldry(k)*gasvmr(iplon,k1,2))  ! n2o
712                   colamt(k,5)=max(temcol(k), coldry(k)*gasvmr(iplon,k1,3))  ! ch4
713     !             colamt(k,6)=max(f_zero,    coldry(k)*gasvmr(iplon,k1,4))  ! o2 - not used
714     !             colamt(k,7)=max(f_zero,    coldry(k)*gasvmr(iplon,k1,5))  ! co - not used
715                 enddo
716               else
717                 do k = 1, NLAY
718                   colamt(k,4) = f_zero     ! n2o
719                   colamt(k,5) = f_zero     ! ch4
720     !             colamt(k,6) = f_zero     ! o2 - not used
721     !             colamt(k,7) = f_zero     ! co - not used
722                 enddo
723               endif
724     
725               if (icfclw == 1) then
726                 do k = 1, NLAY
727                   k1 = NLP1 - k
728                   wx(k,1) = max( f_zero, coldry(k)*gasvmr(iplon,k1,9) )   ! ccl4
729                   wx(k,2) = max( f_zero, coldry(k)*gasvmr(iplon,k1,6) )   ! cf11
730                   wx(k,3) = max( f_zero, coldry(k)*gasvmr(iplon,k1,7) )   ! cf12
731                   wx(k,4) = max( f_zero, coldry(k)*gasvmr(iplon,k1,8) )   ! cf22
732                 enddo
733               else
734                 wx(:,:) = f_zero
735               endif
736     
737     ! mji - for variable diffusivity angle, sum moist atmos and water over column
738               tem1 = f_zero
739               tem2 = f_zero
740               do k = 1, NLAY
741                 tem1 = tem1 + coldry(k) + colamt(k,1)
742                 tem2 = tem2 + colamt(k,1)
743               enddo
744     ! mji - calculate column precipitable water and variable diffusivity angle
745     !         tem3 = tem2 / (amdw * tem1)
746     !         pwvcm = tem3 * (1.0e3 * pint(iplon,NLP1)) / (1.0e2 * con_g)
747               pwvcm = 10.0 * pint(iplon,NLP1) * tem2 /(amdw * tem1* con_g)
748     
749             else                        ! input from sfc to toa
750     
751               if (irgaslw == 1) then
752                 do k = 1, NLAY
753                   colamt(k,4)=max(temcol(k), coldry(k)*gasvmr(iplon,k,2))  ! n2o
754                   colamt(k,5)=max(temcol(k), coldry(k)*gasvmr(iplon,k,3))  ! ch4
755     !             colamt(k,6)=max(f_zero,    coldry(k)*gasvmr(iplon,k,4))  ! o2 - not used
756     !             colamt(k,7)=max(f_zero,    coldry(k)*gasvmr(iplon,k,5))  ! co - not used
757                 enddo
758               else
759                 do k = 1, NLAY
760                   colamt(k,4) = f_zero     ! n2o
761                   colamt(k,5) = f_zero     ! ch4
762     !             colamt(k,6) = f_zero     ! o2 - not used
763     !             colamt(k,7) = f_zero     ! co - not used
764                 enddo
765               endif
766     
767               if (icfclw == 1) then
768                 do k = 1, NLAY
769                   wx(k,1) = max( f_zero, coldry(k)*gasvmr(iplon,k,9) )   ! ccl4
770                   wx(k,2) = max( f_zero, coldry(k)*gasvmr(iplon,k,6) )   ! cf11
771                   wx(k,3) = max( f_zero, coldry(k)*gasvmr(iplon,k,7) )   ! cf12
772                   wx(k,4) = max( f_zero, coldry(k)*gasvmr(iplon,k,8) )   ! cf22
773                 enddo
774               else
775                 wx(:,:) = f_zero
776               endif
777     
778     ! mji - for variable diffusivity angle, sum moist atmos and water over column
779               tem1 = f_zero
780               tem2 = f_zero
781               do k = 1, NLAY
782                 tem1 = tem1 + coldry(k) + colamt(k,1)
783                 tem2 = tem2 + colamt(k,1)
784               enddo
785     ! mji - calculate column precipitable water and variable diffusivity angle
786     !         tem3 = tem2 / (amdw * tem1)
787     !         pwvcm = tem3 * (1.0e3 * pint(iplon,1)) / (1.0e2 * con_g)
788               pwvcm = 10.0 * pint(iplon,1) * tem2 /(amdw * tem1* con_g)
789     
790             endif                       ! if_iflip
791     
792             do j = 1, NBANDS
793               if (j==1 .or. j==4 .or. j==10) then
794                 secdiff(j) = 1.66
795               else
796                 secdiff(j) = a0(j) + a1(j) * exp( a2(j)*pwvcm )
797               endif
798             enddo
799             if (pwvcm < 1.0) secdiff(6) = 1.80
800             if (pwvcm > 7.1) secdiff(7) = 1.50
801     ! mji
802     
803             do k = 1, NLAY
804     !  ...  using e = 1334.2 cm-1.
805               tem1 = co2fac * coldry(k)
806               co2mult(k) = (colamt(k,2) - tem1) * 272.63                    &
807          &               * exp(-1919.4/tavel(k)) / (8.7604e-4*tavel(k))
808               forfac(k)  = pavel(k)*stpfac / (tavel(k)*(1.0 + h2ovmr(k)))
809             enddo
810             
811     !     if (lprnt .and. ipt == iplon) then
812     !     print *,'  coldry',coldry
813     !     print *,' wx(*,1) ',(wx(k,1),k=1,NLAY)
814     !     print *,' wx(*,2) ',(wx(k,2),k=1,NLAY)
815     !     print *,' wx(*,3) ',(wx(k,3),k=1,NLAY)
816     !     print *,' wx(*,4) ',(wx(k,4),k=1,NLAY)
817     !     print *,' iplon ',iplon
818     !     print *,'  pavel ',pavel
819     !     print *,'  delp ',delp
820     !     print *,'  tavel ',tavel
821     !     print *,'  tz ',tz
822     !     print *,' h2ovmr ',h2ovmr
823     !     print *,' o3vmr ',o3vmr
824     !     endif
825     
826     !  ---  calculate cloud optical properties
827     
828             call cldprop                                                    &
829     !  ---  inputs:
830          &     ( cldfrac, cwp1, cip1, rew1, rei1, cda1, cda2, cda3, cda4,   &
831          &       NLAY, NLP1,                                                &
832     !  ---  output:
833          &       taucloud                                                   &
834          &     )
835     
836     !     if (lprnt) then
837     !     print *,' after cldprop'
838     !     print *,' cwp1',cwp1
839     !     print *,' cip1',cip1
840     !     print *,' rew1',rew1
841     !     print *,' rei1',rei1
842     !     print *,' taucl',cda1
843     !     print *,' cldfrac',cldfrac
844     !     print *,' taucloud',taucloud
845     !     endif
846     
847     !  ---  calculate information needed by the radiative transfer routine
848     !       that is specific to this atmosphere, especially some of the 
849     !       coefficients and indices needed to compute the optical depths
850     !       by interpolating data from stored reference atmospheres. 
851     
852             laytrop = 0
853             layswtch= 0
854             laylow  = 0
855     
856             do k = 1, NLAY
857     
858     !  ---  find the two reference pressures on either side of the
859     !       layer pressure.  store them in jp and jp1.  store in fp the
860     !       fraction of the difference (in ln(pressure)) between these
861     !       two values that the layer pressure lies.
862     
863               plog(k) = log(pavel(k))
864               jp(k)= max(1, min(58, int(36.0 - 5.0*(plog(k)+0.04)) ))
865               jp1  = jp(k) + 1
866     !  ---  limit pressure extrapolation at the top
867               fp   = max(f_zero, min(1.0, 5.0*(preflog(jp(k))-plog(k)) ))
868     !org      fp   = 5.0 * (preflog(jp(k)) - plog(k))
869     
870     !  ---  determine, for each reference pressure (jp and jp1), which
871     !       reference temperature (these are different for each
872     !       reference pressure) is nearest the layer temperature but does
873     !       not exceed it.  store these indices in jt and jt1, resp.
874     !       store in ft (resp. ft1) the fraction of the way between jt
875     !       (jt1) and the next highest reference temperature that the
876     !       layer temperature falls.
877     
878               tem1 = (tavel(k) - tref(jp(k))) / 15.0
879               tem2 = (tavel(k) - tref(jp1  )) / 15.0
880               jt (k) = max(1, min(4, int(3.0 + tem1) ))
881               jt1(k) = max(1, min(4, int(3.0 + tem2) ))
882     !  ---  restrict extrapolation ranges by limiting abs(det t) < 37.5 deg
883               ft  = max(-0.5, min(1.5, tem1 - float(jt (k) - 3) ))
884               ft1 = max(-0.5, min(1.5, tem2 - float(jt1(k) - 3) ))
885     !org      ft  = tem1 - float(jt (k) - 3)
886     !org      ft1 = tem2 - float(jt1(k) - 3)
887     
888     !  ---  we have now isolated the layer ln pressure and temperature,
889     !       between two reference pressures and two reference temperatures
890     !       (for each reference pressure).  we multiply the pressure
891     !       fraction fp with the appropriate temperature fractions to get
892     !       the factors that will be needed for the interpolation that yields
893     !       the optical depths (performed in routines taugbn for band n).
894     
895               fac10(k) = (1.0 - fp) * ft
896               fac00(k) = (1.0 - fp) * (1.0 - ft)
897               fac11(k) = fp * ft1
898               fac01(k) = fp * (1.0 - ft1)
899     
900             enddo
901     
902     !  ---  if the pressure is less than ~100mb, perform a different
903     !       set of species interpolations.
904     
905             do k = 1, NLAY
906               if (plog(k) > 4.56) then
907                 laytrop =  laytrop + 1
908     
909     !  ---  for one band, the "switch" occurs at ~300 mb.
910                 if (plog(k) >= 5.76) layswtch = layswtch + 1
911                 if (plog(k) >= 6.62) laylow   = laylow   + 1
912     
913     !  ---  set up factors (tem1) needed to separately include the water
914     !       vapor self-continuum in the calculation of absorption
915     !       coefficient.
916     
917                 tem1 = (tavel(k) - 188.0) / 7.2
918     
919                 selffac(k) = h2ovmr(k) * forfac(k)
920                 indself(k) = min(9, max(1, int(tem1)-7 ))
921                 selffrac(k)= tem1 - float(indself(k) + 7)
922     
923               else
924     
925                 selffac(k) = f_zero
926                 indself(k) = 0
927                 selffrac(k)= f_zero
928     
929               endif
930             enddo
931     
932     !  ---  set laylow for profiles with surface pressure less than 750mb.
933             if (laylow == 0) laylow = 1
934     
935     !     if (lprnt) then
936     !      print *,'laytrop,layswtch,laylow',laytrop,layswtch,laylow
937     !      print *,'colh2o',(colamt(k,1),k=1,NLAY)
938     !      print *,'colco2',(colamt(k,2),k=1,NLAY)
939     !      print *,'colo3', (colamt(k,3),k=1,NLAY)
940     !      print *,'coln2o',(colamt(k,4),k=1,NLAY)
941     !      print *,'colch4',(colamt(k,5),k=1,NLAY)
942     !      print *,'co2mult',co2mult
943     !      print *,'fac00',fac00
944     !      print *,'fac01',fac01
945     !      print *,'fac10',fac10
946     !      print *,'fac11',fac11
947     !      print *,'jp',jp
948     !      print *,'jt',jt
949     !      print *,'jt1',jt1
950     !      print *,'selffac',selffac
951     !      print *,'selffrac',selffrac
952     !      print *,'indself',indself
953     !      print *,'forfac',forfac
954     !     endif
955     
956             call taumol                                                     &
957     !  ---  inputs:
958          &     ( laytrop,layswtch,laylow,h2ovmr,colamt,wx,co2mult,          &
959          &       fac00,fac01,fac10,fac11,jp,jt,jt1,selffac,selffrac,        &
960          &       indself,forfac,secdiff,tauaer, NLAY,                       &
961     !  ---  outputs:
962          &       itr, pfrac                                                 &
963          &     )
964     
965     !     if (lprnt) then
966     !     print *,' after taumol'
967     !     do k=1,NLAY
968     !       write(6,123) k
969     !123    format(' k =',i3,5x,'PFRAC')
970     !       write(6,122) (pfrac(j,k),j=1,NGPTLW)
971     !122    format(10e14.7)
972     !       write(6,124) k
973     !124    format(' k =',i3,5x,'ITR')
974     !       write(6,125) (itr(j,k),j=1,NGPTLW)
975     !125    format(10i10)
976     !     enddo
977     !     endif
978     
979     !  ---  call the radiative transfer routine.
980     
981             if (iovrlw == 0) then
982     
983               call rtrn                                                     &
984     !  ---  inputs:
985          &     ( tavel,tz,delp,semiss,cldfrac,taucloud,pfrac,               &
986          &       secdiff, stemp, itr, NLAY, NLP1,                           &
987     !  ---  outputs:
988          &       totuflux,totdflux,htr, totuclfl,totdclfl,htrcl, htrb       &
989          &     )
990     
991             else
992     
993               call rtrnmr                                                   &
994     !  ---  inputs:
995          &     ( tavel,tz,delp,semiss,cldfrac,taucloud,pfrac,               &
996          &       secdiff, stemp, itr, NLAY, NLP1,                           &
997     !  ---  outputs:
998          &       totuflux,totdflux,htr, totuclfl,totdclfl,htrcl, htrb       &
999          &     )
1000     
1001             endif   ! end if_iovrlw_block
1002     
1003     
1004     !  ---  output total-sky and clear-sky fluxes and heating rates.
1005     
1006             topflx(iplon)%upfxc = totuflux(NLAY)
1007             topflx(iplon)%upfx0 = totuclfl(NLAY)
1008     
1009             sfcflx(iplon)%upfxc = totuflux(0)
1010             sfcflx(iplon)%upfx0 = totuclfl(0)
1011             sfcflx(iplon)%dnfxc = totdflux(0)
1012             sfcflx(iplon)%dnfx0 = totdclfl(0)
1013     
1014             if (iflip == 0) then        ! output from toa to sfc
1015     
1016     !! ---  optional fluxes
1017               if ( lflxprf ) then
1018                 do k = 0, NLAY
1019                   k1 = NLP1 - k
1020                   flxprf(iplon,k1)%upfxc = totuflux(k)
1021                   flxprf(iplon,k1)%dnfxc = totdflux(k)
1022                   flxprf(iplon,k1)%upfx0 = totuclfl(k)
1023                   flxprf(iplon,k1)%dnfx0 = totdclfl(k)
1024                 enddo
1025               endif
1026     
1027               do k = 1, NLAY
1028                 k1 = NLP1 - k
1029                 hlwc(iplon,k1) = htr(k)
1030               enddo
1031     
1032     !! ---  optional clear sky heating rate
1033               if ( lhlw0 ) then
1034                 do k = 1, NLAY
1035                   k1 = NLP1 - k
1036                   hlw0(iplon,k1) = htrcl(k)
1037                 enddo
1038               endif
1039     
1040     !! ---  optional spectral band heating rate
1041               if ( lhlwb ) then
1042                 do j = 1, NBANDS
1043                 do k = 1, NLAY
1044                   k1 = NLP1 - k
1045                   hlwb(iplon,k1,j) = htrb(k,j)
1046                 enddo
1047                 enddo
1048               endif
1049     
1050             else                        ! output from sfc to toa
1051     
1052     !! ---  optional fluxes
1053               if ( lflxprf ) then
1054                 do k = 0, NLAY
1055                   flxprf(iplon,k+1)%upfxc = totuflux(k)
1056                   flxprf(iplon,k+1)%dnfxc = totdflux(k)
1057                   flxprf(iplon,k+1)%upfx0 = totuclfl(k)
1058                   flxprf(iplon,k+1)%dnfx0 = totdclfl(k)
1059                 enddo
1060               endif
1061     
1062               do k = 1, NLAY
1063                 hlwc(iplon,k) = htr(k)
1064               enddo
1065     
1066     !! ---  optional clear sky heating rate
1067               if ( lhlw0 ) then
1068                 do k = 1, NLAY
1069                   hlw0(iplon,k) = htrcl(k)
1070                 enddo
1071               endif
1072     
1073     !! ---  optional spectral band heating rate
1074               if ( lhlwb ) then
1075                 do j = 1, NBANDS
1076                 do k = 1, NLAY
1077                   hlwb(iplon,k,j) = htrb(k,j)
1078                 enddo
1079                 enddo
1080               endif
1081     
1082             endif                       ! if_iflip
1083     
1084           enddo  lab_do_iplon
1085     
1086           return
1087     !...................................
1088           end subroutine lwrad
1089     !-----------------------------------
1090     
1091     
1092     
1093     !-----------------------------------
1094           subroutine rlwinit                                                &
1095     !...................................
1096     
1097     !  ---  inputs:
1098          &     ( icwp, me, NLAY, iovr, isubc )
1099     !  ---  outputs: (none)
1100     
1101     !  *******************************************************************  !
1102     !                                                                       !
1103     !  rrtm longwave radiative transfer model                               !
1104     !  atmospheric and environmental research, inc., cambridge, ma          !
1105     !                                                                       !
1106     !                                                                       !
1107     !  original version:       michael j. iacono; july, 1998                !
1108     !  revision for ncar ccm:  michael j. iacono; september, 1998           !
1109     !                                                                       !
1110     !  this subroutine performs calculations necessary for the initialization
1111     !  of the lw model, rrtm.  lookup tables are computed for use in the lw !
1112     !  radiative transfer, and input absorption coefficient data for each   !
1113     !  spectral band are reduced from 256 g-points to 140 for use in rrtm.  !
1114     !                                                                       !
1115     !  *******************************************************************  !
1116     !                                                                       !
1117     ! definitions:                                                          !
1118     !     arrays for 5000-point look-up tables:                             !
1119     !     tau  - clear-sky optical depth (used in cloudy radiative transfer)!
1120     !     tf     tau transition function; i.e. the transition of the planck !
1121     !            function from that for the mean layer temperature to that  !
1122     !            for the layer boundary temperature as a function of optical!
1123     !            depth. the "linear in tau" method is used to make the table!
1124     !     trans- transmittance                                              !
1125     !                                                                       !
1126     !  *******************************************************************  !
1127     !                                                                       !
1128     !  inputs:                                                              !
1129     !    icwp     -  flag of cloud schemes used by model                    !
1130     !                =0: diagnostic scheme gives cloud tau, omiga, and g    !
1131     !                =1: prognostic scheme gives cloud liq/ice path, etc.   !
1132     !    me       - print control for parallel process                      !
1133     !    NLAY     - number of vertical layers                               !
1134     !    iovr     - cloud overlapping control flag                          !
1135     !                =0: random overlapping clouds                          !
1136     !                =1: maximum/random overlapping clouds                  !
1137     !    isubc    - mcica sub-column cloud approximation control flag       !
1138     !                *** not used in this version of code                   !
1139     !                                                                       !
1140     !  outputs: (none)                                                      !
1141     !                                                                       !
1142     !  control flags in module "module_radlw_cntr_para":                    !
1143     !     ilwrate - heating rate unit selections                            !
1144     !               =1: output in k/day                                     !
1145     !               =2: output in k/second                                  !
1146     !     iaerlw  - control flag for aerosols                               !
1147     !               =0: do not include aerosol effect                       !
1148     !               >0: include aerosol effect                              !
1149     !     irgaslw - control flag for rare gases (ch4,n2o,o2, etc.)          !
1150     !               =0: do not include rare gases                           !
1151     !               =1: include all rare gases                              !
1152     !     icfclw  - control flag for cfc gases                              !
1153     !               =0: do not include cfc gases                            !
1154     !               =1: include all cfc gases                               !
1155     !     iflagliq- cloud optical properties contrl flag                    !
1156     !               =0: input cloud opt depth from diagnostic scheme        !
1157     !               >0: input cwp,cip, and other cloud content parameters   !
1158     !                                                                       !
1159     !  *******************************************************************  !
1160     !
1161           implicit none
1162     !
1163     !  ---  inputs:
1164           integer, intent(in) :: icwp, me, NLAY, iovr, isubc
1165     
1166     !  ---  outputs: none
1167     
1168     !  ---  locals:
1169           real (kind=kind_phys) :: tfn, fp, rtfp, pival, explimit
1170           integer               :: i
1171     !
1172     !===> ... begin here
1173     !
1174           iovrlw  = iovr     ! assign module variable of overlap flag
1175     
1176           if ( iovrlw<0 .or. iovrlw>1 ) then
1177             print *,'  *** Error in specification of cloud overlap flag',   &
1178          &          ' IOVRLW=',iovrlw,' in RLWINIT !!'
1179             stop
1180           endif
1181     
1182           if (me == 0) then
1183             print *,' - Using AER Longwave Radiation, Version: ', VTAGLW
1184     
1185             if (iaerlw > 0) then
1186               print *,'   --- Using input aerosol parameters for LW'
1187             else
1188               print *,'   --- Aerosol effect is NOT included in LW, all'    &
1189          &           ,' internal aerosol parameters are reset to zeros'
1190             endif
1191     
1192             if (irgaslw == 1) then
1193               print *,'   --- Include rare gases N2O, CH4, O2, absorptions',&
1194          &            ' in LW'
1195             else
1196               print *,'   --- Rare gases effect is NOT included in LW'
1197             endif
1198     
1199             if (icfclw == 1) then
1200               print *,'   --- Include CFC gases absorptions in LW'
1201             else
1202               print *,'   --- CFC gases effect is NOT included in LW'
1203             endif
1204     
1205             if ( isubc == 0 ) then
1206               print *,'   --- Using standard grid average clouds, no sub-', &
1207          &            'column clouds approximation'
1208             else
1209               print *,'   --- Sub-column cloud scheme is not available in', &
1210          &            ' this version of code.  Using standard grid',        &
1211          &            ' average of clouds'
1212             endif
1213           endif
1214     
1215     !  --- ...  check cloud flags for consistency
1216     
1217           if ((icwp == 0 .and. iflagliq /= 0) .or.                          &
1218          &    (icwp == 1 .and. iflagliq == 0)) then
1219             print *, ' *** Model cloud scheme inconsistent with LW',        &
1220          &           ' radiation cloud radiative property setup !!'
1221             stop
1222           endif
1223     
1224     !  --- ...  setup default surface emissivity for each band here
1225     
1226           semiss0(:) = 1.0
1227     
1228     !  --- ...  setup constant factors for flux and heating rate
1229     !           the 1.0e-2 is to convert pressure from mb to N/m**2
1230     
1231           pival = 2.0*asin(1.0)
1232           fluxfac = pival * 2.0d4
1233     !     fluxfac = 62831.85307179586                   ! = 2 * pi * 1.0e4
1234     
1235           if (ilwrate == 1) then
1236     !       heatfac = con_g * 86400. * 1.0e-2 / con_cp  !   (in k/day)
1237             heatfac = con_g * 864.0 / con_cp            !   (in k/day)
1238           else
1239             heatfac = con_g * 1.0e-2 / con_cp           !   (in k/second)
1240           endif
1241     
1242     !  --- ...  compute lookup tables for transmittance, tau transition
1243     !           function, and clear sky tau (for the cloudy sky radiative
1244     !           transfer).  tau is computed as a function of the tau 
1245     !           transition function, transmittance is calculated as a 
1246     !           function of tau, and the tau transition function is 
1247     !           calculated using the linear in tau formulation at values of
1248     !           tau above 0.01.  tf is approximated as tau/6 for tau < 0.01.
1249     !           all tables are computed at intervals of 0.001.  the inverse
1250     !           of the constant used in the pade approximation to the tau
1251     !           transition function is set to b.
1252     
1253           tau  (0) = f_zero
1254           tf   (0) = f_zero
1255           trans(0) = 1.0
1256     
1257           tau  (N5000) = 1.e10
1258           tf   (N5000) = 1.0
1259           trans(N5000) = f_zero
1260     
1261           explimit = aint( -log(tiny(trans(0))) )
1262     
1263           do i = 1, N5000-1
1264              tfn = real(i, kind_phys) / real(N5000-i, kind_phys)
1265              tau  (i) = bpade * tfn
1266              if (tau(i) >= explimit) then
1267                trans(i) = f_zero
1268              else
1269                trans(i) = exp(-tau(i))
1270              endif
1271     
1272              if (tau(i) < 0.1) then
1273                 tf(i) = tau(i) / 6.0
1274              else
1275                 tf(i) = 1. - 2.*( (1./tau(i)) - (trans(i)/(1.-trans(i))) )
1276              endif
1277           enddo
1278     
1279     !  --- ...  calculate lookup tables for functions needed in routine
1280     !           taumol (taugb2)
1281     
1282           corr1(0) = 1.0
1283           corr2(0) = 1.0
1284     
1285           corr1(N200) = 1.0
1286           corr2(N200) = 1.0
1287     
1288           do i = 1, N200-1
1289              fp = 0.005 * float(i)
1290              rtfp = sqrt(fp)
1291              corr1(i) = rtfp / fp
1292              corr2(i) = (1.0 - rtfp) / (1.0 - fp)
1293           enddo
1294     
1295     !...................................
1296           end subroutine rlwinit
1297     !-----------------------------------
1298     
1299     
1300     
1301     !-----------------------------------
1302           subroutine cldprop                                                &
1303     !...................................
1304     
1305     !  ---  inputs:
1306          &     ( cldfrac,cliqp,cicep,reliq,reice,cdat1,cdat2,cdat3,cdat4,   &
1307          &       NLAY, NLP1,                                                &
1308     !  ---  output:
1309          &       taucloud                                                   &
1310          &     )
1311     
1312     !  *******************************************************************  !
1313     !                                                                       !
1314     !    purpose:  compute the cloud optical depth(s) for each cloudy layer.!
1315     !                                                                       !
1316     !  *******************************************************************  !
1317     !                                                                       !
1318     !  inputs:                                                              !
1319     !     cldfrac - layer cloud fraction                                 L  !
1320     !        - - -  for iflagliq > 0  (prognostic cloud sckeme)  - - -      !
1321     !     cliqp   - layer cloud liquid water path  (g/m**2)              L  !
1322     !     reliq   - effective radius for water cloud (micron)            L  !
1323     !     cicep   - layer cloud ice water path  (g/m**2)                 L  !
1324     !     reice   - effective radius for ice cloud (micron)              L  !
1325     !     cdat1   - layer rain drop water path  (g/m**2)                 L  !
1326     !     cdat2   - effective radius for rain drop (microm)              L  !
1327     !     cdat3   - layer snow flake water path (g/m**2)                 L  !
1328     !               (if use fu's formula it needs to be normalized by       !
1329     !                snow density (g/m**3/1.0e6) to get unit of micron)     !
1330     !     cdat4   - effective radius for snow flakes (micron)            L  !
1331     !        - - -  for iflagliq = 0  (diagnostic cloud sckeme)  - - -      !
1332     !     cdat1   - input cloud optical depth                            L  !
1333     !     cdat2   - optional use                                         L  !
1334     !     cdat3   - optional use                                         L  !
1335     !     cdat4   - optional use                                         L  !
1336     !     cliqp   - not used                                             L  !
1337     !     reliq   - not used                                             L  !
1338     !     cicep   - not used                                             L  !
1339     !     reice   - not used                                             L  !
1340     !                                                                       !
1341     !     NLAY/NLP1-vertical layer/level numbers                         1  !
1342     !                                                                       !
1343     !    explanation of the method for each value of iflagliq, and iflagice.!
1344     !    set up in module "module_radlw_cntr_para"                          !
1345     !                                                                       !
1346     !     iflagliq=0 and =1 do not distingish being liquid and ice clouds.  !
1347     !     iflagliq=2 and =3 does distinguish between liquid and ice clouds, !
1348     !                  and requires further user input (iflagice) to specify!
1349     !                  the method to be used to compute the aborption due to!
1350     !                  liquid and ice parts.                                !
1351     !  ...................................................................  !
1352     !                                                                       !
1353     !     iflagliq=0:  for each cloudy layer, the cloud fraction and (gray) !
1354     !                  optical depth are input.                             !
1355     !     iflagliq=1:  for each cloudy layer, the cloud fraction and cloud  !
1356     !                  water path (g/m2) are input.  using clp only. the    !
1357     !                  (gray) cloud optical depth is computed as in ccm2.   !
1358     !     iflagliq=2:  the optical depths due to water clouds are computed  !
1359     !                  as in ccm3.                                          !
1360     !     iflagliq=3:  the water droplet effective radius (microns) is input!
1361     !                  and the opt depths due to water clouds are computed  !
1362     !                  as in hu and stamnes, j., clim., 6, 728-742, (1993). !
1363     !                  the values for absorption coefficients appropriate for
1364     !                  the spectral bands in rrtm have been obtained for a  !
1365     !                  range of effective radii by an averaging procedure   !
1366     !                  based on the work of j. pinto (private communication).
1367     !                  linear interpolation is used to get the absorption   !
1368     !                  coefficients for the input effective radius.         !
1369     !                                                                       !
1370     !     iflagice=0:  the cloud ice path (g/m2) and ice effective radius   !
1371     !                  (microns) are input and the optical depths due to ice!
1372     !                  clouds are computed as in ccm3.                      !
1373     !     iflagice=1:  the cloud ice path (g/m2) and ice effective radius   !
1374     !                  (microns) are input and the optical depths due to ice!
1375     !                  clouds are computed as in ebert and curry, jgr, 97,  !
1376     !                  3831-3836 (1992).  the spectral regions in this work !
1377     !                  have been matched with the spectral bands in rrtm to !
1378     !                  as great an extent as possible:                      !
1379     !                     e&c 1      ib = 5      rrtm bands 9-16            !
1380     !                     e&c 2      ib = 4      rrtm bands 6-8             !
1381     !                     e&c 3      ib = 3      rrtm bands 3-5             !
1382     !                     e&c 4      ib = 2      rrtm band 2                !
1383     !                     e&c 5      ib = 1      rrtm band 1                !
1384     !     iflagice=2:  the cloud ice path (g/m2) and ice effective radius   !
1385     !                  (microns) are input and the optical depths due to ice!
1386     !                  clouds are computed as in streamer (reference: j. key,
1387     !                  streamer user's guide, technical report 96-01,       !
1388     !                  department of geography, boston university, 85 pp.   !
1389     !                  (1996)).  the values of absorption coefficients      !
1390     !                  appropriate for the spectral bands of rrtm were      !
1391     !                  obtained by an averaging procedure based on the work !
1392     !                  of j. pinto (private communication).                 !
1393     !                                                                       !
1394     !  outputs:                                                             !
1395     !     taucloud - cloud optical depth                        NBANDS*L    !
1396     !                                                                       !
1397     !  *******************************************************************  !
1398     !
1399           use module_radlw_cldprlw
1400     
1401           implicit none
1402     
1403     !  ---  inputs:
1404           integer, intent(in) :: NLAY, NLP1
1405     
1406           real (kind=kind_phys), dimension(0:), intent(in) :: cldfrac
1407     
1408           real (kind=kind_phys), dimension(:),  intent(in) :: cliqp, cicep, &
1409          &       reliq, reice, cdat1, cdat2, cdat3, cdat4
1410     
1411     !  ---  outputs:
1412           real (kind=kind_phys), dimension(:,:), intent(out) :: taucloud
1413     
1414     !  ---  locals:
1415           real (kind=kind_phys) :: cliq, cice, radliq, radice, factor, fint
1416           real (kind=kind_phys) :: taurain, tausnow
1417           integer               :: j, k, index
1418     
1419     !
1420     !===> ... begin here
1421     !
1422           do k = 1, NLAY
1423             do j = 1, NBANDS
1424               taucloud(j,k) = f_zero
1425             enddo
1426           enddo
1427     
1428           lab_do_k : do k = 1, NLAY
1429     
1430             lab_if_cld : if (cldfrac(k) > eps) then
1431     
1432     !  ---  ice clouds and water clouds combined.
1433               lab_if_liq : if (iflagliq == 0) then
1434     
1435                 do j = 1, NBANDS
1436                   taucloud(j,k) = cdat1(k)
1437                 enddo
1438     
1439               elseif (iflagliq == 1) then  lab_if_liq
1440     
1441                 taurain = absrain * cdat1(k)                 ! ncar formula
1442                 tausnow = abssnow0 * cdat3(k)                ! ncar formula
1443     !           tausnow = abssnow1 * cdat3(k) / cdat4(k)     ! fu's formula
1444     
1445     !           taucloud(1,k) = absliq1 * (cliqp(k) + cicep(k))
1446     !           taucloud(1,k) = absliq1 * cliqp(k)
1447                 taucloud(1,k) = absliq1*cliqp(k) + taurain + tausnow
1448                 do j = 2,NBANDS
1449                   taucloud(j,k) = taucloud(1,k)
1450                 enddo
1451     
1452     !  ---  separate treatement of ice clouds and water clouds.
1453               else  lab_if_liq
1454     
1455                 taurain = absrain * cdat1(k)                 ! ncar formula
1456                 tausnow = abssnow0 * cdat3(k)                ! ncar formula
1457     !           tausnow = abssnow1 * cdat3(k) / cdat4(k)     ! fu's formula
1458     
1459                 cliq = cliqp(k)
1460                 cice = cicep(k)
1461                 radliq = max(2.5e0, min(60.0e0, real(reliq(k)) ))
1462                 radice = reice(k)
1463     !           radice = max(13.e0, min(130.e0, real(reice(k)) ))
1464     
1465     !  ---  calculation of absorption coefficients due to liquid clouds.
1466                 if (cliq == f_zero) then
1467                   do j = 1, NBANDS
1468                     abscoliq(j) = f_zero
1469                   enddo
1470                 elseif (iflagliq == 2) then
1471                   abscoliq(1) = cliq * absliq2
1472                   do j = 2, NBANDS
1473                     abscoliq(j) = abscoliq(1)
1474                   enddo
1475                 elseif (iflagliq == 3) then
1476                   factor = radliq - 1.5
1477                   index = min(57, int(factor))
1478                   fint = factor - index
1479                   do j = 1, NBANDS
1480                     abscoliq(j) = cliq * (absliq3(index,j) + fint *         &
1481          &                    (absliq3(index+1,j) - (absliq3(index,j))))
1482                   enddo
1483                 endif
1484     
1485     !  ---  calculation of absorption coefficients due to ice clouds.
1486                 if (cice == f_zero) then
1487                   do j = 1, NBANDS
1488                     abscoice(j) = f_zero
1489                   enddo
1490                 elseif (iflagice == 0) then
1491                   abscoice(1) = cice * (absice0(1) + absice0(2)/radice)
1492                   do j = 2, NBANDS
1493                     abscoice(j) = abscoice(1)
1494                   enddo
1495                 elseif (iflagice == 1) then
1496                   do j = 1, NBANDS
1497                     index = ipat(j)
1498                     abscoice(j) = cice * (absice1(1,index)                  &
1499          &                      + absice1(2,index)/radice)
1500                   enddo
1501                 elseif (iflagice == 2) then
1502                   factor = (radice - 10.0) / 3.0
1503                   index = min(39, int(factor))
1504                   fint = factor - index
1505                   do j = 1, NBANDS
1506                     abscoice(j) = cice * (absice2(index,j) + fint *         &
1507          &                    (absice2(index+1,j) - (absice2(index,j))))
1508                   enddo
1509                 endif
1510     
1511                 do j = 1, NBANDS
1512     !             taucloud(j,k) = abscoice(j) + abscoliq(j)
1513                   taucloud(j,k) = abscoice(j) + abscoliq(j)                 &
1514          &                      + taurain + tausnow
1515                 enddo
1516     
1517               endif  lab_if_liq
1518     
1519             endif  lab_if_cld
1520     
1521           enddo  lab_do_k
1522     
1523           return
1524     !...................................
1525           end subroutine cldprop
1526     !-----------------------------------
1527     
1528     
1529     
1530     !-----------------------------------
1531           subroutine rtrn                                                   &
1532     !...................................
1533     
1534     !  ---  inputs:
1535          &     ( tavel,tz,delp,semiss,cldfrac,taucloud,pfrac,               &
1536          &       secdiff, stemp, itr, NLAY, NLP1,                           &
1537     !  ---  outputs:
1538          &       totuflux,totdflux,htr, totuclfl,totdclfl,htrcl, htrb       &
1539          &     )
1540     
1541     !  *******************************************************************  !
1542     !                                                                       !
1543     !  rrtm longwave radiative transfer model                               !
1544     !  atmospheric and environmental research, inc., cambridge, ma          !
1545     !                                                                       !
1546     !  original version:       e. j. mlawer, et al.                         !
1547     !  revision for ncar ccm:  michael j. iacono; september, 1998           !
1548     !  revision to use variable diffusivity angle instead of the original   !
1549     !     fixed value of 1.66: m. j. iacono apr 2004                        !
1550     !                                                                       !
1551     !  this program calculates the upward fluxes, downward fluxes, and      !
1552     !  heating rates for an arbitrary clear or cloudy atmosphere.  the input!
1553     !  to this program is the atmospheric profile, all planck function      !
1554     !  information, and the cloud fraction by layer.  a variable diffusivity!
1555     !  angle (secdiff) is used forthe angle integration.  bands 2-3 and 5-9 !
1556     !  use a value of secdiff that varies from 1.50 to 1.80 as a function of!
1557     !  the column water vapor, and other bands use a value of 1.66.  the    !
1558     !  gaussian weight appropriate to this angle (wtnum=0.5) is applied     !
1559     !  here.  note that use of a single angle for the flux integration      !
1560     !  can cause errors of 1 to 4 w/m2 within cloudy layers.                !
1561     !                                                                       !
1562     !  *******************************************************************  !
1563     !
1564           use module_radlw_avplank
1565     !
1566           implicit none
1567     
1568     !  ---  inputs:
1569           integer, intent(in)  ::  NLAY, NLP1
1570     
1571           integer, intent(in)  ::  itr(:,:)
1572     
1573           real (kind=kind_phys), dimension(0:), intent(in) :: tz, cldfrac
1574     
1575           real (kind=kind_phys), dimension(:),  intent(in) :: tavel, delp,  &
1576          &       semiss, secdiff
1577     
1578           real (kind=kind_phys), dimension(:,:),intent(in) :: taucloud,     &
1579          &       pfrac
1580     
1581           real (kind=kind_phys), intent(in) :: stemp
1582     
1583     !  ---  outputs:
1584           real (kind=kind_phys), dimension(:),  intent(out) :: htr, htrcl
1585           real (kind=kind_phys), dimension(:,:),intent(out) :: htrb
1586     
1587           real (kind=kind_phys), dimension(0:), intent(out) ::              &
1588          &       totuflux, totdflux, totuclfl, totdclfl
1589     
1590     !  ---  locals:
1591           real (kind=kind_phys), dimension(NGPTLW,NLAY)   :: gassrcu,       &
1592          &       cldsrcu, trans0
1593           real (kind=kind_phys), dimension(NGPTLW,0:NLAY) :: bglev
1594           real (kind=kind_phys), dimension(NGPTLW)        :: radclru,       &
1595          &       radclrd, radtotu, radtotd, bgsfc
1596           real (kind=kind_phys), dimension(NBANDS,0:NLAY) :: plvl,          &
1597          &       totufxsb, totdfxsb
1598           real (kind=kind_phys), dimension(NBANDS,NLAY)   :: play, odcld,   &
1599          &       trncld, efcfr1
1600           real (kind=kind_phys), dimension(NBANDS)        :: plksfc
1601           real (kind=kind_phys), dimension(0:NLAY)        :: fnet, fnetc
1602     
1603           real (kind=kind_phys) :: totdrad, clrdrad, toturad, clrurad
1604           real (kind=kind_phys) :: delbgup, delbgdn, bglay, tau0, tauc,     &
1605          &       transc, cldsrcd, gassrcd, factot, odsm, tem1, tem2
1606     
1607           integer :: j, k, ind, inb, itm1, itm2, jtm1, jtm2
1608     
1609     !  ====================  defination of variables  ====================  !
1610     !                                                                       !
1611     !  input variables:                                                     !
1612     !    tavel   (NLAY)       ! layer temperatures (k)                      !
1613     !    tz      (0:NLAY)     ! level (interface) temperature (k)           !
1614     !    delp    (NLAY)       ! layer pressure thickness (mb)               !
1615     !    semiss  (NBANDS)     ! surface emissivities for each band          !
1616     !    cldfrac (0:NLP1)     ! layer cloud fraction (padded at 2 ends)     !
1617     !    taucloud(NBANDS,NLAY)! layer cloud optical depth                   !
1618     !    pfrac   (NGPTLW,NLAY)! planck fractions                            !
1619     !    secdiff(NBANDS)      ! variable diffusivity angle defined as an    !
1620     !                           exponential function of the column water    !
1621     !                           amount in bands 2-3 and 5-9. this reduces   !
1622     !                           the bias of several w/m2 in downward surface!
1623     !                           flux in high water profiles caused by using !
1624     !                           the constant diffusivity angle of 1.66.(mji)!
1625     !    stemp                ! surface ground temperature (k)              !
1626     !    itr     (NGPTLW,NLAY)! integer look-up table index                 !
1627     !    NLAY/NLP1            ! number of model layers/levels               !
1628     !                                                                       !
1629     !  constants or shared variables:                                       !
1630     !    NGPTLW               ! total number of g-point subintervals        !
1631     !    NBANDS               ! number of longwave spectral bands           !
1632     !    wtnum                ! weight for radiance to flux conversion      !
1633     !    bpade                ! pade constant                               !
1634     !    tau                  ! clear sky optical depth look-up table       !
1635     !    tf                   ! tau transition function look-up table       !
1636     !    trans                ! clear sky transmittance look-up table       !
1637     !                                                                       !
1638     !  output variables:                                                    !
1639     !    totuflux(0:NLAY)     ! upward longwave flux (w/m2)                 !
1640     !    totdflux(0:NLAY)     ! downward longwave flux (w/m2)               !
1641     !    htr     (NLAY)       ! longwave heating rate (k/day)               !
1642     !    totuclfl(0:NLAY)     ! clear sky upward longwave flux (w/m2)       !
1643     !    totdclfl(0:NLAY)     ! clear sky downward longwave flux (w/m2)     !
1644     !    htrcl   (NLAY)       ! clear sky longwave heating rate (k/day)     !
1645     !    htrb    (NLAY,NBANDS)! spectral band lw heating rate (k/day)       !
1646     !                                                                       !
1647     !  local variables:                                                     !
1648     !    odcld   (NBANDS,NLAY)! cloud optical depth                         !
1649     !    trncld  (NBANDS,NLAY)! cloud transmittance                         !
1650     !    efcfr1  (NBANDS,NLAY)! effective clear  sky fraction               !
1651     !    radtotu (NGPTLW)     ! upward radiance                             !
1652     !    radtotd (NGPTLW)     ! downward radiance                           !
1653     !    radclru (NGPTLW)     ! clear sky upward radiance                   !
1654     !    radclrd (NGPTLW)     ! clear sky downward radiance                 !
1655     !    toturad              ! spectrally summed upward radiance           !
1656     !    totdrad              ! spectrally summed downward radiance         !
1657     !    clrurad              ! spectrally summed clear sky upward radiance !
1658     !    clrdrad              ! spectrally summed clear sky dnward radiance !
1659     !    totufxsb(NBANDS,NLAY)! spectral band upward longwave flux (w/m2)   !
1660     !    totdfxsb(NBANDS,NLAY)! spectral band downward longwave flux (w/m2) !
1661     !                                                                       !
1662     !    fnet    (0:NLAY)     ! net longwave flux (w/m2)                    !
1663     !    fnetc   (0:NLAY)     ! clear sky net longwave flux (w/m2)          !
1664     !                                                                       !
1665     !  =====================    end of definitions    ====================  !
1666     
1667     !
1668     !===> ... begin here
1669     !
1670     !  --- ... calculate the integrated planck functions at the surface
1671     
1672           itm1 = min(NPLNK, max(1, int(stemp-159.0) ))
1673           itm2 = min(NPLNK, itm1+1)
1674           tem1 = stemp - int(stemp)
1675           do j = 1, NBANDS
1676             plksfc(j) = delwave(j) * ( totplnk(itm1,j)                      &
1677          &           + tem1 * (totplnk(itm2,j) - totplnk(itm1,j)) )
1678           enddo
1679     
1680           do j = 1, NGPTLW
1681             inb = ngb(j)                 ! band index
1682             bgsfc(j) = pfrac(j,1) * plksfc(inb)
1683           enddo
1684     
1685     !  --- ... calculate the integrated planck functions at the level and
1686     !          layer temperatures.
1687     
1688           itm1 = min(NPLNK, max(1, int(tz(0)-159.0) ))
1689           itm2 = min(NPLNK, itm1+1)
1690           tem1 = tz(0) - int(tz(0))
1691           do j = 1, NBANDS
1692             plvl(j,0) = delwave(j) * ( totplnk(itm1,j)                      &
1693          &            + tem1 * (totplnk(itm2,j) - totplnk(itm1,j)) )
1694           enddo
1695     
1696           do k = 1, NLAY
1697             itm1 = min(NPLNK, max(1, int(tz(k)   -159.0) ))
1698             itm2 = min(NPLNK, itm1+1)
1699             jtm1 = min(NPLNK, max(1, int(tavel(k)-159.0) ))
1700             jtm2 = min(NPLNK, jtm1+1)
1701     
1702             tem1 = tz(k)    - int(tz(k))
1703             tem2 = tavel(k) - int(tavel(k))
1704     
1705             do j = 1, NBANDS
1706               plvl(j,k) = delwave(j) * ( totplnk(itm1,j)                    &
1707          &              + tem1 * (totplnk(itm2,j) - totplnk(itm1,j)) )
1708               play(j,k) = delwave(j) * ( totplnk(jtm1,j)                    &
1709          &              + tem2 * (totplnk(jtm2,j) - totplnk(jtm1,j)) )
1710     
1711     !  --- ... cloudy sky optical depth and absorptivity.
1712     ! mji     odcld(j,k)  = secang * taucloud(j,k)
1713               odcld(j,k)  = secdiff(j) * taucloud(j,k)
1714               trncld(j,k) = exp( -odcld(j,k) )
1715               efcfr1(j,k) = 1.0 - cldfrac(k) + trncld(j,k)*cldfrac(k)
1716             enddo
1717     
1718             do j = 1, NGPTLW
1719               inb = ngb(j)                 ! band index
1720               bglev(j,k-1) = pfrac(j,k) * plvl(inb,k-1)
1721             enddo
1722           enddo
1723     
1724     !  --- ...  initialize for radiative transfer.
1725     
1726           if ( lhlwb ) then
1727             do k = 0, NLAY
1728               do j = 1, NBANDS
1729                 totufxsb(j,k) = f_zero
1730                 totdfxsb(j,k) = f_zero
1731               enddo
1732             enddo
1733           endif
1734     
1735           do j = 1, NGPTLW
1736              inb = ngb(j)                 ! band index
1737              radclrd(j) = f_zero
1738              radtotd(j) = f_zero
1739              bglev(j,NLAY) = pfrac(j,NLAY) * plvl(inb,NLAY)
1740           enddo
1741     
1742     !===> ...  downward radiative transfer
1743     !          totdrad holds summed radiance for total sky stream
1744     !          clrdrad holds summed radiance for clear sky stream
1745     
1746           do k = NLAY, 1, -1
1747     
1748             totdrad = f_zero
1749             clrdrad = f_zero
1750     
1751             if (cldfrac(k) > eps) then
1752     !  --- ... cloudy layer
1753     
1754               do j = 1, NGPTLW
1755     !  --- ... get lookup table index
1756                 ind = itr(j,k)
1757                 inb = ngb(j)                 ! band index
1758     
1759     !  --- ... get clear sky transmittance from lookup table
1760                 tau0 = tf(ind)
1761                 trans0(j,k) = trans(ind)
1762                 transc      = trans0(j,k) * trncld(inb,k)
1763     
1764     !  --- ... add clear sky and cloud optical depths
1765                 odsm = tau(ind) + odcld(inb,k)
1766                 tauc = odsm / (bpade + odsm)
1767     
1768                 bglay = pfrac(j,k) * play(inb,k)
1769                 delbgup = bglev(j,k) - bglay
1770                 tem1 = bglay + tau0*delbgup
1771                 tem2 = bglay + tauc*delbgup
1772                 gassrcu(j,k) = tem1 - trans0(j,k)*tem1
1773                 cldsrcu(j,k) = tem2 - transc     *tem2
1774     
1775                 delbgdn = bglev(j,k-1) - bglay
1776                 tem1 = bglay + tau0*delbgdn
1777                 tem2 = bglay + tauc*delbgdn
1778                 gassrcd = tem1 - trans0(j,k)*tem1
1779                 cldsrcd = tem2 - transc     *tem2
1780     
1781     !  --- ... total sky radiance
1782                 radtotd(j) = radtotd(j)*trans0(j,k)*efcfr1(inb,k)           &
1783          &                 + gassrcd + cldfrac(k)*(cldsrcd - gassrcd)
1784                 totdrad = totdrad + radtotd(j)
1785     
1786     !  --- ... clear sky radiance
1787                 radclrd(j) = radclrd(j)*trans0(j,k) + gassrcd
1788                 clrdrad = clrdrad + radclrd(j)
1789               enddo
1790     
1791             else
1792     
1793     !  --- ... clear layer
1794     
1795               do j = 1, NGPTLW
1796                 ind = itr(j,k)
1797                 inb = ngb(j)                 ! band index
1798     
1799     !  --- ... get clear sky transmittance from lookup table
1800                 tau0 = tf(ind)
1801                 trans0(j,k) = trans(ind)
1802     
1803                 bglay = pfrac(j,k) * play(inb,k)
1804     
1805                 delbgup = bglev(j,k) - bglay
1806                 tem1 = bglay + tau0*delbgup
1807                 gassrcu(j,k) = tem1 - trans0(j,k)*tem1
1808     !           cldsrcu(j,k) = 0.0
1809     
1810                 delbgdn  = bglev(j,k-1) - bglay
1811                 tem2 = bglay + tau0*delbgdn
1812                 gassrcd = tem2 - trans0(j,k)*tem2
1813     
1814     !  --- ... total sky radiance
1815                 radtotd(j) = radtotd(j)*trans0(j,k) + gassrcd
1816                 totdrad = totdrad + radtotd(j)
1817     
1818     !  --- ... clear sky radiance
1819                 radclrd(j) = radclrd(j)*trans0(j,k) + gassrcd
1820                 clrdrad = clrdrad + radclrd(j)
1821               enddo
1822     
1823             endif
1824     
1825             totdflux(k-1) = totdrad
1826             totdclfl(k-1) = clrdrad
1827     
1828     !  --- ... total sky radiance for each of the spectral bands
1829             if ( lhlwb ) then
1830               do j = 1, NGPTLW
1831                 inb = ngb(j)                 ! band index
1832                 totdfxsb(inb,k-1) = totdfxsb(inb,k-1) + radtotd(j)
1833               enddo
1834     
1835               totdfxsb(:,NLAY) = f_zero
1836             endif
1837     
1838           enddo   ! end do_k_loop
1839     
1840           totdflux(NLAY) = f_zero
1841           totdclfl(NLAY) = f_zero
1842     
1843     !  --- ...  spectral emissivity & reflectance
1844     !           include the contribution of spectrally varying longwave
1845     !           emissivity and reflection from the surface to the upward
1846     !           radiative transfer.
1847     !    note: spectral and lambertian reflection are identical for the one
1848     !           angle flux integration used here.
1849     
1850           toturad = f_zero
1851           clrurad = f_zero
1852     
1853           do j = 1, NGPTLW
1854             inb = ngb(j)                 ! band index
1855             tem1 = 1.0 - semiss(inb)
1856     !       tem2 = bglev(j,0) * semiss(inb)
1857             tem2 = bgsfc(j) * semiss(inb)
1858     
1859     !  --- ... total sky radiance
1860             radtotu(j) = tem2 + tem1 * radtotd(j)
1861             toturad = toturad + radtotu(j)
1862     
1863     !  --- ... clear sky radiance
1864             radclru(j) = tem2 + tem1 * radclrd(j)
1865             clrurad = clrurad + radclru(j)
1866           enddo
1867     
1868           totuflux(0) = toturad
1869           totuclfl(0) = clrurad
1870     
1871     !  --- ... total sky radiance for each of the spectral bands
1872           if ( lhlwb ) then
1873             do j = 1, NGPTLW
1874               inb = ngb(j)                 ! band index
1875               totufxsb(inb,0) = totufxsb(inb,0) + radtotu(j)
1876             enddo
1877           endif
1878     
1879     !     print *,' toturad(0)=',totuflux(0)
1880     !     print *,' clrurad(0)=',totuclfl(0)
1881     
1882     !===> ...  upward radiative transfer
1883     !          toturad holds summed radiance for total sky stream
1884     !          clrurad holds summed radiance for clear sky stream
1885     
1886           do k = 1, NLAY
1887     
1888             toturad = f_zero
1889             clrurad = f_zero
1890     
1891     !  --- ... check flag for cloud in current layer
1892             if (cldfrac(k) > eps) then
1893     
1894     !  --- ... cloudy layers
1895     
1896               do j = 1, NGPTLW
1897                 inb   = ngb(j)                 ! band index
1898     
1899     !  --- ... total sky radiance
1900                 radtotu(j) = radtotu(j)*trans0(j,k)*efcfr1(inb,k)           &
1901          &         + gassrcu(j,k) + cldfrac(k)*(cldsrcu(j,k)-gassrcu(j,k))
1902                 toturad = toturad + radtotu(j)
1903     
1904     !  --- ... clear sky radiance
1905                 radclru(j) = radclru(j)*trans0(j,k) + gassrcu(j,k)
1906                 clrurad = clrurad + radclru(j)
1907               enddo
1908     
1909             else
1910     
1911     !  --- ... clear layer
1912     
1913               do j = 1, NGPTLW
1914     
1915     !  --- ... total sky radiance
1916                 radtotu(j) = radtotu(j)*trans0(j,k) + gassrcu(j,k)
1917                 toturad = toturad + radtotu(j)
1918     
1919     !  --- ... clear sky radiance
1920                 radclru(j) = radclru(j)*trans0(j,k) + gassrcu(j,k)
1921                 clrurad = clrurad + radclru(j)
1922               enddo
1923     
1924             endif
1925     
1926             totuflux(k) = toturad
1927             totuclfl(k) = clrurad
1928     
1929     !  --- ... total sky radiance for each of the spectral bands
1930             if ( lhlwb ) then
1931               do j = 1, NGPTLW
1932                 inb = ngb(j)                 ! band index
1933                 totufxsb(inb,k) = totufxsb(inb,k) + radtotu(j)
1934               enddo
1935             endif
1936     
1937           enddo
1938     
1939     !===> ...  convert radiances to fluxes and heating rates for total sky.
1940     !          calculates clear sky surface and toa values.  to compute clear
1941     !          sky profiles, uncomment relevant lines below.
1942     
1943           factot = wtnum * fluxfac
1944     
1945           totuflux(0) = totuflux(0) * factot
1946           totdflux(0) = totdflux(0) * factot
1947           totuclfl(0) = totuclfl(0) * factot
1948           totdclfl(0) = totdclfl(0) * factot
1949           fnet(0) = totuflux(0) - totdflux(0)
1950     
1951           do k = 1, NLAY
1952             totuflux(k) = totuflux(k) * factot
1953             totdflux(k) = totdflux(k) * factot
1954     
1955             totuclfl(k) = totuclfl(k) * factot
1956             totdclfl(k) = totdclfl(k) * factot
1957     
1958             fnet(k) = totuflux(k) - totdflux(k)
1959             htr (k) = heatfac * (fnet(k-1) - fnet(k)) / delp(k)
1960           enddo
1961     
1962     !! --- ...  optional clear sky heating rates
1963           if ( lhlw0 ) then
1964             fnetc(0) = totuclfl(0) - totdclfl(0)
1965     
1966             do k = 1, NLAY
1967               fnetc(k) = totuclfl(k) - totdclfl(k)
1968               htrcl(k) = heatfac * (fnetc(k-1) - fnetc(k)) / delp(k)
1969             enddo
1970           endif
1971     
1972     !! --- ...  optional spectral band heating rates
1973           if ( lhlwb ) then
1974             do j = 1, NBANDS
1975               totufxsb(j,0) = totufxsb(j,0) * factot
1976               totdfxsb(j,0) = totdfxsb(j,0) * factot
1977               fnet(0) = totufxsb(j,0) - totdfxsb(j,0)
1978     
1979               do k = 1, NLAY
1980                 totufxsb(j,k) = totufxsb(j,k) * factot
1981                 totdfxsb(j,k) = totdfxsb(j,k) * factot
1982                 fnet(k) = totufxsb(j,k) - totdfxsb(j,k)
1983                 htrb(k,j) = heatfac * (fnet(k-1) - fnet(k)) / delp(k)
1984               enddo
1985             enddo
1986           endif
1987     
1988           return
1989     !...................................
1990           end subroutine rtrn
1991     !-----------------------------------
1992     
1993     
1994     
1995     !-----------------------------------
1996           subroutine rtrnmr                                                 &
1997     !...................................
1998     
1999     !  ---  inputs:
2000          &     ( tavel,tz,delp,semiss,cldfrac,taucloud,pfrac,               &
2001          &       secdiff, stemp, itr, NLAY, NLP1,                           &
2002     !  ---  outputs:
2003          &       totuflux,totdflux,htr, totuclfl,totdclfl,htrcl, htrb       &
2004          &     )
2005     
2006     !  *******************************************************************  !
2007     !                                                                       !
2008     !  rrtm longwave radiative transfer model                               !
2009     !  atmospheric and environmental research, inc., cambridge, ma          !
2010     !                                                                       !
2011     !  original version:       e. j. mlawer, et al.                         !
2012     !  revision for ncar ccm:  michael j. iacono; september, 1998           !
2013     !  revision to use variable diffusivity angle instead of the original   !
2014     !     fixed value of 1.66: m. j. iacono apr 2004                        !
2015     !                                                                       !
2016     !  this program calculates the upward fluxes, downward fluxes, and      !
2017     !  heating rates for an arbitrary clear or cloudy atmosphere.  the input!
2018     !  to this program is the atmospheric profile, all planck function      !
2019     !  information, and the cloud fraction by layer.  a variable diffusivity!
2020     !  angle (secdiff) is used forthe angle integration.  bands 2-3 and 5-9 !
2021     !  use a value of secdiff that varies from 1.50 to 1.80 as a function of!
2022     !  the column water vapor, and other bands use a value of 1.66.  the    !
2023     !  gaussian weight appropriate to this angle (wtnum=0.5) is applied     !
2024     !  here.  note that use of a single angle for the flux integration      !
2025     !  can cause errors of 1 to 4 w/m2 within cloudy layers.                !
2026     !  this routine computes a generalized maximum/random cloud overlap.    !
2027     !  adjacent cloud layers are treated with maximum overlap in which up   !
2028     !  to two previous layers of cloud information is considered.  non-     !
2029     !  adjacent groups of clouds are treated with random overlap.           !
2030     !                                                                       !
2031     !  *******************************************************************  !
2032     !
2033           use module_radlw_avplank
2034     !
2035           implicit none
2036     
2037     !  ---  inputs:
2038           integer, intent(in)  ::  NLAY, NLP1
2039     
2040           integer, intent(in)  ::  itr(:,:)
2041     
2042           real (kind=kind_phys), dimension(0:), intent(in) :: tz, cldfrac
2043     
2044           real (kind=kind_phys), dimension(:),  intent(in) :: tavel, delp,  &
2045          &       semiss, secdiff
2046     
2047           real (kind=kind_phys), dimension(:,:),intent(in) :: taucloud,     &
2048          &       pfrac
2049     
2050           real (kind=kind_phys), intent(in) :: stemp
2051     
2052     !  ---  outputs:
2053           real (kind=kind_phys), dimension(:),  intent(out) :: htr, htrcl
2054           real (kind=kind_phys), dimension(:,:),intent(out) :: htrb
2055     
2056           real (kind=kind_phys), dimension(0:), intent(out) ::              &
2057          &       totuflux, totdflux, totuclfl, totdclfl
2058     
2059     !  ---  locals:
2060     !  dimensions for radiative transfer
2061           real (kind=kind_phys), dimension(NGPTLW,NLAY)   :: gassrcu,       &
2062          &       cldsrcu, trans0, transc
2063           real (kind=kind_phys), dimension(NGPTLW,0:NLAY) :: bglev
2064           real (kind=kind_phys), dimension(NGPTLW)        :: radclru,       &
2065          &       radclrd, radtotu, radtotd, bgsfc
2066           real (kind=kind_phys), dimension(NBANDS,0:NLAY) :: plvl,          &
2067          &       totufxsb, totdfxsb
2068           real (kind=kind_phys), dimension(NBANDS,NLAY)   :: play,          &
2069          &       odcld, trncld
2070           real (kind=kind_phys), dimension(NBANDS)        :: plksfc
2071           real (kind=kind_phys), dimension(0:NLAY)        :: fnet, fnetc
2072     
2073           real (kind=kind_phys) :: totdrad, clrdrad, toturad, clrurad
2074           real (kind=kind_phys) :: delbgup, delbgdn, bglay, tau0, tauc,     &
2075          &       cldsrcd, gassrcd, factot, odsm, tem1, tem2
2076     
2077           integer :: j, k, ind, inb, itm1, itm2, jtm1, jtm2
2078     
2079     !  dimensions for cloud overlap adjustment
2080           real (kind=kind_phys), dimension(NGPTLW) ::  clrradu, cldradu,    &
2081          &       clrradd, cldradd, rad
2082           real (kind=kind_phys), dimension(1:NLP1) ::  faccld1u, faccld2u,  &
2083          &       facclr1u, facclr2u, faccmb1u, faccmb2u
2084           real (kind=kind_phys), dimension(0:NLAY) ::  faccld1d, faccld2d,  &
2085          &       facclr1d, facclr2d, faccmb1d, faccmb2d
2086           real (kind=kind_phys) :: fmax, fmin, rat1, rat2, radmod, cldsrc
2087     
2088           logical :: istcldu(NLAY), istcldd(NLAY)
2089     
2090     !  ====================  defination of variables  ====================  !
2091     !                                                                       !
2092     !  input variables:                                                     !
2093     !    tavel   (NLAY)       ! layer temperatures (k)                      !
2094     !    tZ      (0:NLAY)     ! level (interface) temperatures (k)          !
2095     !    delp    (NLAY)       ! layer pressure thickness (mb)               !
2096     !    semiss  (NBANDS)     ! surface emissivities for each band          !
2097     !    cldfrac (0:NLP1)     ! layer cloud fraction (padded at 2 ends)     !
2098     !    taucloud(NBANDS,NLAY)! layer cloud optical depth                   !
2099     !    pfrac   (NGPTLW,NLAY)! planck fractions                            !
2100     !    secdiff(NBANDS)      ! variable diffusivity angle defined as an    !
2101     !                           exponential function of the column water    !
2102     !                           amount in bands 2-3 and 5-9. this reduces   !
2103     !                           the bias of several w/m2 in downward surface!
2104     !                           flux in high water profiles caused by using !
2105     !                           the constant diffusivity angle of 1.66.(mji)!
2106     !    stemp   (NBANDS)     ! surface ground temperature (k)              !
2107     !    itr     (NGPTLW,NLAY)! integer look-up table index                 !
2108     !    NLAY/NLP1            ! number of model layers/levels               !
2109     !                                                                       !
2110     !  constants or shared variables:                                       !
2111     !    NGPTLW               ! total number of g-point subintervals        !
2112     !    NBANDS               ! number of longwave spectral bands           !
2113     !    wtnum                ! weight for radiance to flux conversion      !
2114     !    bpade                ! pade constant                               !
2115     !    tau                  ! clear sky optical depth look-up table       !
2116     !    tf                   ! tau transition function look-up table       !
2117     !    trans                ! clear sky transmittance look-up table       !
2118     !                                                                       !
2119     !  output variables:                                                    !
2120     !    totuflux(0:NLAY)     ! upward longwave flux (w/m2)                 !
2121     !    totdflux(0:NLAY)     ! downward longwave flux (w/m2)               !
2122     !    htr     (NLAY)       ! longwave heating rate (k/d or k/s)          !
2123     !    totuclfl(0:NLAY)     ! clear sky upward longwave flux (w/m2)       !
2124     !    totdclfl(0:NLAY)     ! clear sky downward longwave flux (w/m2)     !
2125     !    htrcl   (NLAY)       ! clear sky longwave heating rate (k/d or k/s)!
2126     !    htrb    (NLAY,NBANDS)! spectral band lw heating rate (k/d or k/s)  !
2127     !                                                                       !
2128     !  local variables:                                                     !
2129     !    odcld   (NBANDS,NLAY)! cloud optical depth                         !
2130     !    trncld  (NBANDS,NLAY)! cloud transmittance                         !
2131     !    radtotu (NGPTLW)     ! upward radiance                             !
2132     !    radtotd (NGPTLW)     ! downward radiance                           !
2133     !    radclru (NGPTLW)     ! clear sky upward radiance                   !
2134     !    radclrd (NGPTLW)     ! clear sky downward radiance                 !
2135     !    toturad              ! spectrally summed upward radiance           !
2136     !    totdrad              ! spectrally summed downward radiance         !
2137     !    clrurad              ! spectrally summed clear sky upward radiance !
2138     !    clrdrad              ! spectrally summed clear sky dnward radiance !
2139     !    totufxsb(NBANDS,NLAY)! spectral band upward longwave flux (w/m2)   !
2140     !    totdfxsb(NBANDS,NLAY)! spectral band downward longwave flux (w/m2) !
2141     !                                                                       !
2142     !    fnet    (0:NLAY)     ! net longwave flux (w/m2)                    !
2143     !    fnetc   (0:NLAY)     ! clear sky net longwave flux (w/m2)          !
2144     !                                                                       !
2145     !  =====================    end of definitions    ====================  !
2146     !
2147     
2148     !
2149     !===> ... begin here
2150     !
2151           do k = 1, NLP1
2152             faccld1u(k) = f_zero
2153             faccld2u(k) = f_zero
2154             facclr1u(k) = f_zero
2155             facclr2u(k) = f_zero
2156             faccmb1u(k) = f_zero
2157             faccmb2u(k) = f_zero
2158           enddo
2159     
2160           istcldu(1) = cldfrac(1) > eps
2161           rat1 = f_zero
2162           rat2 = f_zero
2163     
2164           do k = 1, NLAY-1
2165     
2166             istcldu(k+1) = cldfrac(k+1)>eps .and. cldfrac(k)<=eps
2167     
2168             if (cldfrac(k) > eps) then
2169     !  --- ... maximum/random cloud overlap
2170     
2171               if (cldfrac(k+1) >= cldfrac(k)) then
2172                 if (istcldu(k)) then
2173                   if (cldfrac(k) < 1.0) then
2174                     facclr2u(k+1) = (cldfrac(k+1) - cldfrac(k))             &
2175          &                        / (1.0 - cldfrac(k))
2176                   endif
2177                   facclr2u(k) = f_zero
2178                   faccld2u(k) = f_zero
2179                 else
2180                   fmax = max(cldfrac(k), cldfrac(k-1))
2181                   if (cldfrac(k+1) > fmax) then
2182                     facclr1u(k+1) = rat2
2183                     facclr2u(k+1) = (cldfrac(k+1) - fmax)/(1.0 - fmax)
2184                   elseif (cldfrac(k+1) < fmax) then
2185                     facclr1u(k+1) = (cldfrac(k+1) - cldfrac(k))             &
2186          &                        / (cldfrac(k-1) - cldfrac(k))
2187                   else
2188                     facclr1u(k+1) = rat2
2189                   endif
2190                 endif
2191     
2192                 if (facclr1u(k+1)>f_zero .or. facclr2u(k+1)>f_zero) then
2193                   rat1 = 1.0
2194                   rat2 = f_zero
2195                 else
2196                   rat1 = f_zero
2197                   rat2 = f_zero
2198                 endif
2199               else
2200                 if (istcldu(k)) then
2201                   faccld2u(k+1) = (cldfrac(k) - cldfrac(k+1)) /  cldfrac(k)
2202                   facclr2u(k) = f_zero
2203                   faccld2u(k) = f_zero
2204                 else
2205                   fmin = min(cldfrac(k), cldfrac(k-1))
2206                   if (cldfrac(k+1) <= fmin) then
2207                     faccld1u(k+1) = rat1
2208                     faccld2u(k+1) = (fmin - cldfrac(k+1)) / fmin
2209                   else
2210                     faccld1u(k+1) = (cldfrac(k) - cldfrac(k+1))             &
2211          &                        / (cldfrac(k) - fmin)
2212                   endif
2213                 endif
2214     
2215                 if (faccld1u(k+1)>f_zero .or. faccld2u(k+1)>f_zero) then
2216                   rat1 = f_zero
2217                   rat2 = 1.0
2218                 else
2219                   rat1 = f_zero
2220                   rat2 = f_zero
2221                 endif
2222               endif
2223     
2224               faccmb1u(k+1) = facclr1u(k+1) * faccld2u(k) * cldfrac(k-1)
2225               faccmb2u(k+1) = faccld1u(k+1) * facclr2u(k)                   &
2226          &                  * (1.0 - cldfrac(k-1))
2227             endif
2228     
2229           enddo
2230     
2231           do k = 0, NLAY
2232             faccld1d(k) = f_zero
2233             faccld2d(k) = f_zero
2234             facclr1d(k) = f_zero
2235             facclr2d(k) = f_zero
2236             faccmb1d(k) = f_zero
2237             faccmb2d(k) = f_zero
2238           enddo
2239     
2240           istcldd(NLAY) = cldfrac(NLAY) > eps
2241           rat1 = f_zero
2242           rat2 = f_zero
2243     
2244           do k = NLAY, 2, -1
2245     
2246             istcldd(k-1) = cldfrac(k-1) > eps .and. cldfrac(k)<=eps
2247     
2248             if (cldfrac(k) > eps) then
2249     
2250               if (cldfrac(k-1) >= cldfrac(k)) then
2251                 if (istcldd(k)) then
2252                   if (cldfrac(k) < 1.0) then
2253                     facclr2d(k-1) = (cldfrac(k-1) - cldfrac(k))             &
2254          &                        / (1.0 - cldfrac(k))
2255                   endif
2256                   facclr2d(k) = f_zero
2257                   faccld2d(k) = f_zero
2258                 else
2259                   fmax = max(cldfrac(k), cldfrac(k+1))
2260     
2261                   if (cldfrac(k-1) > fmax) then
2262                     facclr1d(k-1) = rat2
2263                     facclr2d(k-1) = (cldfrac(k-1) - fmax)/(1.0 - fmax)
2264                   elseif (cldfrac(k-1) < fmax) then
2265                     facclr1d(k-1) = (cldfrac(k-1) - cldfrac(k))             &
2266          &                        / (cldfrac(k+1) - cldfrac(k))
2267                   else
2268                     facclr1d(k-1) = rat2
2269                   endif
2270                 endif
2271     
2272                 if (facclr1d(k-1)>f_zero .or. facclr2d(k-1)>f_zero) then
2273                   rat1 = 1.0
2274                   rat2 = f_zero
2275                 else
2276                   rat1 = f_zero
2277                   rat2 = f_zero
2278                 endif
2279               else
2280                 if (istcldd(k)) then
2281                   faccld2d(k-1) = (cldfrac(k) - cldfrac(k-1)) / cldfrac(k)
2282                   facclr2d(k) = f_zero
2283                   faccld2d(k) = f_zero
2284                 else
2285                   fmin = min(cldfrac(k), cldfrac(k+1))
2286     
2287                   if (cldfrac(k-1) <= fmin) then
2288                     faccld1d(k-1) = rat1
2289                     faccld2d(k-1) = (fmin - cldfrac(k-1)) / fmin
2290                   else
2291                     faccld1d(k-1) = (cldfrac(k) - cldfrac(k-1))             &
2292          &                        / (cldfrac(k) - fmin)
2293                   endif
2294                 endif
2295     
2296                 if (faccld1d(k-1)>f_zero .or. faccld2d(k-1)>f_zero) then
2297                   rat1 = f_zero
2298                   rat2 = 1.0
2299                 else
2300                   rat1 = f_zero
2301                   rat2 = f_zero
2302                 endif
2303               endif
2304     
2305               faccmb1d(k-1) = facclr1d(k-1) * faccld2d(k) * cldfrac(k+1)
2306               faccmb2d(k-1) = faccld1d(k-1) * facclr2d(k)                   &
2307          &                  * (1.0 - cldfrac(k+1))
2308     
2309             endif
2310     
2311           enddo
2312     
2313     !  --- ... calculate the integrated planck functions at the surface
2314     
2315           itm1 = min(NPLNK, max(1, int(stemp-159.0) ))
2316           itm2 = min(NPLNK, itm1+1)
2317           tem1 = stemp - int(stemp)
2318           do j = 1, NBANDS
2319             plksfc(j) = delwave(j) * ( totplnk(itm1,j)                      &
2320          &            + tem1 * (totplnk(itm2,j) - totplnk(itm1,j)) )
2321           enddo
2322     
2323           do j = 1, NGPTLW
2324             inb = ngb(j)                 ! band index
2325             bgsfc(j) = pfrac(j,1) * plksfc(inb)
2326           enddo
2327     
2328     !  --- ... calculate the integrated planck functions at the level and
2329     !          layer temperatures.
2330     
2331           itm1 = min(NPLNK, max(1, int(tz(0)-159.0) ))
2332           itm2 = min(NPLNK, itm1+1)
2333           tem1 = tz(0) - int(tz(0))
2334           do j = 1, NBANDS
2335             plvl(j,0) = delwave(j) * ( totplnk(itm1,j)                      &
2336          &            + tem1 * (totplnk(itm2,j) - totplnk(itm1,j)) )
2337           enddo
2338     
2339           do k = 1, NLAY
2340             itm1 = min(NPLNK, max(1, int(tz(k)   -159.0) ))
2341             itm2 = min(NPLNK, itm1+1)
2342             jtm1 = min(NPLNK, max(1, int(tavel(k)-159.0) ))
2343             jtm2 = min(NPLNK, jtm1+1)
2344     
2345             tem1 = tz(k)    - int(tz(k))
2346             tem2 = tavel(k) - int(tavel(k))
2347     
2348             do j = 1, NBANDS
2349               plvl(j,k) = delwave(j) * ( totplnk(itm1,j)                    &
2350          &              + tem1 * (totplnk(itm2,j) - totplnk(itm1,j)) )
2351               play(j,k) = delwave(j) * ( totplnk(jtm1,j)                    &
2352          &              + tem2 * (totplnk(jtm2,j) - totplnk(jtm1,j)) )
2353     
2354     !  --- ... cloudy sky optical depth and absorptivity.
2355     ! mji     odcld(j,k)  = secang * taucloud(j,k)
2356               odcld(j,k)  = secdiff(j) * taucloud(j,k)
2357               trncld(j,k) = exp( -odcld(j,k) )
2358             enddo
2359     
2360             do j = 1, NGPTLW
2361               inb = ngb(j)                 ! band index
2362               bglev(j,k-1) = pfrac(j,k) * plvl(inb,k-1)
2363             enddo
2364           enddo
2365     
2366     !  --- ...  initialize for radiative transfer.
2367     
2368           if ( lhlwb ) then
2369             do k = 0, NLAY
2370               do j = 1, NBANDS
2371                 totufxsb(j,k) = f_zero
2372                 totdfxsb(j,k) = f_zero
2373               enddo
2374             enddo
2375           endif
2376     
2377           do j = 1, NGPTLW
2378             inb = ngb(j)                 ! band index
2379             radclrd(j) = f_zero
2380             radtotd(j) = f_zero
2381             bglev (j,NLAY) = pfrac(j,NLAY) * plvl(inb,NLAY)
2382           enddo
2383     
2384     !===> ...  downward radiative transfer
2385     !          totdrad holds summed radiance for total sky stream
2386     !          clrdrad holds summed radiance for clear sky stream
2387     
2388           do k = NLAY, 1, -1
2389     
2390             totdrad = f_zero
2391             clrdrad = f_zero
2392     
2393             if (istcldd(k)) then
2394               do j = 1, NGPTLW
2395                 cldradd(j) = cldfrac(k) * radtotd(j)
2396                 clrradd(j) = radtotd(j) - cldradd(j)
2397                 rad    (j) = f_zero
2398               enddo
2399             endif
2400     
2401             if (cldfrac(k) > eps) then
2402     !  --- ... cloudy layer
2403     
2404               do j = 1, NGPTLW
2405     !  --- ... get lookup table index
2406                 ind = itr(j,k)
2407                 inb = ngb(j)                 ! band index
2408     
2409     !  --- ...  get tf from lookup table
2410                 tau0 = tf(ind)
2411                 trans0(j,k) = trans(ind)
2412                 transc(j,k) = trans(ind) * trncld(inb,k)
2413     
2414     !  --- ... add clear sky and cloud optical depths
2415                 odsm = tau(ind) + odcld(inb,k)
2416                 tauc = odsm / (bpade + odsm)
2417     
2418                 bglay = pfrac(j,k) * play(inb,k)
2419                 delbgup = bglev(j,k) - bglay
2420                 tem1 = bglay + tau0*delbgup
2421                 tem2 = bglay + tauc*delbgup
2422                 gassrcu(j,k) = tem1 - trans0(j,k)*tem1
2423                 cldsrcu(j,k) = tem2 - transc(j,k)*tem2
2424     
2425                 delbgdn = bglev(j,k-1) - bglay
2426                 tem1 = bglay + tau0*delbgdn
2427                 tem2 = bglay + tauc*delbgdn
2428                 gassrcd = tem1 - trans0(j,k)*tem1
2429                 cldsrcd = tem2 - transc(j,k)*tem2
2430     
2431                 cldradd(j) = cldradd(j)*transc(j,k) + cldfrac(k)*cldsrcd
2432                 clrradd(j) = clrradd(j)*trans0(j,k)                         &
2433          &                                    + (1.0-cldfrac(k))*gassrcd
2434     
2435     !  --- ... total sky radiance
2436                 radtotd(j) = cldradd(j) + clrradd(j)
2437                 totdrad = totdrad + radtotd(j)
2438     
2439     !  --- ... clear sky radiance
2440                 radclrd(j) = radclrd(j)*trans0(j,k) + gassrcd
2441                 clrdrad = clrdrad + radclrd(j)
2442     
2443                 radmod = rad(j)                                             &
2444          &        * (facclr1d(k-1)*trans0(j,k) + faccld1d(k-1)*transc(j,k)) &
2445          &         - faccmb1d(k-1)*gassrcd + faccmb2d(k-1)*cldsrcd
2446     
2447                 rad(j) = -radmod + facclr2d(k-1)*(clrradd(j) + radmod)      &
2448          &                       - faccld2d(k-1)*(cldradd(j) - radmod)
2449                 cldradd(j) = cldradd(j) + rad(j)
2450                 clrradd(j) = clrradd(j) - rad(j)
2451     
2452               enddo
2453     
2454             else
2455     
2456     !  --- ... clear layer
2457     
2458               do j = 1, NGPTLW
2459                 ind = itr(j,k)
2460                 inb = ngb(j)                 ! band index
2461     
2462     !  --- ... get tf from lookup table
2463                 tau0 = tf(ind)
2464                 trans0(j,k) = trans(ind)
2465                 transc(j,k) = f_zero
2466     
2467                 bglay = pfrac(j,k) * play(inb,k)
2468     
2469                 delbgup = bglev(j,k) - bglay
2470                 tem1 = bglay + tau0*delbgup
2471                 gassrcu(j,k) = tem1 - trans0(j,k)*tem1
2472     !           cldsrcu(j,k) = 0.0
2473     
2474                 delbgdn  = bglev(j,k-1) - bglay
2475                 tem2 = bglay + tau0*delbgdn
2476                 gassrcd = tem2 - trans0(j,k)*tem2
2477     
2478     !  --- ... total sky radiance
2479                 radtotd(j) = radtotd(j)*trans0(j,k) + gassrcd
2480                 totdrad = totdrad + radtotd(j)
2481     
2482     !  --- ... clear sky radiance
2483                 radclrd(j) = radclrd(j)*trans0(j,k) + gassrcd
2484                 clrdrad = clrdrad + radclrd(j)
2485               enddo
2486     
2487             endif
2488     
2489             totdflux(k-1) = totdrad
2490             totdclfl(k-1) = clrdrad
2491     
2492     !  --- ... total sky radiance for each of the spectral bands
2493             if ( lhlwb ) then
2494               do j = 1, NGPTLW
2495                 inb = ngb(j)                 ! band index
2496                 totdfxsb(inb,k-1) = totdfxsb(inb,k-1) + radtotd(j)
2497               enddo
2498     
2499               totdfxsb(:,NLAY) = f_zero
2500             endif
2501     
2502           enddo   ! end do_k_loop
2503     
2504           totdflux(NLAY) = f_zero
2505           totdclfl(NLAY) = f_zero
2506     
2507     !  --- ...  spectral emissivity & reflectance
2508     !           include the contribution of spectrally varying longwave
2509     !           emissivity and reflection from the surface to the upward
2510     !           radiative transfer.
2511     !    note: spectral and lambertian reflection are identical for the one
2512     !           angle flux integration used here.
2513     
2514           toturad = f_zero
2515           clrurad = f_zero
2516     
2517           do j = 1, NGPTLW
2518             inb = ngb(j)                 ! band index
2519             tem1 = 1.0 - semiss(inb)
2520     !       tem2 = bglev(j,0) * semiss(inb)
2521             tem2 = bgsfc(j) * semiss(inb)
2522     
2523     !  --- ... total sky radiance
2524             radtotu(j) = tem2 + tem1 * radtotd(j)
2525             toturad = toturad + radtotu(j)
2526     
2527     !  --- ... clear sky radiance
2528             radclru(j) = tem2 + tem1 * radclrd(j)
2529             clrurad = clrurad + radclru(j)
2530           enddo
2531     
2532           totuflux(0) = toturad
2533           totuclfl(0) = clrurad
2534     
2535     !  --- ... total sky radiance for each of the spectral bands
2536           if ( lhlwb ) then
2537             do j = 1, NGPTLW
2538               inb = ngb(j)                 ! band index
2539               totufxsb(inb,0) = totufxsb(inb,0) + radtotu(j)
2540             enddo
2541           endif
2542     
2543     !===> ...   upward radiative transfer
2544     !           toturad holds the summed radiance for total sky stream
2545     !           clrurad holds the summed radiance for clear sky stream
2546     
2547           do k = 1, NLAY
2548     
2549             toturad = f_zero
2550             clrurad = f_zero
2551     
2552             if (istcldu(k)) then
2553               do j = 1, NGPTLW
2554                 cldradu(j) = radtotu(j) * cldfrac(k)
2555                 clrradu(j) = radtotu(j) - cldradu(j)
2556                 rad(j)     = f_zero
2557               enddo
2558             endif
2559     
2560     !  --- ... check flag for cloud in current layer
2561             if (cldfrac(k) > eps) then
2562     
2563     !  --- ... cloudy layers
2564     
2565               do j = 1, NGPTLW
2566                 cldradu(j) = cldradu(j)*transc(j,k)+cldfrac(k)*cldsrcu(j,k)
2567                 clrradu(j) = clrradu(j)*trans0(j,k)                         &
2568          &                                + (1.0 - cldfrac(k))*gassrcu(j,k)
2569     
2570     !  --- ... total sky radiance
2571                 radtotu(j) = cldradu(j) + clrradu(j)
2572                 toturad = toturad + radtotu(j)
2573     
2574     !  --- ... clear sky radiance
2575                 radclru(j) = radclru(j)*trans0(j,k) + gassrcu(j,k)
2576                 clrurad = clrurad + radclru(j)
2577     
2578                 radmod = rad(j)                                             &
2579          &        * (facclr1u(k+1)*trans0(j,k) + faccld1u(k+1)*transc(j,k)) &
2580          &        - faccmb1u(k+1)*gassrcu(j,k) + faccmb2u(k+1)*cldsrcu(j,k)
2581     
2582                 rad(j) = -radmod + facclr2u(k+1)*(clrradu(j) + radmod)      &
2583          &                       - faccld2u(k+1)*(cldradu(j) - radmod)
2584                 cldradu(j) = cldradu(j) + rad(j)
2585                 clrradu(j) = clrradu(j) - rad(j)
2586               enddo
2587     
2588             else
2589     
2590     !  --- ... clear layer
2591     
2592               do j = 1, NGPTLW
2593     
2594     !  --- ... total sky radiance
2595                 radtotu(j) = radtotu(j)*trans0(j,k) + gassrcu(j,k)
2596                 toturad = toturad + radtotu(j)
2597     
2598     !  --- ... clear sky radiance
2599     !          upward clear and total sky streams must remain separate
2600     !          because surface reflectance is different for each.
2601                 radclru(j) = radclru(j)*trans0(j,k) + gassrcu(j,k)
2602                 clrurad = clrurad + radclru(j)
2603               enddo
2604     
2605             endif
2606     
2607             totuflux(k) = toturad
2608             totuclfl(k) = clrurad
2609     
2610     !  --- ... total sky radiance for each of the spectral bands
2611             if ( lhlwb ) then
2612               do j = 1, NGPTLW
2613                 inb = ngb(j)                 ! band index
2614                 totufxsb(inb,k) = totufxsb(inb,k) + radtotu(j)
2615               enddo
2616             endif
2617     
2618           enddo
2619     
2620     !===> ...  convert radiances to fluxes and heating rates for total sky.
2621     !          calculates clear sky surface and toa values. to compute clear
2622     !          sky profiles, uncomment relevant lines below.
2623     
2624           factot = fluxfac * wtnum
2625     
2626           totuflux(0) = totuflux(0) * factot
2627           totdflux(0) = totdflux(0) * factot
2628           totuclfl(0) = totuclfl(0) * factot
2629           totdclfl(0) = totdclfl(0) * factot
2630           fnet(0) = totuflux(0) - totdflux(0)
2631     
2632           do k = 1, NLAY
2633             totuflux(k) = totuflux(k) * factot
2634             totdflux(k) = totdflux(k) * factot
2635     
2636             totuclfl(k) = totuclfl(k) * factot
2637             totdclfl(k) = totdclfl(k) * factot
2638     
2639             fnet(k) = totuflux(k) - totdflux(k)
2640             htr (k) = heatfac * (fnet(k-1) - fnet(k)) / delp(k)
2641           enddo
2642     
2643     !! --- ...  optional clear sky heating rates
2644           if ( lhlw0 ) then
2645             fnetc(0) = totuclfl(0) - totdclfl(0)
2646     
2647             do k = 1, NLAY
2648               fnetc(k) = totuclfl(k) - totdclfl(k)
2649               htrcl(k) = heatfac * (fnetc(k-1) - fnetc(k)) / delp(k)
2650             enddo
2651           endif
2652     
2653     !! --- ...  optional spectral band heating rates
2654           if ( lhlwb ) then
2655             do j = 1, NBANDS
2656               totufxsb(j,0) = totufxsb(j,0) * factot
2657               totdfxsb(j,0) = totdfxsb(j,0) * factot
2658               fnet(0) = totufxsb(j,0) - totdfxsb(j,0)
2659     
2660               do k = 1, NLAY
2661                 totufxsb(j,k) = totufxsb(j,k) * factot
2662                 totdfxsb(j,k) = totdfxsb(j,k) * factot
2663                 fnet(k) = totufxsb(j,k) - totdfxsb(j,k)
2664                 htrb(k,j) = heatfac * (fnet(k-1) - fnet(k)) / delp(k)
2665               enddo
2666             enddo
2667           endif
2668     
2669           return
2670     !...................................
2671           end subroutine rtrnmr
2672     !-----------------------------------
2673     
2674     
2675     
2676     
2677     !-----------------------------------
2678           subroutine taumol                                                 &
2679     !...................................
2680     !  ---  inputs:
2681          &     ( laytrop,layswtch,laylow,h2ovmr,colamt,wx,co2mult,          &
2682          &       fac00,fac01,fac10,fac11,jp,jt,jt1,selffac,selffrac,        &
2683          &       indself,forfac,secdiff,tauaer, NLAY,                       &
2684     !  ---  outputs:
2685          &       itr, pfrac                                                 &
2686          &     )
2687     
2688     !  ************    original subprogram description    ***************  *
2689     !                                                                      *
2690     !                  optical depths developed for the                    *
2691     !                                                                      *
2692     !                rapid radiative transfer model (rrtm)                 *
2693     !                                                                      *
2694     !            atmospheric and environmental research, inc.              *
2695     !                        840 memorial drive                            *
2696     !                        cambridge, ma 02139                           *
2697     !                                                                      *
2698     !                           eli j. mlawer                              *
2699     !                         steven j. taubman                            *
2700     !                                                                      *
2701     !                       email:  mlawer@aer.com                         *
2702     !                                                                      *
2703     !        the authors wish to acknowledge the contributions of the      *
2704     !        following people:  patrick d. brown, michael j. iacono,       *
2705     !        ronald e. farren, luke chen, robert bergstrom.                *
2706     !                                                                      *
2707     !  revision for ncar ccm:  michael j. iacono; september, 1998          *
2708     !                                                                      *
2709     !     taumol                                                           *
2710     !                                                                      *
2711     !     this file contains the subroutines taugbn (where n goes from     *
2712     !     1 to 16).  taugbn calculates the optical depths and planck       *
2713     !     fractions per g-value and layer for band n.                      *
2714     !                                                                      *
2715     !  output:  transmittance look-up table index (unitless)               *
2716     !           fractions needed to compute planck functions at every layer*
2717     !           and g-value                                                *
2718     !                                                                      !
2719     !  description:                                                        !
2720     !     NG##        - number of g-values in band ## (##=01-16)           !
2721     !     nspa(iband) - for the lower atmosphere, the number of reference  !
2722     !                   atmospheres that are stored for band iband per     !
2723     !                   pressure level and temperature.  each of these     !
2724     !                   atmospheres has different relative amounts of the  !
2725     !                   key species for the band (i.e. different binary    !
2726     !                   species parameters).                               !
2727     !     nspb(iband) - same for upper atmosphere                          !
2728     !     oneminus    - since problems are caused in some cases by         !
2729     !                   interpolation parameters equal to or greater than  !
2730     !                   1, for these cases these parameters are set to this!
2731     !                   value, slightly < 1.                               !
2732     !     laytrop, layswtch,laylow                                         !
2733     !                 - layer at which switch is made from one combination !
2734     !                   of key species to another                          !
2735     !     h2ovmr(NLAY)- layer h2o volume mixing ratio (vmr)                !
2736     !     colamt(NLAY,MAXGAS)                                              !
2737     !                 - column amounts of water vapor,carbon dioxide,      !
2738     !                   ozone, nitrous oxide, methane, o2, co, respectively!
2739     !                   (molecules/cm**2)                                  !
2740     !     co2mult(NLAY)-for bands in which carbon dioxide is implemented   !
2741     !                   as a trace species, this is the factor used to     !
2742     !                   multiply the band's average co2 absorption         !
2743     !                   coefficient to get the added contribution to the   !
2744     !                   optical depth relative to 355 ppm.                 !
2745     !     facij (NLAY)- for layer lay, these are factors that are needed to!
2746     !                   compute the interpolation factors that multiply the!
2747     !                   appropriate reference k-values.  a value of 0 (1)  !
2748     !                   for i,j indicates that the corresponding factor    !
2749     !                   multiplies reference k-value for the lower (higher)!
2750     !                   of the two appropriate temperatures, and altitudes,!
2751     !                   respectively.                                      !
2752     !     jp (NLAY)   - the index of the lower (in altitude) of the two    !
2753     !                   appropriate reference pressure levels needed for   !
2754     !                   interpolation.                                     !
2755     !     jt, jt1(NLAY)-the indices of the lower of the two appropriate    !
2756     !                   reference temperatures needed for interpolation    !
2757     !                   (for pressure levels jp and jp+1, respectively)    !
2758     !     selffac(NLAY)-scale factor needed to water vapor self-continuum, !
2759     !                   equals (water vapor density)/(atmospheric density  !
2760     !                   at 296k and 1013 mb)                               !
2761     !     selffrac(NLAY)factor needed for temperature interpolation of     !
2762     !                   reference water vapor self-continuum data          !
2763     !     indself(NLAY)-index of the lower of the two appropriate reference!
2764     !                   temperatures needed for the self-continuum         !
2765     !                   interpolation                                      !
2766     !     secdiff(NBANDS)- variable diffusivity angle defined as an        !
2767     !                   exponential function of the column water amount in !
2768     !                   bands 2-3 and 5-9. this reduces the bias of several!
2769     !                   w/m2 in downward surface flux in high water        !
2770     !                   profiles caused by using the constant diffusivity  !
2771     !                   angle of 1.66.(mji)                                !
2772     !     tauaer(NBANDS,NLAY)- aerosols optical depth                      !
2773     !                                                                      !
2774     !  data input                                                          !
2775     !     absa(nspa(nn),5,13,NGnn), absb(nspb(nn),5,13:59,NGnn),           !
2776     !     selfref(10,NGnn)                                                 !
2777     !                      (note:  nn is the band number)                  !
2778     !                                                                      !
2779     !     absa    - k-values for low reference atmospheres (no water vapor !
2780     !               self-continuum) (units: cm**2/molecule)                !
2781     !     absb    - k-values for high reference atmospheres (all sources)  !
2782     !               (units: cm**2/molecule)                                !
2783     !     selfref - k-values for water vapor self-continuum for reference  !
2784     !               atmospheres (used below laytrop)                       !
2785     !               (units: cm**2/molecule)                                !
2786     !                                                                      !
2787     !  ******************************************************************  !
2788     !
2789           implicit none
2790     !
2791     !  ---  inputs:
2792           integer, intent(in) :: laytrop, layswtch, laylow, NLAY
2793     
2794           integer, dimension(:), intent(in) :: jp, jt, jt1, indself
2795     
2796           real (kind=kind_phys), dimension(:),  intent(in) :: h2ovmr,       &
2797          &       co2mult, fac00, fac01, fac10, fac11, selffac, selffrac,    &
2798          &       forfac, secdiff
2799     
2800           real (kind=kind_phys), dimension(:,:),intent(in) :: colamt, wx,   &
2801          &       tauaer
2802     
2803     !  ---  outputs:
2804           real (kind=kind_phys), dimension(:,:), intent(out) :: pfrac
2805     
2806           integer,               dimension(:,:), intent(out) :: itr
2807     
2808     !  ---  locals:
2809           real (kind=kind_phys) :: taug(NGPTLW,NLAY), tem1, tem2
2810           integer :: j, k, ja, jb, kk, id0(NLAY,NBANDS), id1(NLAY,NBANDS),  &
2811          &           inb
2812     !
2813     !===> ... begin here
2814     !
2815           do j = 1, NBANDS
2816             ja = 13
2817             jb = 12
2818             kk = laytrop
2819             if (j == 8) then
2820               ja = 7
2821               jb = 6
2822               kk = layswtch
2823             endif
2824     
2825             do k = 1, kk
2826               id0(k,j) = ((jp(k)-1) *5 + jt (k) - 1) * nspa(j)
2827               id1(k,j) = ( jp(k)    *5 + jt1(k) - 1) * nspa(j)
2828             enddo
2829             do k = kk+1, NLAY
2830               id0(k,j) = ((jp(k)-ja)*5 + jt (k) - 1) * nspb(j)
2831               id1(k,j) = ((jp(k)-jb)*5 + jt1(k) - 1) * nspb(j)
2832             enddo
2833           enddo
2834     
2835           call taugb01
2836           call taugb02
2837           call taugb03
2838           call taugb04
2839           call taugb05
2840           call taugb06
2841           call taugb07
2842           call taugb08
2843           call taugb09
2844           call taugb10
2845           call taugb11
2846           call taugb12
2847           call taugb13
2848           call taugb14
2849           call taugb15
2850           call taugb16
2851     
2852     ! mji do k = 1, NLAY
2853     !       do j = 1, NGPTLW
2854     !         tem1 = max( f_zero, secang*taug(j,k) )
2855     !         tem2 = tem1 / (bpade + tem1)
2856     !         itr(j,k) = 5.0e3 * tem2 + 0.5
2857     !       enddo
2858     ! mji enddo
2859     
2860           do j = 1, NGPTLW
2861             inb = ngb(j)
2862     
2863             do k = 1, NLAY
2864     !         tem1 = max( f_zero, secdiff(inb)*taug(j,k) )
2865               tem1 = max( f_zero, secdiff(inb)*(taug(j,k)+tauaer(inb,k)) )
2866               tem2 = tem1 / (bpade + tem1)
2867               itr(j,k) = max(0, min(N5000, int(5.0e3*tem2+0.5) ))
2868             enddo
2869           enddo
2870     
2871     
2872     ! =================
2873           contains
2874     ! =================
2875     
2876     !-----------------------------------
2877           subroutine taugb01
2878     !...................................
2879     
2880     !  ------------------------------------------------------------------  !
2881     !     written by eli j. mlawer, atmospheric & environmental research.  !
2882     !     revised by michael j.iacono, atmospheric & environmental research!
2883     !                                                                      !
2884     !     band 1:  10-250 cm-1 (low - h2o; high - h2o)                     !
2885     !                                                                      !
2886     !     compute the optical depth by interpolating in ln(pressure) and   !
2887     !     temperature.  below laytrop, the water vapor self-continuum      !
2888     !     is interpolated (in temperature) separately.                     !
2889     !  ------------------------------------------------------------------  !
2890     !
2891           use module_radlw_kgb01
2892     !
2893           implicit none
2894     !
2895           integer :: j, k, ind01, ind02, ind11, ind12, inds
2896     
2897           do k = 1, laytrop
2898             ind01 = id0(k,1) + 1
2899             ind02 = ind01 + 1
2900             ind11 = min(MSA01, id1(k,1) + 1 )
2901             ind12 = min(MSA01, ind11 + 1)
2902             inds = indself(k)
2903     
2904             do j = 1, NG01
2905               taug(j,k) = colamt(k,1)                                       &
2906          &         * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) +    &
2907          &             fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) +    &
2908          &             selffac(k)*( selfref(inds,j) + selffrac(k) *         &
2909          &            (selfref(inds+1,j) - selfref(inds,j)) ) +             &
2910          &             forfac(k)*forref(j) )
2911     
2912               pfrac(j,k) = fracrefa(j)
2913             enddo
2914           enddo
2915     
2916           do k = laytrop+1, NLAY
2917             ind01 = id0(k,1) + 1
2918             ind02 = ind01 + 1
2919             ind11 = id1(k,1) + 1
2920             ind12 = ind11 + 1
2921     
2922             do j = 1, NG01
2923               taug(j,k) = colamt(k,1)                                       &
2924          &         * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) +    &
2925          &             fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) +    &
2926          &             forfac(k) * forref(j) )
2927     
2928               pfrac(j,k) = fracrefb(j)
2929             enddo
2930           enddo
2931     
2932           return
2933     !...................................
2934           end subroutine taugb01
2935     !-----------------------------------
2936     
2937     
2938     !-----------------------------------
2939           subroutine taugb02
2940     !...................................
2941     
2942     !     band 2:  250-500 cm-1 (low - h2o; high - h2o)
2943     !
2944           use module_radlw_kgb02
2945     !
2946           implicit none
2947     !
2948           integer :: i, j, k, ind01, ind02, ind11, ind12, inds, ifrac, ifp
2949     !
2950           real (kind=kind_phys) :: fc00, fc01, fc10, fc11, h2oparam,        &
2951          &      fracint, one
2952           data  one / 1.0 /
2953     
2954     !     compute the optical depth by interpolating in ln(pressure) and 
2955     !     temperature.  below laytrop, the water vapor self-continuum is 
2956     !     interpolated (in temperature) separately.
2957     
2958           do k = 1, laytrop
2959             h2oparam = h2ovmr(k) / (h2ovmr(k) + 0.002)
2960     
2961             ifrac = 13
2962     !!changed b_do_ifrac : do i = 2, 12
2963             lab_do_ifrac : do i = 2, 13
2964               if (h2oparam >= refparam(i)) then
2965                 ifrac = i
2966                 exit lab_do_ifrac
2967               endif
2968             enddo lab_do_ifrac
2969     
2970             fracint = max(-one, min(one, (h2oparam - refparam(ifrac))       &
2971          &          / (refparam(ifrac-1) - refparam(ifrac)) ))
2972     
2973             ifp = max( 0, int(2.e2*(fac11(k)+fac01(k))+0.5) )
2974             fc00 = fac00(k) * corr2(ifp)
2975             fc10 = fac10(k) * corr2(ifp)
2976             fc01 = fac01(k) * corr1(ifp)
2977             fc11 = fac11(k) * corr1(ifp)
2978     
2979             ind01 = id0(k,2) + 1
2980             ind02 = ind01 + 1
2981             ind11 = min(MSA02, id1(k,2) + 1 )
2982             ind12 = min(MSA02, ind11 + 1 )
2983             inds = indself(k)
2984     
2985             do j = 1, NG02
2986               taug(NS02+j,k) = colamt(k,1)                                  &
2987          &         * ( fc00*absa(ind01,j) + fc10*absa(ind02,j) +            &
2988          &             fc01*absa(ind11,j) + fc11*absa(ind12,j) +            &
2989          &             selffac(k)*(selfref(inds,j) + selffrac(k) *          &
2990          &            ( selfref(inds+1,j) - selfref(inds,j)) ) +            &
2991          &             forfac(k) * forref(j) )
2992     
2993               pfrac(NS02+j,k) = fracrefa(j,ifrac) + fracint                 &
2994          &         * (fracrefa(j,ifrac-1) - fracrefa(j,ifrac))
2995             enddo
2996           enddo
2997     
2998           do k = laytrop+1, NLAY
2999             ifp = max( 0, int(2.e2*(fac11(k)+fac01(k))+0.5) )
3000             fc00 = fac00(k) * corr2(ifp)
3001             fc10 = fac10(k) * corr2(ifp)
3002             fc01 = fac01(k) * corr1(ifp)
3003             fc11 = fac11(k) * corr1(ifp)
3004     
3005             ind01 = id0(k,2) + 1
3006             ind02 = ind01 + 1
3007             ind11 = id1(k,2) + 1
3008             ind12 = ind11 + 1
3009     
3010             do j = 1, NG02
3011               taug(NS02+j,k) = colamt(k,1)                                  &
3012          &         * ( fc00*absb(ind01,j) + fc10*absb(ind02,j) +            &
3013          &             fc01*absb(ind11,j) + fc11*absb(ind12,j) +            &
3014          &             forfac(k) * forref(j) )
3015     
3016               pfrac(NS02+j,k) = fracrefb(j)
3017             enddo
3018           enddo
3019     
3020           return
3021     !...................................
3022           end subroutine taugb02
3023     !-----------------------------------
3024     
3025     
3026     !-----------------------------------
3027           subroutine taugb03
3028     !...................................
3029     
3030     !     band 3:  500-630 cm-1 (low - h2o,co2; high - h2o,co2)
3031     !
3032           use module_radlw_kgb03
3033     !
3034           implicit none
3035     !
3036           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
3037          &           ind14, inds, js, ns, jp0, jp1
3038     !
3039           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
3040          &      fac011, fac101, fac111, strrat, speccomb, specmult,         &
3041          &      fs, fs1, fp, ratio, n2omult, tem0, tem1, tem2
3042     
3043           strrat = 1.19268
3044     
3045     !     compute the optical depth by interpolating in ln(pressure), 
3046     !     temperature, and appropriate species.  below laytrop, the water
3047     !     vapor self-continuum is interpolated (in temperature) separately.  
3048     
3049           do k = 1, laytrop
3050             speccomb = colamt(k,1) + strrat*colamt(k,2)
3051             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
3052     
3053             js = 1 + int(specmult)
3054             fs = specmult - int(specmult)
3055             if (js == 8) then
3056               if (fs >= 0.9) then
3057                 js = 9
3058                 fs = 10.0 * (fs - 0.9)
3059               else
3060                 fs = fs / 0.9
3061               endif
3062             endif
3063     
3064             fs1 = 1.0 - fs
3065             fac000 = fs1 * fac00(k)
3066             fac010 = fs1 * fac10(k)
3067             fac100 = fs  * fac00(k)
3068             fac110 = fs  * fac10(k)
3069             fac001 = fs1 * fac01(k)
3070             fac011 = fs1 * fac11(k)
3071             fac101 = fs  * fac01(k)
3072             fac111 = fs  * fac11(k)
3073     
3074             ind01 = id0(k,3) + js
3075             ind02 = ind01 + 1
3076             ind03 = ind01 + 10
3077             ind04 = ind01 + 11
3078             ind11 = min(MSA03, id1(k,3) + js )
3079             ind12 = min(MSA03, ind11 + 1 )
3080             ind13 = min(MSA03, ind11 + 10)
3081             ind14 = min(MSA03, ind11 + 11)
3082             inds = indself(k)
3083     
3084             jp0 = jp(k)
3085             jp1 = jp0 + 1
3086             ns = js + int(fs + 0.5)
3087             if (ns == 10) then
3088               tem1 = n2oref(jp0) / h2oref(jp0)
3089               tem2 = n2oref(jp1) / h2oref(jp1)
3090             else
3091               tem0 = (1.0 - etaref(ns)) / strrat
3092               tem1 = tem0 * n2oref(jp0) / co2ref(jp0)
3093               tem2 = tem0 * n2oref(jp1) / co2ref(jp1)
3094             endif
3095             ratio = tem1 + (fac01(k) + fac11(k)) * (tem2 - tem1)
3096             n2omult = colamt(k,4) - speccomb*ratio
3097     
3098             do j = 1, NG03
3099               taug(NS03+j,k) = speccomb                                     &
3100          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
3101          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
3102          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
3103          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
3104          &         + colamt(k,1) * ( selffac(k)*(selfref(inds,j) +          &
3105          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) ) +  &
3106          &             forfac(k) * forref(j) )                              &
3107          &         + n2omult * absn2oa(j)
3108     
3109               pfrac(NS03+j,k) = fracrefa(j,js)                              &
3110          &          + fs * (fracrefa(j,js+1) - fracrefa(j,js))
3111             enddo
3112           enddo
3113     
3114           do k = laytrop+1, NLAY
3115             speccomb = colamt(k,1) + strrat*colamt(k,2)
3116             specmult = 4.0 * min(oneminus, colamt(k,1)/speccomb)
3117     
3118             js = 1 + int(specmult)
3119             fs = specmult - int(specmult)
3120     
3121             fs1 = 1.0 - fs
3122             fac000 = fs1 * fac00(k)
3123             fac010 = fs1 * fac10(k)
3124             fac100 = fs  * fac00(k)
3125             fac110 = fs  * fac10(k)
3126             fac001 = fs1 * fac01(k)
3127             fac011 = fs1 * fac11(k)
3128             fac101 = fs  * fac01(k)
3129             fac111 = fs  * fac11(k)
3130     
3131             ind01 = id0(k,3) + js
3132             ind02 = ind01 + 1
3133             ind03 = ind01 + 5
3134             ind04 = ind01 + 6
3135             ind11 = id1(k,3) + js
3136             ind12 = ind11 + 1
3137             ind13 = ind11 + 5
3138             ind14 = ind11 + 6
3139     
3140             jp0 = jp(k)
3141             jp1 = jp0 + 1
3142             ns = js + int(fs + 0.5)
3143             if (ns == 5) then
3144               tem1 = n2oref(jp0) / h2oref(jp0)
3145               tem2 = n2oref(jp1) / h2oref(jp1)
3146             else
3147               tem0 = (1.0 - etaref(ns)) / strrat
3148               tem1 = tem0 * n2oref(jp0) / co2ref(jp0)
3149               tem2 = tem0 * n2oref(jp1) / co2ref(jp1)
3150             endif
3151             ratio = tem1 + (fac01(k) + fac11(k)) * (tem2 - tem1)
3152             n2omult = colamt(k,4) - speccomb * ratio
3153     
3154             do j = 1, NG03
3155               taug(NS03+j,k) = speccomb                                     &
3156          &         * ( fac000*absb(ind01,j) + fac100*absb(ind02,j) +        &
3157          &             fac010*absb(ind03,j) + fac110*absb(ind04,j) +        &
3158          &             fac001*absb(ind11,j) + fac101*absb(ind12,j) +        &
3159          &             fac011*absb(ind13,j) + fac111*absb(ind14,j) )        &
3160          &         + colamt(k,1) * forfac(k)*forref(j)                      &
3161          &         + n2omult * absn2ob(j)
3162     
3163               pfrac(NS03+j,k) = fracrefb(j,js)                              &
3164          &          + fs * (fracrefb(j,js+1) - fracrefb(j,js))
3165             enddo
3166           enddo
3167     
3168           return
3169     !...................................
3170           end subroutine taugb03
3171     !-----------------------------------
3172     
3173     
3174     !-----------------------------------
3175           subroutine taugb04
3176     !...................................
3177     
3178     !     band 4:  630-700 cm-1 (low - h2o,co2; high - o3,co2)
3179     !
3180           use module_radlw_kgb04
3181     !
3182           implicit none
3183     !
3184           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
3185          &           ind14, inds, js, ns
3186     !
3187           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
3188          &      fac011, fac101, fac111, strrat1, strrat2, speccomb,         &
3189          &      specmult, fs, fs1
3190     
3191           strrat1 = 850.577
3192           strrat2 = 35.7416
3193     
3194     !     compute the optical depth by interpolating in ln(pressure), 
3195     !     temperature, and appropriate species.  below laytrop, the water
3196     !     vapor self-continuum is interpolated (in temperature) separately.  
3197     
3198           do k = 1, laytrop
3199             speccomb = colamt(k,1) + strrat1*colamt(k,2)
3200             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
3201     
3202             js = 1 + int(specmult)
3203             fs = specmult - int(specmult)
3204     
3205             fs1 = 1.0 - fs
3206             fac000 = fs1 * fac00(k)
3207             fac010 = fs1 * fac10(k)
3208             fac100 = fs  * fac00(k)
3209             fac110 = fs  * fac10(k)
3210             fac001 = fs1 * fac01(k)
3211             fac011 = fs1 * fac11(k)
3212             fac101 = fs  * fac01(k)
3213             fac111 = fs  * fac11(k)
3214     
3215             ind01 = id0(k,4) + js
3216             ind02 = ind01 + 1
3217             ind03 = ind01 + 9
3218             ind04 = ind01 + 10
3219             ind11 = min(MSA04, id1(k,4) + js )
3220             ind12 = min(MSA04, ind11 + 1 )
3221             ind13 = min(MSA04, ind11 + 9 )
3222             ind14 = min(MSA04, ind11 + 10)
3223             inds = indself(k)
3224     
3225             do j = 1, NG04
3226               taug(NS04+j,k) = speccomb                                     &
3227          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
3228          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
3229          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
3230          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
3231          &         + colamt(k,1) * selffac(k)*( selfref(inds,j) +           &
3232          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )
3233     
3234               pfrac(NS04+j,k) = fracrefa(j,js)                              &
3235          &         + fs * (fracrefa(j,js+1) - fracrefa(j,js))
3236             enddo
3237           enddo
3238     
3239           do k = laytrop+1, NLAY
3240             speccomb = colamt(k,3) + strrat2*colamt(k,2)
3241             specmult = 4.0 * min(oneminus, colamt(k,3)/speccomb)
3242     
3243             js = 1 + int(specmult)
3244             fs = specmult - int(specmult)
3245             if (js > 1) then
3246               js = js + 1
3247             elseif (fs >= 0.0024) then
3248               js = 2
3249               fs = (fs - 0.0024) / 0.9976
3250             else
3251               js = 1
3252               fs = fs / 0.0024
3253             endif
3254     
3255             fs1 = 1.0 - fs
3256             fac000 = fs1 * fac00(k)
3257             fac010 = fs1 * fac10(k)
3258             fac100 = fs  * fac00(k)
3259             fac110 = fs  * fac10(k)
3260             fac001 = fs1 * fac01(k)
3261             fac011 = fs1 * fac11(k)
3262             fac101 = fs  * fac01(k)
3263             fac111 = fs  * fac11(k)
3264     
3265             ind01 = id0(k,4) + js
3266             ind02 = ind01 + 1
3267             ind03 = ind01 + 6
3268             ind04 = ind01 + 7
3269             ind11 = id1(k,4) + js
3270             ind12 = ind11 + 1
3271             ind13 = ind11 + 6
3272             ind14 = ind11 + 7
3273     
3274             do j = 1, NG04
3275               taug(NS04+j,k) = speccomb                                     &
3276          &         * ( fac000*absb(ind01,j) + fac100*absb(ind02,j) +        &
3277          &             fac010*absb(ind03,j) + fac110*absb(ind04,j) +        &
3278          &             fac001*absb(ind11,j) + fac101*absb(ind12,j) +        &
3279          &             fac011*absb(ind13,j) + fac111*absb(ind14,j) )
3280     
3281               pfrac(NS04+j,k) = fracrefb(j,js)                              &
3282          &          + fs * (fracrefb(j,js+1) - fracrefb(j,js))
3283             enddo
3284     
3285     !     empirical modification to code to improve stratospheric cooling rates
3286     !     for co2. revised to apply weighting for g-point reduction in this band.
3287     !     from mike iacono, aer, april 01, 2003.
3288     
3289             taug(NS04+8, k) = taug(NS04+8, k) * 0.92
3290             taug(NS04+9, k) = taug(NS04+9, k) * 0.88
3291             taug(NS04+10,k) = taug(NS04+10,k) * 1.07
3292             taug(NS04+11,k) = taug(NS04+11,k) * 1.10
3293             taug(NS04+12,k) = taug(NS04+12,k) * 0.99
3294             taug(NS04+13,k) = taug(NS04+13,k) * 0.88
3295             taug(NS04+14,k) = taug(NS04+14,k) * 0.943
3296     
3297           enddo
3298     
3299           return
3300     !...................................
3301           end subroutine taugb04
3302     !-----------------------------------
3303     
3304     
3305     !-----------------------------------
3306           subroutine taugb05
3307     !...................................
3308     
3309     !     band 5:  700-820 cm-1 (low - h2o,co2; high - o3,co2)
3310     !
3311           use module_radlw_kgb05
3312     !
3313           implicit none
3314     !
3315           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
3316          &           ind14, inds, js, ns
3317     !
3318           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
3319          &      fac011, fac101, fac111, strrat1, strrat2, speccomb,         &
3320          &      specmult, fs, fs1
3321     
3322           strrat1 = 90.4894
3323           strrat2 = 0.900502
3324     
3325     !     compute the optical depth by interpolating in ln(pressure), 
3326     !     temperature, and appropriate species.  below laytrop, the water
3327     !     vapor self-continuum is interpolated (in temperature) separately.  
3328     
3329           do k = 1, laytrop
3330             speccomb = colamt(k,1) + strrat1*colamt(k,2)
3331             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
3332     
3333             js = 1 + int(specmult)
3334             fs = specmult - int(specmult)
3335     
3336             fs1 = 1.0 - fs
3337             fac000 = fs1 * fac00(k)
3338             fac010 = fs1 * fac10(k)
3339             fac100 = fs  * fac00(k)
3340             fac110 = fs  * fac10(k)
3341             fac001 = fs1 * fac01(k)
3342             fac011 = fs1 * fac11(k)
3343             fac101 = fs  * fac01(k)
3344             fac111 = fs  * fac11(k)
3345     
3346             ind01 = id0(k,5) + js
3347             ind02 = ind01 + 1
3348             ind03 = ind01 + 9
3349             ind04 = ind01 + 10
3350             ind11 = min(MSA05, id1(k,5) + js )
3351             ind12 = min(MSA05, ind11 + 1 )
3352             ind13 = min(MSA05, ind11 + 9 )
3353             ind14 = min(MSA05, ind11 + 10)
3354             inds = indself(k)
3355     
3356             do j = 1, NG05
3357               taug(NS05+j,k) = speccomb                                     &
3358          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
3359          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
3360          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
3361          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
3362          &         + colamt(k,1) * selffac(k)*( selfref(inds,j) +           &
3363          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )    &
3364          &         + wx(k,1) * ccl4(j)
3365     
3366               pfrac(NS05+j,k) = fracrefa(j,js)                              &
3367          &          + fs * (fracrefa(j,js+1) - fracrefa(j,js))
3368             enddo
3369           enddo
3370     
3371           do k = laytrop+1, NLAY
3372             speccomb = colamt(k,3) + strrat2*colamt(k,2)
3373             specmult = 4.0 * min(oneminus, colamt(k,3)/speccomb)
3374     
3375             js = 1 + int(specmult)
3376             fs = specmult - int(specmult)
3377     
3378             fs1 = 1.0 - fs
3379             fac000 = fs1 * fac00(k)
3380             fac010 = fs1 * fac10(k)
3381             fac100 = fs  * fac00(k)
3382             fac110 = fs  * fac10(k)
3383             fac001 = fs1 * fac01(k)
3384             fac011 = fs1 * fac11(k)
3385             fac101 = fs  * fac01(k)
3386             fac111 = fs  * fac11(k)
3387     
3388             ind01 = id0(k,5) + js
3389             ind02 = ind01 + 1
3390             ind03 = ind01 + 5
3391             ind04 = ind01 + 6
3392             ind11 = id1(k,5) + js
3393             ind12 = ind11 + 1
3394             ind13 = ind11 + 5
3395             ind14 = ind11 + 6
3396     
3397             do j = 1, NG05
3398               taug(NS05+j,k) = speccomb                                     &
3399          &         * ( fac000*absb(ind01,j) + fac100*absb(ind02,j) +        &
3400          &             fac010*absb(ind03,j) + fac110*absb(ind04,j) +        &
3401          &             fac001*absb(ind11,j) + fac101*absb(ind12,j) +        &
3402          &             fac011*absb(ind13,j) + fac111*absb(ind14,j) )        &
3403          &         + wx(k,1) * ccl4(j)
3404     
3405               pfrac(NS05+j,k) = fracrefb(j,js)                              &
3406          &          + fs * (fracrefb(j,js+1) - fracrefb(j,js))
3407             enddo
3408           enddo
3409     
3410           return
3411     !...................................
3412           end subroutine taugb05
3413     !-----------------------------------
3414     
3415     
3416     !-----------------------------------
3417           subroutine taugb06
3418     !...................................
3419     
3420     !     band 6:  820-980 cm-1 (low - h2o; high - nothing)
3421     !
3422           use module_radlw_kgb06
3423     !
3424           implicit none
3425     !
3426           integer :: j, k, ind01, ind02, ind11, ind12, inds, js, ns
3427     
3428     !     compute the optical depth by interpolating in ln(pressure) and
3429     !     temperature. the water vapor self-continuum is interpolated
3430     !     (in temperature) separately.  
3431     
3432           do k = 1, laytrop
3433             ind01 = id0(k,6) + 1
3434             ind02 = ind01 + 1
3435             ind11 = min(MSA06, id1(k,6) + 1 )
3436             ind12 = min(MSA06, ind11 + 1 )
3437             inds = indself(k)
3438     
3439             do j = 1, NG06
3440               taug(NS06+j,k) = colamt(k,1)                                  &
3441          &         * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) +    &
3442          &             fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) +    &
3443          &             selffac(k) *( selfref(inds,j) +                      &
3444          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) ) )  &
3445          &         + wx(k,2) * cfc11adj(j) + wx(k,3) * cfc12(j)             &
3446          &         + co2mult(k) * absco2(j)
3447     
3448               pfrac(NS06+j,k) = fracrefa(j)
3449             enddo
3450           enddo
3451     
3452     !     nothing important goes on above laytrop in this band.
3453           do k = laytrop+1, NLAY
3454             do j = 1, NG06
3455               taug(NS06+j,k) = wx(k,2) * cfc11adj(j) + wx(k,3) * cfc12(j)
3456     
3457               pfrac(NS06+j,k) = fracrefa(j)
3458             enddo
3459           enddo
3460     
3461           return
3462     !...................................
3463           end subroutine taugb06
3464     !-----------------------------------
3465     
3466     
3467     !-----------------------------------
3468           subroutine taugb07
3469     !...................................
3470     
3471     !     band 7:  980-1080 cm-1 (low - h2o,o3; high - o3)
3472     !
3473           use module_radlw_kgb07
3474     !
3475           implicit none
3476     !
3477           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
3478          &           ind14, inds, js
3479     !
3480           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
3481          &      fac011, fac101, fac111, strrat, speccomb, specmult, fs, fs1
3482     
3483           strrat = 8.21104e4
3484     
3485     !     compute the optical depth by interpolating in ln(pressure), 
3486     !     temperature, and appropriate species.  below laytrop, the water
3487     !     vapor self-continuum is interpolated (in temperature) separately.  
3488     
3489           do k = 1, laytrop
3490             speccomb = colamt(k,1) + strrat*colamt(k,3)
3491             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
3492     
3493             js = 1 + int(specmult)
3494             fs = specmult - int(specmult)
3495     
3496             fs1 = 1.0 - fs
3497             fac000 = fs1 * fac00(k)
3498             fac010 = fs1 * fac10(k)
3499             fac100 = fs  * fac00(k)
3500             fac110 = fs  * fac10(k)
3501             fac001 = fs1 * fac01(k)
3502             fac011 = fs1 * fac11(k)
3503             fac101 = fs  * fac01(k)
3504             fac111 = fs  * fac11(k)
3505     
3506             ind01 = id0(k,7) + js
3507             ind02 = ind01 + 1
3508             ind03 = ind01 + 9
3509             ind04 = ind01 + 10
3510             ind11 = min(MSA07, id1(k,7) + js )
3511             ind12 = min(MSA07, ind11 + 1 )
3512             ind13 = min(MSA07, ind11 + 9 )
3513             ind14 = min(MSA07, ind11 + 10)
3514             inds = indself(k)
3515     
3516             do j = 1, NG07
3517               taug(NS07+j,k) = speccomb                                     &
3518          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
3519          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
3520          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
3521          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
3522          &         + colamt(k,1) * selffac(k)*( selfref(inds,j) +           &
3523          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )    &
3524          &         + co2mult(k) * absco2(j)
3525     
3526               pfrac(NS07+j,k) = fracrefa(j,js)                              &
3527          &              + fs * (fracrefa(j,js+1) - fracrefa(j,js))
3528             enddo
3529           enddo
3530     
3531           do k = laytrop+1, NLAY
3532     
3533             ind01 = id0(k,7) + 1
3534             ind02 = ind01 + 1
3535             ind11 = id1(k,7) + 1
3536             ind12 = ind11 + 1
3537     
3538             do j = 1, NG07
3539               taug(NS07+j,k) = colamt(k,3)                                  &
3540          &         * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) +    &
3541          &             fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )    &
3542          &         + co2mult(k) * absco2(j)
3543     
3544               pfrac(NS07+j,k) = fracrefb(j)
3545             enddo
3546     
3547     !     empirical modification to code to improve stratospheric cooling rates
3548     !     for o3. revised to apply weighting for g-point reduction in this band.
3549     !     from mike iacono, aer, april 01, 2003.
3550     
3551             taug(NS07+6, k) = taug(NS07+6, k) * 0.92
3552             taug(NS07+7, k) = taug(NS07+7, k) * 0.88
3553             taug(NS07+8, k) = taug(NS07+8, k) * 1.07
3554             taug(NS07+9, k) = taug(NS07+9, k) * 1.10
3555             taug(NS07+10,k) = taug(NS07+10,k) * 0.99
3556             taug(NS07+11,k) = taug(NS07+11,k) * 0.855
3557     
3558           enddo
3559     
3560           return
3561     !...................................
3562           end subroutine taugb07
3563     !-----------------------------------
3564     
3565     
3566     !-----------------------------------
3567           subroutine taugb08
3568     !...................................
3569     
3570     !     band 8:  1080-1180 cm-1 (low (i.e.>~300mb) - h2o; high - o3)
3571     !
3572           use module_radlw_kgb08
3573     !
3574           implicit none
3575     !
3576           integer :: j, k, ind01, ind02, ind11, ind12, inds, jp0, jp1
3577     !
3578           real (kind=kind_phys) :: ratio, n2omult, tem1, tem2
3579     
3580     !     compute the optical depth by interpolating in ln(pressure) and 
3581     !     temperature.  
3582     
3583           do k = 1, layswtch
3584             ind01 = id0(k,8) + 1
3585             ind02 = ind01 + 1
3586             ind11 = min(MSA08, id1(k,8) + 1 )
3587             ind12 = min(MSA08, ind11 + 1 )
3588             inds = indself(k)
3589     
3590             jp0 = jp(k)
3591             jp1 = jp0 + 1
3592             tem1 = n2oref(jp0) / h2oref(jp0)
3593             tem2 = n2oref(jp1) / h2oref(jp1)
3594             ratio = tem1 + (fac01(k) + fac11(k)) * (tem2 - tem1)
3595             n2omult = colamt(k,4) - colamt(K,1)*ratio
3596     
3597             do j = 1, NG08
3598               taug(NS08+j,k) = colamt(k,1)                                  &
3599          &         * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) +    &
3600          &             fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) +    &
3601          &             selffac(k)*( selfref(inds,j) +                       &
3602          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) ) )  &
3603          &         + wx(k,3) * cfc12(j) + wx(k,4) * cfc22adj(j)             &
3604          &         + co2mult(k) * absco2a(j) + n2omult * absn2oa(j)
3605     
3606               pfrac(NS08+j,k) = fracrefa(j)
3607             enddo
3608           enddo
3609     
3610           do k = layswtch+1, NLAY
3611             ind01 = id0(k,8) + 1
3612             ind02 = ind01 + 1
3613             ind11 = id1(k,8) + 1
3614             ind12 = ind11 + 1
3615     
3616             jp0 = jp(k)
3617             jp1 = jp0 + 1
3618             tem1 = n2oref(jp0) / o3ref(jp0)
3619             tem2 = n2oref(jp1) / o3ref(jp1)
3620             ratio = tem1 + (fac01(k) + fac11(k)) * (tem2 - tem1)
3621             n2omult = colamt(k,4) - colamt(k,3) * ratio
3622     
3623             do j = 1, NG08
3624               taug(NS08+j,k) = colamt(k,3)                                  &
3625          &         * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) +    &
3626          &             fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )    &
3627          &         + wx(k,3) * cfc12(j) + wx(k,4) * cfc22adj(j)             &
3628          &         + co2mult(k) * absco2b(j) + n2omult * absn2ob(j)
3629     
3630               pfrac(NS08+j,k) = fracrefb(j)
3631             enddo
3632           enddo
3633     
3634           return
3635     !...................................
3636           end subroutine taugb08
3637     !-----------------------------------
3638     
3639     
3640     !-----------------------------------
3641           subroutine taugb09
3642     !...................................
3643     
3644     !     band 9:  1180-1390 cm-1 (low - h2o,ch4; high - ch4)
3645     !
3646           use module_radlw_kgb09
3647     !
3648           implicit none
3649     !
3650           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
3651          &           ind14, inds, js, ns, jfrac, ioff, jp0, jp1
3652     !
3653           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
3654          &      fac011, fac101, fac111, strrat, speccomb, specmult, fs, fs1,&
3655          &      ffrac, ratio, n2omult, tem0, tem1, tem2
3656     
3657           strrat = 21.6282
3658           ioff = 0
3659     
3660     !     compute the optical depth by interpolating in ln(pressure), 
3661     !     temperature, and appropriate species.  below laytrop, the water
3662     !     vapor self-continuum is interpolated (in temperature) separately.  
3663     
3664           do k = 1, laytrop
3665             if (k == laylow) ioff = NG09
3666             if (k == layswtch) ioff = 2 * NG09
3667     
3668             speccomb = colamt(k,1) + strrat*colamt(k,5)
3669             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
3670     
3671             js = 1 + int(specmult)
3672             jfrac = js
3673             fs = specmult - int(specmult)
3674             ffrac = fs
3675             if (js == 8) then
3676               if (fs <= 0.68) then
3677                 fs = fs / 0.68
3678               elseif (fs <= 0.92) then
3679                 js = js + 1
3680                 fs = (fs - 0.68) / 0.24
3681               else
3682                 js = js + 2
3683                 fs = (fs - 0.92) / 0.08
3684               endif
3685             elseif (js == 9) then
3686               js = 10
3687               fs = 1.0
3688               jfrac = 8
3689               ffrac = 1.0
3690             endif
3691     
3692             fs1 = 1.0 - fs
3693             fac000 = fs1 * fac00(k)
3694             fac010 = fs1 * fac10(k)
3695             fac100 = fs  * fac00(k)
3696             fac110 = fs  * fac10(k)
3697             fac001 = fs1 * fac01(k)
3698             fac011 = fs1 * fac11(k)
3699             fac101 = fs  * fac01(k)
3700             fac111 = fs  * fac11(k)
3701     
3702             ind01 = id0(k,9) + js
3703             ind02 = ind01 + 1
3704             ind03 = ind01 + 11
3705             ind04 = ind01 + 12
3706             ind11 = min(MSA09, id1(k,9) + js )
3707             ind12 = min(MSA09, ind11 + 1 )
3708             ind13 = min(MSA09, ind11 + 11)
3709             ind14 = min(MSA09, ind11 + 12)
3710             inds = indself(k)
3711     
3712             jp0 = jp(k)
3713             jp1 = jp0 + 1
3714             ns = js + int(fs + 0.5)
3715             tem0 = (1.0 - etaref(ns)) / strrat
3716             if (ns == 11) then
3717               tem1 = n2oref(jp0) / h2oref(jp0)
3718               tem2 = n2oref(jp1) / h2oref(jp1)
3719             else
3720               tem1 = tem0 * n2oref(jp0) / ch4ref(jp0)
3721               tem2 = tem0 * n2oref(jp1) / ch4ref(jp1)
3722             endif
3723             ratio = tem1 + (fac01(k) + fac11(k)) * (tem2 - tem1)
3724             n2omult = colamt(k,4) - speccomb*ratio
3725     
3726             do j = 1, NG09
3727               taug(NS09+j,k) = speccomb                                     &
3728          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
3729          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
3730          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
3731          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
3732          &         + colamt(k,1) * selffac(k)*( selfref(inds,j) +           &
3733          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )    &
3734          &         + n2omult * absn2o(j+ioff)
3735     
3736               pfrac(NS09+j,k) = fracrefa(j,jfrac)                           &
3737          &          + ffrac * (fracrefa(j,jfrac+1)-fracrefa(j,jfrac))
3738             enddo
3739           enddo
3740     
3741           do k = laytrop+1, NLAY
3742             ind01 = id0(k,9) + 1
3743             ind02 = ind01 + 1
3744             ind11 = id1(k,9) + 1
3745             ind12 = ind11 + 1
3746     
3747             do j = 1, NG09
3748               taug(NS09+j,k) = colamt(k,5)                                  &
3749          &         * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) +    &
3750          &             fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) )
3751     
3752               pfrac(NS09+j,k) = fracrefb(j)
3753             enddo
3754           enddo
3755     
3756           return
3757     !...................................
3758           end subroutine taugb09
3759     !-----------------------------------
3760     
3761     
3762     !-----------------------------------
3763           subroutine taugb10
3764     !...................................
3765     
3766     !     band 10:  1390-1480 cm-1 (low - h2o; high - h2o)
3767     !
3768           use module_radlw_kgb10
3769     !
3770           implicit none
3771     !
3772           integer :: j, k, ind01, ind02, ind11, ind12
3773     !
3774     !     compute the optical depth by interpolating in ln(pressure) and 
3775     !     temperature.  
3776     
3777           do k = 1, laytrop
3778             ind01 = id0(k,10) + 1
3779             ind02 = ind01 + 1
3780             ind11 = min(MSA10, id1(k,10) + 1 )
3781             ind12 = min(MSA10, ind11 + 1 )
3782     
3783             do j = 1, NG10
3784               taug(NS10+j,k) = colamt(k,1)                                  &
3785          &         * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) +    &
3786          &             fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) ) 
3787     
3788               pfrac(NS10+j,k) = fracrefa(j)
3789             enddo
3790           enddo
3791     
3792           do k = laytrop+1, NLAY
3793             ind01 = id0(k,10) + 1
3794             ind02 = ind01 + 1
3795             ind11 = id1(k,10) + 1
3796             ind12 = ind11 + 1
3797     
3798             do j = 1, NG10
3799               taug(NS10+j,k) = colamt(k,1)                                  &
3800          &         * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) +    &
3801          &             fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) 
3802     
3803               pfrac(NS10+j,k) = fracrefb(j)
3804             enddo
3805           enddo
3806     
3807           return
3808     !...................................
3809           end subroutine taugb10
3810     !-----------------------------------
3811     
3812     
3813     !-----------------------------------
3814           subroutine taugb11
3815     !...................................
3816     
3817     !     band 11:  1480-1800 cm-1 (low - h2o; high - h2o)
3818     !
3819           use module_radlw_kgb11
3820     !
3821           implicit none
3822     !
3823           integer :: j, k, ind01, ind02, ind11, ind12, inds
3824     !
3825     
3826     !     compute the optical depth by interpolating in ln(pressure) and 
3827     !     temperature.  below laytrop, the water vapor self-continuum 
3828     !     is interpolated (in temperature) separately.  
3829     
3830           do k = 1, laytrop
3831             ind01 = id0(k,11) + 1
3832             ind02 = ind01 + 1
3833             ind11 = min(MSA11, id1(k,11) + 1 )
3834             ind12 = min(MSA11, ind11 + 1 )
3835             inds = indself(k)
3836     
3837             do j = 1, NG11
3838               taug(NS11+j,k) = colamt(k,1)                                  &
3839          &         * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) +    &
3840          &             fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) +    &
3841          &             selffac(k)*( selfref(inds,j) +  selffrac(k) *        &
3842          &            (selfref(inds+1,j)-selfref(inds,j)) ) )
3843     
3844               pfrac(NS11+j,k) = fracrefa(j)
3845             enddo
3846           enddo
3847     
3848           do k = laytrop+1, NLAY
3849             ind01 = id0(k,11) + 1
3850             ind02 = ind01 + 1
3851             ind11 = id1(k,11) + 1
3852             ind12 = ind11 + 1
3853     
3854             do j = 1, NG11
3855               taug(NS11+j,k) = colamt(k,1)                                  &
3856          &         * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) +    &
3857          &             fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) 
3858     
3859               pfrac(NS11+j,k) = fracrefb(j)
3860             enddo
3861           enddo
3862     
3863           return
3864     !...................................
3865           end subroutine taugb11
3866     !-----------------------------------
3867     
3868     
3869     !-----------------------------------
3870           subroutine taugb12
3871     !...................................
3872     
3873     !     band 12:  1800-2080 cm-1 (low - h2o,co2; high - nothing)
3874     !
3875           use module_radlw_kgb12
3876     !
3877           implicit none
3878     !
3879           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
3880          &           ind14, inds, js
3881     !
3882           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
3883          &      fac011, fac101, fac111, strrat, speccomb, specmult, fs, fs1
3884     
3885           strrat = 0.009736757
3886     
3887     !     compute the optical depth by interpolating in ln(pressure), 
3888     !     temperature, and appropriate species.  below laytrop, the water
3889     !     vapor self-continuum is interpolated (in temperature) separately.  
3890     
3891           do k = 1, laytrop
3892             speccomb = colamt(k,1) + strrat*colamt(k,2)
3893             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
3894     
3895             js = 1 + int(specmult)
3896             fs = specmult - int(specmult)
3897     
3898             fs1 = 1.0 - fs
3899             fac000 = fs1 * fac00(k)
3900             fac010 = fs1 * fac10(k)
3901             fac100 = fs  * fac00(k)
3902             fac110 = fs  * fac10(k)
3903             fac001 = fs1 * fac01(k)
3904             fac011 = fs1 * fac11(k)
3905             fac101 = fs  * fac01(k)
3906             fac111 = fs  * fac11(k)
3907     
3908             ind01 = id0(k,12) + js
3909             ind02 = ind01 + 1
3910             ind03 = ind01 + 9
3911             ind04 = ind01 + 10
3912             ind11 = min(MSA12, id1(k,12) + js )
3913             ind12 = min(MSA12, ind11 + 1 )
3914             ind13 = min(MSA12, ind11 + 9 )
3915             ind14 = min(MSA12, ind11 + 10)
3916             inds = indself(k)
3917     
3918             do j = 1, NG12
3919               taug(NS12+j,k) = speccomb                                     &
3920          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
3921          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
3922          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
3923          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
3924          &         + colamt(k,1)*selffac(k)*( selfref(inds,j) +             &
3925          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )
3926     
3927               pfrac(NS12+j,k) = fracrefa(j,js)                              &
3928          &          + fs * (fracrefa(j,js+1) - fracrefa(j,js))
3929             enddo
3930           enddo
3931     
3932           do k = laytrop+1, NLAY
3933             do j = 1, NG12
3934               taug(NS12+j,k) = f_zero
3935               pfrac(NS12+j,k) = f_zero
3936             enddo
3937           enddo
3938     
3939           return
3940     !...................................
3941           end subroutine taugb12
3942     !-----------------------------------
3943     
3944     
3945     !-----------------------------------
3946           subroutine taugb13
3947     !...................................
3948     
3949     !     band 13:  2080-2250 cm-1 (low - h2o,n2o; high - nothing)
3950     !
3951           use module_radlw_kgb13
3952     !
3953           implicit none
3954     !
3955           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
3956          &           ind14, inds, js
3957     !
3958           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
3959          &      fac011, fac101, fac111, strrat, speccomb, specmult, fs, fs1
3960     
3961           strrat = 16658.87
3962     
3963     !     compute the optical depth by interpolating in ln(pressure), 
3964     !     temperature, and appropriate species.  below laytrop, the water
3965     !     vapor self-continuum is interpolated (in temperature) separately.  
3966     
3967           do k = 1, laytrop
3968             speccomb = colamt(k,1) + strrat*colamt(k,4)
3969             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
3970     
3971             js = 1 + int(specmult)
3972             fs = specmult - int(specmult)
3973     
3974             fs1 = 1.0 - fs
3975             fac000 = fs1 * fac00(k)
3976             fac010 = fs1 * fac10(k)
3977             fac100 = fs  * fac00(k)
3978             fac110 = fs  * fac10(k)
3979             fac001 = fs1 * fac01(k)
3980             fac011 = fs1 * fac11(k)
3981             fac101 = fs  * fac01(k)
3982             fac111 = fs  * fac11(k)
3983     
3984             ind01 = id0(k,13) + js
3985             ind02 = ind01 + 1
3986             ind03 = ind01 + 9
3987             ind04 = ind01 + 10
3988             ind11 = min(MSA13, id1(k,13) + js )
3989             ind12 = min(MSA13, ind11 + 1 )
3990             ind13 = min(MSA13, ind11 + 9 )
3991             ind14 = min(MSA13, ind11 + 10)
3992             inds = indself(k)
3993     
3994             do j = 1, NG13
3995               taug(NS13+j,k) = speccomb                                     &
3996          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
3997          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
3998          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
3999          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
4000          &         + colamt(k,1)*selffac(k)*( selfref(inds,j) +             &
4001          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )
4002     
4003               pfrac(NS13+j,k) = fracrefa(j,js)                              &
4004          &          + fs * (fracrefa(j,js+1) - fracrefa(j,js))
4005             enddo
4006           enddo
4007     
4008           do k = laytrop+1, NLAY
4009             do j = 1, NG13
4010               taug(NS13+j,k) = f_zero
4011               pfrac(NS13+j,k) = f_zero
4012             enddo
4013           enddo
4014     
4015           return
4016     !...................................
4017           end subroutine taugb13
4018     !-----------------------------------
4019     
4020     
4021     !-----------------------------------
4022           subroutine taugb14
4023     !...................................
4024     
4025     !     band 14:  2250-2380 cm-1 (low - co2; high - co2)
4026     !
4027           use module_radlw_kgb14
4028     !
4029           implicit none
4030     !
4031           integer :: j, k, ind01, ind02, ind11, ind12, inds
4032     !
4033     
4034     !     compute the optical depth by interpolating in ln(pressure) and 
4035     !     temperature.  below laytrop, the water vapor self-continuum 
4036     !     is interpolated (in temperature) separately.  
4037     
4038           do k = 1, laytrop
4039             ind01 = id0(k,14) + 1
4040             ind02 = ind01 + 1
4041             ind11 = min(MSA14, id1(k,14) + 1 )
4042             ind12 = min(MSA14, ind11 + 1 )
4043             inds = indself(k)
4044     
4045             do j = 1, NG14
4046               taug(NS14+j,k) = colamt(k,2)                                  &
4047          &         * ( fac00(k)*absa(ind01,j) + fac10(k)*absa(ind02,j) +    &
4048          &             fac01(k)*absa(ind11,j) + fac11(k)*absa(ind12,j) )    &
4049          &         + colamt(k,1)*selffac(k)*( selfref(inds,j) +             &
4050          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )
4051     
4052               pfrac(NS14+j,k) = fracrefa(j)
4053             enddo
4054           enddo
4055     
4056           do k = laytrop+1, NLAY
4057             ind01 = id0(k,14) + 1
4058             ind02 = ind01 + 1
4059             ind11 = id1(k,14) + 1
4060             ind12 = ind11 + 1
4061     
4062             do j = 1, NG14
4063               taug(NS14+j,k) = colamt(k,2)                                  &
4064          &         * ( fac00(k)*absb(ind01,j) + fac10(k)*absb(ind02,j) +    &
4065          &             fac01(k)*absb(ind11,j) + fac11(k)*absb(ind12,j) ) 
4066     
4067               pfrac(NS14+j,k) = fracrefb(j)
4068             enddo
4069           enddo
4070     
4071           return
4072     !...................................
4073           end subroutine taugb14
4074     !-----------------------------------
4075     
4076     
4077     !-----------------------------------
4078           subroutine taugb15
4079     !...................................
4080     
4081     !     band 15:  2380-2600 cm-1 (low - n2o,co2; high - nothing)
4082     !
4083           use module_radlw_kgb15
4084     !
4085           implicit none
4086     !
4087           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
4088          &           ind14, inds, js
4089     !
4090           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
4091          &      fac011, fac101, fac111, strrat, speccomb, specmult, fs, fs1
4092     
4093           strrat = 0.2883201
4094     
4095     !     compute the optical depth by interpolating in ln(pressure), 
4096     !     temperature, and appropriate species.  below laytrop, the water
4097     !     vapor self-continuum is interpolated (in temperature) separately.  
4098     
4099           do k = 1, laytrop
4100             speccomb = colamt(k,4) + strrat*colamt(k,2)
4101             specmult = 8.0 * min(oneminus, colamt(k,4)/speccomb)
4102     
4103             js = 1 + int(specmult)
4104             fs = specmult - int(specmult)
4105     
4106             fs1 = 1.0 - fs
4107             fac000 = fs1 * fac00(k)
4108             fac010 = fs1 * fac10(k)
4109             fac100 = fs  * fac00(k)
4110             fac110 = fs  * fac10(k)
4111             fac001 = fs1 * fac01(k)
4112             fac011 = fs1 * fac11(k)
4113             fac101 = fs  * fac01(k)
4114             fac111 = fs  * fac11(k)
4115     
4116             ind01 = id0(k,15) + js
4117             ind02 = ind01 + 1
4118             ind03 = ind01 + 9
4119             ind04 = ind01 + 10
4120             ind11 = min(MSA15, id1(k,15) + js )
4121             ind12 = min(MSA15, ind11 + 1 )
4122             ind13 = min(MSA15, ind11 + 9 )
4123             ind14 = min(MSA15, ind11 + 10)
4124             inds = indself(k)
4125     
4126             do j = 1, NG15
4127               taug(NS15+j,k) = speccomb                                     &
4128          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
4129          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
4130          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
4131          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
4132          &         + colamt(k,1)*selffac(k)*( selfref(inds,j) +             &
4133          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )
4134     
4135               pfrac(NS15+j,k) = fracrefa(j,js)                              &
4136          &          + fs * (fracrefa(j,js+1) - fracrefa(j,js))
4137             enddo
4138           enddo
4139     
4140           do k = laytrop+1, NLAY
4141             do j = 1, NG15
4142               taug(NS15+j,k) = f_zero
4143               pfrac(NS15+j,k) = f_zero
4144             enddo
4145           enddo
4146     
4147           return
4148     !...................................
4149           end subroutine taugb15
4150     !-----------------------------------
4151     
4152     
4153     !-----------------------------------
4154           subroutine taugb16
4155     !...................................
4156     
4157     !     band 16:  2600-3000 cm-1 (low - h2o,ch4; high - nothing)
4158     !
4159           use module_radlw_kgb16
4160     !
4161           implicit none
4162     !
4163           integer :: j, k, ind01, ind02, ind03, ind04, ind11, ind12, ind13, &
4164          &           ind14, inds, js
4165     !
4166           real (kind=kind_phys) :: fac000, fac010, fac100, fac110, fac001,  &
4167          &      fac011, fac101, fac111, strrat, speccomb, specmult, fs, fs1
4168     
4169           strrat = 830.411
4170     
4171     !     compute the optical depth by interpolating in ln(pressure), 
4172     !     temperature, and appropriate species.  below laytrop, the water
4173     !     vapor self-continuum is interpolated (in temperature) separately.  
4174     
4175           do k = 1, laytrop
4176             speccomb = colamt(k,1) + strrat*colamt(k,5)
4177             specmult = 8.0 * min(oneminus, colamt(k,1)/speccomb)
4178     
4179             js = 1 + int(specmult)
4180             fs = specmult - int(specmult)
4181     
4182             fs1 = 1.0 - fs
4183             fac000 = fs1 * fac00(k)
4184             fac010 = fs1 * fac10(k)
4185             fac100 = fs  * fac00(k)
4186             fac110 = fs  * fac10(k)
4187             fac001 = fs1 * fac01(k)
4188             fac011 = fs1 * fac11(k)
4189             fac101 = fs  * fac01(k)
4190             fac111 = fs  * fac11(k)
4191     
4192             ind01 = id0(k,16) + js
4193             ind02 = ind01 + 1
4194             ind03 = ind01 + 9
4195             ind04 = ind01 + 10
4196             ind11 = min(MSA16, id1(k,16) + js )
4197             ind12 = min(MSA16, ind11 + 1 )
4198             ind13 = min(MSA16, ind11 + 9 )
4199             ind14 = min(MSA16, ind11 + 10)
4200             inds = indself(k)
4201     
4202             do j = 1, NG16
4203               taug(NS16+j,k) = speccomb                                     &
4204          &         * ( fac000*absa(ind01,j) + fac100*absa(ind02,j) +        &
4205          &             fac010*absa(ind03,j) + fac110*absa(ind04,j) +        &
4206          &             fac001*absa(ind11,j) + fac101*absa(ind12,j) +        &
4207          &             fac011*absa(ind13,j) + fac111*absa(ind14,j) )        &
4208          &         + colamt(k,1)*selffac(k)*( selfref(inds,j) +             &
4209          &             selffrac(k)*(selfref(inds+1,j)-selfref(inds,j)) )
4210     
4211               pfrac(NS16+j,k) = fracrefa(j,js)                              &
4212          &          + fs * (fracrefa(j,js+1) - fracrefa(j,js))
4213             enddo
4214           enddo
4215     
4216           do k = laytrop+1, NLAY
4217             do j = 1, NG16
4218               taug(NS16+j,k) = f_zero
4219               pfrac(NS16+j,k) = f_zero
4220             enddo
4221           enddo
4222     
4223           return
4224     !...................................
4225           end subroutine taugb16
4226     !-----------------------------------
4227     
4228     
4229     !...................................
4230           end subroutine taumol
4231     !-----------------------------------
4232     
4233     
4234     !
4235     !........................................!
4236           end module module_radlw_main       !
4237     !========================================!
4238     
4239