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