source: palm/trunk/LIB/rrtmg/rrtmg_sw_setcoef.f90 @ 4318

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

Added support for RRTMG radiation code

File size: 15.8 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_setcoef
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 : mxmol
22      use rrsw_ref, only : pref, preflog, tref
23      use rrsw_vsn, only : hvrset, hnamset
24
25      implicit none
26
27      contains
28
29!----------------------------------------------------------------------------
30      subroutine setcoef_sw(nlayers, pavel, tavel, pz, tz, tbound, coldry, wkl, &
31                            laytrop, layswtch, laylow, jp, jt, jt1, &
32                            co2mult, colch4, colco2, colh2o, colmol, coln2o, &
33                            colo2, colo3, fac00, fac01, fac10, fac11, &
34                            selffac, selffrac, indself, forfac, forfrac, indfor)
35!----------------------------------------------------------------------------
36!
37! Purpose:  For a given atmosphere, calculate the indices and
38! fractions related to the pressure and temperature interpolations.
39
40! Modifications:
41! Original: J. Delamere, AER, Inc. (version 2.5, 02/04/01)
42! Revised: Rewritten and adapted to ECMWF F90, JJMorcrette 030224
43! Revised: For uniform rrtmg formatting, MJIacono, Jul 2006
44
45! ------ Declarations -------
46
47! ----- Input -----
48      integer(kind=im), intent(in) :: nlayers         ! total number of layers
49
50      real(kind=rb), intent(in) :: pavel(:)           ! layer pressures (mb)
51                                                      !    Dimensions: (nlayers)
52      real(kind=rb), intent(in) :: tavel(:)           ! layer temperatures (K)
53                                                      !    Dimensions: (nlayers)
54      real(kind=rb), intent(in) :: pz(0:)             ! level (interface) pressures (hPa, mb)
55                                                      !    Dimensions: (0:nlayers)
56      real(kind=rb), intent(in) :: tz(0:)             ! level (interface) temperatures (K)
57                                                      !    Dimensions: (0:nlayers)
58      real(kind=rb), intent(in) :: tbound             ! surface temperature (K)
59      real(kind=rb), intent(in) :: coldry(:)          ! dry air column density (mol/cm2)
60                                                      !    Dimensions: (nlayers)
61      real(kind=rb), intent(in) :: wkl(:,:)           ! molecular amounts (mol/cm-2)
62                                                      !    Dimensions: (mxmol,nlayers)
63
64! ----- Output -----
65      integer(kind=im), intent(out) :: laytrop        ! tropopause layer index
66      integer(kind=im), intent(out) :: layswtch       !
67      integer(kind=im), intent(out) :: laylow         !
68
69      integer(kind=im), intent(out) :: jp(:)          !
70                                                      !    Dimensions: (nlayers)
71      integer(kind=im), intent(out) :: jt(:)          !
72                                                      !    Dimensions: (nlayers)
73      integer(kind=im), intent(out) :: jt1(:)         !
74                                                      !    Dimensions: (nlayers)
75
76      real(kind=rb), intent(out) :: colh2o(:)         ! column amount (h2o)
77                                                      !    Dimensions: (nlayers)
78      real(kind=rb), intent(out) :: colco2(:)         ! column amount (co2)
79                                                      !    Dimensions: (nlayers)
80      real(kind=rb), intent(out) :: colo3(:)          ! column amount (o3)
81                                                      !    Dimensions: (nlayers)
82      real(kind=rb), intent(out) :: coln2o(:)         ! column amount (n2o)
83                                                      !    Dimensions: (nlayers)
84      real(kind=rb), intent(out) :: colch4(:)         ! column amount (ch4)
85                                                      !    Dimensions: (nlayers)
86      real(kind=rb), intent(out) :: colo2(:)          ! column amount (o2)
87                                                      !    Dimensions: (nlayers)
88      real(kind=rb), intent(out) :: colmol(:)         !
89                                                      !    Dimensions: (nlayers)
90      real(kind=rb), intent(out) :: co2mult(:)        !
91                                                      !    Dimensions: (nlayers)
92
93      integer(kind=im), intent(out) :: indself(:)
94                                                      !    Dimensions: (nlayers)
95      integer(kind=im), intent(out) :: indfor(:)
96                                                      !    Dimensions: (nlayers)
97      real(kind=rb), intent(out) :: selffac(:)
98                                                      !    Dimensions: (nlayers)
99      real(kind=rb), intent(out) :: selffrac(:)
100                                                      !    Dimensions: (nlayers)
101      real(kind=rb), intent(out) :: forfac(:)
102                                                      !    Dimensions: (nlayers)
103      real(kind=rb), intent(out) :: forfrac(:)
104                                                      !    Dimensions: (nlayers)
105
106      real(kind=rb), intent(out) :: &                 !
107                         fac00(:), fac01(:), &        !    Dimensions: (nlayers)
108                         fac10(:), fac11(:) 
109
110! ----- Local -----
111
112      integer(kind=im) :: indbound
113      integer(kind=im) :: indlev0
114      integer(kind=im) :: lay
115      integer(kind=im) :: jp1
116
117      real(kind=rb) :: stpfac
118      real(kind=rb) :: tbndfrac
119      real(kind=rb) :: t0frac
120      real(kind=rb) :: plog
121      real(kind=rb) :: fp
122      real(kind=rb) :: ft
123      real(kind=rb) :: ft1
124      real(kind=rb) :: water
125      real(kind=rb) :: scalefac
126      real(kind=rb) :: factor
127      real(kind=rb) :: co2reg
128      real(kind=rb) :: compfp
129
130
131! Initializations
132      stpfac = 296._rb/1013._rb
133
134      indbound = tbound - 159._rb
135      tbndfrac = tbound - int(tbound)
136      indlev0  = tz(0) - 159._rb
137      t0frac   = tz(0) - int(tz(0))
138
139      laytrop  = 0
140      layswtch = 0
141      laylow   = 0
142
143! Begin layer loop
144      do lay = 1, nlayers
145! Find the two reference pressures on either side of the
146! layer pressure.  Store them in JP and JP1.  Store in FP the
147! fraction of the difference (in ln(pressure)) between these
148! two values that the layer pressure lies.
149
150         plog = log(pavel(lay))
151         jp(lay) = int(36._rb - 5*(plog+0.04_rb))
152         if (jp(lay) .lt. 1) then
153            jp(lay) = 1
154         elseif (jp(lay) .gt. 58) then
155            jp(lay) = 58
156         endif
157         jp1 = jp(lay) + 1
158         fp = 5._rb * (preflog(jp(lay)) - plog)
159
160! Determine, for each reference pressure (JP and JP1), which
161! reference temperature (these are different for each 
162! reference pressure) is nearest the layer temperature but does
163! not exceed it.  Store these indices in JT and JT1, resp.
164! Store in FT (resp. FT1) the fraction of the way between JT
165! (JT1) and the next highest reference temperature that the
166! layer temperature falls.
167
168         jt(lay) = int(3._rb + (tavel(lay)-tref(jp(lay)))/15._rb)
169         if (jt(lay) .lt. 1) then
170            jt(lay) = 1
171         elseif (jt(lay) .gt. 4) then
172            jt(lay) = 4
173         endif
174         ft = ((tavel(lay)-tref(jp(lay)))/15._rb) - real((jt(lay)-3),kind=rb)
175         jt1(lay) = int(3._rb + (tavel(lay)-tref(jp1))/15._rb)
176         if (jt1(lay) .lt. 1) then
177            jt1(lay) = 1
178         elseif (jt1(lay) .gt. 4) then
179            jt1(lay) = 4
180         endif
181         ft1 = ((tavel(lay)-tref(jp1))/15._rb) - real((jt1(lay)-3),kind=rb)
182
183         water = wkl(1,lay)/coldry(lay)
184         scalefac = pavel(lay) * stpfac / tavel(lay)
185
186! If the pressure is less than ~100mb, perform a different
187! set of species interpolations.
188
189         if (plog .le. 4.56_rb) go to 5300
190         laytrop =  laytrop + 1
191         if (plog .ge. 6.62_rb) laylow = laylow + 1
192
193! Set up factors needed to separately include the water vapor
194! foreign-continuum in the calculation of absorption coefficient.
195
196         forfac(lay) = scalefac / (1.+water)
197         factor = (332.0_rb-tavel(lay))/36.0_rb
198         indfor(lay) = min(2, max(1, int(factor)))
199         forfrac(lay) = factor - real(indfor(lay),kind=rb)
200
201! Set up factors needed to separately include the water vapor
202! self-continuum in the calculation of absorption coefficient.
203
204         selffac(lay) = water * forfac(lay)
205         factor = (tavel(lay)-188.0_rb)/7.2_rb
206         indself(lay) = min(9, max(1, int(factor)-7))
207         selffrac(lay) = factor - real((indself(lay) + 7),kind=rb)
208
209! Calculate needed column amounts.
210
211         colh2o(lay) = 1.e-20_rb * wkl(1,lay)
212         colco2(lay) = 1.e-20_rb * wkl(2,lay)
213         colo3(lay) = 1.e-20_rb * wkl(3,lay)
214!           colo3(lay) = 0._rb
215!           colo3(lay) = colo3(lay)/1.16_rb
216         coln2o(lay) = 1.e-20_rb * wkl(4,lay)
217         colch4(lay) = 1.e-20_rb * wkl(6,lay)
218         colo2(lay) = 1.e-20_rb * wkl(7,lay)
219         colmol(lay) = 1.e-20_rb * coldry(lay) + colh2o(lay)
220!           colco2(lay) = 0._rb
221!           colo3(lay) = 0._rb
222!           coln2o(lay) = 0._rb
223!           colch4(lay) = 0._rb
224!           colo2(lay) = 0._rb
225!           colmol(lay) = 0._rb
226         if (colco2(lay) .eq. 0._rb) colco2(lay) = 1.e-32_rb * coldry(lay)
227         if (coln2o(lay) .eq. 0._rb) coln2o(lay) = 1.e-32_rb * coldry(lay)
228         if (colch4(lay) .eq. 0._rb) colch4(lay) = 1.e-32_rb * coldry(lay)
229         if (colo2(lay) .eq. 0._rb) colo2(lay) = 1.e-32_rb * coldry(lay)
230! Using E = 1334.2 cm-1.
231         co2reg = 3.55e-24_rb * coldry(lay)
232         co2mult(lay)= (colco2(lay) - co2reg) * &
233               272.63_rb*exp(-1919.4_rb/tavel(lay))/(8.7604e-4_rb*tavel(lay))
234         goto 5400
235
236! Above laytrop.
237 5300    continue
238
239! Set up factors needed to separately include the water vapor
240! foreign-continuum in the calculation of absorption coefficient.
241
242         forfac(lay) = scalefac / (1.+water)
243         factor = (tavel(lay)-188.0_rb)/36.0_rb
244         indfor(lay) = 3
245         forfrac(lay) = factor - 1.0_rb
246
247! Calculate needed column amounts.
248
249         colh2o(lay) = 1.e-20_rb * wkl(1,lay)
250         colco2(lay) = 1.e-20_rb * wkl(2,lay)
251         colo3(lay)  = 1.e-20_rb * wkl(3,lay)
252         coln2o(lay) = 1.e-20_rb * wkl(4,lay)
253         colch4(lay) = 1.e-20_rb * wkl(6,lay)
254         colo2(lay)  = 1.e-20_rb * wkl(7,lay)
255         colmol(lay) = 1.e-20_rb * coldry(lay) + colh2o(lay)
256         if (colco2(lay) .eq. 0._rb) colco2(lay) = 1.e-32_rb * coldry(lay)
257         if (coln2o(lay) .eq. 0._rb) coln2o(lay) = 1.e-32_rb * coldry(lay)
258         if (colch4(lay) .eq. 0._rb) colch4(lay) = 1.e-32_rb * coldry(lay)
259         if (colo2(lay)  .eq. 0._rb) colo2(lay)  = 1.e-32_rb * coldry(lay)
260         co2reg = 3.55e-24_rb * coldry(lay)
261         co2mult(lay)= (colco2(lay) - co2reg) * &
262               272.63_rb*exp(-1919.4_rb/tavel(lay))/(8.7604e-4_rb*tavel(lay))
263
264         selffac(lay) = 0._rb
265         selffrac(lay)= 0._rb
266         indself(lay) = 0
267
268 5400    continue
269
270! We have now isolated the layer ln pressure and temperature,
271! between two reference pressures and two reference temperatures
272! (for each reference pressure).  We multiply the pressure
273! fraction FP with the appropriate temperature fractions to get
274! the factors that will be needed for the interpolation that yields
275! the optical depths (performed in routines TAUGBn for band n).
276
277         compfp = 1._rb - fp
278         fac10(lay) = compfp * ft
279         fac00(lay) = compfp * (1._rb - ft)
280         fac11(lay) = fp * ft1
281         fac01(lay) = fp * (1._rb - ft1)
282
283! End layer loop
284      enddo
285
286      end subroutine setcoef_sw
287
288!***************************************************************************
289      subroutine swatmref
290!***************************************************************************
291
292      save
293 
294! These pressures are chosen such that the ln of the first pressure
295! has only a few non-zero digits (i.e. ln(PREF(1)) = 6.96000) and
296! each subsequent ln(pressure) differs from the previous one by 0.2.
297
298      pref(:) = (/ &
299          1.05363e+03_rb,8.62642e+02_rb,7.06272e+02_rb,5.78246e+02_rb,4.73428e+02_rb, &
300          3.87610e+02_rb,3.17348e+02_rb,2.59823e+02_rb,2.12725e+02_rb,1.74164e+02_rb, &
301          1.42594e+02_rb,1.16746e+02_rb,9.55835e+01_rb,7.82571e+01_rb,6.40715e+01_rb, &
302          5.24573e+01_rb,4.29484e+01_rb,3.51632e+01_rb,2.87892e+01_rb,2.35706e+01_rb, &
303          1.92980e+01_rb,1.57998e+01_rb,1.29358e+01_rb,1.05910e+01_rb,8.67114e+00_rb, &
304          7.09933e+00_rb,5.81244e+00_rb,4.75882e+00_rb,3.89619e+00_rb,3.18993e+00_rb, &
305          2.61170e+00_rb,2.13828e+00_rb,1.75067e+00_rb,1.43333e+00_rb,1.17351e+00_rb, &
306          9.60789e-01_rb,7.86628e-01_rb,6.44036e-01_rb,5.27292e-01_rb,4.31710e-01_rb, &
307          3.53455e-01_rb,2.89384e-01_rb,2.36928e-01_rb,1.93980e-01_rb,1.58817e-01_rb, &
308          1.30029e-01_rb,1.06458e-01_rb,8.71608e-02_rb,7.13612e-02_rb,5.84256e-02_rb, &
309          4.78349e-02_rb,3.91639e-02_rb,3.20647e-02_rb,2.62523e-02_rb,2.14936e-02_rb, &
310          1.75975e-02_rb,1.44076e-02_rb,1.17959e-02_rb,9.65769e-03_rb /)
311
312      preflog(:) = (/ &
313           6.9600e+00_rb, 6.7600e+00_rb, 6.5600e+00_rb, 6.3600e+00_rb, 6.1600e+00_rb, &
314           5.9600e+00_rb, 5.7600e+00_rb, 5.5600e+00_rb, 5.3600e+00_rb, 5.1600e+00_rb, &
315           4.9600e+00_rb, 4.7600e+00_rb, 4.5600e+00_rb, 4.3600e+00_rb, 4.1600e+00_rb, &
316           3.9600e+00_rb, 3.7600e+00_rb, 3.5600e+00_rb, 3.3600e+00_rb, 3.1600e+00_rb, &
317           2.9600e+00_rb, 2.7600e+00_rb, 2.5600e+00_rb, 2.3600e+00_rb, 2.1600e+00_rb, &
318           1.9600e+00_rb, 1.7600e+00_rb, 1.5600e+00_rb, 1.3600e+00_rb, 1.1600e+00_rb, &
319           9.6000e-01_rb, 7.6000e-01_rb, 5.6000e-01_rb, 3.6000e-01_rb, 1.6000e-01_rb, &
320          -4.0000e-02_rb,-2.4000e-01_rb,-4.4000e-01_rb,-6.4000e-01_rb,-8.4000e-01_rb, &
321          -1.0400e+00_rb,-1.2400e+00_rb,-1.4400e+00_rb,-1.6400e+00_rb,-1.8400e+00_rb, &
322          -2.0400e+00_rb,-2.2400e+00_rb,-2.4400e+00_rb,-2.6400e+00_rb,-2.8400e+00_rb, &
323          -3.0400e+00_rb,-3.2400e+00_rb,-3.4400e+00_rb,-3.6400e+00_rb,-3.8400e+00_rb, &
324          -4.0400e+00_rb,-4.2400e+00_rb,-4.4400e+00_rb,-4.6400e+00_rb /)
325
326! These are the temperatures associated with the respective
327! pressures for the MLS standard atmosphere.
328
329      tref(:) = (/ &
330           2.9420e+02_rb, 2.8799e+02_rb, 2.7894e+02_rb, 2.6925e+02_rb, 2.5983e+02_rb, &
331           2.5017e+02_rb, 2.4077e+02_rb, 2.3179e+02_rb, 2.2306e+02_rb, 2.1578e+02_rb, &
332           2.1570e+02_rb, 2.1570e+02_rb, 2.1570e+02_rb, 2.1706e+02_rb, 2.1858e+02_rb, &
333           2.2018e+02_rb, 2.2174e+02_rb, 2.2328e+02_rb, 2.2479e+02_rb, 2.2655e+02_rb, &
334           2.2834e+02_rb, 2.3113e+02_rb, 2.3401e+02_rb, 2.3703e+02_rb, 2.4022e+02_rb, &
335           2.4371e+02_rb, 2.4726e+02_rb, 2.5085e+02_rb, 2.5457e+02_rb, 2.5832e+02_rb, &
336           2.6216e+02_rb, 2.6606e+02_rb, 2.6999e+02_rb, 2.7340e+02_rb, 2.7536e+02_rb, &
337           2.7568e+02_rb, 2.7372e+02_rb, 2.7163e+02_rb, 2.6955e+02_rb, 2.6593e+02_rb, &
338           2.6211e+02_rb, 2.5828e+02_rb, 2.5360e+02_rb, 2.4854e+02_rb, 2.4348e+02_rb, & 
339           2.3809e+02_rb, 2.3206e+02_rb, 2.2603e+02_rb, 2.2000e+02_rb, 2.1435e+02_rb, &
340           2.0887e+02_rb, 2.0340e+02_rb, 1.9792e+02_rb, 1.9290e+02_rb, 1.8809e+02_rb, &
341           1.8329e+02_rb, 1.7849e+02_rb, 1.7394e+02_rb, 1.7212e+02_rb /)
342
343      end subroutine swatmref
344
345      end module rrtmg_sw_setcoef
346
347
Note: See TracBrowser for help on using the repository browser.