source: palm/tags/release-5.0/LIB/rrtmg/rrtmg_sw_reftra.f90 @ 2977

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

Added support for RRTMG radiation code

File size: 11.0 KB
Line 
1!     path:      $Source$
2!     author:    $Author: mike $
3!     revision:  $Revision: 11661 $
4!     created:   $Date: 2009-05-22 18:22:22 -0400 (Fri, 22 May 2009) $
5
6      module rrtmg_sw_reftra
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 rrsw_tbl, only : tblint, bpade, od_lo, exp_tbl
22      use rrsw_vsn, only : hvrrft, hnamrft
23
24      implicit none
25
26      contains
27
28! --------------------------------------------------------------------
29      subroutine reftra_sw(nlayers, lrtchk, pgg, prmuz, ptau, pw, &
30                           pref, prefd, ptra, ptrad)
31! --------------------------------------------------------------------
32 
33! Purpose: computes the reflectivity and transmissivity of a clear or
34!   cloudy layer using a choice of various approximations.
35!
36! Interface:  *rrtmg_sw_reftra* is called by *rrtmg_sw_spcvrt*
37!
38! Description:
39! explicit arguments :
40! --------------------
41! inputs
42! ------
43!      lrtchk  = .t. for all layers in clear profile
44!      lrtchk  = .t. for cloudy layers in cloud profile
45!              = .f. for clear layers in cloud profile
46!      pgg     = assymetry factor
47!      prmuz   = cosine solar zenith angle
48!      ptau    = optical thickness
49!      pw      = single scattering albedo
50!
51! outputs
52! -------
53!      pref    : collimated beam reflectivity
54!      prefd   : diffuse beam reflectivity
55!      ptra    : collimated beam transmissivity
56!      ptrad   : diffuse beam transmissivity
57!
58!
59! Method:
60! -------
61!      standard delta-eddington, p.i.f.m., or d.o.m. layer calculations.
62!      kmodts  = 1 eddington (joseph et al., 1976)
63!              = 2 pifm (zdunkowski et al., 1980)
64!              = 3 discrete ordinates (liou, 1973)
65!
66!
67! Modifications:
68! --------------
69! Original: J-JMorcrette, ECMWF, Feb 2003
70! Revised for F90 reformatting: MJIacono, AER, Jul 2006
71! Revised to add exponential lookup table: MJIacono, AER, Aug 2007
72!
73! ------------------------------------------------------------------
74
75! ------- Declarations ------
76
77! ------- Input -------
78
79      integer(kind=im), intent(in) :: nlayers
80
81      logical, intent(in) :: lrtchk(:)                         ! Logical flag for reflectivity and
82                                                               ! and transmissivity calculation;
83                                                               !   Dimensions: (nlayers)
84
85      real(kind=rb), intent(in) :: pgg(:)                      ! asymmetry parameter
86                                                               !   Dimensions: (nlayers)
87      real(kind=rb), intent(in) :: ptau(:)                     ! optical depth
88                                                               !   Dimensions: (nlayers)
89      real(kind=rb), intent(in) :: pw(:)                       ! single scattering albedo
90                                                               !   Dimensions: (nlayers)
91      real(kind=rb), intent(in) :: prmuz                       ! cosine of solar zenith angle
92
93! ------- Output -------
94
95      real(kind=rb), intent(inout) :: pref(:)                  ! direct beam reflectivity
96                                                               !   Dimensions: (nlayers+1)
97      real(kind=rb), intent(inout) :: prefd(:)                 ! diffuse beam reflectivity
98                                                               !   Dimensions: (nlayers+1)
99      real(kind=rb), intent(inout) :: ptra(:)                  ! direct beam transmissivity
100                                                               !   Dimensions: (nlayers+1)
101      real(kind=rb), intent(inout) :: ptrad(:)                 ! diffuse beam transmissivity
102                                                               !   Dimensions: (nlayers+1)
103
104! ------- Local -------
105
106      integer(kind=im) :: jk, jl, kmodts
107      integer(kind=im) :: itind
108
109      real(kind=rb) :: tblind
110      real(kind=rb) :: za, za1, za2
111      real(kind=rb) :: zbeta, zdend, zdenr, zdent
112      real(kind=rb) :: ze1, ze2, zem1, zem2, zemm, zep1, zep2
113      real(kind=rb) :: zg, zg3, zgamma1, zgamma2, zgamma3, zgamma4, zgt
114      real(kind=rb) :: zr1, zr2, zr3, zr4, zr5
115      real(kind=rb) :: zrk, zrk2, zrkg, zrm1, zrp, zrp1, zrpp
116      real(kind=rb) :: zsr3, zt1, zt2, zt3, zt4, zt5, zto1
117      real(kind=rb) :: zw, zwcrit, zwo
118
119      real(kind=rb), parameter :: eps = 1.e-08_rb
120
121!     ------------------------------------------------------------------
122
123! Initialize
124
125      hvrrft = '$Revision: 11661 $'
126
127      zsr3=sqrt(3._rb)
128      zwcrit=0.9999995_rb
129      kmodts=2
130
131      do jk=1, nlayers
132         if (.not.lrtchk(jk)) then
133            pref(jk) =0._rb
134            ptra(jk) =1._rb
135            prefd(jk)=0._rb
136            ptrad(jk)=1._rb
137         else
138            zto1=ptau(jk)
139            zw  =pw(jk)
140            zg  =pgg(jk) 
141
142! General two-stream expressions
143
144            zg3= 3._rb * zg
145            if (kmodts == 1) then
146               zgamma1= (7._rb - zw * (4._rb + zg3)) * 0.25_rb
147               zgamma2=-(1._rb - zw * (4._rb - zg3)) * 0.25_rb
148               zgamma3= (2._rb - zg3 * prmuz ) * 0.25_rb
149            else if (kmodts == 2) then 
150               zgamma1= (8._rb - zw * (5._rb + zg3)) * 0.25_rb
151               zgamma2=  3._rb *(zw * (1._rb - zg )) * 0.25_rb
152               zgamma3= (2._rb - zg3 * prmuz ) * 0.25_rb
153            else if (kmodts == 3) then 
154               zgamma1= zsr3 * (2._rb - zw * (1._rb + zg)) * 0.5_rb
155               zgamma2= zsr3 * zw * (1._rb - zg ) * 0.5_rb
156               zgamma3= (1._rb - zsr3 * zg * prmuz ) * 0.5_rb
157            end if
158            zgamma4= 1._rb - zgamma3
159   
160! Recompute original s.s.a. to test for conservative solution
161
162            zwo= zw / (1._rb - (1._rb - zw) * (zg / (1._rb - zg))**2)
163   
164            if (zwo >= zwcrit) then
165! Conservative scattering
166
167               za  = zgamma1 * prmuz 
168               za1 = za - zgamma3
169               zgt = zgamma1 * zto1
170       
171! Homogeneous reflectance and transmittance,
172! collimated beam
173
174               ze1 = min ( zto1 / prmuz , 500._rb)
175!               ze2 = exp( -ze1 )
176
177! Use exponential lookup table for transmittance, or expansion of
178! exponential for low tau
179               if (ze1 .le. od_lo) then
180                  ze2 = 1._rb - ze1 + 0.5_rb * ze1 * ze1
181               else
182                  tblind = ze1 / (bpade + ze1)
183                  itind = tblint * tblind + 0.5_rb
184                  ze2 = exp_tbl(itind)
185               endif
186!
187
188               pref(jk) = (zgt - za1 * (1._rb - ze2)) / (1._rb + zgt)
189               ptra(jk) = 1._rb - pref(jk)
190
191! isotropic incidence
192
193               prefd(jk) = zgt / (1._rb + zgt)
194               ptrad(jk) = 1._rb - prefd(jk)       
195
196! This is applied for consistency between total (delta-scaled) and direct (unscaled)
197! calculations at very low optical depths (tau < 1.e-4) when the exponential lookup
198! table returns a transmittance of 1.0.
199               if (ze2 .eq. 1.0_rb) then
200                  pref(jk) = 0.0_rb
201                  ptra(jk) = 1.0_rb
202                  prefd(jk) = 0.0_rb
203                  ptrad(jk) = 1.0_rb
204               endif
205
206            else
207! Non-conservative scattering
208
209               za1 = zgamma1 * zgamma4 + zgamma2 * zgamma3
210               za2 = zgamma1 * zgamma3 + zgamma2 * zgamma4
211               zrk = sqrt ( zgamma1**2 - zgamma2**2)
212               zrp = zrk * prmuz               
213               zrp1 = 1._rb + zrp
214               zrm1 = 1._rb - zrp
215               zrk2 = 2._rb * zrk
216               zrpp = 1._rb - zrp*zrp
217               zrkg = zrk + zgamma1
218               zr1  = zrm1 * (za2 + zrk * zgamma3)
219               zr2  = zrp1 * (za2 - zrk * zgamma3)
220               zr3  = zrk2 * (zgamma3 - za2 * prmuz )
221               zr4  = zrpp * zrkg
222               zr5  = zrpp * (zrk - zgamma1)
223               zt1  = zrp1 * (za1 + zrk * zgamma4)
224               zt2  = zrm1 * (za1 - zrk * zgamma4)
225               zt3  = zrk2 * (zgamma4 + za1 * prmuz )
226               zt4  = zr4
227               zt5  = zr5
228
229! mji - reformulated code to avoid potential floating point exceptions
230!               zbeta = - zr5 / zr4
231               zbeta = (zgamma1 - zrk) / zrkg
232!!
233       
234! Homogeneous reflectance and transmittance
235
236               ze1 = min ( zrk * zto1, 500._rb)
237               ze2 = min ( zto1 / prmuz , 500._rb)
238!
239! Original
240!              zep1 = exp( ze1 )
241!              zem1 = exp(-ze1 )
242!              zep2 = exp( ze2 )
243!              zem2 = exp(-ze2 )
244!
245! Revised original, to reduce exponentials
246!              zep1 = exp( ze1 )
247!              zem1 = 1._rb / zep1
248!              zep2 = exp( ze2 )
249!              zem2 = 1._rb / zep2
250!
251! Use exponential lookup table for transmittance, or expansion of
252! exponential for low tau
253               if (ze1 .le. od_lo) then
254                  zem1 = 1._rb - ze1 + 0.5_rb * ze1 * ze1
255                  zep1 = 1._rb / zem1
256               else
257                  tblind = ze1 / (bpade + ze1)
258                  itind = tblint * tblind + 0.5_rb
259                  zem1 = exp_tbl(itind)
260                  zep1 = 1._rb / zem1
261               endif
262
263               if (ze2 .le. od_lo) then
264                  zem2 = 1._rb - ze2 + 0.5_rb * ze2 * ze2
265                  zep2 = 1._rb / zem2
266               else
267                  tblind = ze2 / (bpade + ze2)
268                  itind = tblint * tblind + 0.5_rb
269                  zem2 = exp_tbl(itind)
270                  zep2 = 1._rb / zem2
271               endif
272
273! collimated beam
274
275! mji - reformulated code to avoid potential floating point exceptions
276!               zdenr = zr4*zep1 + zr5*zem1
277!               pref(jk) = zw * (zr1*zep1 - zr2*zem1 - zr3*zem2) / zdenr
278!               zdent = zt4*zep1 + zt5*zem1
279!               ptra(jk) = zem2 - zem2 * zw * (zt1*zep1 - zt2*zem1 - zt3*zep2) / zdent
280
281               zdenr = zr4*zep1 + zr5*zem1
282               zdent = zt4*zep1 + zt5*zem1
283               if (zdenr .ge. -eps .and. zdenr .le. eps) then
284                  pref(jk) = eps
285                  ptra(jk) = zem2
286               else
287                  pref(jk) = zw * (zr1*zep1 - zr2*zem1 - zr3*zem2) / zdenr
288                  ptra(jk) = zem2 - zem2 * zw * (zt1*zep1 - zt2*zem1 - zt3*zep2) / zdent
289               endif
290!!
291
292! diffuse beam
293
294               zemm = zem1*zem1
295               zdend = 1._rb / ( (1._rb - zbeta*zemm ) * zrkg)
296               prefd(jk) =  zgamma2 * (1._rb - zemm) * zdend
297               ptrad(jk) =  zrk2*zem1*zdend
298
299            endif
300
301         endif         
302
303      enddo   
304
305      end subroutine reftra_sw
306
307      end module rrtmg_sw_reftra
308
Note: See TracBrowser for help on using the repository browser.