Line data Source code
1 :
2 : module radlw
3 : !-----------------------------------------------------------------------
4 : !
5 : ! Purpose: Longwave radiation calculations.
6 : !
7 : !-----------------------------------------------------------------------
8 : use shr_kind_mod, only: r8 => shr_kind_r8
9 : use ppgrid, only: pcols, pver, pverp
10 : use cam_abortutils, only: endrun
11 : use scamMod, only: single_column, scm_crm_mode
12 : use radconstants, only: nlwbands
13 : implicit none
14 :
15 : private
16 : save
17 :
18 : ! Public methods
19 :
20 : public ::&
21 : radlw_init, &! initialize constants
22 : radclwmx ! driver for longwave radiation code
23 :
24 : ! Private module data
25 :
26 : real(r8) :: gravit_cgs ! Acceleration of gravity (cm/s**2)
27 : real(r8) :: stebol_cgs ! Stefan-Boltzmann constant (CGS)
28 :
29 : !===============================================================================
30 : CONTAINS
31 : !===============================================================================
32 :
33 749232 : subroutine radclwmx(lchnk ,ncol ,doabsems , &
34 : lwupcgs ,tnm ,qnm ,o3 , &
35 : pmid ,pint ,pmln ,piln , &
36 : n2o ,ch4 ,cfc11 ,cfc12 , &
37 : cld ,emis ,pmxrgn ,nmxrgn ,qrl,qrlc, &
38 : flns ,flnt ,flnsc ,flntc ,flwds , &
39 : fldsc ,flut ,flutc , &
40 : fnl ,fcnl ,co2mmr, odap_aer)
41 : !-----------------------------------------------------------------------
42 : !
43 : ! Purpose:
44 : ! Compute longwave radiation heating rates and boundary fluxes
45 : !
46 : ! Method:
47 : ! Uses broad band absorptivity/emissivity method to compute clear sky;
48 : ! assumes randomly overlapped clouds with variable cloud emissivity to
49 : ! include effects of clouds.
50 : !
51 : ! Computes clear sky absorptivity/emissivity at lower frequency (in
52 : ! general) than the model radiation frequency; uses previously computed
53 : ! and stored values for efficiency
54 : !
55 : ! Note: This subroutine contains vertical indexing which proceeds
56 : ! from bottom to top rather than the top to bottom indexing
57 : ! used in the rest of the model.
58 : !
59 : ! Author: B. Collins
60 : !
61 : !-----------------------------------------------------------------------
62 : use radae, only: nbands, radems, radabs, radtpl, abstot_3d, &
63 : absnxt_3d, emstot_3d, ntoplw, radoz2, trcpth
64 : use cam_history, only: outfld
65 : use quicksort, only: quick_sort
66 : use radconstants, only: nlwbands
67 : use phys_control, only: cam_physpkg_is
68 : use ref_pres, only: trop_cloud_top_lev
69 :
70 : integer, parameter :: pverp2=pver+2, pverp3=pver+3, pverp4=pver+4
71 : real(r8), parameter :: cldmin = 1.0e-80_r8
72 :
73 : !------------------------------Arguments--------------------------------
74 : !
75 : ! Input arguments
76 : !
77 : integer, intent(in) :: lchnk ! chunk identifier
78 : integer, intent(in) :: ncol ! number of atmospheric columns
79 : logical, intent(in) :: doabsems ! True => abs/emiss calculation this timestep
80 :
81 : ! maximally overlapped region.
82 : ! 0->pmxrgn(i,1) is range of pmid for
83 : ! 1st region, pmxrgn(i,1)->pmxrgn(i,2) for
84 : ! 2nd region, etc
85 : integer, intent(in) :: nmxrgn(pcols) ! Number of maximally overlapped regions
86 :
87 : real(r8), intent(in) :: pmxrgn(pcols,pverp) ! Maximum values of pmid for each
88 : real(r8), intent(in) :: lwupcgs(pcols) ! Longwave up flux in CGS units
89 : !
90 : ! Input arguments which are only passed to other routines
91 : !
92 : real(r8), intent(in) :: tnm(pcols,pver) ! Level temperature
93 : real(r8), intent(in) :: qnm(pcols,pver) ! Level moisture field
94 : real(r8), intent(in) :: o3(pcols,pver) ! ozone mass mixing ratio
95 : real(r8), intent(in) :: pmid(pcols,pver) ! Level pressure
96 : real(r8), intent(in) :: pint(pcols,pverp) ! Model interface pressure
97 : real(r8), intent(in) :: pmln(pcols,pver) ! Ln(pmid)
98 : real(r8), intent(in) :: piln(pcols,pverp) ! Ln(pint)
99 : real(r8), intent(in) :: co2mmr(pcols) ! co2 column mean mass mixing ratio
100 : real(r8), intent(in) :: n2o(pcols,pver) ! nitrous oxide mass mixing ratio
101 : real(r8), intent(in) :: ch4(pcols,pver) ! methane mass mixing ratio
102 : real(r8), intent(in) :: cfc11(pcols,pver) ! cfc11 mass mixing ratio
103 : real(r8), intent(in) :: cfc12(pcols,pver) ! cfc12 mass mixing ratio
104 : real(r8), intent(in) :: cld(pcols,pver) ! Cloud cover
105 : real(r8), intent(in) :: emis(pcols,pver) ! Cloud emissivity
106 :
107 : ! [fraction] absorbtion optical depth, cumulative from top
108 : real(r8), intent(in) :: odap_aer(pcols,pver,nlwbands)
109 :
110 :
111 : !
112 : ! Output arguments
113 : !
114 : real(r8), intent(out) :: qrl (pcols,pver) ! Longwave heating rate
115 : real(r8), intent(out) :: qrlc(pcols,pver) ! Clearsky longwave heating rate
116 : real(r8), intent(out) :: flns(pcols) ! Surface cooling flux
117 : real(r8), intent(out) :: flnt(pcols) ! Net outgoing flux
118 : real(r8), intent(out) :: flut(pcols) ! Upward flux at top of model
119 : real(r8), intent(out) :: flnsc(pcols) ! Clear sky surface cooing
120 : real(r8), intent(out) :: flntc(pcols) ! Net clear sky outgoing flux
121 : real(r8), intent(out) :: flutc(pcols) ! Upward clear-sky flux at top of model
122 : real(r8), intent(out) :: flwds(pcols) ! Down longwave flux at surface
123 : real(r8), intent(out) :: fldsc(pcols) ! Down clear-sky longwave flux at surface
124 : real(r8),intent(out) :: fcnl(pcols,pverp) ! clear sky net flux at interfaces
125 : real(r8),intent(out) :: fnl(pcols,pverp) ! net flux at interfaces
126 : !
127 : !---------------------------Local variables-----------------------------
128 : !
129 : ! Implicit save here is fine since this should have the same value for
130 : ! the entire run.
131 : integer :: ntopcld = 2 ! Lowest layer without clouds.
132 : ! Shouldn't this be turned off by default?
133 :
134 : integer i ! Longitude index
135 : integer ilon ! Longitude index
136 : integer ii ! Longitude index
137 : integer iimx ! Longitude index (max overlap)
138 : integer k ! Level index
139 : integer k1 ! Level index
140 : integer k2 ! Level index
141 : integer k3 ! Level index
142 : integer km ! Level index
143 : integer km1 ! Level index
144 : integer km3 ! Level index
145 : integer km4 ! Level index
146 : integer irgn ! Index for max-overlap regions
147 : integer l ! Index for clouds to overlap
148 : integer l1 ! Index for clouds to overlap
149 : integer n ! Counter
150 :
151 : !
152 : real(r8) :: plco2(pcols,pverp) ! Path length co2
153 : real(r8) :: plh2o(pcols,pverp) ! Path length h2o
154 : real(r8) tmp(pcols) ! Temporary workspace
155 : real(r8) tmp2(pcols) ! Temporary workspace
156 : real(r8) tmp3(0:pverp) ! Temporary workspace
157 : real(r8) tmp4 ! Temporary workspace
158 : real(r8) tfdl ! Temporary workspace
159 : real(r8) tful ! Temporary workspace
160 : real(r8) absbt(pcols) ! Downward emission at model top
161 : real(r8) plol(pcols,pverp) ! O3 pressure wghted path length
162 : real(r8) plos(pcols,pverp) ! O3 path length
163 : real(r8) co2em(pcols,pverp) ! Layer co2 normalized planck funct. derivative
164 : real(r8) co2eml(pcols,pver) ! Interface co2 normalized planck funct. deriv.
165 : real(r8) delt(pcols) ! Diff t**4 mid layer to top interface
166 : real(r8) delt1(pcols) ! Diff t**4 lower intrfc to mid layer
167 : real(r8) bk1(pcols) ! Absrptvty for vertical quadrature
168 : real(r8) bk2(pcols) ! Absrptvty for vertical quadrature
169 : real(r8) cldp(pcols,pverp) ! Cloud cover with extra layer
170 : real(r8) ful(pcols,pverp) ! Total upwards longwave flux
171 : real(r8) fsul(pcols,pverp) ! Clear sky upwards longwave flux
172 : real(r8) fdl(pcols,pverp) ! Total downwards longwave flux
173 : real(r8) fsdl(pcols,pverp) ! Clear sky downwards longwv flux
174 : real(r8) fclb4(pcols,-1:pver) ! Sig t**4 for cld bottom interfc
175 : real(r8) fclt4(pcols,0:pver) ! Sig t**4 for cloud top interfc
176 : real(r8) s(pcols,pverp,pverp) ! Flx integral sum
177 : real(r8) tplnka(pcols,pverp) ! Planck fnctn temperature
178 : real(r8) s2c(pcols,pverp) ! H2o cont amount
179 : real(r8) tcg(pcols,pverp) ! H2o-mass-wgted temp. (Curtis-Godson approx.)
180 : real(r8) w(pcols,pverp) ! H2o path
181 : real(r8) tplnke(pcols) ! Planck fnctn temperature
182 : real(r8) h2otr(pcols,pverp) ! H2o trnmsn for o3 overlap
183 : real(r8) co2t(pcols,pverp) ! Prs wghted temperature path
184 : real(r8) tint(pcols,pverp) ! Interface temperature
185 : real(r8) tint4(pcols,pverp) ! Interface temperature**4
186 : real(r8) tlayr(pcols,pverp) ! Level temperature
187 : real(r8) tlayr4(pcols,pverp) ! Level temperature**4
188 : real(r8) plh2ob(nbands,pcols,pverp)! Pressure weighted h2o path with
189 : ! Hulst-Curtis-Godson temp. factor
190 : ! for H2O bands
191 : real(r8) wb(nbands,pcols,pverp) ! H2o path length with
192 : ! Hulst-Curtis-Godson temp. factor
193 : ! for H2O bands
194 :
195 : real(r8) cld0 ! previous cloud amt (for max overlap)
196 : real(r8) cld1 ! next cloud amt (for max overlap)
197 : real(r8) emx(0:pverp) ! Emissivity factors (max overlap)
198 : real(r8) emx0 ! Emissivity factors for BCs (max overlap)
199 : real(r8) trans ! 1 - emis
200 : real(r8) asort(pver) ! 1 - cloud amounts to be sorted for max ovrlp.
201 : real(r8) atmp ! Temporary storage for sort when nxs = 2
202 : real(r8) maxcld(pcols) ! Maximum cloud at any layer
203 :
204 : integer indx(pcols) ! index vector of gathered array values
205 : integer indxmx(pcols+1,pverp)! index vector of gathered array values
206 : ! integer indxmx(pcols,pverp)! index vector of gathered array values
207 : ! (max overlap)
208 : integer nrgn(pcols) ! Number of max overlap regions at longitude
209 : integer npts ! number of values satisfying some criterion
210 : integer ncolmx(pverp) ! number of columns with clds in region
211 : integer kx1(pcols,pverp) ! Level index for top of max-overlap region
212 : integer kx2(pcols,0:pverp)! Level index for bottom of max-overlap region
213 : integer kxs(0:pverp,pcols,pverp)! Level indices for cld layers sorted by cld()
214 : ! in descending order
215 : integer nxs(pcols,pverp) ! Number of cloudy layers between kx1 and kx2
216 : integer nxsk ! Number of cloudy layers between (kx1/kx2)&k
217 : integer ksort(0:pverp) ! Level indices of cloud amounts to be sorted
218 : ! for max ovrlp. calculation
219 : integer ktmp ! Temporary storage for sort when nxs = 2
220 :
221 : !
222 : ! Pointer variables to 3d structures
223 : !
224 749232 : real(r8), pointer :: abstot(:,:,:)
225 749232 : real(r8), pointer :: absnxt(:,:,:)
226 749232 : real(r8), pointer :: emstot(:,:)
227 :
228 : ! [fraction] Total transmission between interfaces k1 and k2
229 : real(r8) :: aer_trn_ttl(pcols,pverp,pverp,nlwbands)
230 : !
231 : ! Trace gas variables
232 : !
233 : real(r8) ucfc11(pcols,pverp) ! CFC11 path length
234 : real(r8) ucfc12(pcols,pverp) ! CFC12 path length
235 : real(r8) un2o0(pcols,pverp) ! N2O path length
236 : real(r8) un2o1(pcols,pverp) ! N2O path length (hot band)
237 : real(r8) uch4(pcols,pverp) ! CH4 path length
238 : real(r8) uco211(pcols,pverp) ! CO2 9.4 micron band path length
239 : real(r8) uco212(pcols,pverp) ! CO2 9.4 micron band path length
240 : real(r8) uco213(pcols,pverp) ! CO2 9.4 micron band path length
241 : real(r8) uco221(pcols,pverp) ! CO2 10.4 micron band path length
242 : real(r8) uco222(pcols,pverp) ! CO2 10.4 micron band path length
243 : real(r8) uco223(pcols,pverp) ! CO2 10.4 micron band path length
244 : real(r8) bn2o0(pcols,pverp) ! pressure factor for n2o
245 : real(r8) bn2o1(pcols,pverp) ! pressure factor for n2o
246 : real(r8) bch4(pcols,pverp) ! pressure factor for ch4
247 : real(r8) uptype(pcols,pverp) ! p-type continuum path length
248 : real(r8) abplnk1(14,pcols,pverp) ! non-nearest layer Plack factor
249 : real(r8) abplnk2(14,pcols,pverp) ! nearest layer factor
250 : !
251 : !
252 : !-----------------------------------------------------------------------
253 : !
254 : !
255 : ! Set pointer variables
256 : !
257 749232 : abstot => abstot_3d(:,:,:,lchnk)
258 749232 : absnxt => absnxt_3d(:,:,:,lchnk)
259 749232 : emstot => emstot_3d(:,:,lchnk)
260 :
261 : !
262 : ! Calculate some temperatures needed to derive absorptivity and
263 : ! emissivity, as well as some h2o path lengths
264 : !
265 : call radtpl(ncol , &
266 : tnm ,lwupcgs ,qnm ,pint ,plco2 ,plh2o , &
267 : tplnka ,s2c ,tcg ,w ,tplnke , &
268 : tint ,tint4 ,tlayr ,tlayr4 ,pmln , &
269 749232 : piln ,plh2ob ,wb ,co2mmr)
270 :
271 :
272 749232 : if (doabsems) then
273 : !
274 : ! Compute ozone path lengths at frequency of a/e calculation.
275 : !
276 68112 : call radoz2(ncol, o3, pint, plol, plos)
277 : !
278 : ! Compute trace gas path lengths
279 : !
280 : call trcpth(ncol , &
281 : tnm ,pint ,cfc11 ,cfc12 ,n2o , &
282 : ch4 ,qnm ,ucfc11 ,ucfc12 ,un2o0 , &
283 : un2o1 ,uch4 ,uco211 ,uco212 ,uco213 , &
284 : uco221 ,uco222 ,uco223 ,bn2o0 ,bn2o1 , &
285 68112 : bch4 ,uptype ,co2mmr)
286 : !
287 : ! Compute transmission through aerosols from (absorption) optical depths
288 : !
289 68112 : call aer_trans_from_od(ncol, odap_aer, aer_trn_ttl)
290 : !
291 : ! Compute total emissivity:
292 : !
293 : call radems(lchnk ,ncol , &
294 : s2c ,tcg ,w ,tplnke ,plh2o , &
295 : pint ,plco2 ,tint ,tint4 ,tlayr , &
296 : tlayr4 ,plol ,plos ,ucfc11 ,ucfc12 , &
297 : un2o0 ,un2o1 ,uch4 ,uco211 ,uco212 , &
298 : uco213 ,uco221 ,uco222 ,uco223 ,uptype , &
299 : bn2o0 ,bn2o1 ,bch4 ,co2em ,co2eml , &
300 : co2t ,h2otr ,abplnk1 ,abplnk2 ,emstot , &
301 : plh2ob ,wb , &
302 68112 : aer_trn_ttl, co2mmr)
303 : !
304 : ! Compute total absorptivity:
305 : !
306 :
307 : call radabs(lchnk ,ncol , &
308 : pmid ,pint ,co2em ,co2eml ,tplnka , &
309 : s2c ,tcg ,w ,h2otr ,plco2 , &
310 : plh2o ,co2t ,tint ,tlayr ,plol , &
311 : plos ,pmln ,piln ,ucfc11 ,ucfc12 , &
312 : un2o0 ,un2o1 ,uch4 ,uco211 ,uco212 , &
313 : uco213 ,uco221 ,uco222 ,uco223 ,uptype , &
314 : bn2o0 ,bn2o1 ,bch4 ,abplnk1 ,abplnk2 , &
315 : abstot ,absnxt ,plh2ob ,wb , &
316 68112 : odap_aer,aer_trn_ttl, co2mmr)
317 : end if
318 : !
319 : ! Compute sums used in integrals (all longitude points)
320 : !
321 : ! Definition of bk1 & bk2 depends on finite differencing. for
322 : ! trapezoidal rule bk1=bk2. trapezoidal rule applied for nonadjacent
323 : ! layers only.
324 : !
325 : ! delt=t**4 in layer above current sigma level km.
326 : ! delt1=t**4 in layer below current sigma level km.
327 : !
328 12510432 : do i=1,ncol
329 11761200 : delt(i) = tint4(i,pver) - tlayr4(i,pverp)
330 11761200 : delt1(i) = tlayr4(i,pverp) - tint4(i,pverp)
331 11761200 : s(i,pverp,pverp) = stebol_cgs*(delt1(i)*absnxt(i,pver,1) + delt (i)*absnxt(i,pver,4))
332 12510432 : s(i,pver,pverp) = stebol_cgs*(delt (i)*absnxt(i,pver,2) + delt1(i)*absnxt(i,pver,3))
333 : end do
334 19480032 : do k=ntoplw,pver-1
335 313510032 : do i=1,ncol
336 294030000 : bk2(i) = (abstot_3d(i,k,pver,lchnk) + abstot_3d(i,k,pverp,lchnk))*0.5_r8
337 294030000 : bk1(i) = bk2(i)
338 312760800 : s(i,k,pverp) = stebol_cgs*(bk2(i)*delt(i) + bk1(i)*delt1(i))
339 : end do
340 : end do
341 : !
342 : ! All k, km>1
343 : !
344 19480032 : do km=pver,ntoplw+1,-1
345 312760800 : do i=1,ncol
346 294030000 : delt(i) = tint4(i,km-1) - tlayr4(i,km)
347 312760800 : delt1(i) = tlayr4(i,km) - tint4(i,km)
348 : end do
349 : !CSD$ PARALLEL DO PRIVATE( i, k, bk1, bk2 )
350 525211632 : do k=pverp,ntoplw,-1
351 505731600 : if (k == km) then
352 312760800 : do i=1,ncol
353 294030000 : bk2(i) = absnxt(i,km-1,4)
354 312760800 : bk1(i) = absnxt(i,km-1,1)
355 : end do
356 487000800 : else if (k == km-1) then
357 312760800 : do i=1,ncol
358 294030000 : bk2(i) = absnxt(i,km-1,2)
359 312760800 : bk1(i) = absnxt(i,km-1,3)
360 : end do
361 : else
362 7819020000 : do i=1,ncol
363 7350750000 : bk2(i) = (abstot_3d(i,k,km-1,lchnk) + abstot_3d(i,k,km,lchnk))*0.5_r8
364 7819020000 : bk1(i) = bk2(i)
365 : end do
366 : end if
367 8463272400 : do i=1,ncol
368 8444541600 : s(i,k,km) = s(i,k,km+1) + stebol_cgs*(bk2(i)*delt(i) + bk1(i)*delt1(i))
369 : end do
370 : end do
371 : !CSD$ END PARALLEL DO
372 : end do
373 : !
374 : ! Computation of clear sky fluxes always set first level of fsul
375 : !
376 12510432 : do i=1,ncol
377 12510432 : fsul(i,pverp) = lwupcgs(i)
378 : end do
379 : !
380 : ! Downward clear sky fluxes store intermediate quantities in down flux
381 : ! Initialize fluxes to clear sky values.
382 : !
383 12510432 : do i=1,ncol
384 11761200 : tmp(i) = fsul(i,pverp) - stebol_cgs*tint4(i,pverp)
385 11761200 : fsul(i,ntoplw) = fsul(i,pverp) - abstot_3d(i,ntoplw,pverp,lchnk)*tmp(i) + s(i,ntoplw,ntoplw+1)
386 12510432 : fsdl(i,ntoplw) = stebol_cgs*(tplnke(i)**4)*emstot(i,ntoplw)
387 : end do
388 : !
389 : ! fsdl(i,pverp) assumes isothermal layer
390 : !
391 19480032 : do k=ntoplw+1,pver
392 313510032 : do i=1,ncol
393 294030000 : fsul(i,k) = fsul(i,pverp) - abstot_3d(i,k,pverp,lchnk)*tmp(i) + s(i,k,k+1)
394 312760800 : fsdl(i,k) = stebol_cgs*(tplnke(i)**4)*emstot(i,k) - (s(i,k,ntoplw+1) - s(i,k,k+1))
395 : end do
396 : end do
397 : !
398 : ! Store the downward emission from level 1 = total gas emission * sigma
399 : ! t**4. fsdl does not yet include all terms
400 : !
401 12510432 : do i=1,ncol
402 11761200 : absbt(i) = stebol_cgs*(tplnke(i)**4)*emstot(i,pverp)
403 12510432 : fsdl(i,pverp) = absbt(i) - s(i,pverp,ntoplw+1)
404 : end do
405 :
406 20978496 : do k = ntoplw,pverp
407 338530896 : do i = 1,ncol
408 337781664 : fcnl(i,k) = fsul(i,k) - fsdl(i,k)
409 : end do
410 : end do
411 : !
412 : !----------------------------------------------------------------------
413 : ! Modifications for clouds -- max/random overlap assumption
414 : !
415 : ! The column is divided into sets of adjacent layers, called regions,
416 : ! in which the clouds are maximally overlapped. The clouds are
417 : ! randomly overlapped between different regions. The number of
418 : ! regions in a column is set by nmxrgn, and the range of pressures
419 : ! included in each region is set by pmxrgn. The max/random overlap
420 : ! can be written in terms of the solutions of random overlap with
421 : ! cloud amounts = 1. The random overlap assumption is equivalent to
422 : ! setting the flux boundary conditions (BCs) at the edges of each region
423 : ! equal to the mean all-sky flux at those boundaries. Since the
424 : ! emissivity array for propogating BCs is only computed for the
425 : ! TOA BC, the flux BCs elsewhere in the atmosphere have to be formulated
426 : ! in terms of solutions to the random overlap equations. This is done
427 : ! by writing the flux BCs as the sum of a clear-sky flux and emission
428 : ! from a cloud outside the region weighted by an emissivity. This
429 : ! emissivity is determined from the location of the cloud and the
430 : ! flux BC.
431 : !
432 : ! Copy cloud amounts to buffer with extra layer (needed for overlap logic)
433 : !
434 :
435 749232 : ntopcld = max(ntopcld, trop_cloud_top_lev)
436 :
437 25770096 : cldp(:ncol,1:ntopcld) = 0.0_r8
438 300999600 : cldp(:ncol,ntopcld+1:pver) = cld(:ncol,ntopcld+1:pver)
439 12510432 : cldp(:ncol,pverp) = 0.0_r8
440 : !
441 : !
442 : ! Select only those locations where there are no clouds
443 : ! (maximum cloud fraction <= 1.e-3 treated as clear)
444 : ! Set all-sky fluxes to clear-sky values.
445 : !
446 749232 : maxcld(1:ncol) = maxval(cldp(1:ncol,ntoplw:pver),dim=2)
447 :
448 749232 : npts = 0
449 12510432 : do i=1,ncol
450 12510432 : if (maxcld(i) < cldmin) then
451 316626 : npts = npts + 1
452 316626 : indx(npts) = i
453 : end if
454 : end do
455 :
456 1065858 : do ii = 1, npts
457 316626 : i = indx(ii)
458 9614760 : do k = ntoplw, pverp
459 8548902 : fdl(i,k) = fsdl(i,k)
460 8865528 : ful(i,k) = fsul(i,k)
461 : end do
462 : end do
463 : !
464 : ! Select only those locations where there are clouds
465 : !
466 : npts = 0
467 12510432 : do i=1,ncol
468 12510432 : if (maxcld(i) >= cldmin) then
469 11444574 : npts = npts + 1
470 11444574 : indx(npts) = i
471 : end if
472 : end do
473 :
474 : !
475 : ! Initialize all-sky fluxes. fdl(i,1) & ful(i,pverp) are boundary conditions
476 : !
477 12193806 : do ii = 1, npts
478 11444574 : i = indx(ii)
479 11444574 : fdl(i,ntoplw) = fsdl(i,ntoplw)
480 11444574 : fdl(i,pverp) = 0.0_r8
481 11444574 : ful(i,ntoplw) = 0.0_r8
482 11444574 : ful(i,pverp) = fsul(i,pverp)
483 297558924 : do k = ntoplw+1, pver
484 286114350 : fdl(i,k) = 0.0_r8
485 297558924 : ful(i,k) = 0.0_r8
486 : end do
487 : !
488 : ! Initialize Planck emission from layer boundaries
489 : !
490 309003498 : do k = ntoplw, pver
491 297558924 : fclt4(i,k-1) = stebol_cgs*tint4(i,k)
492 309003498 : fclb4(i,k-1) = stebol_cgs*tint4(i,k+1)
493 : enddo
494 11444574 : fclb4(i,ntoplw-2) = stebol_cgs*tint4(i,ntoplw)
495 11444574 : fclt4(i,pver) = stebol_cgs*tint4(i,pverp)
496 : !
497 : ! Initialize indices for layers to be max-overlapped
498 : !
499 40753317 : do irgn = 0, nmxrgn(i)
500 40753317 : kx2(i,irgn) = ntoplw-1
501 : end do
502 12193806 : nrgn(i) = 0
503 : end do
504 :
505 : !----------------------------------------------------------------------
506 : ! INDEX CALCULATIONS FOR MAX OVERLAP
507 :
508 : !CSD$ PARALLEL DO PRIVATE( ii, ilon, irgn, n, k1, k2, k, indxmx, ncolmx, iimx, i, ksort ) &
509 : !CSD$ PRIVATE( asort, ktmp, atmp, km1, km4, k3, emx0, nxsk, emx, cld0 ) &
510 : !CSD$ PRIVATE( tmp4, l1, tmp3, tfdl, l, cld1, trans, km3, tful )
511 12193806 : do ii = 1, npts
512 11444574 : ilon = indx(ii)
513 :
514 : !
515 : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
516 : !
517 29308743 : do irgn = 1, nmxrgn(ilon)
518 : !
519 : ! Calculate min/max layer indices inside region.
520 : !
521 17864169 : n = 0
522 17864169 : if (kx2(ilon,irgn-1) < pver) then
523 17864169 : nrgn(ilon) = irgn
524 17864169 : k1 = kx2(ilon,irgn-1)+1
525 17864169 : kx1(ilon,irgn) = k1
526 17864169 : kx2(ilon,irgn) = k1 - 1
527 89740248 : do k2 = pver, k1, -1
528 89740248 : if (pmid(ilon,k2) <= pmxrgn(ilon,irgn)) then
529 17864169 : kx2(ilon,irgn) = k2
530 17864169 : exit
531 : end if
532 : end do
533 : !
534 : ! Identify columns with clouds in the given region.
535 : !
536 217832889 : do k = k1, k2
537 217832889 : if (cldp(ilon,k) >= cldmin) then
538 17864169 : n = n+1
539 17864169 : indxmx(n,irgn) = ilon
540 17864169 : exit
541 : endif
542 : end do
543 : endif
544 17864169 : ncolmx(irgn) = n
545 : !
546 : ! Dummy value for handling clear-sky regions
547 : !
548 17864169 : indxmx(ncolmx(irgn)+1,irgn) = ncol+1
549 : !
550 : ! Outer loop over columns with clouds in the max-overlap region
551 : !
552 47172912 : do iimx = 1, ncolmx(irgn)
553 17864169 : i = indxmx(iimx,irgn)
554 : !
555 : ! Sort cloud areas and corresponding level indices.
556 : !
557 17864169 : n = 0
558 315423093 : do k = kx1(i,irgn),kx2(i,irgn)
559 315423093 : if (cldp(i,k) >= cldmin) then
560 78579288 : n = n+1
561 78579288 : ksort(n) = k
562 : !
563 : ! We need indices for clouds in order of largest to smallest, so
564 : ! sort 1-cld in ascending order
565 : !
566 78579288 : asort(n) = 1.0_r8-cldp(i,k)
567 : end if
568 : end do
569 17864169 : nxs(i,irgn) = n
570 : !
571 : ! If nxs(i,irgn) eq 1, no need to sort.
572 : ! If nxs(i,irgn) eq 2, sort by swapping if necessary
573 : ! If nxs(i,irgn) ge 3, sort using local sort routine
574 : !
575 17864169 : if (nxs(i,irgn) == 2) then
576 4105690 : if (asort(2) < asort(1)) then
577 3130963 : ktmp = ksort(1)
578 3130963 : ksort(1) = ksort(2)
579 3130963 : ksort(2) = ktmp
580 :
581 3130963 : atmp = asort(1)
582 3130963 : asort(1) = asort(2)
583 3130963 : asort(2) = atmp
584 : endif
585 13758479 : else if (nxs(i,irgn) >= 3) then
586 11268485 : call quick_sort(asort(1:nxs(i,irgn)),ksort(1:nxs(i,irgn)))
587 : endif
588 :
589 114307626 : do l = 1, nxs(i,irgn)
590 96443457 : kxs(l,i,irgn) = ksort(l)
591 : end do
592 : !
593 : ! End loop over longitude i for fluxes
594 : !
595 : end do
596 : !
597 : ! End loop over regions irgn for max-overlap
598 : !
599 : end do
600 : !
601 : !----------------------------------------------------------------------
602 : ! DOWNWARD FLUXES:
603 : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
604 : !
605 29308743 : do irgn = 1, nmxrgn(ilon)
606 : !
607 : ! Compute clear-sky fluxes for regions without clouds
608 : !
609 17864169 : iimx = 1
610 17864169 : if (ilon < indxmx(iimx,irgn) .and. irgn <= nrgn(ilon)) then
611 : !
612 : ! Calculate emissivity so that downward flux at upper boundary of region
613 : ! can be cast in form of solution for downward flux from cloud above
614 : ! that boundary. Then solutions for fluxes at other levels take form of
615 : ! random overlap expressions. Try to locate "cloud" as close as possible
616 : ! to TOA such that the "cloud" pseudo-emissivity is between 0 and 1.
617 : !
618 0 : k1 = kx1(ilon,irgn)
619 0 : do km1 = ntoplw-2, k1-2
620 0 : km4 = km1+3
621 0 : k2 = k1
622 0 : k3 = k2+1
623 0 : tmp(ilon) = s(ilon,k2,min(k3,pverp))*min(1,pverp2-k3)
624 : emx0 = (fdl(ilon,k1)-fsdl(ilon,k1))/ &
625 0 : ((fclb4(ilon,km1)-s(ilon,k2,km4)+tmp(ilon))- fsdl(ilon,k1))
626 0 : if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
627 : end do
628 0 : km1 = min(km1,k1-2)
629 0 : do k2 = kx1(ilon,irgn)+1, kx2(ilon,irgn)+1
630 0 : k3 = k2+1
631 0 : tmp(ilon) = s(ilon,k2,min(k3,pverp))*min(1,pverp2-k3)
632 : fdl(ilon,k2) = (1.0_r8-emx0)*fsdl(ilon,k2) + &
633 0 : emx0*(fclb4(ilon,km1)-s(ilon,k2,km4)+tmp(ilon))
634 : end do
635 : else if (ilon==indxmx(iimx,irgn) .and. iimx<=ncolmx(irgn)) then
636 17864169 : iimx = iimx+1
637 : end if
638 : !
639 : ! Outer loop over columns with clouds in the max-overlap region
640 : !
641 47172912 : do iimx = 1, ncolmx(irgn)
642 17864169 : i = indxmx(iimx,irgn)
643 :
644 : !
645 : ! Calculate emissivity so that downward flux at upper boundary of region
646 : ! can be cast in form of solution for downward flux from cloud above that
647 : ! boundary. Then solutions for fluxes at other levels take form of
648 : ! random overlap expressions. Try to locate "cloud" as close as possible
649 : ! to TOA such that the "cloud" pseudo-emissivity is between 0 and 1.
650 : !
651 17864169 : k1 = kx1(i,irgn)
652 26075955 : do km1 = ntoplw-2,k1-2
653 26074827 : km4 = km1+3
654 26074827 : k2 = k1
655 26074827 : k3 = k2 + 1
656 26074827 : tmp(i) = s(i,k2,min(k3,pverp))*min(1,pverp2-k3)
657 26074827 : tmp2(i) = s(i,k2,min(km4,pverp))*min(1,pverp2-km4)
658 26074827 : emx0 = (fdl(i,k1)-fsdl(i,k1))/((fclb4(i,km1)-tmp2(i)+tmp(i))-fsdl(i,k1))
659 26075955 : if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
660 : end do
661 17864169 : km1 = min(km1,k1-2)
662 17864169 : ksort(0) = km1 + 1
663 : !
664 : ! Loop to calculate fluxes at level k
665 : !
666 17864169 : nxsk = 0
667 333287262 : do k = kx1(i,irgn), kx2(i,irgn)
668 : !
669 : ! Identify clouds (largest to smallest area) between kx1 and k
670 : ! Since nxsk will increase with increasing k up to nxs(i,irgn), once
671 : ! nxsk == nxs(i,irgn) then use the list constructed for previous k
672 : !
673 297558924 : if (nxsk < nxs(i,irgn)) then
674 : nxsk = 0
675 1710707644 : do l = 1, nxs(i,irgn)
676 1432159636 : k1 = kxs(l,i,irgn)
677 1710707644 : if (k >= k1) then
678 304166143 : nxsk = nxsk + 1
679 304166143 : ksort(nxsk) = k1
680 : endif
681 : end do
682 : endif
683 : !
684 : ! Dummy value of index to insure computation of cloud amt is valid for l=nxsk+1
685 : !
686 297558924 : ksort(nxsk+1) = pverp
687 : !
688 : ! Initialize iterated emissivity factors
689 : !
690 674792305 : do l = 1, nxsk
691 674792305 : emx(l) = emis(i,ksort(l))
692 : end do
693 : !
694 : ! Initialize iterated emissivity factor for bnd. condition at upper interface
695 : !
696 297558924 : emx(0) = emx0
697 : !
698 : ! Initialize previous cloud amounts
699 : !
700 297558924 : cld0 = 1.0_r8
701 : !
702 : ! Indices for flux calculations
703 : !
704 297558924 : k2 = k+1
705 297558924 : k3 = k2+1
706 297558924 : tmp4 = s(i,k2,min(k3,pverp))*min(1,pverp2-k3)
707 : !
708 : ! Special case nxsk == 0
709 : !
710 297558924 : if ( nxsk == 0 ) then
711 199968720 : fdl(i,k2) = fdl(i,k2)+fsdl(i,k2)
712 199968720 : if ( emx0 /= 0.0_r8 ) then
713 29492369 : km1 = ksort(0)-1
714 29492369 : km4 = km1+3
715 : fdl(i,k2) = fdl(i,k2)+emx0* &
716 29492369 : (fclb4(i,km1)-(s(i,k2,min(km4,pverp))*min(1,pverp2-km4))+tmp4-fsdl(i,k2))
717 : end if
718 : cycle
719 : end if ! nxsk == 0
720 :
721 : !
722 : ! Loop over number of cloud levels inside region (biggest to smallest cld area)
723 : !
724 572413789 : do l1 = 0, nxsk
725 474823585 : km1 = ksort(l1)-1
726 474823585 : km4 = km1+3
727 572413789 : tmp3(l1) = fclb4(i,km1)-(s(i,k2,min(km4,pverp))*min(1,pverp2-km4))+tmp4-fsdl(i,k2)
728 : end do
729 :
730 : tfdl = 0.0_r8
731 :
732 572413789 : do l = 1, nxsk+1
733 : !
734 : ! Calculate downward fluxes
735 : !
736 474823585 : cld1 = cldp(i,ksort(l))*min(1,nxsk+1-l)
737 474823585 : if (cld0 /= cld1) then
738 459449201 : tfdl = tfdl+(cld0-cld1)*fsdl(i,k2)
739 2249198454 : do l1 = 0, l - 1
740 2249198454 : tfdl = tfdl+(cld0-cld1)*emx(l1)*tmp3(l1)
741 : end do
742 : endif
743 474823585 : cld0 = cld1
744 : !
745 : ! Multiply emissivity factors by current cloud transmissivity
746 : !
747 572413789 : if (l <= nxsk) then
748 377233381 : k1 = ksort(l)
749 377233381 : trans = 1.0_r8-emis(i,k1)
750 : !
751 : ! Ideally the upper bound on l1 would be l-1, but the sort routine
752 : ! scrambles the order of layers with identical cloud amounts
753 : !
754 3052200275 : do l1 = 0, nxsk
755 3052200275 : if (ksort(l1) < k1) then
756 1337483447 : emx(l1) = emx(l1)*trans
757 : endif
758 : end do
759 : end if
760 : !
761 : ! End loop over number l of cloud levels
762 : !
763 : end do
764 115454373 : fdl(i,k2) = tfdl
765 : !
766 : ! End loop over level k for fluxes
767 : !
768 : end do
769 : !
770 : ! End loop over longitude i for fluxes
771 : !
772 : end do
773 : !
774 : ! End loop over regions irgn for max-overlap
775 : !
776 : end do
777 :
778 : !
779 : !----------------------------------------------------------------------
780 : ! UPWARD FLUXES:
781 : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
782 : !
783 30057975 : do irgn = nmxrgn(ilon), 1, -1
784 : !
785 : ! Compute clear-sky fluxes for regions without clouds
786 : !
787 17864169 : iimx = 1
788 17864169 : if (ilon < indxmx(iimx,irgn) .and. irgn <= nrgn(ilon)) then
789 : !
790 : ! Calculate emissivity so that upward flux at lower boundary of region
791 : ! can be cast in form of solution for upward flux from cloud below that
792 : ! boundary. Then solutions for fluxes at other levels take form of
793 : ! random overlap expressions. Try to locate "cloud" as close as possible
794 : ! to surface such that the "cloud" pseudo-emissivity is between 0 and 1.
795 : ! Include allowance for surface emissivity (both numerator and denominator
796 : ! equal 1)
797 : !
798 0 : k1 = kx2(ilon,irgn)+1
799 0 : if (k1 < pverp) then
800 0 : do km1 = pver-1,kx2(ilon,irgn),-1
801 0 : km3 = km1+2
802 0 : k2 = k1
803 0 : k3 = k2+1
804 0 : tmp(ilon) = s(ilon,k2,min(km3,pverp))* min(1,pverp2-km3)
805 : emx0 = (ful(ilon,k1)-fsul(ilon,k1))/ &
806 0 : ((fclt4(ilon,km1)+s(ilon,k2,k3)-tmp(ilon))- fsul(ilon,k1))
807 0 : if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
808 : end do
809 0 : km1 = max(km1,kx2(ilon,irgn))
810 : else
811 0 : km1 = k1-1
812 0 : km3 = km1+2
813 0 : emx0 = 1.0_r8
814 : endif
815 :
816 0 : do k2 = kx1(ilon,irgn), kx2(ilon,irgn)
817 0 : k3 = k2+1
818 : !
819 : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
820 : !
821 0 : tmp(ilon) = s(ilon,k2,min(km3,pverp))* min(1,pverp2-km3)
822 : ful(ilon,k2) =(1.0_r8-emx0)*fsul(ilon,k2) + emx0* &
823 0 : (fclt4(ilon,km1)+s(ilon,k2,k3)-tmp(ilon))
824 : end do
825 17864169 : else if (ilon==indxmx(iimx,irgn) .and. iimx<=ncolmx(irgn)) then
826 17864169 : iimx = iimx+1
827 : end if
828 : !
829 : ! Outer loop over columns with clouds in the max-overlap region
830 : !
831 47172912 : do iimx = 1, ncolmx(irgn)
832 17864169 : i = indxmx(iimx,irgn)
833 :
834 : !
835 : ! Calculate emissivity so that upward flux at lower boundary of region
836 : ! can be cast in form of solution for upward flux from cloud at that
837 : ! boundary. Then solutions for fluxes at other levels take form of
838 : ! random overlap expressions. Try to locate "cloud" as close as possible
839 : ! to surface such that the "cloud" pseudo-emissivity is between 0 and 1.
840 : ! Include allowance for surface emissivity (both numerator and denominator
841 : ! equal 1)
842 : !
843 17864169 : k1 = kx2(i,irgn)+1
844 17864169 : if (k1 < pverp) then
845 19704636 : do km1 = pver-1,kx2(i,irgn),-1
846 19702327 : km3 = km1+2
847 19702327 : k2 = k1
848 19702327 : k3 = k2+1
849 19702327 : tmp(i) = s(i,k2,min(km3,pverp))*min(1,pverp2-km3)
850 19702327 : emx0 = (ful(i,k1)-fsul(i,k1))/((fclt4(i,km1)+s(i,k2,k3)-tmp(i))-fsul(i,k1))
851 19704636 : if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
852 : end do
853 6419595 : km1 = max(km1,kx2(i,irgn))
854 : else
855 : emx0 = 1.0_r8
856 : km1 = k1-1
857 : endif
858 17864169 : ksort(0) = km1 + 1
859 :
860 : !
861 : ! Loop to calculate fluxes at level k
862 : !
863 17864169 : nxsk = 0
864 333287262 : do k = kx2(i,irgn), kx1(i,irgn), -1
865 : !
866 : ! Identify clouds (largest to smallest area) between k and kx2
867 : ! Since nxsk will increase with decreasing k up to nxs(i,irgn), once
868 : ! nxsk == nxs(i,irgn) then use the list constructed for previous k
869 : !
870 297558924 : if (nxsk < nxs(i,irgn)) then
871 : nxsk = 0
872 700410440 : do l = 1, nxs(i,irgn)
873 602820236 : k1 = kxs(l,i,irgn)
874 700410440 : if (k <= k1) then
875 304166143 : nxsk = nxsk + 1
876 304166143 : ksort(nxsk) = k1
877 : endif
878 : end do
879 : endif
880 : !
881 : ! Dummy value of index to insure computation of cloud amt is valid for l=nxsk+1
882 : !
883 297558924 : ksort(nxsk+1) = pverp
884 : !
885 : ! Initialize iterated emissivity factors
886 : !
887 1504131705 : do l = 1, nxsk
888 1504131705 : emx(l) = emis(i,ksort(l))
889 : end do
890 : !
891 : ! Initialize iterated emissivity factor for bnd. condition at lower interface
892 : !
893 297558924 : emx(0) = emx0
894 : !
895 : ! Initialize previous cloud amounts
896 : !
897 297558924 : cld0 = 1.0_r8
898 : !
899 : ! Indices for flux calculations
900 : !
901 297558924 : k2 = k
902 297558924 : k3 = k2+1
903 : !
904 : ! Special case nxsk == 0
905 : !
906 297558924 : if ( nxsk == 0 ) then
907 19010916 : ful(i,k2) = ful(i,k2)+fsul(i,k2)
908 19010916 : if ( emx0 /= 0.0_r8 ) then
909 19010916 : km1 = ksort(0)-1
910 19010916 : km3 = km1+2
911 : !
912 : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
913 : !
914 : ful(i,k2) = ful(i,k2)+emx0* &
915 19010916 : (fclt4(i,km1)+s(i,k2,k3)-(s(i,k2,min(km3,pverp))*min(1,pverp2-km3))-fsul(i,k2))
916 :
917 : end if
918 : cycle
919 : end if
920 : !
921 : ! Loop over number of cloud levels inside region (biggest to smallest cld area)
922 : !
923 1763668797 : do l1 = 0, nxsk
924 1485120789 : km1 = ksort(l1)-1
925 1485120789 : km3 = km1+2
926 : !
927 : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
928 : !
929 1763668797 : tmp3(l1) = fclt4(i,km1)+s(i,k2,k3)-(s(i,k2,min(km3,pverp))*min(1,pverp2-km3))- fsul(i,k2)
930 : end do
931 :
932 : tful = 0.0_r8
933 :
934 1763668797 : do l = 1, nxsk+1
935 : !
936 : ! Calculate upward fluxes
937 : !
938 1485120789 : cld1 = cldp(i,ksort(l))*min(1,nxsk+1-l)
939 1485120789 : if (cld0 /= cld1) then
940 1449201158 : tful = tful+(cld0-cld1)*fsul(i,k2)
941 7420388250 : do l1 = 0, l - 1
942 7420388250 : tful = tful+(cld0-cld1)*emx(l1)*tmp3(l1)
943 : end do
944 : endif
945 1485120789 : cld0 = cld1
946 : !
947 : ! Multiply emissivity factors by current cloud transmissivity
948 : !
949 1763668797 : if (l <= nxsk) then
950 1206572781 : k1 = ksort(l)
951 1206572781 : trans = 1.0_r8-emis(i,k1)
952 : !
953 : ! Ideally the upper bound on l1 would be l-1, but the sort routine
954 : ! scrambles the order of layers with identical cloud amounts
955 : !
956 10275808737 : do l1 = 0, nxsk
957 10275808737 : if (ksort(l1) > k1) then
958 4534617978 : emx(l1) = emx(l1)*trans
959 : endif
960 : end do
961 : end if
962 : !
963 : ! End loop over number l of cloud levels
964 : !
965 : end do
966 296412177 : ful(i,k2) = tful
967 : !
968 : ! End loop over level k for fluxes
969 : !
970 : end do
971 : !
972 : ! End loop over longitude i for fluxes
973 : !
974 : end do
975 : !
976 : ! End loop over regions irgn for max-overlap
977 : !
978 : end do
979 : !
980 : ! End outermost longitude loop
981 : !
982 : end do
983 : !CSD$ END PARALLEL DO
984 : !
985 : ! End cloud modification loops
986 : !
987 : !----------------------------------------------------------------------
988 : ! All longitudes: store history tape quantities
989 : !
990 12510432 : do i=1,ncol
991 11761200 : flwds(i) = fdl (i,pverp )
992 11761200 : fldsc(i) = fsdl(i,pverp )
993 11761200 : flns(i) = ful (i,pverp ) - fdl (i,pverp )
994 11761200 : flnsc(i) = fsul(i,pverp ) - fsdl(i,pverp )
995 11761200 : flnt(i) = ful (i,ntoplw) - fdl (i,ntoplw)
996 11761200 : flntc(i) = fsul(i,ntoplw) - fsdl(i,ntoplw)
997 11761200 : flut(i) = ful (i,ntoplw)
998 12510432 : flutc(i) = fsul(i,ntoplw)
999 : end do
1000 :
1001 749232 : if (single_column.and.scm_crm_mode) then
1002 0 : call outfld('FUL ',ful*1.e-3_r8,pcols,lchnk)
1003 0 : call outfld('FDL ',fdl*1.e-3_r8,pcols,lchnk)
1004 0 : call outfld('FULC ',fsul*1.e-3_r8,pcols,lchnk)
1005 0 : call outfld('FDLC ',fsdl*1.e-3_r8,pcols,lchnk)
1006 : endif
1007 :
1008 20978496 : do k = ntoplw,pverp
1009 338530896 : do i = 1,ncol
1010 337781664 : fnl(i,k) = ful(i,k) - fdl(i,k)
1011 : end do
1012 : end do
1013 : !
1014 : ! Computation of longwave heating (J/kg/s)
1015 : !
1016 20229264 : do k=ntoplw,pver
1017 326020464 : do i=1,ncol
1018 917373600 : qrl(i,k) = (ful(i,k) - fdl(i,k) - ful(i,k+1) + fdl(i,k+1))* &
1019 917373600 : 1.E-4_r8*gravit_cgs/((pint(i,k) - pint(i,k+1)))
1020 : qrlc(i,k) = (fsul(i,k) - fsdl(i,k) - fsul(i,k+1) + fsdl(i,k+1))* &
1021 325271232 : 1.E-4_r8*gravit_cgs/((pint(i,k) - pint(i,k+1)))
1022 : end do
1023 : end do
1024 : ! Return 0 above solution domain
1025 749232 : if ( ntoplw > 1 )then
1026 0 : qrl(:ncol,:ntoplw-1) = 0._r8
1027 0 : qrlc(:ncol,:ntoplw-1) = 0._r8
1028 : end if
1029 : !
1030 749232 : return
1031 1498464 : end subroutine radclwmx
1032 :
1033 :
1034 :
1035 : !-------------------------------------------------------------------------------
1036 :
1037 1536 : subroutine radlw_init(gravit, stebol)
1038 : !-----------------------------------------------------------------------
1039 : !
1040 : ! Purpose:
1041 : ! Initialize various constants for radiation scheme; note that
1042 : ! the radiation scheme uses cgs units.
1043 : !-----------------------------------------------------------------------
1044 : !
1045 : ! Input arguments
1046 : !
1047 : real(r8), intent(in) :: gravit ! Acceleration of gravity (MKS)
1048 : real(r8), intent(in) :: stebol ! Stefan-Boltzmann's constant (MKS)
1049 : !
1050 : !-----------------------------------------------------------------------
1051 : !
1052 : ! Set general radiation consts; convert to cgs units where appropriate:
1053 : !
1054 1536 : gravit_cgs = 100._r8*gravit
1055 1536 : stebol_cgs = 1.e3_r8*stebol
1056 :
1057 749232 : end subroutine radlw_init
1058 :
1059 : !-------------------------------------------------------------------------------
1060 68112 : subroutine aer_trans_from_od(ncol, odap_aer, aer_trn_ttl)
1061 : use radconstants, only: nlwbands
1062 : use ppgrid, only: pcols, pver, pverp
1063 : integer, intent(in) :: ncol
1064 : ! [fraction] absorption optical depth, per layer
1065 : real(r8), intent(in) :: odap_aer(pcols,pver,nlwbands)
1066 : ! [fraction] Total transmission between interfaces k1 and k2
1067 : real(r8), intent(out) :: aer_trn_ttl(pcols,pverp,pverp,nlwbands)
1068 : ! [fraction] absorption optical depth, cumulative from top
1069 : real(r8) :: odap_aer_ttl(pcols,pverp,nlwbands)
1070 :
1071 : integer i, k1, k2, bnd_idx ! column iterator, level iterators, band iterator
1072 : ! odap_aer_ttl is cumulative total optical depth from top of atmosphere
1073 68112 : odap_aer_ttl = 0._r8
1074 544896 : do bnd_idx=1,nlwbands
1075 12941280 : do k1=1,pver
1076 207467568 : do i=1,ncol
1077 206990784 : odap_aer_ttl(i,k1+1,bnd_idx) = odap_aer_ttl(i,k1,bnd_idx) + odap_aer(i,k1,bnd_idx)
1078 : end do
1079 : end do
1080 : end do
1081 :
1082 : ! compute transmission from top of atmosphere to level
1083 : ! where angular dependence of optical depth has been
1084 : ! integrated (giving factor 1.666666)
1085 1907136 : do k1=1,pverp
1086 220750992 : aer_trn_ttl(1:pcols,k1,k1,:)=1._r8
1087 : enddo
1088 : aer_trn_ttl(1:ncol,1,2:pverp,1:nlwbands) = &
1089 207535680 : exp(-1.66_r8 * odap_aer_ttl(1:ncol,2:pverp,1:nlwbands) )
1090 :
1091 : ! compute transmission between a given layer (k1) and lower layers.
1092 1770912 : do k1=2,pver
1093 23907312 : do k2=k1+1,pverp
1094 : aer_trn_ttl(1:ncol,k1,k2,1:nlwbands) = &
1095 22136400 : aer_trn_ttl(1:ncol,1,k2,1:nlwbands) / &
1096 2633360400 : aer_trn_ttl(1:ncol,1,k1,1:nlwbands)
1097 : end do
1098 : end do
1099 :
1100 : ! transmission from k1 to k2 is same as transmission from k2 to k1.
1101 1839024 : do k1=2,pverp
1102 25746336 : do k2=1,k1-1
1103 2820053808 : aer_trn_ttl(1:ncol,k1,k2,1:nlwbands)=aer_trn_ttl(1:ncol,k2,k1,1:nlwbands)
1104 : end do
1105 : end do
1106 :
1107 68112 : return
1108 : end subroutine aer_trans_from_od
1109 :
1110 : end module radlw
|
