source: palm/trunk/LIB/rrtmg/rrtmg_lw_rtrn.f90 @ 3052

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

Added support for RRTMG radiation code

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