source: palm/trunk/LIB/rrtmg/rrtmg_lw_rtrnmr.f90 @ 4476

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

Added support for RRTMG radiation code

File size: 34.8 KB
Line 
1!     path:      $Source: /storm/rc1/cvsroot/rc/rrtmg_lw/src/rrtmg_lw_rtrnmr.f90,v $
2!     author:    $Author: mike $
3!     revision:  $Revision: 1.7 $
4!     created:   $Date: 2009/11/12 20:52:26 $
5!
6      module rrtmg_lw_rtrnmr
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: hvrrtx, hnamrtx
26
27      implicit none
28
29      contains
30
31!-----------------------------------------------------------------------------
32      subroutine rtrnmr(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 maximum-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      real(kind=rb) :: secdiff(nbndlw)                 ! secant of diffusivity angle
141      real(kind=rb) :: a0(nbndlw),a1(nbndlw),a2(nbndlw)! diffusivity angle adjustment coefficients
142      real(kind=rb) :: wtdiff, rec_6
143      real(kind=rb) :: transcld, radld, radclrd, plfrac, blay, dplankup, dplankdn
144      real(kind=rb) :: odepth, odtot, odepth_rec, odtot_rec, gassrc, ttot
145      real(kind=rb) :: tblind, tfactot, bbd, bbdtot, tfacgas, transc, tausfac
146      real(kind=rb) :: rad0, reflect, radlu, radclru
147
148      real(kind=rb) :: duflux_dt(0:nlayers)
149      real(kind=rb) :: duclfl_dt(0:nlayers)
150      real(kind=rb) :: d_urad_dt(0:nlayers)
151      real(kind=rb) :: d_clrurad_dt(0:nlayers)
152      real(kind=rb) :: d_rad0_dt, d_radlu_dt, d_radclru_dt
153
154      integer(kind=im) :: icldlyr(nlayers)             ! flag for cloud in layer
155      integer(kind=im) :: ibnd, ib, iband, lay, lev, l ! loop indices
156      integer(kind=im) :: igc                          ! g-point interval counter
157      integer(kind=im) :: iclddn                       ! flag for cloud in down path
158      integer(kind=im) :: ittot, itgas, itr            ! lookup table indices
159      integer(kind=im) :: ipat(16,0:2)
160
161! Declarations for cloud overlap adjustment
162      real(kind=rb) :: faccld1(nlayers+1),faccld2(nlayers+1)
163      real(kind=rb) :: facclr1(nlayers+1),facclr2(nlayers+1)
164      real(kind=rb) :: faccmb1(nlayers+1),faccmb2(nlayers+1)
165      real(kind=rb) :: faccld1d(0:nlayers),faccld2d(0:nlayers)
166      real(kind=rb) :: facclr1d(0:nlayers),facclr2d(0:nlayers)
167      real(kind=rb) :: faccmb1d(0:nlayers),faccmb2d(0:nlayers)
168
169      real(kind=rb) :: fmax, fmin, rat1, rat2
170      real(kind=rb) :: clrradd, cldradd, clrradu, cldradu, oldclr, oldcld
171      real(kind=rb) :: rad, cldsrc, radmod
172
173      integer(kind=im) :: istcld(nlayers+1),istcldd(0:nlayers)
174
175! ------- Definitions -------
176! input
177!    nlayers                      ! number of model layers
178!    ngptlw                       ! total number of g-point subintervals
179!    nbndlw                       ! number of longwave spectral bands
180!    ncbands                      ! number of spectral bands for clouds
181!    secdiff                      ! diffusivity angle
182!    wtdiff                       ! weight for radiance to flux conversion
183!    pavel                        ! layer pressures (mb)
184!    pz                           ! level (interface) pressures (mb)
185!    tavel                        ! layer temperatures (k)
186!    tz                           ! level (interface) temperatures(mb)
187!    tbound                       ! surface temperature (k)
188!    cldfrac                      ! layer cloud fraction
189!    taucloud                     ! layer cloud optical depth
190!    itr                          ! integer look-up table index
191!    icldlyr                      ! flag for cloudy layers
192!    iclddn                       ! flag for cloud in column at any layer
193!    semiss                       ! surface emissivities for each band
194!    reflect                      ! surface reflectance
195!    bpade                        ! 1/(pade constant)
196!    tau_tbl                      ! clear sky optical depth look-up table
197!    exp_tbl                      ! exponential look-up table for transmittance
198!    tfn_tbl                      ! tau transition function look-up table
199
200! local
201!    atrans                       ! gaseous absorptivity
202!    abscld                       ! cloud absorptivity
203!    atot                         ! combined gaseous and cloud absorptivity
204!    odclr                        ! clear sky (gaseous) optical depth
205!    odcld                        ! cloud optical depth
206!    odtot                        ! optical depth of gas and cloud
207!    tfacgas                      ! gas-only pade factor, used for planck fn
208!    tfactot                      ! gas and cloud pade factor, used for planck fn
209!    bbdgas                       ! gas-only planck function for downward rt
210!    bbugas                       ! gas-only planck function for upward rt
211!    bbdtot                       ! gas and cloud planck function for downward rt
212!    bbutot                       ! gas and cloud planck function for upward calc.
213!    gassrc                       ! source radiance due to gas only
214!    efclfrac                     ! effective cloud fraction
215!    radlu                        ! spectrally summed upward radiance
216!    radclru                      ! spectrally summed clear sky upward radiance
217!    urad                         ! upward radiance by layer
218!    clrurad                      ! clear sky upward radiance by layer
219!    radld                        ! spectrally summed downward radiance
220!    radclrd                      ! spectrally summed clear sky downward radiance
221!    drad                         ! downward radiance by layer
222!    clrdrad                      ! clear sky downward radiance by layer
223!    d_radlu_dt                   ! spectrally summed upward radiance
224!    d_radclru_dt                 ! spectrally summed clear sky upward radiance
225!    d_urad_dt                    ! upward radiance by layer
226!    d_clrurad_dt                 ! clear sky upward radiance by layer
227
228! output
229!    totuflux                     ! upward longwave flux (w/m2)
230!    totdflux                     ! downward longwave flux (w/m2)
231!    fnet                         ! net longwave flux (w/m2)
232!    htr                          ! longwave heating rate (k/day)
233!    totuclfl                     ! clear sky upward longwave flux (w/m2)
234!    totdclfl                     ! clear sky downward longwave flux (w/m2)
235!    fnetc                        ! clear sky net longwave flux (w/m2)
236!    htrc                         ! clear sky longwave heating rate (k/day)
237!    dtotuflux_dt                 ! change in upward longwave flux (w/m2/k)
238!                                 ! with respect to surface temperature
239!    dtotuclfl_dt                 ! change in clear sky upward longwave flux (w/m2/k)
240!                                 ! with respect to surface temperature
241
242
243! These arrays indicate the spectral 'region' (used in the
244! calculation of ice cloud optical depths) corresponding
245! to each spectral band.  See cldprop.f for more details.
246      data ipat /1,1,1,1,1,1,1,1,1, 1, 1, 1, 1, 1, 1, 1, &
247                 1,2,3,3,3,4,4,4,5, 5, 5, 5, 5, 5, 5, 5, &
248                 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16/
249
250! This secant and weight corresponds to the standard diffusivity
251! angle.  This initial value is redefined below for some bands.
252      data wtdiff /0.5_rb/
253      data rec_6 /0.166667_rb/
254
255! Reset diffusivity angle for Bands 2-3 and 5-9 to vary (between 1.50
256! and 1.80) as a function of total column water vapor.  The function
257! has been defined to minimize flux and cooling rate errors in these bands
258! over a wide range of precipitable water values.
259      data a0 / 1.66_rb,  1.55_rb,  1.58_rb,  1.66_rb, &
260                1.54_rb, 1.454_rb,  1.89_rb,  1.33_rb, &
261               1.668_rb,  1.66_rb,  1.66_rb,  1.66_rb, &
262                1.66_rb,  1.66_rb,  1.66_rb,  1.66_rb /
263      data a1 / 0.00_rb,  0.25_rb,  0.22_rb,  0.00_rb, &
264                0.13_rb, 0.446_rb, -0.10_rb,  0.40_rb, &
265              -0.006_rb,  0.00_rb,  0.00_rb,  0.00_rb, &
266                0.00_rb,  0.00_rb,  0.00_rb,  0.00_rb /
267      data a2 / 0.00_rb, -12.0_rb, -11.7_rb,  0.00_rb, &
268               -0.72_rb,-0.243_rb,  0.19_rb,-0.062_rb, &
269               0.414_rb,  0.00_rb,  0.00_rb,  0.00_rb, &
270                0.00_rb,  0.00_rb,  0.00_rb,  0.00_rb /
271
272      do ibnd = 1,nbndlw
273         if (ibnd.eq.1 .or. ibnd.eq.4 .or. ibnd.ge.10) then
274           secdiff(ibnd) = 1.66_rb
275         else
276           secdiff(ibnd) = a0(ibnd) + a1(ibnd)*exp(a2(ibnd)*pwvcm)
277           if (secdiff(ibnd) .gt. 1.80_rb) secdiff(ibnd) = 1.80_rb
278           if (secdiff(ibnd) .lt. 1.50_rb) secdiff(ibnd) = 1.50_rb
279         endif
280      enddo
281
282      hvrrtx = '$Revision: 1.7 $'
283
284      urad(0) = 0.0_rb
285      drad(0) = 0.0_rb
286      totuflux(0) = 0.0_rb
287      totdflux(0) = 0.0_rb
288      clrurad(0) = 0.0_rb
289      clrdrad(0) = 0.0_rb
290      totuclfl(0) = 0.0_rb
291      totdclfl(0) = 0.0_rb
292      if (idrv .eq. 1) then
293         d_urad_dt(0) = 0.0_rb
294         d_clrurad_dt(0) = 0.0_rb
295         dtotuflux_dt(0) = 0.0_rb
296         dtotuclfl_dt(0) = 0.0_rb
297      endif
298
299      do lay = 1, nlayers
300         urad(lay) = 0.0_rb
301         drad(lay) = 0.0_rb
302         totuflux(lay) = 0.0_rb
303         totdflux(lay) = 0.0_rb
304         clrurad(lay) = 0.0_rb
305         clrdrad(lay) = 0.0_rb
306         totuclfl(lay) = 0.0_rb
307         totdclfl(lay) = 0.0_rb
308         if (idrv .eq. 1) then
309            d_urad_dt(lay) = 0.0_rb
310            d_clrurad_dt(lay) = 0.0_rb
311            dtotuflux_dt(lay) = 0.0_rb
312            dtotuclfl_dt(lay) = 0.0_rb
313         endif
314
315         do ib = 1, ncbands
316            if (cldfrac(lay) .ge. 1.e-6_rb) then
317               odcld(lay,ib) = secdiff(ib) * taucloud(lay,ib)
318               icldlyr(lay) = 1
319            else
320               odcld(lay,ib) = 0.0_rb
321               icldlyr(lay) = 0
322            endif
323         enddo
324      enddo
325
326! Maximum/Random cloud overlap parameter
327
328      istcld(1) = 1
329      istcldd(nlayers) = 1
330      do lev = 1, nlayers
331
332         if (icldlyr(lev).eq.1) then
333! Maximum/random cloud overlap
334            istcld(lev+1) = 0
335            if (lev .eq. nlayers) then
336               faccld1(lev+1) = 0._rb
337               faccld2(lev+1) = 0._rb
338               facclr1(lev+1) = 0._rb
339               facclr2(lev+1) = 0._rb
340               faccmb1(lev+1) = 0._rb
341               faccmb2(lev+1) = 0._rb
342            elseif (cldfrac(lev+1) .ge. cldfrac(lev)) then
343               faccld1(lev+1) = 0._rb
344               faccld2(lev+1) = 0._rb
345               if (istcld(lev) .eq. 1) then
346                  facclr1(lev+1) = 0._rb
347                  facclr2(lev+1) = 0._rb
348                  if (cldfrac(lev) .lt. 1._rb) facclr2(lev+1) = &
349                     (cldfrac(lev+1)-cldfrac(lev))/(1._rb-cldfrac(lev))
350                  facclr2(lev) = 0._rb
351                  faccld2(lev) = 0._rb
352               else
353                  fmax = max(cldfrac(lev),cldfrac(lev-1))
354                  if (cldfrac(lev+1) .gt. fmax) then
355                     facclr1(lev+1) = rat2
356                     facclr2(lev+1) = (cldfrac(lev+1)-fmax)/(1._rb-fmax)
357                  elseif (cldfrac(lev+1) .lt. fmax) then
358                     facclr1(lev+1) = (cldfrac(lev+1)-cldfrac(lev))/ &
359                        (cldfrac(lev-1)-cldfrac(lev))
360                     facclr2(lev+1) = 0._rb
361                  else
362                     facclr1(lev+1) = rat2
363                     facclr2(lev+1) = 0._rb
364                  endif
365               endif
366               if (facclr1(lev+1).gt.0._rb .or. facclr2(lev+1).gt.0._rb) then
367                  rat1 = 1._rb
368                  rat2 = 0._rb
369               else
370                  rat1 = 0._rb
371                  rat2 = 0._rb
372               endif
373            else
374               facclr1(lev+1) = 0._rb
375               facclr2(lev+1) = 0._rb
376               if (istcld(lev) .eq. 1) then
377                  faccld1(lev+1) = 0._rb
378                  faccld2(lev+1) = (cldfrac(lev)-cldfrac(lev+1))/cldfrac(lev)
379
380                  facclr2(lev) = 0._rb
381                  faccld2(lev) = 0._rb
382               else
383                  fmin = min(cldfrac(lev),cldfrac(lev-1))
384                  if (cldfrac(lev+1) .le. fmin) then
385                     faccld1(lev+1) = rat1
386                     faccld2(lev+1) = (fmin-cldfrac(lev+1))/fmin
387                  else
388                     faccld1(lev+1) = (cldfrac(lev)-cldfrac(lev+1))/(cldfrac(lev)-fmin)
389                     faccld2(lev+1) = 0._rb
390                  endif
391               endif
392               if (faccld1(lev+1).gt.0._rb .or. faccld2(lev+1).gt.0._rb) then
393                  rat1 = 0._rb
394                  rat2 = 1._rb
395               else
396                  rat1 = 0._rb
397                  rat2 = 0._rb
398               endif
399            endif
400            faccmb1(lev+1) = facclr1(lev+1) * faccld2(lev) * cldfrac(lev-1) 
401            faccmb2(lev+1) = faccld1(lev+1) * facclr2(lev) * (1._rb - cldfrac(lev-1)) 
402         else
403            istcld(lev+1) = 1
404         endif
405      enddo
406
407      do lev = nlayers, 1, -1
408         if (icldlyr(lev).eq.1) then
409            istcldd(lev-1) = 0
410            if (lev .eq. 1) then
411               faccld1d(lev-1) = 0._rb
412               faccld2d(lev-1) = 0._rb
413               facclr1d(lev-1) = 0._rb
414               facclr2d(lev-1) = 0._rb
415               faccmb1d(lev-1) = 0._rb
416               faccmb2d(lev-1) = 0._rb
417            elseif (cldfrac(lev-1) .ge. cldfrac(lev)) then
418               faccld1d(lev-1) = 0._rb
419               faccld2d(lev-1) = 0._rb
420               if (istcldd(lev) .eq. 1) then
421                  facclr1d(lev-1) = 0._rb
422                  facclr2d(lev-1) = 0._rb
423                  if (cldfrac(lev) .lt. 1._rb) facclr2d(lev-1) = &
424                     (cldfrac(lev-1)-cldfrac(lev))/(1._rb-cldfrac(lev))
425                  facclr2d(lev) = 0._rb
426                  faccld2d(lev) = 0._rb
427               else
428                  fmax = max(cldfrac(lev),cldfrac(lev+1))
429                  if (cldfrac(lev-1) .gt. fmax) then
430                     facclr1d(lev-1) = rat2
431                     facclr2d(lev-1) = (cldfrac(lev-1)-fmax)/(1._rb-fmax)
432                  elseif (cldfrac(lev-1) .lt. fmax) then
433                     facclr1d(lev-1) = (cldfrac(lev-1)-cldfrac(lev))/ &
434                        (cldfrac(lev+1)-cldfrac(lev))
435                     facclr2d(lev-1) = 0.
436                  else
437                     facclr1d(lev-1) = rat2
438                     facclr2d(lev-1) = 0._rb
439                  endif
440               endif
441               if (facclr1d(lev-1).gt.0._rb .or. facclr2d(lev-1).gt.0._rb)then
442                  rat1 = 1._rb
443                  rat2 = 0._rb
444               else
445                  rat1 = 0._rb
446                  rat2 = 0._rb
447               endif
448            else
449               facclr1d(lev-1) = 0._rb
450               facclr2d(lev-1) = 0._rb
451               if (istcldd(lev) .eq. 1) then
452                  faccld1d(lev-1) = 0._rb
453                  faccld2d(lev-1) = (cldfrac(lev)-cldfrac(lev-1))/cldfrac(lev)
454                  facclr2d(lev) = 0._rb
455                  faccld2d(lev) = 0._rb
456               else
457                  fmin = min(cldfrac(lev),cldfrac(lev+1))
458                  if (cldfrac(lev-1) .le. fmin) then
459                     faccld1d(lev-1) = rat1
460                     faccld2d(lev-1) = (fmin-cldfrac(lev-1))/fmin
461                  else
462                     faccld1d(lev-1) = (cldfrac(lev)-cldfrac(lev-1))/(cldfrac(lev)-fmin)
463                     faccld2d(lev-1) = 0._rb
464                  endif
465               endif
466               if (faccld1d(lev-1).gt.0._rb .or. faccld2d(lev-1).gt.0._rb)then
467                  rat1 = 0._rb
468                  rat2 = 1._rb
469               else
470                  rat1 = 0._rb
471                  rat2 = 0._rb
472               endif
473            endif
474            faccmb1d(lev-1) = facclr1d(lev-1) * faccld2d(lev) * cldfrac(lev+1) 
475            faccmb2d(lev-1) = faccld1d(lev-1) * facclr2d(lev) * (1._rb - cldfrac(lev+1))
476         else
477            istcldd(lev-1) = 1
478         endif
479      enddo
480
481      igc = 1
482! Loop over frequency bands.
483      do iband = istart, iend
484
485! Reinitialize g-point counter for each band if output for each band is requested.
486         if (iout.gt.0.and.iband.ge.2) igc = ngs(iband-1)+1
487         if (ncbands .eq. 1) then
488            ib = ipat(iband,0)
489         elseif (ncbands .eq.  5) then
490            ib = ipat(iband,1)
491         elseif (ncbands .eq. 16) then
492            ib = ipat(iband,2)
493         endif
494
495! Loop over g-channels.
496 1000    continue
497
498! Radiative transfer starts here.
499         radld = 0._rb
500         radclrd = 0._rb
501         iclddn = 0
502
503! Downward radiative transfer loop. 
504
505         do lev = nlayers, 1, -1
506               plfrac = fracs(lev,igc)
507               blay = planklay(lev,iband)
508               dplankup = planklev(lev,iband) - blay
509               dplankdn = planklev(lev-1,iband) - blay
510               odepth = secdiff(iband) * taut(lev,igc)
511
512               if (odepth .lt. 0.0_rb) odepth = 0.0_rb
513! Cloudy layer
514               if (icldlyr(lev).eq.1) then
515                  iclddn = 1
516                  odtot = odepth + odcld(lev,ib)
517                  if (odtot .lt. 0.06_rb) then
518                     atrans(lev) = odepth - 0.5_rb*odepth*odepth
519                     odepth_rec = rec_6*odepth
520                     gassrc = plfrac*(blay+dplankdn*odepth_rec)*atrans(lev)
521
522                     atot(lev) =  odtot - 0.5_rb*odtot*odtot
523                     odtot_rec = rec_6*odtot
524                     bbdtot =  plfrac * (blay+dplankdn*odtot_rec)
525                     bbd = plfrac*(blay+dplankdn*odepth_rec)
526                 
527                     bbugas(lev) =  plfrac * (blay+dplankup*odepth_rec)
528                     bbutot(lev) =  plfrac * (blay+dplankup*odtot_rec)
529                  elseif (odepth .le. 0.06_rb) then
530                     atrans(lev) = odepth - 0.5_rb*odepth*odepth
531                     odepth_rec = rec_6*odepth
532                     gassrc = plfrac*(blay+dplankdn*odepth_rec)*atrans(lev)
533
534                     odtot = odepth + odcld(lev,ib)
535                     tblind = odtot/(bpade+odtot)
536                     ittot = tblint*tblind + 0.5_rb
537                     tfactot = tfn_tbl(ittot)
538                     bbdtot = plfrac * (blay + tfactot*dplankdn)
539                     bbd = plfrac*(blay+dplankdn*odepth_rec)
540                     atot(lev) = 1._rb - exp_tbl(ittot)
541
542                     bbugas(lev) = plfrac * (blay + dplankup*odepth_rec)
543                     bbutot(lev) = plfrac * (blay + tfactot * dplankup)
544                  else
545                     tblind = odepth/(bpade+odepth)
546                     itgas = tblint*tblind+0.5_rb
547                     odepth = tau_tbl(itgas)
548                     atrans(lev) = 1._rb - exp_tbl(itgas)
549                     tfacgas = tfn_tbl(itgas)
550                     gassrc = atrans(lev) * plfrac * (blay + tfacgas*dplankdn)
551
552                     odtot = odepth + odcld(lev,ib)
553                     tblind = odtot/(bpade+odtot)
554                     ittot = tblint*tblind + 0.5_rb
555                     tfactot = tfn_tbl(ittot)
556                     bbdtot = plfrac * (blay + tfactot*dplankdn)
557                     bbd = plfrac*(blay+tfacgas*dplankdn)
558                     atot(lev) = 1._rb - exp_tbl(ittot)
559
560                     bbugas(lev) = plfrac * (blay + tfacgas * dplankup)
561                     bbutot(lev) = plfrac * (blay + tfactot * dplankup)
562                  endif
563
564                  if (istcldd(lev) .eq. 1) then
565                     cldradd = cldfrac(lev) * radld
566                     clrradd = radld - cldradd
567                     oldcld = cldradd
568                     oldclr = clrradd
569                     rad = 0._rb
570                  endif
571                  ttot = 1._rb - atot(lev)
572                  cldsrc = bbdtot * atot(lev)
573                  cldradd = cldradd * ttot + cldfrac(lev) * cldsrc
574                  clrradd = clrradd * (1._rb-atrans(lev)) + (1._rb - cldfrac(lev)) * gassrc
575                  radld = cldradd + clrradd
576                  drad(lev-1) = drad(lev-1) + radld
577
578                  radmod = rad * &
579                       (facclr1d(lev-1) * (1.-atrans(lev)) + &
580                       faccld1d(lev-1) *  ttot) - &
581                       faccmb1d(lev-1) * gassrc + &
582                       faccmb2d(lev-1) * cldsrc
583
584                  oldcld = cldradd - radmod
585                  oldclr = clrradd + radmod
586                  rad = -radmod + facclr2d(lev-1)*oldclr - faccld2d(lev-1)*oldcld
587                  cldradd = cldradd + rad
588                  clrradd = clrradd - rad
589! Clear layer
590               else
591                  if (odepth .le. 0.06_rb) then
592                     atrans(lev) = odepth-0.5_rb*odepth*odepth
593                     odepth = rec_6*odepth
594                     bbd = plfrac*(blay+dplankdn*odepth)
595                     bbugas(lev) = plfrac*(blay+dplankup*odepth)
596                  else
597                     tblind = odepth/(bpade+odepth)
598                     itr = tblint*tblind+0.5_rb
599                     transc = exp_tbl(itr)
600                     atrans(lev) = 1._rb-transc
601                     tausfac = tfn_tbl(itr)
602                     bbd = plfrac*(blay+tausfac*dplankdn)
603                     bbugas(lev) = plfrac * (blay + tausfac * dplankup)
604                  endif   
605                  radld = radld + (bbd-radld)*atrans(lev)
606                  drad(lev-1) = drad(lev-1) + radld
607                endif
608!  Set clear sky stream to total sky stream as long as layers
609!  remain clear.  Streams diverge when a cloud is reached (iclddn=1),
610!  and clear sky stream must be computed separately from that point.
611                 if (iclddn.eq.1) then
612                     radclrd = radclrd + (bbd-radclrd) * atrans(lev) 
613                     clrdrad(lev-1) = clrdrad(lev-1) + radclrd
614                  else
615                     radclrd = radld
616                     clrdrad(lev-1) = drad(lev-1)
617                  endif
618            enddo
619
620! Spectral emissivity & reflectance
621!  Include the contribution of spectrally varying longwave emissivity
622!  and reflection from the surface to the upward radiative transfer.
623!  Note: Spectral and Lambertian reflection are identical for the
624!  diffusivity angle flux integration used here.
625!  Note: The emissivity is applied to plankbnd and dplankbnd_dt when
626!  they are defined in subroutine setcoef.
627
628         rad0 = fracs(1,igc) * plankbnd(iband)
629         if (idrv .eq. 1) then
630            d_rad0_dt = fracs(1,igc) * dplankbnd_dt(iband)
631         endif
632
633!  Add in reflection of surface downward radiance.
634         reflect = 1._rb - semiss(iband)
635         radlu = rad0 + reflect * radld
636         radclru = rad0 + reflect * radclrd
637
638! Upward radiative transfer loop.
639
640         urad(0) = urad(0) + radlu
641         clrurad(0) = clrurad(0) + radclru
642         if (idrv .eq. 1) then
643            d_radlu_dt = d_rad0_dt
644            d_urad_dt(0) = d_urad_dt(0) + d_radlu_dt
645            d_radclru_dt = d_rad0_dt
646            d_clrurad_dt(0) = d_clrurad_dt(0) + d_radclru_dt
647         endif
648
649         do lev = 1, nlayers
650! Cloudy layer
651            if (icldlyr(lev) .eq. 1) then
652               gassrc = bbugas(lev) * atrans(lev)
653               if (istcld(lev) .eq. 1) then
654                  cldradu = cldfrac(lev) * radlu
655                  clrradu = radlu - cldradu
656                  oldcld = cldradu
657                  oldclr = clrradu
658                  rad = 0._rb
659               endif
660               ttot = 1._rb - atot(lev)
661               cldsrc = bbutot(lev) * atot(lev)
662               cldradu = cldradu * ttot + cldfrac(lev) * cldsrc
663               clrradu = clrradu * (1.0_rb-atrans(lev)) + (1._rb - cldfrac(lev)) * gassrc
664! Total sky radiance
665               radlu = cldradu + clrradu
666               urad(lev) = urad(lev) + radlu
667               radmod = rad * &
668                   (facclr1(lev+1)*(1.0_rb-atrans(lev))+ &
669                   faccld1(lev+1) *  ttot) - &
670                   faccmb1(lev+1) * gassrc + &
671                   faccmb2(lev+1) * cldsrc
672               oldcld = cldradu - radmod
673               oldclr = clrradu + radmod
674               rad = -radmod + facclr2(lev+1)*oldclr - faccld2(lev+1)*oldcld
675               cldradu = cldradu + rad
676               clrradu = clrradu - rad
677               if (idrv .eq. 1) then
678                  d_radlu_dt = d_radlu_dt * cldfrac(lev) * (1.0_rb - atot(lev)) + &
679                         d_radlu_dt * (1.0_rb - cldfrac(lev)) * (1.0_rb - atrans(lev))
680                  d_urad_dt(lev) = d_urad_dt(lev) + d_radlu_dt
681               endif
682! Clear layer
683            else
684               radlu = radlu + (bbugas(lev)-radlu)*atrans(lev)
685               urad(lev) = urad(lev) + radlu
686               if (idrv .eq. 1) then
687                  d_radlu_dt = d_radlu_dt * (1.0_rb - atrans(lev))
688                  d_urad_dt(lev) = d_urad_dt(lev) + d_radlu_dt
689               endif
690            endif
691!  Set clear sky stream to total sky stream as long as all layers
692!  are clear (iclddn=0).  Streams must be calculated separately at
693!  all layers when a cloud is present (iclddn=1), because surface
694!  reflectance is different for each stream.
695               if (iclddn.eq.1) then
696                  radclru = radclru + (bbugas(lev)-radclru)*atrans(lev) 
697                  clrurad(lev) = clrurad(lev) + radclru
698               else
699                  radclru = radlu
700                  clrurad(lev) = urad(lev)
701               endif
702               if (idrv .eq. 1) then
703                  if (iclddn.eq.1) then
704                     d_radclru_dt = d_radclru_dt * (1.0_rb - atrans(lev))
705                     d_clrurad_dt(lev) = d_clrurad_dt(lev) + d_radclru_dt
706                  else
707                     d_radclru_dt = d_radlu_dt
708                     d_clrurad_dt(lev) = d_urad_dt(lev)
709                  endif
710               endif
711         enddo
712
713! Increment g-point counter
714         igc = igc + 1
715! Return to continue radiative transfer for all g-channels in present band
716         if (igc .le. ngs(iband)) go to 1000
717
718! Process longwave output from band.
719! Calculate upward, downward, and net flux.
720         do lev = nlayers, 0, -1
721            uflux(lev) = urad(lev)*wtdiff
722            dflux(lev) = drad(lev)*wtdiff
723            urad(lev) = 0.0_rb
724            drad(lev) = 0.0_rb
725            totuflux(lev) = totuflux(lev) + uflux(lev) * delwave(iband)
726            totdflux(lev) = totdflux(lev) + dflux(lev) * delwave(iband)
727            uclfl(lev) = clrurad(lev)*wtdiff
728            dclfl(lev) = clrdrad(lev)*wtdiff
729            clrurad(lev) = 0.0_rb
730            clrdrad(lev) = 0.0_rb
731            totuclfl(lev) = totuclfl(lev) + uclfl(lev) * delwave(iband)
732            totdclfl(lev) = totdclfl(lev) + dclfl(lev) * delwave(iband)
733         enddo
734
735! Calculate total change in upward flux wrt surface temperature
736         if (idrv .eq. 1) then
737            do lev = nlayers, 0, -1
738               duflux_dt(lev) = d_urad_dt(lev) * wtdiff
739               d_urad_dt(lev) = 0.0_rb
740               dtotuflux_dt(lev) = dtotuflux_dt(lev) + duflux_dt(lev) * delwave(iband) * fluxfac
741
742               duclfl_dt(lev) = d_clrurad_dt(lev) * wtdiff
743               d_clrurad_dt(lev) = 0.0_rb
744               dtotuclfl_dt(lev) = dtotuclfl_dt(lev) + duclfl_dt(lev) * delwave(iband) * fluxfac
745            enddo
746         endif
747
748! End spectral band loop
749      enddo
750
751! Calculate fluxes at surface
752      totuflux(0) = totuflux(0) * fluxfac
753      totdflux(0) = totdflux(0) * fluxfac
754      fnet(0) = totuflux(0) - totdflux(0)
755
756      totuclfl(0) = totuclfl(0) * fluxfac
757      totdclfl(0) = totdclfl(0) * fluxfac
758      fnetc(0) = totuclfl(0) - totdclfl(0)
759
760! Calculate fluxes at model levels
761      do lev = 1, nlayers
762         totuflux(lev) = totuflux(lev) * fluxfac
763         totdflux(lev) = totdflux(lev) * fluxfac
764         fnet(lev) = totuflux(lev) - totdflux(lev)
765         totuclfl(lev) = totuclfl(lev) * fluxfac
766         totdclfl(lev) = totdclfl(lev) * fluxfac
767         fnetc(lev) = totuclfl(lev) - totdclfl(lev)
768         l = lev - 1
769
770! Calculate heating rates at model layers
771         htr(l)=heatfac*(fnet(l)-fnet(lev))/(pz(l)-pz(lev)) 
772         htrc(l)=heatfac*(fnetc(l)-fnetc(lev))/(pz(l)-pz(lev)) 
773      enddo
774
775! Set heating rate to zero in top layer
776      htr(nlayers) = 0.0_rb
777      htrc(nlayers) = 0.0_rb
778
779      end subroutine rtrnmr
780
781      end module rrtmg_lw_rtrnmr
782
Note: See TracBrowser for help on using the repository browser.