source: palm/trunk/LIB/rrtmg/rrtmg_lw_rtrnmc.f90 @ 1585

Last change on this file since 1585 was 1585, checked in by maronga, 6 years ago

Added support for RRTMG radiation code

File size: 26.9 KB
Line 
1!     path:      $Source: /storm/rc1/cvsroot/rc/rrtmg_lw/src/rrtmg_lw_rtrnmc.f90,v $
2!     author:    $Author: mike $
3!     revision:  $Revision: 1.7 $
4!     created:   $Date: 2009/11/12 20:52:25 $
5!
6      module rrtmg_lw_rtrnmc
7
8!  --------------------------------------------------------------------------
9! |                                                                          |
10! |  Copyright 2002-2009, Atmospheric & Environmental Research, Inc. (AER).  |
11! |  This software may be used, copied, or redistributed as long as it is    |
12! |  not sold and this copyright notice is reproduced on each copy made.     |
13! |  This model is provided as is without any express or implied warranties. |
14! |                       (http://www.rtweb.aer.com/)                        |
15! |                                                                          |
16!  --------------------------------------------------------------------------
17
18! --------- Modules ----------
19
20      use parkind, only : im => kind_im, rb => kind_rb
21      use parrrtm, only : mg, nbndlw, ngptlw
22      use rrlw_con, only: fluxfac, heatfac
23      use rrlw_wvn, only: delwave, ngb, ngs
24      use rrlw_tbl, only: tblint, bpade, tau_tbl, exp_tbl, tfn_tbl
25      use rrlw_vsn, only: hvrrtc, hnamrtc
26
27      implicit none
28
29      contains
30
31!-----------------------------------------------------------------------------
32      subroutine rtrnmc(nlayers, istart, iend, iout, pz, semiss, ncbands, &
33                        cldfmc, taucmc, planklay, planklev, plankbnd, &
34                        pwvcm, fracs, taut, &
35                        totuflux, totdflux, fnet, htr, &
36                        totuclfl, totdclfl, fnetc, htrc, &
37                        idrv, dplankbnd_dt, dtotuflux_dt, dtotuclfl_dt )
38!-----------------------------------------------------------------------------
39!
40!  Original version:   E. J. Mlawer, et al. RRTM_V3.0
41!  Revision for GCMs:  Michael J. Iacono; October, 2002
42!  Revision for F90:  Michael J. Iacono; June, 2006
43!  Revision for dFdT option: M. J. Iacono and E. J. Mlawer, November 2009
44!
45!  This program calculates the upward fluxes, downward fluxes, and
46!  heating rates for an arbitrary clear or cloudy atmosphere.  The input
47!  to this program is the atmospheric profile, all Planck function
48!  information, and the cloud fraction by layer.  A variable diffusivity
49!  angle (SECDIFF) is used for the angle integration.  Bands 2-3 and 5-9
50!  use a value for SECDIFF that varies from 1.50 to 1.80 as a function of
51!  the column water vapor, and other bands use a value of 1.66.  The Gaussian
52!  weight appropriate to this angle (WTDIFF=0.5) is applied here.  Note that
53!  use of the emissivity angle for the flux integration can cause errors of
54!  1 to 4 W/m2 within cloudy layers. 
55!  Clouds are treated with the McICA stochastic approach and maximum-random
56!  cloud overlap.
57!  This subroutine also provides the optional capability to calculate
58!  the derivative of upward flux respect to surface temperature using
59!  the pre-tabulated derivative of the Planck function with respect to
60!  temperature integrated over each spectral band.
61!***************************************************************************
62
63! ------- Declarations -------
64
65! ----- Input -----
66      integer(kind=im), intent(in) :: nlayers         ! total number of layers
67      integer(kind=im), intent(in) :: istart          ! beginning band of calculation
68      integer(kind=im), intent(in) :: iend            ! ending band of calculation
69      integer(kind=im), intent(in) :: iout            ! output option flag
70
71! Atmosphere
72      real(kind=rb), intent(in) :: pz(0:)             ! level (interface) pressures (hPa, mb)
73                                                      !    Dimensions: (0:nlayers)
74      real(kind=rb), intent(in) :: pwvcm              ! precipitable water vapor (cm)
75      real(kind=rb), intent(in) :: semiss(:)          ! lw surface emissivity
76                                                      !    Dimensions: (nbndlw)
77      real(kind=rb), intent(in) :: planklay(:,:)      !
78                                                      !    Dimensions: (nlayers,nbndlw)
79      real(kind=rb), intent(in) :: planklev(0:,:)     !
80                                                      !    Dimensions: (0:nlayers,nbndlw)
81      real(kind=rb), intent(in) :: plankbnd(:)        !
82                                                      !    Dimensions: (nbndlw)
83      real(kind=rb), intent(in) :: fracs(:,:)         !
84                                                      !    Dimensions: (nlayers,ngptw)
85      real(kind=rb), intent(in) :: taut(:,:)          ! gaseous + aerosol optical depths
86                                                      !    Dimensions: (nlayers,ngptlw)
87
88! Clouds
89      integer(kind=im), intent(in) :: ncbands         ! number of cloud spectral bands
90      real(kind=rb), intent(in) :: cldfmc(:,:)        ! layer cloud fraction [mcica]
91                                                      !    Dimensions: (ngptlw,nlayers)
92      real(kind=rb), intent(in) :: taucmc(:,:)        ! layer cloud optical depth [mcica]
93                                                      !    Dimensions: (ngptlw,nlayers)
94      integer(kind=im), intent(in) :: idrv            ! flag for calculation of dF/dt from
95                                                      ! Planck derivative [0=off, 1=on]
96      real(kind=rb), intent(in) :: dplankbnd_dt(:)    ! derivative of Planck function wrt temp
97                                                      !    Dimensions: (nbndlw)
98
99! ----- Output -----
100      real(kind=rb), intent(out) :: totuflux(0:)      ! upward longwave flux (w/m2)
101                                                      !    Dimensions: (0:nlayers)
102      real(kind=rb), intent(out) :: totdflux(0:)      ! downward longwave flux (w/m2)
103                                                      !    Dimensions: (0:nlayers)
104      real(kind=rb), intent(out) :: fnet(0:)          ! net longwave flux (w/m2)
105                                                      !    Dimensions: (0:nlayers)
106      real(kind=rb), intent(out) :: htr(0:)           ! longwave heating rate (k/day)
107                                                      !    Dimensions: (0:nlayers)
108      real(kind=rb), intent(out) :: totuclfl(0:)      ! clear sky upward longwave flux (w/m2)
109                                                      !    Dimensions: (0:nlayers)
110      real(kind=rb), intent(out) :: totdclfl(0:)      ! clear sky downward longwave flux (w/m2)
111                                                      !    Dimensions: (0:nlayers)
112      real(kind=rb), intent(out) :: fnetc(0:)         ! clear sky net longwave flux (w/m2)
113                                                      !    Dimensions: (0:nlayers)
114      real(kind=rb), intent(out) :: htrc(0:)          ! clear sky longwave heating rate (k/day)
115                                                      !    Dimensions: (0:nlayers)
116      real(kind=rb), intent(out) :: dtotuflux_dt(0:)  ! change in upward longwave flux (w/m2/k)
117                                                      ! with respect to surface temperature
118                                                      !    Dimensions: (0:nlayers)
119      real(kind=rb), intent(out) :: dtotuclfl_dt(0:)  ! change in upward longwave flux (w/m2/k)
120                                                      ! with respect to surface temperature
121                                                      !    Dimensions: (0:nlayers)
122
123! ----- Local -----
124! Declarations for radiative transfer
125      real(kind=rb) :: abscld(nlayers,ngptlw)
126      real(kind=rb) :: atot(nlayers)
127      real(kind=rb) :: atrans(nlayers)
128      real(kind=rb) :: bbugas(nlayers)
129      real(kind=rb) :: bbutot(nlayers)
130      real(kind=rb) :: clrurad(0:nlayers)
131      real(kind=rb) :: clrdrad(0:nlayers)
132      real(kind=rb) :: efclfrac(nlayers,ngptlw)
133      real(kind=rb) :: uflux(0:nlayers)
134      real(kind=rb) :: dflux(0:nlayers)
135      real(kind=rb) :: urad(0:nlayers)
136      real(kind=rb) :: drad(0:nlayers)
137      real(kind=rb) :: uclfl(0:nlayers)
138      real(kind=rb) :: dclfl(0:nlayers)
139      real(kind=rb) :: odcld(nlayers,ngptlw)
140
141
142      real(kind=rb) :: secdiff(nbndlw)                 ! secant of diffusivity angle
143      real(kind=rb) :: a0(nbndlw),a1(nbndlw),a2(nbndlw)! diffusivity angle adjustment coefficients
144      real(kind=rb) :: wtdiff, rec_6
145      real(kind=rb) :: transcld, radld, radclrd, plfrac, blay, dplankup, dplankdn
146      real(kind=rb) :: odepth, odtot, odepth_rec, odtot_rec, gassrc
147      real(kind=rb) :: tblind, tfactot, bbd, bbdtot, tfacgas, transc, tausfac
148      real(kind=rb) :: rad0, reflect, radlu, radclru
149
150      real(kind=rb) :: duflux_dt(0:nlayers)
151      real(kind=rb) :: duclfl_dt(0:nlayers)
152      real(kind=rb) :: d_urad_dt(0:nlayers)
153      real(kind=rb) :: d_clrurad_dt(0:nlayers)
154      real(kind=rb) :: d_rad0_dt, d_radlu_dt, d_radclru_dt
155
156      integer(kind=im) :: icldlyr(nlayers)                  ! flag for cloud in layer
157      integer(kind=im) :: ibnd, ib, iband, lay, lev, l, ig  ! loop indices
158      integer(kind=im) :: igc                               ! g-point interval counter
159      integer(kind=im) :: iclddn                            ! flag for cloud in down path
160      integer(kind=im) :: ittot, itgas, itr                 ! lookup table indices
161
162! ------- Definitions -------
163! input
164!    nlayers                      ! number of model layers
165!    ngptlw                       ! total number of g-point subintervals
166!    nbndlw                       ! number of longwave spectral bands
167!    ncbands                      ! number of spectral bands for clouds
168!    secdiff                      ! diffusivity angle
169!    wtdiff                       ! weight for radiance to flux conversion
170!    pavel                        ! layer pressures (mb)
171!    pz                           ! level (interface) pressures (mb)
172!    tavel                        ! layer temperatures (k)
173!    tz                           ! level (interface) temperatures(mb)
174!    tbound                       ! surface temperature (k)
175!    cldfrac                      ! layer cloud fraction
176!    taucloud                     ! layer cloud optical depth
177!    itr                          ! integer look-up table index
178!    icldlyr                      ! flag for cloudy layers
179!    iclddn                       ! flag for cloud in column at any layer
180!    semiss                       ! surface emissivities for each band
181!    reflect                      ! surface reflectance
182!    bpade                        ! 1/(pade constant)
183!    tau_tbl                      ! clear sky optical depth look-up table
184!    exp_tbl                      ! exponential look-up table for transmittance
185!    tfn_tbl                      ! tau transition function look-up table
186
187! local
188!    atrans                       ! gaseous absorptivity
189!    abscld                       ! cloud absorptivity
190!    atot                         ! combined gaseous and cloud absorptivity
191!    odclr                        ! clear sky (gaseous) optical depth
192!    odcld                        ! cloud optical depth
193!    odtot                        ! optical depth of gas and cloud
194!    tfacgas                      ! gas-only pade factor, used for planck fn
195!    tfactot                      ! gas and cloud pade factor, used for planck fn
196!    bbdgas                       ! gas-only planck function for downward rt
197!    bbugas                       ! gas-only planck function for upward rt
198!    bbdtot                       ! gas and cloud planck function for downward rt
199!    bbutot                       ! gas and cloud planck function for upward calc.
200!    gassrc                       ! source radiance due to gas only
201!    efclfrac                     ! effective cloud fraction
202!    radlu                        ! spectrally summed upward radiance
203!    radclru                      ! spectrally summed clear sky upward radiance
204!    urad                         ! upward radiance by layer
205!    clrurad                      ! clear sky upward radiance by layer
206!    radld                        ! spectrally summed downward radiance
207!    radclrd                      ! spectrally summed clear sky downward radiance
208!    drad                         ! downward radiance by layer
209!    clrdrad                      ! clear sky downward radiance by layer
210!    d_radlu_dt                   ! spectrally summed upward radiance
211!    d_radclru_dt                 ! spectrally summed clear sky upward radiance
212!    d_urad_dt                    ! upward radiance by layer
213!    d_clrurad_dt                 ! clear sky upward radiance by layer
214
215! output
216!    totuflux                     ! upward longwave flux (w/m2)
217!    totdflux                     ! downward longwave flux (w/m2)
218!    fnet                         ! net longwave flux (w/m2)
219!    htr                          ! longwave heating rate (k/day)
220!    totuclfl                     ! clear sky upward longwave flux (w/m2)
221!    totdclfl                     ! clear sky downward longwave flux (w/m2)
222!    fnetc                        ! clear sky net longwave flux (w/m2)
223!    htrc                         ! clear sky longwave heating rate (k/day)
224!    dtotuflux_dt                 ! change in upward longwave flux (w/m2/k)
225!                                 ! with respect to surface temperature
226!    dtotuclfl_dt                 ! change in clear sky upward longwave flux (w/m2/k)
227!                                 ! with respect to surface temperature
228
229
230! This secant and weight corresponds to the standard diffusivity
231! angle.  This initial value is redefined below for some bands.
232      data wtdiff /0.5_rb/
233      data rec_6 /0.166667_rb/
234
235! Reset diffusivity angle for Bands 2-3 and 5-9 to vary (between 1.50
236! and 1.80) as a function of total column water vapor.  The function
237! has been defined to minimize flux and cooling rate errors in these bands
238! over a wide range of precipitable water values.
239      data a0 / 1.66_rb,  1.55_rb,  1.58_rb,  1.66_rb, &
240                1.54_rb, 1.454_rb,  1.89_rb,  1.33_rb, &
241               1.668_rb,  1.66_rb,  1.66_rb,  1.66_rb, &
242                1.66_rb,  1.66_rb,  1.66_rb,  1.66_rb /
243      data a1 / 0.00_rb,  0.25_rb,  0.22_rb,  0.00_rb, &
244                0.13_rb, 0.446_rb, -0.10_rb,  0.40_rb, &
245              -0.006_rb,  0.00_rb,  0.00_rb,  0.00_rb, &
246                0.00_rb,  0.00_rb,  0.00_rb,  0.00_rb /
247      data a2 / 0.00_rb, -12.0_rb, -11.7_rb,  0.00_rb, &
248               -0.72_rb,-0.243_rb,  0.19_rb,-0.062_rb, &
249               0.414_rb,  0.00_rb,  0.00_rb,  0.00_rb, &
250                0.00_rb,  0.00_rb,  0.00_rb,  0.00_rb /
251
252      hvrrtc = '$Revision: 1.7 $'
253
254      do ibnd = 1,nbndlw
255         if (ibnd.eq.1 .or. ibnd.eq.4 .or. ibnd.ge.10) then
256           secdiff(ibnd) = 1.66_rb
257         else
258           secdiff(ibnd) = a0(ibnd) + a1(ibnd)*exp(a2(ibnd)*pwvcm)
259           if (secdiff(ibnd) .gt. 1.80_rb) secdiff(ibnd) = 1.80_rb
260           if (secdiff(ibnd) .lt. 1.50_rb) secdiff(ibnd) = 1.50_rb
261         endif
262      enddo
263
264      urad(0) = 0.0_rb
265      drad(0) = 0.0_rb
266      totuflux(0) = 0.0_rb
267      totdflux(0) = 0.0_rb
268      clrurad(0) = 0.0_rb
269      clrdrad(0) = 0.0_rb
270      totuclfl(0) = 0.0_rb
271      totdclfl(0) = 0.0_rb
272      if (idrv .eq. 1) then
273         d_urad_dt(0) = 0.0_rb
274         d_clrurad_dt(0) = 0.0_rb
275         dtotuflux_dt(0) = 0.0_rb
276         dtotuclfl_dt(0) = 0.0_rb
277      endif
278
279      do lay = 1, nlayers
280         urad(lay) = 0.0_rb
281         drad(lay) = 0.0_rb
282         totuflux(lay) = 0.0_rb
283         totdflux(lay) = 0.0_rb
284         clrurad(lay) = 0.0_rb
285         clrdrad(lay) = 0.0_rb
286         totuclfl(lay) = 0.0_rb
287         totdclfl(lay) = 0.0_rb
288         icldlyr(lay) = 0
289         if (idrv .eq. 1) then
290            d_urad_dt(lay) = 0.0_rb
291            d_clrurad_dt(lay) = 0.0_rb
292            dtotuflux_dt(lay) = 0.0_rb
293            dtotuclfl_dt(lay) = 0.0_rb
294         endif
295
296! Change to band loop?
297         do ig = 1, ngptlw
298            if (cldfmc(ig,lay) .eq. 1._rb) then
299               ib = ngb(ig)
300               odcld(lay,ig) = secdiff(ib) * taucmc(ig,lay)
301               transcld = exp(-odcld(lay,ig))
302               abscld(lay,ig) = 1._rb - transcld
303               efclfrac(lay,ig) = abscld(lay,ig) * cldfmc(ig,lay)
304               icldlyr(lay) = 1
305            else
306               odcld(lay,ig) = 0.0_rb
307               abscld(lay,ig) = 0.0_rb
308               efclfrac(lay,ig) = 0.0_rb
309            endif
310         enddo
311
312      enddo
313
314      igc = 1
315! Loop over frequency bands.
316      do iband = istart, iend
317
318! Reinitialize g-point counter for each band if output for each band is requested.
319         if (iout.gt.0.and.iband.ge.2) igc = ngs(iband-1)+1
320
321! Loop over g-channels.
322 1000    continue
323
324! Radiative transfer starts here.
325         radld = 0._rb
326         radclrd = 0._rb
327         iclddn = 0
328
329! Downward radiative transfer loop. 
330
331         do lev = nlayers, 1, -1
332               plfrac = fracs(lev,igc)
333               blay = planklay(lev,iband)
334               dplankup = planklev(lev,iband) - blay
335               dplankdn = planklev(lev-1,iband) - blay
336               odepth = secdiff(iband) * taut(lev,igc)
337               if (odepth .lt. 0.0_rb) odepth = 0.0_rb
338!  Cloudy layer
339               if (icldlyr(lev).eq.1) then
340                  iclddn = 1
341                  odtot = odepth + odcld(lev,igc)
342                  if (odtot .lt. 0.06_rb) then
343                     atrans(lev) = odepth - 0.5_rb*odepth*odepth
344                     odepth_rec = rec_6*odepth
345                     gassrc = plfrac*(blay+dplankdn*odepth_rec)*atrans(lev)
346
347                     atot(lev) =  odtot - 0.5_rb*odtot*odtot
348                     odtot_rec = rec_6*odtot
349                     bbdtot =  plfrac * (blay+dplankdn*odtot_rec)
350                     bbd = plfrac*(blay+dplankdn*odepth_rec)
351                     radld = radld - radld * (atrans(lev) + &
352                         efclfrac(lev,igc) * (1. - atrans(lev))) + &
353                         gassrc + cldfmc(igc,lev) * &
354                         (bbdtot * atot(lev) - gassrc)
355                     drad(lev-1) = drad(lev-1) + radld
356                 
357                     bbugas(lev) =  plfrac * (blay+dplankup*odepth_rec)
358                     bbutot(lev) =  plfrac * (blay+dplankup*odtot_rec)
359
360                  elseif (odepth .le. 0.06_rb) then
361                     atrans(lev) = odepth - 0.5_rb*odepth*odepth
362                     odepth_rec = rec_6*odepth
363                     gassrc = plfrac*(blay+dplankdn*odepth_rec)*atrans(lev)
364
365                     odtot = odepth + odcld(lev,igc)
366                     tblind = odtot/(bpade+odtot)
367                     ittot = tblint*tblind + 0.5_rb
368                     tfactot = tfn_tbl(ittot)
369                     bbdtot = plfrac * (blay + tfactot*dplankdn)
370                     bbd = plfrac*(blay+dplankdn*odepth_rec)
371                     atot(lev) = 1. - exp_tbl(ittot)
372
373                     radld = radld - radld * (atrans(lev) + &
374                         efclfrac(lev,igc) * (1._rb - atrans(lev))) + &
375                         gassrc + cldfmc(igc,lev) * &
376                         (bbdtot * atot(lev) - gassrc)
377                     drad(lev-1) = drad(lev-1) + radld
378
379                     bbugas(lev) = plfrac * (blay + dplankup*odepth_rec)
380                     bbutot(lev) = plfrac * (blay + tfactot * dplankup)
381
382                  else
383
384                     tblind = odepth/(bpade+odepth)
385                     itgas = tblint*tblind+0.5_rb
386                     odepth = tau_tbl(itgas)
387                     atrans(lev) = 1._rb - exp_tbl(itgas)
388                     tfacgas = tfn_tbl(itgas)
389                     gassrc = atrans(lev) * plfrac * (blay + tfacgas*dplankdn)
390
391                     odtot = odepth + odcld(lev,igc)
392                     tblind = odtot/(bpade+odtot)
393                     ittot = tblint*tblind + 0.5_rb
394                     tfactot = tfn_tbl(ittot)
395                     bbdtot = plfrac * (blay + tfactot*dplankdn)
396                     bbd = plfrac*(blay+tfacgas*dplankdn)
397                     atot(lev) = 1._rb - exp_tbl(ittot)
398
399                  radld = radld - radld * (atrans(lev) + &
400                    efclfrac(lev,igc) * (1._rb - atrans(lev))) + &
401                    gassrc + cldfmc(igc,lev) * &
402                    (bbdtot * atot(lev) - gassrc)
403                  drad(lev-1) = drad(lev-1) + radld
404                  bbugas(lev) = plfrac * (blay + tfacgas * dplankup)
405                  bbutot(lev) = plfrac * (blay + tfactot * dplankup)
406                  endif
407!  Clear layer
408               else
409                  if (odepth .le. 0.06_rb) then
410                     atrans(lev) = odepth-0.5_rb*odepth*odepth
411                     odepth = rec_6*odepth
412                     bbd = plfrac*(blay+dplankdn*odepth)
413                     bbugas(lev) = plfrac*(blay+dplankup*odepth)
414                  else
415                     tblind = odepth/(bpade+odepth)
416                     itr = tblint*tblind+0.5_rb
417                     transc = exp_tbl(itr)
418                     atrans(lev) = 1._rb-transc
419                     tausfac = tfn_tbl(itr)
420                     bbd = plfrac*(blay+tausfac*dplankdn)
421                     bbugas(lev) = plfrac * (blay + tausfac * dplankup)
422                  endif   
423                  radld = radld + (bbd-radld)*atrans(lev)
424                  drad(lev-1) = drad(lev-1) + radld
425               endif
426!  Set clear sky stream to total sky stream as long as layers
427!  remain clear.  Streams diverge when a cloud is reached (iclddn=1),
428!  and clear sky stream must be computed separately from that point.
429                  if (iclddn.eq.1) then
430                     radclrd = radclrd + (bbd-radclrd) * atrans(lev) 
431                     clrdrad(lev-1) = clrdrad(lev-1) + radclrd
432                  else
433                     radclrd = radld
434                     clrdrad(lev-1) = drad(lev-1)
435                  endif
436            enddo
437
438! Spectral emissivity & reflectance
439!  Include the contribution of spectrally varying longwave emissivity
440!  and reflection from the surface to the upward radiative transfer.
441!  Note: Spectral and Lambertian reflection are identical for the
442!  diffusivity angle flux integration used here.
443!  Note: The emissivity is applied to plankbnd and dplankbnd_dt when
444!  they are defined in subroutine setcoef.
445
446         rad0 = fracs(1,igc) * plankbnd(iband)
447         if (idrv .eq. 1) then
448            d_rad0_dt = fracs(1,igc) * dplankbnd_dt(iband)
449         endif
450
451!  Add in specular reflection of surface downward radiance.
452         reflect = 1._rb - semiss(iband)
453         radlu = rad0 + reflect * radld
454         radclru = rad0 + reflect * radclrd
455
456
457! Upward radiative transfer loop.
458         urad(0) = urad(0) + radlu
459         clrurad(0) = clrurad(0) + radclru
460         if (idrv .eq. 1) then
461            d_radlu_dt = d_rad0_dt
462            d_urad_dt(0) = d_urad_dt(0) + d_radlu_dt
463            d_radclru_dt = d_rad0_dt
464            d_clrurad_dt(0) = d_clrurad_dt(0) + d_radclru_dt
465         endif
466
467         do lev = 1, nlayers
468!  Cloudy layer
469            if (icldlyr(lev) .eq. 1) then
470               gassrc = bbugas(lev) * atrans(lev)
471               radlu = radlu - radlu * (atrans(lev) + &
472                   efclfrac(lev,igc) * (1._rb - atrans(lev))) + &
473                   gassrc + cldfmc(igc,lev) * &
474                   (bbutot(lev) * atot(lev) - gassrc)
475               urad(lev) = urad(lev) + radlu
476               if (idrv .eq. 1) then
477                  d_radlu_dt = d_radlu_dt * cldfmc(igc,lev) * (1.0_rb - atot(lev)) + &
478                         d_radlu_dt * (1.0_rb - cldfmc(igc,lev)) * (1.0_rb - atrans(lev))
479                  d_urad_dt(lev) = d_urad_dt(lev) + d_radlu_dt
480               endif
481!  Clear layer
482            else
483               radlu = radlu + (bbugas(lev)-radlu)*atrans(lev)
484               urad(lev) = urad(lev) + radlu
485               if (idrv .eq. 1) then
486                  d_radlu_dt = d_radlu_dt * (1.0_rb - atrans(lev))
487                  d_urad_dt(lev) = d_urad_dt(lev) + d_radlu_dt
488               endif
489            endif
490!  Set clear sky stream to total sky stream as long as all layers
491!  are clear (iclddn=0).  Streams must be calculated separately at
492!  all layers when a cloud is present (ICLDDN=1), because surface
493!  reflectance is different for each stream.
494               if (iclddn.eq.1) then
495                  radclru = radclru + (bbugas(lev)-radclru)*atrans(lev) 
496                  clrurad(lev) = clrurad(lev) + radclru
497               else
498                  radclru = radlu
499                  clrurad(lev) = urad(lev)
500               endif
501               if (idrv .eq. 1) then
502                  if (iclddn.eq.1) then
503                     d_radclru_dt = d_radclru_dt * (1.0_rb - atrans(lev))
504                     d_clrurad_dt(lev) = d_clrurad_dt(lev) + d_radclru_dt
505                  else
506                     d_radclru_dt = d_radlu_dt
507                     d_clrurad_dt(lev) = d_urad_dt(lev)
508                  endif
509               endif
510         enddo
511
512! Increment g-point counter
513         igc = igc + 1
514! Return to continue radiative transfer for all g-channels in present band
515         if (igc .le. ngs(iband)) go to 1000
516
517! Process longwave output from band for total and clear streams.
518! Calculate upward, downward, and net flux.
519         do lev = nlayers, 0, -1
520            uflux(lev) = urad(lev)*wtdiff
521            dflux(lev) = drad(lev)*wtdiff
522            urad(lev) = 0.0_rb
523            drad(lev) = 0.0_rb
524            totuflux(lev) = totuflux(lev) + uflux(lev) * delwave(iband)
525            totdflux(lev) = totdflux(lev) + dflux(lev) * delwave(iband)
526            uclfl(lev) = clrurad(lev)*wtdiff
527            dclfl(lev) = clrdrad(lev)*wtdiff
528            clrurad(lev) = 0.0_rb
529            clrdrad(lev) = 0.0_rb
530            totuclfl(lev) = totuclfl(lev) + uclfl(lev) * delwave(iband)
531            totdclfl(lev) = totdclfl(lev) + dclfl(lev) * delwave(iband)
532         enddo
533
534! Calculate total change in upward flux wrt surface temperature
535         if (idrv .eq. 1) then
536            do lev = nlayers, 0, -1
537               duflux_dt(lev) = d_urad_dt(lev) * wtdiff
538               d_urad_dt(lev) = 0.0_rb
539               dtotuflux_dt(lev) = dtotuflux_dt(lev) + duflux_dt(lev) * delwave(iband) * fluxfac
540               duclfl_dt(lev) = d_clrurad_dt(lev) * wtdiff
541               d_clrurad_dt(lev) = 0.0_rb
542               dtotuclfl_dt(lev) = dtotuclfl_dt(lev) + duclfl_dt(lev) * delwave(iband) * fluxfac
543            enddo
544         endif
545
546! End spectral band loop
547      enddo
548
549! Calculate fluxes at surface
550      totuflux(0) = totuflux(0) * fluxfac
551      totdflux(0) = totdflux(0) * fluxfac
552      fnet(0) = totuflux(0) - totdflux(0)
553      totuclfl(0) = totuclfl(0) * fluxfac
554      totdclfl(0) = totdclfl(0) * fluxfac
555      fnetc(0) = totuclfl(0) - totdclfl(0)
556
557! Calculate fluxes at model levels
558      do lev = 1, nlayers
559         totuflux(lev) = totuflux(lev) * fluxfac
560         totdflux(lev) = totdflux(lev) * fluxfac
561         fnet(lev) = totuflux(lev) - totdflux(lev)
562         totuclfl(lev) = totuclfl(lev) * fluxfac
563         totdclfl(lev) = totdclfl(lev) * fluxfac
564         fnetc(lev) = totuclfl(lev) - totdclfl(lev)
565         l = lev - 1
566
567! Calculate heating rates at model layers
568         htr(l)=heatfac*(fnet(l)-fnet(lev))/(pz(l)-pz(lev)) 
569         htrc(l)=heatfac*(fnetc(l)-fnetc(lev))/(pz(l)-pz(lev)) 
570      enddo
571
572! Set heating rate to zero in top layer
573      htr(nlayers) = 0.0_rb
574      htrc(nlayers) = 0.0_rb
575
576      end subroutine rtrnmc
577
578      end module rrtmg_lw_rtrnmc
579
Note: See TracBrowser for help on using the repository browser.