source: palm/trunk/LIB/rrtmg/rrtmg_sw_cldprmc.f90 @ 2993

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

Added support for RRTMG radiation code

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