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 33520 : 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 33520 : real(r8), pointer :: abstot(:,:,:)
225 33520 : real(r8), pointer :: absnxt(:,:,:)
226 33520 : 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 33520 : abstot => abstot_3d(:,:,:,lchnk)
258 33520 : absnxt => absnxt_3d(:,:,:,lchnk)
259 33520 : 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 33520 : piln ,plh2ob ,wb ,co2mmr)
270 :
271 :
272 33520 : if (doabsems) then
273 : !
274 : ! Compute ozone path lengths at frequency of a/e calculation.
275 : !
276 3352 : 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 3352 : bch4 ,uptype ,co2mmr)
286 : !
287 : ! Compute transmission through aerosols from (absorption) optical depths
288 : !
289 3352 : 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 3352 : 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 3352 : 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 519520 : do i=1,ncol
329 486000 : delt(i) = tint4(i,pver) - tlayr4(i,pverp)
330 486000 : delt1(i) = tlayr4(i,pverp) - tint4(i,pverp)
331 486000 : s(i,pverp,pverp) = stebol_cgs*(delt1(i)*absnxt(i,pver,1) + delt (i)*absnxt(i,pver,4))
332 519520 : s(i,pver,pverp) = stebol_cgs*(delt (i)*absnxt(i,pver,2) + delt1(i)*absnxt(i,pver,3))
333 : end do
334 871520 : do k=ntoplw,pver-1
335 13021520 : do i=1,ncol
336 12150000 : bk2(i) = (abstot_3d(i,k,pver,lchnk) + abstot_3d(i,k,pverp,lchnk))*0.5_r8
337 12150000 : bk1(i) = bk2(i)
338 12988000 : 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 871520 : do km=pver,ntoplw+1,-1
345 12988000 : do i=1,ncol
346 12150000 : delt(i) = tint4(i,km-1) - tlayr4(i,km)
347 12988000 : delt1(i) = tlayr4(i,km) - tint4(i,km)
348 : end do
349 : !CSD$ PARALLEL DO PRIVATE( i, k, bk1, bk2 )
350 23497520 : do k=pverp,ntoplw,-1
351 22626000 : if (k == km) then
352 12988000 : do i=1,ncol
353 12150000 : bk2(i) = absnxt(i,km-1,4)
354 12988000 : bk1(i) = absnxt(i,km-1,1)
355 : end do
356 21788000 : else if (k == km-1) then
357 12988000 : do i=1,ncol
358 12150000 : bk2(i) = absnxt(i,km-1,2)
359 12988000 : bk1(i) = absnxt(i,km-1,3)
360 : end do
361 : else
362 324700000 : do i=1,ncol
363 303750000 : bk2(i) = (abstot_3d(i,k,km-1,lchnk) + abstot_3d(i,k,km,lchnk))*0.5_r8
364 324700000 : bk1(i) = bk2(i)
365 : end do
366 : end if
367 351514000 : do i=1,ncol
368 350676000 : 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 519520 : do i=1,ncol
377 519520 : 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 519520 : do i=1,ncol
384 486000 : tmp(i) = fsul(i,pverp) - stebol_cgs*tint4(i,pverp)
385 486000 : fsul(i,ntoplw) = fsul(i,pverp) - abstot_3d(i,ntoplw,pverp,lchnk)*tmp(i) + s(i,ntoplw,ntoplw+1)
386 519520 : fsdl(i,ntoplw) = stebol_cgs*(tplnke(i)**4)*emstot(i,ntoplw)
387 : end do
388 : !
389 : ! fsdl(i,pverp) assumes isothermal layer
390 : !
391 871520 : do k=ntoplw+1,pver
392 13021520 : do i=1,ncol
393 12150000 : fsul(i,k) = fsul(i,pverp) - abstot_3d(i,k,pverp,lchnk)*tmp(i) + s(i,k,k+1)
394 12988000 : 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 519520 : do i=1,ncol
402 486000 : absbt(i) = stebol_cgs*(tplnke(i)**4)*emstot(i,pverp)
403 519520 : fsdl(i,pverp) = absbt(i) - s(i,pverp,ntoplw+1)
404 : end do
405 :
406 938560 : do k = ntoplw,pverp
407 14060560 : do i = 1,ncol
408 14027040 : 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 33520 : ntopcld = max(ntopcld, trop_cloud_top_lev)
436 :
437 1072560 : cldp(:ncol,1:ntopcld) = 0.0_r8
438 12502000 : cldp(:ncol,ntopcld+1:pver) = cld(:ncol,ntopcld+1:pver)
439 519520 : 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 33520 : maxcld(1:ncol) = maxval(cldp(1:ncol,ntoplw:pver),dim=2)
447 :
448 33520 : npts = 0
449 519520 : do i=1,ncol
450 519520 : if (maxcld(i) < cldmin) then
451 21840 : npts = npts + 1
452 21840 : indx(npts) = i
453 : end if
454 : end do
455 :
456 55360 : do ii = 1, npts
457 21840 : i = indx(ii)
458 645040 : do k = ntoplw, pverp
459 589680 : fdl(i,k) = fsdl(i,k)
460 611520 : ful(i,k) = fsul(i,k)
461 : end do
462 : end do
463 : !
464 : ! Select only those locations where there are clouds
465 : !
466 33520 : npts = 0
467 519520 : do i=1,ncol
468 519520 : if (maxcld(i) >= cldmin) then
469 464160 : npts = npts + 1
470 464160 : 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 497680 : do ii = 1, npts
478 464160 : i = indx(ii)
479 464160 : fdl(i,ntoplw) = fsdl(i,ntoplw)
480 464160 : fdl(i,pverp) = 0.0_r8
481 464160 : ful(i,ntoplw) = 0.0_r8
482 464160 : ful(i,pverp) = fsul(i,pverp)
483 12068160 : do k = ntoplw+1, pver
484 11604000 : fdl(i,k) = 0.0_r8
485 12068160 : ful(i,k) = 0.0_r8
486 : end do
487 : !
488 : ! Initialize Planck emission from layer boundaries
489 : !
490 12532320 : do k = ntoplw, pver
491 12068160 : fclt4(i,k-1) = stebol_cgs*tint4(i,k)
492 12532320 : fclb4(i,k-1) = stebol_cgs*tint4(i,k+1)
493 : enddo
494 464160 : fclb4(i,ntoplw-2) = stebol_cgs*tint4(i,ntoplw)
495 464160 : fclt4(i,pver) = stebol_cgs*tint4(i,pverp)
496 : !
497 : ! Initialize indices for layers to be max-overlapped
498 : !
499 1695979 : do irgn = 0, nmxrgn(i)
500 1695979 : kx2(i,irgn) = ntoplw-1
501 : end do
502 497680 : 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 497680 : do ii = 1, npts
512 464160 : ilon = indx(ii)
513 :
514 : !
515 : ! Outermost loop over regions (sets of adjacent layers) to be max overlapped
516 : !
517 1231819 : do irgn = 1, nmxrgn(ilon)
518 : !
519 : ! Calculate min/max layer indices inside region.
520 : !
521 767659 : n = 0
522 767659 : if (kx2(ilon,irgn-1) < pver) then
523 767659 : nrgn(ilon) = irgn
524 767659 : k1 = kx2(ilon,irgn-1)+1
525 767659 : kx1(ilon,irgn) = k1
526 767659 : kx2(ilon,irgn) = k1 - 1
527 4077556 : do k2 = pver, k1, -1
528 4077556 : if (pmid(ilon,k2) <= pmxrgn(ilon,irgn)) then
529 767659 : kx2(ilon,irgn) = k2
530 767659 : exit
531 : end if
532 : end do
533 : !
534 : ! Identify columns with clouds in the given region.
535 : !
536 9143397 : do k = k1, k2
537 9143397 : if (cldp(ilon,k) >= cldmin) then
538 767659 : n = n+1
539 767659 : indxmx(n,irgn) = ilon
540 767659 : exit
541 : endif
542 : end do
543 : endif
544 767659 : ncolmx(irgn) = n
545 : !
546 : ! Dummy value for handling clear-sky regions
547 : !
548 767659 : indxmx(ncolmx(irgn)+1,irgn) = ncol+1
549 : !
550 : ! Outer loop over columns with clouds in the max-overlap region
551 : !
552 1999478 : do iimx = 1, ncolmx(irgn)
553 767659 : i = indxmx(iimx,irgn)
554 : !
555 : ! Sort cloud areas and corresponding level indices.
556 : !
557 767659 : n = 0
558 12835819 : do k = kx1(i,irgn),kx2(i,irgn)
559 12835819 : if (cldp(i,k) >= cldmin) then
560 2976647 : n = n+1
561 2976647 : 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 2976647 : asort(n) = 1.0_r8-cldp(i,k)
567 : end if
568 : end do
569 767659 : 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 767659 : if (nxs(i,irgn) == 2) then
576 173989 : if (asort(2) < asort(1)) then
577 118576 : ktmp = ksort(1)
578 118576 : ksort(1) = ksort(2)
579 118576 : ksort(2) = ktmp
580 :
581 118576 : atmp = asort(1)
582 118576 : asort(1) = asort(2)
583 118576 : asort(2) = atmp
584 : endif
585 593670 : else if (nxs(i,irgn) >= 3) then
586 429795 : call quick_sort(asort(1:nxs(i,irgn)),ksort(1:nxs(i,irgn)))
587 : endif
588 :
589 4511965 : do l = 1, nxs(i,irgn)
590 3744306 : 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 1231819 : do irgn = 1, nmxrgn(ilon)
606 : !
607 : ! Compute clear-sky fluxes for regions without clouds
608 : !
609 767659 : iimx = 1
610 767659 : 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 0 : 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 0 : 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 767659 : else if (ilon==indxmx(iimx,irgn) .and. iimx<=ncolmx(irgn)) then
636 767659 : iimx = iimx+1
637 : end if
638 : !
639 : ! Outer loop over columns with clouds in the max-overlap region
640 : !
641 1999478 : do iimx = 1, ncolmx(irgn)
642 767659 : 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 767659 : k1 = kx1(i,irgn)
652 1383063 : do km1 = ntoplw-2,k1-2
653 1380162 : km4 = km1+3
654 1380162 : k2 = k1
655 1380162 : k3 = k2 + 1
656 1380162 : tmp(i) = s(i,k2,min(k3,pverp))*min(1,pverp2-k3)
657 1380162 : tmp2(i) = s(i,k2,min(km4,pverp))*min(1,pverp2-km4)
658 1380162 : emx0 = (fdl(i,k1)-fsdl(i,k1))/((fclb4(i,km1)-tmp2(i)+tmp(i))-fsdl(i,k1))
659 1383063 : if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
660 : end do
661 767659 : km1 = min(km1,k1-2)
662 767659 : ksort(0) = km1 + 1
663 : !
664 : ! Loop to calculate fluxes at level k
665 : !
666 767659 : nxsk = 0
667 13603478 : 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 12068160 : if (nxsk < nxs(i,irgn)) then
674 11352385 : nxsk = 0
675 63403603 : do l = 1, nxs(i,irgn)
676 52051218 : k1 = kxs(l,i,irgn)
677 63403603 : if (k >= k1) then
678 10851794 : nxsk = nxsk + 1
679 10851794 : 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 12068160 : ksort(nxsk+1) = pverp
687 : !
688 : ! Initialize iterated emissivity factors
689 : !
690 25431775 : do l = 1, nxsk
691 25431775 : emx(l) = emis(i,ksort(l))
692 : end do
693 : !
694 : ! Initialize iterated emissivity factor for bnd. condition at upper interface
695 : !
696 12068160 : emx(0) = emx0
697 : !
698 : ! Initialize previous cloud amounts
699 : !
700 12068160 : cld0 = 1.0_r8
701 : !
702 : ! Indices for flux calculations
703 : !
704 12068160 : k2 = k+1
705 12068160 : k3 = k2+1
706 12068160 : tmp4 = s(i,k2,min(k3,pverp))*min(1,pverp2-k3)
707 : !
708 : ! Special case nxsk == 0
709 : !
710 12068160 : if ( nxsk == 0 ) then
711 8375738 : fdl(i,k2) = fdl(i,k2)+fsdl(i,k2)
712 8375738 : if ( emx0 /= 0.0_r8 ) then
713 1365390 : km1 = ksort(0)-1
714 1365390 : km4 = km1+3
715 5461560 : fdl(i,k2) = fdl(i,k2)+emx0* &
716 6826950 : (fclb4(i,km1)-(s(i,k2,min(km4,pverp))*min(1,pverp2-km4))+tmp4-fsdl(i,k2))
717 : end if
718 8375738 : cycle
719 : end if ! nxsk == 0
720 :
721 : !
722 : ! Loop over number of cloud levels inside region (biggest to smallest cld area)
723 : !
724 20748459 : do l1 = 0, nxsk
725 17056037 : km1 = ksort(l1)-1
726 17056037 : km4 = km1+3
727 20748459 : tmp3(l1) = fclb4(i,km1)-(s(i,k2,min(km4,pverp))*min(1,pverp2-km4))+tmp4-fsdl(i,k2)
728 : end do
729 :
730 3692422 : tfdl = 0.0_r8
731 :
732 20748459 : do l = 1, nxsk+1
733 : !
734 : ! Calculate downward fluxes
735 : !
736 17056037 : cld1 = cldp(i,ksort(l))*min(1,nxsk+1-l)
737 17056037 : if (cld0 /= cld1) then
738 16520276 : tfdl = tfdl+(cld0-cld1)*fsdl(i,k2)
739 78478332 : do l1 = 0, l - 1
740 78478332 : tfdl = tfdl+(cld0-cld1)*emx(l1)*tmp3(l1)
741 : end do
742 : endif
743 17056037 : cld0 = cld1
744 : !
745 : ! Multiply emissivity factors by current cloud transmissivity
746 : !
747 20748459 : if (l <= nxsk) then
748 13363615 : k1 = ksort(l)
749 13363615 : 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 105123789 : do l1 = 0, nxsk
755 105123789 : if (ksort(l1) < k1) then
756 45880087 : 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 4460081 : 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 1265339 : do irgn = nmxrgn(ilon), 1, -1
784 : !
785 : ! Compute clear-sky fluxes for regions without clouds
786 : !
787 767659 : iimx = 1
788 767659 : 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 0 : 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 0 : 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 767659 : else if (ilon==indxmx(iimx,irgn) .and. iimx<=ncolmx(irgn)) then
826 767659 : iimx = iimx+1
827 : end if
828 : !
829 : ! Outer loop over columns with clouds in the max-overlap region
830 : !
831 1999478 : do iimx = 1, ncolmx(irgn)
832 767659 : 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 767659 : k1 = kx2(i,irgn)+1
844 767659 : if (k1 < pverp) then
845 815136 : do km1 = pver-1,kx2(i,irgn),-1
846 814237 : km3 = km1+2
847 814237 : k2 = k1
848 814237 : k3 = k2+1
849 814237 : tmp(i) = s(i,k2,min(km3,pverp))*min(1,pverp2-km3)
850 814237 : emx0 = (ful(i,k1)-fsul(i,k1))/((fclt4(i,km1)+s(i,k2,k3)-tmp(i))-fsul(i,k1))
851 815136 : if (emx0 >= 0.0_r8 .and. emx0 <= 1.0_r8) exit
852 : end do
853 303499 : km1 = max(km1,kx2(i,irgn))
854 : else
855 464160 : emx0 = 1.0_r8
856 464160 : km1 = k1-1
857 : endif
858 767659 : ksort(0) = km1 + 1
859 :
860 : !
861 : ! Loop to calculate fluxes at level k
862 : !
863 767659 : nxsk = 0
864 13603478 : 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 12068160 : if (nxsk < nxs(i,irgn)) then
871 3692422 : nxsk = 0
872 24931184 : do l = 1, nxs(i,irgn)
873 21238762 : k1 = kxs(l,i,irgn)
874 24931184 : if (k <= k1) then
875 10851794 : nxsk = nxsk + 1
876 10851794 : 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 12068160 : ksort(nxsk+1) = pverp
884 : !
885 : ! Initialize iterated emissivity factors
886 : !
887 56244231 : do l = 1, nxsk
888 56244231 : emx(l) = emis(i,ksort(l))
889 : end do
890 : !
891 : ! Initialize iterated emissivity factor for bnd. condition at lower interface
892 : !
893 12068160 : emx(0) = emx0
894 : !
895 : ! Initialize previous cloud amounts
896 : !
897 12068160 : cld0 = 1.0_r8
898 : !
899 : ! Indices for flux calculations
900 : !
901 12068160 : k2 = k
902 12068160 : k3 = k2+1
903 : !
904 : ! Special case nxsk == 0
905 : !
906 12068160 : if ( nxsk == 0 ) then
907 715775 : ful(i,k2) = ful(i,k2)+fsul(i,k2)
908 715775 : if ( emx0 /= 0.0_r8 ) then
909 715775 : km1 = ksort(0)-1
910 715775 : km3 = km1+2
911 : !
912 : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
913 : !
914 2863100 : ful(i,k2) = ful(i,k2)+emx0* &
915 3578875 : (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 715775 : cycle
919 : end if
920 : !
921 : ! Loop over number of cloud levels inside region (biggest to smallest cld area)
922 : !
923 66880841 : do l1 = 0, nxsk
924 55528456 : km1 = ksort(l1)-1
925 55528456 : km3 = km1+2
926 : !
927 : ! If km3 == pver+2, one of the s integrals = 0 (integration limits both = p_s)
928 : !
929 66880841 : 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 11352385 : tful = 0.0_r8
933 :
934 66880841 : do l = 1, nxsk+1
935 : !
936 : ! Calculate upward fluxes
937 : !
938 55528456 : cld1 = cldp(i,ksort(l))*min(1,nxsk+1-l)
939 55528456 : if (cld0 /= cld1) then
940 54178102 : tful = tful+(cld0-cld1)*fsul(i,k2)
941 265270114 : do l1 = 0, l - 1
942 265270114 : tful = tful+(cld0-cld1)*emx(l1)*tmp3(l1)
943 : end do
944 : endif
945 55528456 : cld0 = cld1
946 : !
947 : ! Multiply emissivity factors by current cloud transmissivity
948 : !
949 66880841 : if (l <= nxsk) then
950 44176071 : k1 = ksort(l)
951 44176071 : 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 359975901 : do l1 = 0, nxsk
957 359975901 : if (ksort(l1) > k1) then
958 157899915 : 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 12120044 : 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 519520 : do i=1,ncol
991 486000 : flwds(i) = fdl (i,pverp )
992 486000 : fldsc(i) = fsdl(i,pverp )
993 486000 : flns(i) = ful (i,pverp ) - fdl (i,pverp )
994 486000 : flnsc(i) = fsul(i,pverp ) - fsdl(i,pverp )
995 486000 : flnt(i) = ful (i,ntoplw) - fdl (i,ntoplw)
996 486000 : flntc(i) = fsul(i,ntoplw) - fsdl(i,ntoplw)
997 486000 : flut(i) = ful (i,ntoplw)
998 519520 : flutc(i) = fsul(i,ntoplw)
999 : end do
1000 :
1001 33520 : 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 938560 : do k = ntoplw,pverp
1009 14060560 : do i = 1,ncol
1010 14027040 : 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 905040 : do k=ntoplw,pver
1017 13541040 : do i=1,ncol
1018 126360000 : qrl(i,k) = (ful(i,k) - fdl(i,k) - ful(i,k+1) + fdl(i,k+1))* &
1019 138996000 : 1.E-4_r8*gravit_cgs/((pint(i,k) - pint(i,k+1)))
1020 126360000 : qrlc(i,k) = (fsul(i,k) - fsdl(i,k) - fsul(i,k+1) + fsdl(i,k+1))* &
1021 139867520 : 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 33520 : 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 33520 : return
1031 67040 : end subroutine radclwmx
1032 :
1033 :
1034 :
1035 : !-------------------------------------------------------------------------------
1036 :
1037 1024 : 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 1024 : gravit_cgs = 100._r8*gravit
1055 1024 : stebol_cgs = 1.e3_r8*stebol
1056 :
1057 33520 : end subroutine radlw_init
1058 :
1059 : !-------------------------------------------------------------------------------
1060 3352 : 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 3352 : odap_aer_ttl = 0._r8
1074 26816 : do bnd_idx=1,nlwbands
1075 636880 : do k1=1,pver
1076 9478728 : do i=1,ncol
1077 9455264 : 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 93856 : do k1=1,pverp
1086 10863832 : aer_trn_ttl(1:pcols,k1,k1,:)=1._r8
1087 : enddo
1088 0 : aer_trn_ttl(1:ncol,1,2:pverp,1:nlwbands) = &
1089 9482080 : 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 87152 : do k1=2,pver
1093 1176552 : do k2=k1+1,pverp
1094 2178800 : aer_trn_ttl(1:ncol,k1,k2,1:nlwbands) = &
1095 1089400 : aer_trn_ttl(1:ncol,1,k2,1:nlwbands) / &
1096 122632200 : 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 90504 : do k1=2,pverp
1102 1267056 : do k2=1,k1-1
1103 128909768 : 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 3352 : return
1108 : end subroutine aer_trans_from_od
1109 :
1110 : end module radlw
|
