source: palm/trunk/LIB/rrtmg/rrtmg_sw_cldprop.f90 @ 3771

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

Added support for RRTMG radiation code

File size: 17.5 KB
Line 
1!     path:      $Source$
2!     author:    $Author: miacono $
3!     revision:  $Revision: 23308 $
4!     created:   $Date: 2013-12-27 17:23:51 -0500 (Fri, 27 Dec 2013) $
5
6      module rrtmg_sw_cldprop
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 parrrsw, only : nbndsw, jpband, jpb1, jpb2
22      use rrsw_cld, only : extliq1, ssaliq1, asyliq1, &
23                           extice2, ssaice2, asyice2, &
24                           extice3, ssaice3, asyice3, fdlice3, &
25                           abari, bbari, cbari, dbari, ebari, fbari
26      use rrsw_wvn, only : wavenum1, wavenum2
27      use rrsw_vsn, only : hvrcld, hnamcld
28
29      implicit none
30
31      contains
32
33! ----------------------------------------------------------------------------
34      subroutine cldprop_sw(nlayers, inflag, iceflag, liqflag, cldfrac, &
35                            tauc, ssac, asmc, fsfc, ciwp, clwp, rei, rel, &
36                            taucldorig, taucloud, ssacloud, asmcloud)
37! ----------------------------------------------------------------------------
38
39! Purpose: Compute the cloud optical properties for each cloudy layer.
40! Note: Only inflag = 0 and inflag=2/liqflag=1/iceflag=1,2,3 are available;
41! (Hu & Stamnes, Ebert and Curry, Key, and Fu) are implemented.
42
43! ------- Input -------
44
45      integer(kind=im), intent(in) :: nlayers         ! total number of layers
46      integer(kind=im), intent(in) :: inflag          ! see definitions
47      integer(kind=im), intent(in) :: iceflag         ! see definitions
48      integer(kind=im), intent(in) :: liqflag         ! see definitions
49
50      real(kind=rb), intent(in) :: cldfrac(:)         ! cloud fraction
51                                                      !    Dimensions: (nlayers)
52      real(kind=rb), intent(in) :: ciwp(:)            ! cloud ice water path
53                                                      !    Dimensions: (nlayers)
54      real(kind=rb), intent(in) :: clwp(:)            ! cloud liquid water path
55                                                      !    Dimensions: (nlayers)
56      real(kind=rb), intent(in) :: rei(:)             ! cloud ice particle effective size (microns)
57                                                      !    Dimensions: (nlayers)
58                                                      ! specific definition of rei depends on setting of iceflag:
59                                                      ! iceflag = 0: (inactive)
60                                                      !             
61                                                      ! iceflag = 1: ice effective radius, r_ec, (Ebert and Curry, 1992),
62                                                      !              r_ec range is limited to 13.0 to 130.0 microns
63                                                      ! iceflag = 2: ice effective radius, r_k, (Key, Streamer Ref. Manual, 1996)
64                                                      !              r_k range is limited to 5.0 to 131.0 microns
65                                                      ! iceflag = 3: generalized effective size, dge, (Fu, 1996),
66                                                      !              dge range is limited to 5.0 to 140.0 microns
67                                                      !              [dge = 1.0315 * r_ec]
68      real(kind=rb), intent(in) :: rel(:)             ! cloud liquid particle effective radius (microns)
69                                                      !    Dimensions: (nlayers)
70      real(kind=rb), intent(in) :: tauc(:,:)          ! cloud optical depth
71                                                      !    Dimensions: (nbndsw,nlayers)
72      real(kind=rb), intent(in) :: ssac(:,:)          ! single scattering albedo
73                                                      !    Dimensions: (nbndsw,nlayers)
74      real(kind=rb), intent(in) :: asmc(:,:)          ! asymmetry parameter
75                                                      !    Dimensions: (nbndsw,nlayers)
76      real(kind=rb), intent(in) :: fsfc(:,:)          ! forward scattering fraction
77                                                      !    Dimensions: (nbndsw,nlayers)
78
79! ------- Output -------
80
81      real(kind=rb), intent(out) :: taucloud(:,:)     ! cloud optical depth (delta scaled)
82                                                      !    Dimensions: (nlayers,jpband)
83      real(kind=rb), intent(out) :: taucldorig(:,:)   ! cloud optical depth (non-delta scaled)
84                                                      !    Dimensions: (nlayers,jpband)
85      real(kind=rb), intent(out) :: ssacloud(:,:)     ! single scattering albedo (delta scaled)
86                                                      !    Dimensions: (nlayers,jpband)
87      real(kind=rb), intent(out) :: asmcloud(:,:)     ! asymmetry parameter (delta scaled)
88                                                      !    Dimensions: (nlayers,jpband)
89
90! ------- Local -------
91
92!      integer(kind=im) :: ncbands
93      integer(kind=im) :: ib, ib1, ib2, lay, istr, index, icx
94
95      real(kind=rb), parameter :: eps = 1.e-06_rb     ! epsilon
96      real(kind=rb), parameter :: cldmin = 1.e-20_rb  ! minimum value for cloud quantities
97      real(kind=rb) :: cwp                            ! total cloud water path
98      real(kind=rb) :: radliq                         ! cloud liquid droplet radius (microns)
99      real(kind=rb) :: radice                         ! cloud ice effective size (microns)
100      real(kind=rb) :: factor
101      real(kind=rb) :: fint
102      real(kind=rb) :: tauctot(nlayers)               ! band integrated cloud optical depth
103
104      real(kind=rb) :: taucldorig_a, ssacloud_a, taucloud_a, ffp, ffp1, ffpssa
105      real(kind=rb) :: tauiceorig, scatice, ssaice, tauice, tauliqorig, scatliq, ssaliq, tauliq
106
107      real(kind=rb) :: fdelta(jpb1:jpb2)
108      real(kind=rb) :: extcoice(jpb1:jpb2), gice(jpb1:jpb2)
109      real(kind=rb) :: ssacoice(jpb1:jpb2), forwice(jpb1:jpb2)
110      real(kind=rb) :: extcoliq(jpb1:jpb2), gliq(jpb1:jpb2)
111      real(kind=rb) :: ssacoliq(jpb1:jpb2), forwliq(jpb1:jpb2)
112
113! Initialize
114
115      hvrcld = '$Revision: 23308 $'
116
117!      ncbands = 29
118      ib1 = jpb1
119      ib2 = jpb2
120      tauctot(:) = 0._rb
121
122      do lay = 1, nlayers
123         do ib = ib1 , ib2
124            taucldorig(lay,ib) = tauc(ib-15,lay)
125            taucloud(lay,ib) = 0.0_rb
126            ssacloud(lay,ib) = 1.0_rb
127            asmcloud(lay,ib) = 0.0_rb
128            tauctot(lay) = tauctot(lay) + tauc(ib-15,lay)
129         enddo
130      enddo
131
132! Main layer loop
133      do lay = 1, nlayers
134
135         cwp = ciwp(lay) + clwp(lay)
136         if (cldfrac(lay) .ge. cldmin .and. &
137            (cwp .ge. cldmin .or. tauctot(lay) .ge. cldmin)) then
138
139! (inflag=0): Cloud optical properties input directly
140! Cloud optical properties already defined in tauc, ssac, asmc are unscaled;
141! Apply delta-M scaling here
142            if (inflag .eq. 0) then
143
144               do ib = ib1 , ib2
145                  taucldorig_a = tauc(ib-15,lay)
146                  ffp = fsfc(ib-15,lay)
147                  ffp1 = 1.0_rb - ffp
148                  ffpssa = 1.0_rb - ffp * ssac(ib-15,lay)
149                  ssacloud_a = ffp1 * ssac(ib-15,lay) / ffpssa
150                  taucloud_a = ffpssa * taucldorig_a
151
152                  taucldorig(lay,ib) = taucldorig_a
153                  ssacloud(lay,ib) = ssacloud_a
154                  taucloud(lay,ib) = taucloud_a
155                  asmcloud(lay,ib) = (asmc(ib-15,lay) - ffp) / (ffp1)
156               enddo
157
158! (inflag=2): Separate treatement of ice clouds and water clouds.
159            elseif (inflag .eq. 2) then       
160               radice = rei(lay)
161
162! Calculation of absorption coefficients due to ice clouds.
163               if (ciwp(lay) .eq. 0.0_rb) then
164                  do ib = ib1 , ib2
165                     extcoice(ib) = 0.0_rb
166                     ssacoice(ib) = 0.0_rb
167                     gice(ib)     = 0.0_rb
168                     forwice(ib)  = 0.0_rb
169                  enddo
170
171! (iceflag = 1):
172! Note: This option uses Ebert and Curry approach for all particle sizes similar to
173! CAM3 implementation, though this is somewhat ineffective for large ice particles
174               elseif (iceflag .eq. 1) then
175                  if (radice .lt. 13.0_rb .or. radice .gt. 130._rb) stop &
176                     'ICE RADIUS OUT OF BOUNDS'
177                  do ib = ib1, ib2
178                     if (wavenum2(ib) .gt. 1.43e04_rb) then
179                        icx = 1
180                     elseif (wavenum2(ib) .gt. 7.7e03_rb) then
181                        icx = 2
182                     elseif (wavenum2(ib) .gt. 5.3e03_rb) then
183                        icx = 3
184                     elseif (wavenum2(ib) .gt. 4.0e03_rb) then
185                        icx = 4
186                     elseif (wavenum2(ib) .ge. 2.5e03_rb) then
187                        icx = 5
188                     endif
189                     extcoice(ib) = abari(icx) + bbari(icx)/radice
190                     ssacoice(ib) = 1._rb - cbari(icx) - dbari(icx) * radice
191                     gice(ib) = ebari(icx) + fbari(icx) * radice
192
193! Check to ensure upper limit of gice is within physical limits for large particles
194                     if (gice(ib) .ge. 1.0_rb) gice(ib) = 1.0_rb - eps
195                     forwice(ib) = gice(ib)*gice(ib)
196! Check to ensure all calculated quantities are within physical limits.
197                     if (extcoice(ib) .lt. 0.0_rb) stop 'ICE EXTINCTION LESS THAN 0.0'
198                     if (ssacoice(ib) .gt. 1.0_rb) stop 'ICE SSA GRTR THAN 1.0'
199                     if (ssacoice(ib) .lt. 0.0_rb) stop 'ICE SSA LESS THAN 0.0'
200                     if (gice(ib) .gt. 1.0_rb) stop 'ICE ASYM GRTR THAN 1.0'
201                     if (gice(ib) .lt. 0.0_rb) stop 'ICE ASYM LESS THAN 0.0'
202                  enddo
203
204! For iceflag=2 option, ice particle effective radius is limited to 5.0 to 131.0 microns
205
206               elseif (iceflag .eq. 2) then
207                  if (radice .lt. 5.0_rb .or. radice .gt. 131.0_rb) stop 'ICE RADIUS OUT OF BOUNDS'
208                  factor = (radice - 2._rb)/3._rb
209                  index = int(factor)
210                  if (index .eq. 43) index = 42
211                  fint = factor - real(index,kind=rb)
212                  do ib = ib1, ib2
213                     extcoice(ib) = extice2(index,ib) + fint * &
214                                   (extice2(index+1,ib) -  extice2(index,ib))
215                     ssacoice(ib) = ssaice2(index,ib) + fint * &
216                                   (ssaice2(index+1,ib) -  ssaice2(index,ib))
217                     gice(ib) = asyice2(index,ib) + fint * &
218                                   (asyice2(index+1,ib) -  asyice2(index,ib))
219                     forwice(ib) = gice(ib)*gice(ib)
220! Check to ensure all calculated quantities are within physical limits.
221                     if (extcoice(ib) .lt. 0.0_rb) stop 'ICE EXTINCTION LESS THAN 0.0'
222                     if (ssacoice(ib) .gt. 1.0_rb) stop 'ICE SSA GRTR THAN 1.0'
223                     if (ssacoice(ib) .lt. 0.0_rb) stop 'ICE SSA LESS THAN 0.0'
224                     if (gice(ib) .gt. 1.0_rb) stop 'ICE ASYM GRTR THAN 1.0'
225                     if (gice(ib) .lt. 0.0_rb) stop 'ICE ASYM LESS THAN 0.0'
226                  enddo
227
228! For iceflag=3 option, ice particle generalized effective size is limited to 5.0 to 140.0 microns
229
230               elseif (iceflag .eq. 3) then
231                  if (radice .lt. 5.0_rb .or. radice .gt. 140.0_rb) stop 'ICE GENERALIZED EFFECTIVE SIZE OUT OF BOUNDS'
232                  factor = (radice - 2._rb)/3._rb
233                  index = int(factor)
234                  if (index .eq. 46) index = 45
235                  fint = factor - real(index,kind=rb)
236                  do ib = ib1 , ib2
237                     extcoice(ib) = extice3(index,ib) + fint * &
238                                   (extice3(index+1,ib) - extice3(index,ib))
239                     ssacoice(ib) = ssaice3(index,ib) + fint * &
240                                   (ssaice3(index+1,ib) - ssaice3(index,ib))
241                     gice(ib) = asyice3(index,ib) + fint * &
242                               (asyice3(index+1,ib) - asyice3(index,ib))
243                     fdelta(ib) = fdlice3(index,ib) + fint * &
244                                 (fdlice3(index+1,ib) - fdlice3(index,ib))
245                     if (fdelta(ib) .lt. 0.0_rb) stop 'FDELTA LESS THAN 0.0'
246                     if (fdelta(ib) .gt. 1.0_rb) stop 'FDELTA GT THAN 1.0'                     
247                     forwice(ib) = fdelta(ib) + 0.5_rb / ssacoice(ib)
248! See Fu 1996 p. 2067
249                     if (forwice(ib) .gt. gice(ib)) forwice(ib) = gice(ib)
250! Check to ensure all calculated quantities are within physical limits.
251                     if (extcoice(ib) .lt. 0.0_rb) stop 'ICE EXTINCTION LESS THAN 0.0'
252                     if (ssacoice(ib) .gt. 1.0_rb) stop 'ICE SSA GRTR THAN 1.0'
253                     if (ssacoice(ib) .lt. 0.0_rb) stop 'ICE SSA LESS THAN 0.0'
254                     if (gice(ib) .gt. 1.0_rb) stop 'ICE ASYM GRTR THAN 1.0'
255                     if (gice(ib) .lt. 0.0_rb) stop 'ICE ASYM LESS THAN 0.0'
256                  enddo
257
258               endif
259                 
260! Calculation of absorption coefficients due to water clouds.
261                if (clwp(lay) .eq. 0.0_rb) then
262                   do ib = ib1 , ib2
263                      extcoliq(ib) = 0.0_rb
264                      ssacoliq(ib) = 0.0_rb
265                      gliq(ib) = 0.0_rb
266                      forwliq(ib) = 0.0_rb
267                   enddo
268
269                elseif (liqflag .eq. 1) then
270                   radliq = rel(lay)
271                   if (radliq .lt. 2.5_rb .or. radliq .gt. 60._rb) stop &
272                      'LIQUID EFFECTIVE RADIUS OUT OF BOUNDS'
273                   index = int(radliq - 1.5_rb)
274                   if (index .eq. 0) index = 1
275                   if (index .eq. 58) index = 57
276                   fint = radliq - 1.5_rb - real(index,kind=rb)
277                   do ib = ib1 , ib2
278                      extcoliq(ib) = extliq1(index,ib) + fint * &
279                                    (extliq1(index+1,ib) - extliq1(index,ib))
280                      ssacoliq(ib) = ssaliq1(index,ib) + fint * &
281                                    (ssaliq1(index+1,ib) - ssaliq1(index,ib))
282                      if (fint .lt. 0._rb .and. ssacoliq(ib) .gt. 1._rb) &
283                                     ssacoliq(ib) = ssaliq1(index,ib)
284                      gliq(ib) = asyliq1(index,ib) + fint * &
285                                (asyliq1(index+1,ib) - asyliq1(index,ib))
286                      forwliq(ib) = gliq(ib)*gliq(ib)
287! Check to ensure all calculated quantities are within physical limits.
288                      if (extcoliq(ib) .lt. 0.0_rb) stop 'LIQUID EXTINCTION LESS THAN 0.0'
289                      if (ssacoliq(ib) .gt. 1.0_rb) stop 'LIQUID SSA GRTR THAN 1.0'
290                      if (ssacoliq(ib) .lt. 0.0_rb) stop 'LIQUID SSA LESS THAN 0.0'
291                      if (gliq(ib) .gt. 1.0_rb) stop 'LIQUID ASYM GRTR THAN 1.0'
292                      if (gliq(ib) .lt. 0.0_rb) stop 'LIQUID ASYM LESS THAN 0.0'
293                   enddo
294                endif
295
296                do ib = ib1 , ib2
297                   tauliqorig = clwp(lay) * extcoliq(ib)
298                   tauiceorig = ciwp(lay) * extcoice(ib)
299                   taucldorig(lay,ib) = tauliqorig + tauiceorig
300
301                   ssaliq = ssacoliq(ib) * (1.0_rb - forwliq(ib)) / &
302                           (1.0_rb - forwliq(ib) * ssacoliq(ib))
303                   tauliq = (1.0_rb - forwliq(ib) * ssacoliq(ib)) * tauliqorig
304                   ssaice = ssacoice(ib) * (1.0_rb - forwice(ib)) / &
305                           (1.0_rb - forwice(ib) * ssacoice(ib))
306                   tauice = (1.0_rb - forwice(ib) * ssacoice(ib)) * tauiceorig
307
308                   scatliq = ssaliq * tauliq
309                   scatice = ssaice * tauice
310
311                   taucloud(lay,ib) = tauliq + tauice
312
313! Ensure non-zero taucmc and scatice
314                   if (taucloud(lay,ib).eq.0.0_rb) taucloud(lay,ib) = cldmin
315                   if (scatice.eq.0.0_rb) scatice = cldmin
316
317                   ssacloud(lay,ib) = (scatliq + scatice) / taucloud(lay,ib)
318
319                   if (iceflag .eq. 3) then
320! In accordance with the 1996 Fu paper, equation A.3,
321! the moments for ice were calculated depending on whether using spheres
322! or hexagonal ice crystals.
323                      istr = 1
324                      asmcloud(lay,ib) = (1.0_rb/(scatliq+scatice)) * &
325                         (scatliq*(gliq(ib)**istr - forwliq(ib)) / &
326                         (1.0_rb - forwliq(ib)) + scatice * ((gice(ib)-forwice(ib)) / &
327                         (1.0_rb - forwice(ib)))**istr)
328                   else 
329! This code is the standard method for delta-m scaling.
330                      istr = 1
331                      asmcloud(lay,ib) = (scatliq *  &
332                         (gliq(ib)**istr - forwliq(ib)) / &
333                         (1.0_rb - forwliq(ib)) + scatice * (gice(ib)**istr - forwice(ib)) / &
334                         (1.0_rb - forwice(ib)))/(scatliq + scatice)
335                   endif
336
337                enddo
338
339            endif
340
341         endif
342
343! End layer loop
344      enddo
345
346      end subroutine cldprop_sw
347
348      end module rrtmg_sw_cldprop
349
350
Note: See TracBrowser for help on using the repository browser.