Line data Source code
1 : ! path: $Source: /storm/rc1/cvsroot/rc/rrtmg_sw/src/rrtmg_sw_spcvmc.f90,v $
2 : ! author: $Author: mike $
3 : ! revision: $Revision: 1.2 $
4 : ! created: $Date: 2007/08/23 20:40:14 $
5 :
6 : module rrtmg_sw_spcvmc
7 :
8 : ! --------------------------------------------------------------------------
9 : ! | |
10 : ! | Copyright 2002-2007, 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 shr_kind_mod, only: r8 => shr_kind_r8
21 :
22 : use parrrsw, only: nbndsw, ngptsw, mxmol, jpband
23 : use rrsw_tbl, only: tblint, bpade, od_lo, exp_tbl
24 : use rrsw_wvn, only: ngc, ngs
25 : use rrtmg_sw_reftra, only: reftra_sw
26 : use rrtmg_sw_taumol, only: taumol_sw
27 : use rrtmg_sw_vrtqdr, only: vrtqdr_sw
28 :
29 : implicit none
30 :
31 : contains
32 :
33 : ! ---------------------------------------------------------------------------
34 38400 : subroutine spcvmc_sw &
35 : (lchnk, ncol, nlayers, istart, iend, icpr, idelm, iout, &
36 38400 : pavel, tavel, pz, tz, tbound, palbd, palbp, &
37 38400 : pcldfmc, ptaucmc, pasycmc, pomgcmc, ptaormc, &
38 38400 : ptaua, pasya, pomga, prmu0, coldry, wkl, adjflux, &
39 38400 : laytrop, layswtch, laylow, jp, jt, jt1, &
40 38400 : co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, &
41 38400 : fac00, fac01, fac10, fac11, &
42 38400 : selffac, selffrac, indself, forfac, forfrac, indfor, &
43 38400 : pbbfd, pbbfu, pbbcd, pbbcu, puvfd, puvcd, pnifd, pnicd, pnifu, pnicu, &
44 38400 : pbbfddir, pbbcddir, puvfddir, puvcddir, pnifddir, pnicddir, &
45 38400 : pbbfsu, pbbfsd)
46 : ! ---------------------------------------------------------------------------
47 : !
48 : ! Purpose: Contains spectral loop to compute the shortwave radiative fluxes,
49 : ! using the two-stream method of H. Barker and McICA, the Monte-Carlo
50 : ! Independent Column Approximation, for the representation of
51 : ! sub-grid cloud variability (i.e. cloud overlap).
52 : !
53 : ! Interface: *spcvmc_sw* is called from *rrtmg_sw.F90* or rrtmg_sw.1col.F90*
54 : !
55 : ! Method:
56 : ! Adapted from two-stream model of H. Barker;
57 : ! Two-stream model options (selected with kmodts in rrtmg_sw_reftra.F90):
58 : ! 1: Eddington, 2: PIFM, Zdunkowski et al., 3: discret ordinates
59 : !
60 : ! Modifications:
61 : !
62 : ! Original: H. Barker
63 : ! Revision: Merge with RRTMG_SW: J.-J.Morcrette, ECMWF, Feb 2003
64 : ! Revision: Add adjustment for Earth/Sun distance : MJIacono, AER, Oct 2003
65 : ! Revision: Bug fix for use of PALBP and PALBD: MJIacono, AER, Nov 2003
66 : ! Revision: Bug fix to apply delta scaling to clear sky: AER, Dec 2004
67 : ! Revision: Code modified so that delta scaling is not done in cloudy profiles
68 : ! if routine cldprop is used; delta scaling can be applied by swithcing
69 : ! code below if cldprop is not used to get cloud properties.
70 : ! AER, Jan 2005
71 : ! Revision: Modified to use McICA: MJIacono, AER, Nov 2005
72 : ! Revision: Uniform formatting for RRTMG: MJIacono, AER, Jul 2006
73 : ! Revision: Use exponential lookup table for transmittance: MJIacono, AER,
74 : ! Aug 2007
75 : !
76 : ! ------------------------------------------------------------------
77 :
78 : ! ------- Declarations ------
79 :
80 : ! ------- Input -------
81 :
82 : integer, intent(in) :: lchnk
83 :
84 : integer, intent(in) :: nlayers
85 : integer, intent(in) :: istart
86 : integer, intent(in) :: iend
87 : integer, intent(in) :: icpr
88 : integer, intent(in) :: idelm ! delta-m scaling flag
89 : ! [0 = direct and diffuse fluxes are unscaled]
90 : ! [1 = direct and diffuse fluxes are scaled]
91 : integer, intent(in) :: iout
92 : integer, intent(in) :: ncol
93 :
94 : integer, intent(in) :: laytrop(ncol)
95 : integer, intent(in) :: layswtch(ncol)
96 : integer, intent(in) :: laylow(ncol)
97 :
98 : integer, intent(in) :: indfor(ncol,nlayers)
99 : ! Dimensions: (ncol,nlayers)
100 : integer, intent(in) :: indself(ncol,nlayers)
101 : ! Dimensions: (ncol,nlayers)
102 : integer, intent(in) :: jp(ncol,nlayers)
103 : ! Dimensions: (ncol,nlayers)
104 : integer, intent(in) :: jt(ncol,nlayers)
105 : ! Dimensions: (ncol,nlayers)
106 : integer, intent(in) :: jt1(ncol,nlayers)
107 : ! Dimensions: (ncol,nlayers)
108 :
109 : real(kind=r8), intent(in) :: pavel(ncol,nlayers) ! layer pressure (hPa, mb)
110 : ! Dimensions: (ncol,nlayers)
111 : real(kind=r8), intent(in) :: tavel(ncol,nlayers) ! layer temperature (K)
112 : ! Dimensions: (ncol,nlayers)
113 : real(kind=r8), intent(in) :: pz(ncol,0:nlayers) ! level (interface) pressure (hPa, mb)
114 : ! Dimensions: (ncol,0:nlayers)
115 : real(kind=r8), intent(in) :: tz(ncol,0:nlayers) ! level temperatures (hPa, mb)
116 : ! Dimensions: (ncol,0:nlayers)
117 : real(kind=r8), intent(in) :: tbound(ncol) ! surface temperature (K)
118 : real(kind=r8), intent(in) :: wkl(ncol,mxmol,nlayers) ! molecular amounts (mol/cm2)
119 : ! Dimensions: (ncol,mxmol,nlayers)
120 : real(kind=r8), intent(in) :: coldry(ncol,nlayers) ! dry air column density (mol/cm2)
121 : ! Dimensions: (ncol,nlayers)
122 : real(kind=r8), intent(in) :: colmol(ncol,nlayers)
123 : ! Dimensions: (ncol,nlayers)
124 : real(kind=r8), intent(in) :: adjflux(ncol,jpband) ! Earth/Sun distance adjustment
125 : ! Dimensions: (ncol,jpband)
126 :
127 : real(kind=r8), intent(in) :: palbd(ncol,nbndsw) ! surface albedo (diffuse)
128 : ! Dimensions: (ncol,nbndsw)
129 : real(kind=r8), intent(in) :: palbp(ncol,nbndsw) ! surface albedo (direct)
130 : ! Dimensions: (ncol, nbndsw)
131 : real(kind=r8), intent(in) :: prmu0(ncol) ! cosine of solar zenith angle
132 : real(kind=r8), intent(in) :: pcldfmc(ncol,nlayers,ngptsw) ! cloud fraction [mcica]
133 : ! Dimensions: (ncol,nlayers,ngptsw)
134 : real(kind=r8), intent(in) :: ptaucmc(ncol,nlayers,ngptsw) ! cloud optical depth [mcica]
135 : ! Dimensions: (ncol,nlayers,ngptsw)
136 : real(kind=r8), intent(in) :: pasycmc(ncol,nlayers,ngptsw) ! cloud asymmetry parameter [mcica]
137 : ! Dimensions: (ncol,nlayers,ngptsw)
138 : real(kind=r8), intent(in) :: pomgcmc(ncol,nlayers,ngptsw) ! cloud single scattering albedo [mcica]
139 : ! Dimensions: (ncol,nlayers,ngptsw)
140 : real(kind=r8), intent(in) :: ptaormc(ncol,nlayers,ngptsw) ! cloud optical depth, non-delta scaled [mcica]
141 : ! Dimensions: (ncol,nlayers,ngptsw)
142 : real(kind=r8), intent(in) :: ptaua(ncol,nlayers,nbndsw) ! aerosol optical depth
143 : ! Dimensions: (ncol,nlayers,nbndsw)
144 : real(kind=r8), intent(in) :: pasya(ncol,nlayers,nbndsw) ! aerosol asymmetry parameter
145 : ! Dimensions: (ncol,nlayers,nbndsw)
146 : real(kind=r8), intent(in) :: pomga(ncol,nlayers,nbndsw) ! aerosol single scattering albedo
147 : ! Dimensions: (ncol,nlayers,nbndsw)
148 :
149 : real(kind=r8), intent(in) :: colh2o(ncol,nlayers)
150 : ! Dimensions: (ncol,nlayers)
151 : real(kind=r8), intent(in) :: colco2(ncol,nlayers)
152 : ! Dimensions: (ncol,nlayers)
153 : real(kind=r8), intent(in) :: colch4(ncol,nlayers)
154 : ! Dimensions: (ncol,nlayers)
155 : real(kind=r8), intent(in) :: co2mult(ncol,nlayers)
156 : ! Dimensions: (ncol,nlayers)
157 : real(kind=r8), intent(in) :: colo3(ncol,nlayers)
158 : ! Dimensions: (ncol,nlayers)
159 : real(kind=r8), intent(in) :: colo2(ncol,nlayers)
160 : ! Dimensions: (ncol,nlayers)
161 : real(kind=r8), intent(in) :: coln2o(ncol,nlayers)
162 : ! Dimensions: (ncol,nlayers)
163 : real(kind=r8), intent(in) :: forfac(ncol,nlayers)
164 : ! Dimensions: (ncol,nlayers)
165 : real(kind=r8), intent(in) :: forfrac(ncol,nlayers)
166 : ! Dimensions: (ncol,nlayers)
167 : real(kind=r8), intent(in) :: selffac(ncol,nlayers)
168 : ! Dimensions: (ncol,nlayers)
169 : real(kind=r8), intent(in) :: selffrac(ncol,nlayers)
170 : ! Dimensions: (ncol,nlayers)
171 : real(kind=r8), intent(in) :: fac00(ncol,nlayers)
172 : ! Dimensions: (ncol,nlayers)
173 : real(kind=r8), intent(in) :: fac01(ncol,nlayers)
174 : ! Dimensions: (ncol,nlayers)
175 : real(kind=r8), intent(in) :: fac10(ncol,nlayers)
176 : ! Dimensions: (ncol,nlayers)
177 : real(kind=r8), intent(in) :: fac11(ncol,nlayers)
178 : ! Dimensions: (ncol,nlayers)
179 :
180 : ! ------- Output -------
181 : ! All Dimensions: (nlayers+1)
182 : real(kind=r8), intent(out) :: pbbcd(:,:)
183 : real(kind=r8), intent(out) :: pbbcu(:,:)
184 : real(kind=r8), intent(out) :: pbbfd(:,:)
185 : real(kind=r8), intent(out) :: pbbfu(:,:)
186 : real(kind=r8), intent(out) :: pbbfddir(ncol,nlayers+2)
187 : real(kind=r8), intent(out) :: pbbcddir(ncol,nlayers+2)
188 :
189 : real(kind=r8), intent(out) :: puvcd(ncol,nlayers+2)
190 : real(kind=r8), intent(out) :: puvfd(ncol,nlayers+2)
191 : real(kind=r8), intent(out) :: puvcddir(ncol,nlayers+2)
192 : real(kind=r8), intent(out) :: puvfddir(:,:)
193 :
194 : real(kind=r8), intent(out) :: pnicd(ncol,nlayers+2)
195 : real(kind=r8), intent(out) :: pnifd(ncol,nlayers+2)
196 : real(kind=r8), intent(out) :: pnicddir(ncol,nlayers+2)
197 : real(kind=r8), intent(out) :: pnifddir(:,:)
198 :
199 : ! Added for net near-IR flux diagnostic
200 : real(kind=r8), intent(out) :: pnicu(ncol,nlayers+2)
201 : real(kind=r8), intent(out) :: pnifu(ncol,nlayers+2)
202 :
203 : ! Output - inactive ! All Dimensions: (nlayers+1)
204 : ! real(kind=r8), intent(out) :: puvcu(:)
205 : ! real(kind=r8), intent(out) :: puvfu(:)
206 : ! real(kind=r8), intent(out) :: pvscd(:)
207 : ! real(kind=r8), intent(out) :: pvscu(:)
208 : ! real(kind=r8), intent(out) :: pvsfd(:)
209 : ! real(kind=r8), intent(out) :: pvsfu(:)
210 :
211 : real(kind=r8), intent(out) :: pbbfsu(:,:,:) ! shortwave spectral flux up (nswbands,nlayers+1)
212 : real(kind=r8), intent(out) :: pbbfsd(:,:,:) ! shortwave spectral flux down (nswbands,nlayers+1)
213 :
214 :
215 : ! ------- Local -------
216 :
217 76800 : logical :: lrtchkclr(ncol,nlayers),lrtchkcld(ncol,nlayers)
218 :
219 : integer :: klev
220 : integer :: ib1, ib2, ibm, igt, ikl, ikp, ikx
221 : integer :: jb, jg, jl, jk
222 76800 : integer :: iw(ncol), iplon
223 : ! integer, parameter :: nuv = ??
224 : ! integer, parameter :: nvs = ??
225 76800 : integer :: itind(ncol)
226 :
227 76800 : real(kind=r8) :: tblind(ncol), ze1(ncol)
228 76800 : real(kind=r8) :: zclear(ncol), zcloud(ncol)
229 76800 : real(kind=r8) :: zdbt(ncol,nlayers+1), zdbt_nodel(ncol,nlayers+1)
230 76800 : real(kind=r8) :: zgcc(ncol,nlayers), zgco(ncol,nlayers)
231 76800 : real(kind=r8) :: zomcc(ncol,nlayers), zomco(ncol,nlayers)
232 76800 : real(kind=r8) :: zrdnd(ncol,nlayers+1), zrdndc(ncol,nlayers+1)
233 76800 : real(kind=r8) :: zref(ncol,nlayers+1), zrefc(ncol,nlayers+1), zrefo(ncol,nlayers+1)
234 76800 : real(kind=r8) :: zrefd(ncol,nlayers+1), zrefdc(ncol,nlayers+1), zrefdo(ncol,nlayers+1)
235 76800 : real(kind=r8) :: zrup(ncol,nlayers+1), zrupd(ncol,nlayers+1)
236 76800 : real(kind=r8) :: zrupc(ncol,nlayers+1), zrupdc(ncol,nlayers+1)
237 76800 : real(kind=r8) :: ztauc(ncol,nlayers), ztauo(ncol,nlayers)
238 76800 : real(kind=r8) :: ztdbt(ncol,nlayers+1)
239 76800 : real(kind=r8) :: ztra(ncol,nlayers+1), ztrac(ncol,nlayers+1), ztrao(ncol,nlayers+1)
240 76800 : real(kind=r8) :: ztrad(ncol,nlayers+1), ztradc(ncol,nlayers+1), ztrado(ncol,nlayers+1)
241 76800 : real(kind=r8) :: zdbtc(ncol,nlayers+1), ztdbtc(ncol,nlayers+1)
242 76800 : real(kind=r8) :: zincflx(ncol,ngptsw), zdbtc_nodel(ncol,nlayers+1)
243 76800 : real(kind=r8) :: ztdbt_nodel(ncol,nlayers+1), ztdbtc_nodel(ncol,nlayers+1)
244 :
245 76800 : real(kind=r8) :: zdbtmc(ncol), zdbtmo(ncol), zf
246 76800 : real(kind=r8) :: zwf, tauorig(ncol), repclc
247 : ! real(kind=r8) :: zincflux ! inactive
248 :
249 : ! Arrays from rrtmg_sw_taumoln routines
250 :
251 : ! real(kind=r8) :: ztaug(nlayers,16), ztaur(nlayers,16)
252 : ! real(kind=r8) :: zsflxzen(16)
253 76800 : real(kind=r8) :: ztaug(ncol,nlayers,ngptsw), ztaur(ncol,nlayers,ngptsw)
254 76800 : real(kind=r8) :: zsflxzen(ncol,ngptsw)
255 :
256 : ! Arrays from rrtmg_sw_vrtqdr routine
257 :
258 76800 : real(kind=r8) :: zcd(ncol,nlayers+1,ngptsw), zcu(ncol,nlayers+1,ngptsw)
259 76800 : real(kind=r8) :: zfd(ncol,nlayers+1,ngptsw), zfu(ncol,nlayers+1,ngptsw)
260 :
261 : ! Inactive arrays
262 : ! real(kind=r8) :: zbbcd(nlayers+1), zbbcu(nlayers+1)
263 : ! real(kind=r8) :: zbbfd(nlayers+1), zbbfu(nlayers+1)
264 : ! real(kind=r8) :: zbbfddir(nlayers+1), zbbcddir(nlayers+1)
265 : ! ------------------------------------------------------------------
266 :
267 : ! Initializations
268 :
269 38400 : ib1 = istart
270 38400 : ib2 = iend
271 38400 : klev = nlayers
272 :
273 :
274 : ! zincflux = 0.0_r8
275 :
276 38400 : repclc = 1.e-12_r8
277 10782720 : pbbcd(1:ncol,1:klev+1)=0._r8
278 10782720 : pbbcu(1:ncol,1:klev+1)=0._r8
279 10782720 : pbbfd(1:ncol,1:klev+1)=0._r8
280 10782720 : pbbfu(1:ncol,1:klev+1)=0._r8
281 10782720 : pbbcddir(1:ncol,1:klev+1)=0._r8
282 10744320 : pbbfddir(1:ncol,1:klev+1)=0._r8
283 10744320 : puvcd(1:ncol,1:klev+1)=0._r8
284 10744320 : puvfd(1:ncol,1:klev+1)=0._r8
285 10744320 : puvcddir(1:ncol,1:klev+1)=0._r8
286 10782720 : puvfddir(1:ncol,1:klev+1)=0._r8
287 10744320 : pnicd(1:ncol,1:klev+1)=0._r8
288 10744320 : pnifd(1:ncol,1:klev+1)=0._r8
289 10744320 : pnicddir(1:ncol,1:klev+1)=0._r8
290 10782720 : pnifddir(1:ncol,1:klev+1)=0._r8
291 10744320 : pnicu(1:ncol,1:klev+1)=0._r8
292 10744320 : pnifu(1:ncol,1:klev+1)=0._r8
293 142387200 : pbbfsu(:,1:ncol,1:klev+1)= 0._r8
294 142387200 : pbbfsd(:,1:ncol,1:klev+1)=0._r8
295 :
296 :
297 : ! Calculate the optical depths for gaseous absorption and Rayleigh scattering
298 :
299 314880 : do iplon=1,ncol
300 : call taumol_sw(klev, &
301 : colh2o(iplon,:), colco2(iplon,:), colch4(iplon,:),&
302 : colo2(iplon,:), colo3(iplon,:), colmol(iplon,:), &
303 276480 : laytrop(iplon), jp(iplon,:), jt(iplon,:), jt1(iplon,:), &
304 : fac00(iplon,:), fac01(iplon,:), fac10(iplon,:),&
305 : fac11(iplon,:), &
306 : selffac(iplon,:), selffrac(iplon,:),&
307 : indself(iplon,:), forfac(iplon,:), forfrac(iplon,:), indfor(iplon,:), &
308 : zsflxzen(iplon,:), ztaug(iplon,:,:),&
309 591360 : ztaur(iplon,:,:))
310 : enddo
311 :
312 : ! Top of shortwave spectral band loop, jb = 16 -> 29; ibm = 1 -> 14
313 :
314 314880 : jb = ib1-1 ! ???
315 314880 : do iplon=1,ncol
316 314880 : iw(iplon) =0
317 : end do
318 314880 : do iplon=1,ncol
319 : ! Clear-sky
320 : ! TOA direct beam
321 : ! Surface values
322 276480 : ztdbtc(iplon,1)=1.0_r8
323 276480 : ztdbtc_nodel(iplon,1)=1.0_r8
324 276480 : zdbtc(iplon,klev+1) =0.0_r8
325 276480 : ztrac(iplon,klev+1) =0.0_r8
326 276480 : ztradc(iplon,klev+1)=0.0_r8
327 : ! Cloudy-sky
328 : ! Surface values
329 276480 : ztrao(iplon,klev+1) =0.0_r8
330 276480 : ztrado(iplon,klev+1)=0.0_r8
331 : ! Total sky
332 : ! TOA direct beam
333 276480 : ztdbt(iplon,1)=1.0_r8
334 276480 : ztdbt_nodel(iplon,1)=1.0_r8
335 : ! Surface values
336 276480 : zdbt(iplon,klev+1) =0.0_r8
337 276480 : ztra(iplon,klev+1) =0.0_r8
338 314880 : ztrad(iplon,klev+1)=0.0_r8
339 : enddo
340 576000 : do jb = ib1, ib2
341 537600 : ibm = jb-15
342 537600 : igt = ngc(ibm)
343 :
344 : ! Reinitialize g-point counter for each band if output for each band is requested.
345 4408320 : do iplon=1,ncol
346 4408320 : if (iout.gt.0.and.ibm.ge.2) iw(iplon)= ngs(ibm-1)
347 : enddo
348 :
349 :
350 : ! do jk=1,klev+1
351 : ! zbbcd(jk)=0.0_r8
352 : ! zbbcu(jk)=0.0_r8
353 : ! zbbfd(jk)=0.0_r8
354 : ! zbbfu(jk)=0.0_r8
355 : ! enddo
356 :
357 : ! Top of g-point interval loop within each band (iw is cumulative counter)
358 4876800 : do jg = 1,igt
359 35266560 : do iplon=1,ncol
360 30965760 : iw(iplon) = iw(iplon)+1
361 :
362 : ! Apply adjustment for correct Earth/Sun distance and zenith angle to incoming solar flux
363 35266560 : zincflx(iplon,iw(iplon)) = adjflux(iplon,jb) * zsflxzen(iplon,iw(iplon)) * prmu0(iplon)
364 : ! zincflux = zincflux + adjflux(jb) * zsflxzen(iw) * prmu0 ! inactive
365 : enddo
366 :
367 : ! Compute layer reflectances and transmittances for direct and diffuse sources,
368 : ! first clear then cloudy
369 :
370 : ! zrefc(jk) direct albedo for clear
371 : ! zrefo(jk) direct albedo for cloud
372 : ! zrefdc(jk) diffuse albedo for clear
373 : ! zrefdo(jk) diffuse albedo for cloud
374 : ! ztrac(jk) direct transmittance for clear
375 : ! ztrao(jk) direct transmittance for cloudy
376 : ! ztradc(jk) diffuse transmittance for clear
377 : ! ztrado(jk) diffuse transmittance for cloudy
378 : !
379 : ! zref(jk) direct reflectance
380 : ! zrefd(jk) diffuse reflectance
381 : ! ztra(jk) direct transmittance
382 : ! ztrad(jk) diffuse transmittance
383 : !
384 : ! zdbtc(jk) clear direct beam transmittance
385 : ! zdbto(jk) cloudy direct beam transmittance
386 : ! zdbt(jk) layer mean direct beam transmittance
387 : ! ztdbt(jk) total direct beam transmittance at levels
388 35266560 : do iplon=1,ncol
389 30965760 : zrefc(iplon,klev+1) =palbp(iplon,ibm)
390 30965760 : zrefdc(iplon,klev+1)=palbd(iplon,ibm)
391 30965760 : zrupc(iplon,klev+1) =palbp(iplon,ibm)
392 30965760 : zrupdc(iplon,klev+1)=palbd(iplon,ibm)
393 30965760 : zrefo(iplon,klev+1) =palbp(iplon,ibm)
394 30965760 : zrefdo(iplon,klev+1)=palbd(iplon,ibm)
395 30965760 : zref(iplon,klev+1) =palbp(iplon,ibm)
396 30965760 : zrefd(iplon,klev+1)=palbd(iplon,ibm)
397 30965760 : zrup(iplon,klev+1) =palbp(iplon,ibm)
398 35266560 : zrupd(iplon,klev+1)=palbd(iplon,ibm)
399 : enddo
400 : ! Top of layer loop
401 146227200 : do jk=1,klev
402 141926400 : ikl=klev+1-jk
403 1163796480 : do iplon=1,ncol
404 : ! Note: two-stream calculations proceed from top to bottom;
405 : ! RRTMG_SW quantities are given bottom to top and are reversed here
406 :
407 :
408 : ! Set logical flag to do REFTRA calculation
409 : ! Do REFTRA for all clear layers
410 1021870080 : lrtchkclr(iplon,jk)=.true.
411 :
412 : ! Do REFTRA only for cloudy layers in profile, since already done for clear layers
413 1021870080 : lrtchkcld(iplon,jk)=.false.
414 1021870080 : lrtchkcld(iplon,jk)=(pcldfmc(iplon,ikl,iw(iplon)) > repclc)
415 :
416 : ! Clear-sky optical parameters - this section inactive
417 : ! Original
418 : ! ztauc(jk) = ztaur(ikl,iw) + ztaug(ikl,iw)
419 : ! zomcc(jk) = ztaur(ikl,iw) / ztauc(jk)
420 : ! zgcc(jk) = 0.0001_r8
421 : ! Total sky optical parameters
422 : ! ztauo(jk) = ztaur(ikl,iw) + ztaug(ikl,iw) + ptaucmc(ikl,iw)
423 : ! zomco(jk) = ptaucmc(ikl,iw) * pomgcmc(ikl,iw) + ztaur(ikl,iw)
424 : ! zgco (jk) = (ptaucmc(ikl,iw) * pomgcmc(ikl,iw) * pasycmc(ikl,iw) + &
425 : ! ztaur(ikl,iw) * 0.0001_r8) / zomco(jk)
426 : ! zomco(jk) = zomco(jk) / ztauo(jk)
427 :
428 : ! Clear-sky optical parameters including aerosols
429 1021870080 : if (ztaug(iplon,ikl,iw(iplon)) .lt. 0.0_r8) ztaug(iplon,ikl,iw(iplon)) = 0.0_r8
430 :
431 1021870080 : ztauc(iplon,jk) = ztaur(iplon,ikl,iw(iplon)) + ztaug(iplon,ikl,iw(iplon)) + ptaua(iplon,ikl,ibm)
432 1021870080 : zomcc(iplon,jk) = ztaur(iplon,ikl,iw(iplon)) * 1.0_r8 + ptaua(iplon,ikl,ibm) * pomga(iplon,ikl,ibm)
433 1021870080 : zgcc(iplon,jk) = pasya(iplon,ikl,ibm) * pomga(iplon,ikl,ibm) * ptaua(iplon,ikl,ibm) / zomcc(iplon,jk)
434 1163796480 : zomcc(iplon,jk) = zomcc(iplon,jk) / ztauc(iplon,jk)
435 :
436 : ! Pre-delta-scaling clear and cloudy direct beam transmittance (must use 'orig', unscaled cloud OD)
437 : ! \/\/\/ This block of code is only needed for unscaled direct beam calculation
438 : enddo
439 141926400 : if (idelm .eq. 0) then
440 0 : do iplon=1,ncol
441 : !
442 0 : zclear(iplon) = 1.0_r8 - pcldfmc(iplon,ikl,iw(iplon))
443 0 : zcloud(iplon) = pcldfmc(iplon,ikl,iw(iplon))
444 :
445 : ! Clear
446 : ! zdbtmc = exp(-ztauc(jk) / prmu0)
447 :
448 : ! Use exponential lookup table for transmittance, or expansion of exponential for low tau
449 0 : ze1(iplon) = ztauc(iplon,jk) / prmu0(iplon)
450 : enddo
451 0 : do iplon=1,ncol
452 0 : if (ze1(iplon) .le. od_lo) then
453 0 : zdbtmc(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
454 : else
455 0 : tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
456 0 : itind(iplon) = tblint * tblind(iplon) + 0.5_r8
457 0 : zdbtmc(iplon) = exp_tbl(itind(iplon))
458 : endif
459 : enddo
460 0 : do iplon=1,ncol
461 :
462 0 : zdbtc_nodel(iplon,jk) = zdbtmc(iplon)
463 0 : ztdbtc_nodel(iplon,jk+1) = zdbtc_nodel(iplon,jk) * ztdbtc_nodel(iplon,jk)
464 :
465 : ! Clear + Cloud
466 0 : tauorig(iplon) = ztauc(iplon,jk) + ptaormc(iplon,ikl,iw(iplon))
467 : ! zdbtmo = exp(-tauorig / prmu0)
468 :
469 : ! Use exponential lookup table for transmittance, or expansion of exponential for low tau
470 0 : ze1(iplon) = tauorig(iplon) / prmu0(iplon)
471 : enddo
472 0 : do iplon=1,ncol
473 0 : if (ze1(iplon) .le. od_lo) then
474 0 : zdbtmo(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
475 : else
476 0 : tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
477 0 : itind(iplon) = tblint * tblind(iplon) + 0.5_r8
478 0 : zdbtmo(iplon) = exp_tbl(itind(iplon))
479 : endif
480 : enddo
481 0 : do iplon=1,ncol
482 :
483 0 : zdbt_nodel(iplon,jk) = zclear(iplon)*zdbtmc(iplon) + zcloud(iplon)*zdbtmo(iplon)
484 0 : ztdbt_nodel(iplon,jk+1) = zdbt_nodel(iplon,jk) * ztdbt_nodel(iplon,jk)
485 :
486 : enddo
487 : endif
488 1163796480 : do iplon=1,ncol
489 : ! /\/\/\ Above code only needed for unscaled direct beam calculation
490 :
491 :
492 : ! Delta scaling - clear
493 1021870080 : zf = zgcc(iplon,jk) * zgcc(iplon,jk)
494 1021870080 : zwf = zomcc(iplon,jk) * zf
495 1021870080 : ztauc(iplon,jk) = (1.0_r8 - zwf) * ztauc(iplon,jk)
496 1021870080 : zomcc(iplon,jk) = (zomcc(iplon,jk) - zwf) / (1.0_r8 - zwf)
497 1163796480 : zgcc (iplon,jk) = (zgcc(iplon,jk) - zf) / (1.0_r8 - zf)
498 : enddo
499 : ! Total sky optical parameters (cloud properties already delta-scaled)
500 : ! Use this code if cloud properties are derived in rrtmg_sw_cldprop
501 146227200 : if (icpr .ge. 1) then
502 1163796480 : do iplon=1,ncol
503 1021870080 : ztauo(iplon,jk) = ztauc(iplon,jk) + ptaucmc(iplon,ikl,iw(iplon))
504 1021870080 : zomco(iplon,jk) = ztauc(iplon,jk) * zomcc(iplon,jk) + ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon))
505 1021870080 : zgco (iplon,jk) = (ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon)) * pasycmc(iplon,ikl,iw(iplon)) + &
506 1021870080 : ztauc(iplon,jk) * zomcc(iplon,jk) * zgcc(iplon,jk)) / zomco(iplon,jk)
507 1163796480 : zomco(iplon,jk) = zomco(iplon,jk) / ztauo(iplon,jk)
508 : enddo
509 : ! Total sky optical parameters (if cloud properties not delta scaled)
510 : ! Use this code if cloud properties are not derived in rrtmg_sw_cldprop
511 0 : elseif (icpr .eq. 0) then
512 0 : do iplon=1,ncol
513 0 : ztauo(iplon,jk) = ztaur(iplon,ikl,iw(iplon)) + ztaug(iplon,ikl,iw(iplon)) + ptaua(iplon,ikl,ibm) + ptaucmc(iplon,ikl,iw(iplon))
514 0 : zomco(iplon,jk) = ptaua(iplon,ikl,ibm) * pomga(iplon,ikl,ibm) + ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon)) + &
515 0 : ztaur(iplon,ikl,iw(iplon)) * 1.0_r8
516 0 : zgco (iplon,jk) = (ptaucmc(iplon,ikl,iw(iplon)) * pomgcmc(iplon,ikl,iw(iplon)) * pasycmc(iplon,ikl,iw(iplon)) + &
517 0 : ptaua(iplon,ikl,ibm)*pomga(iplon,ikl,ibm)*pasya(iplon,ikl,ibm)) / zomco(iplon,jk)
518 0 : zomco(iplon,jk) = zomco(iplon,jk) / ztauo(iplon,jk)
519 :
520 : ! Delta scaling - clouds
521 : ! Use only if subroutine rrtmg_sw_cldprop is not used to get cloud properties and to apply delta scaling
522 0 : zf = zgco(iplon,jk) * zgco(iplon,jk)
523 0 : zwf = zomco(iplon,jk) * zf
524 0 : ztauo(iplon,jk) = (1._r8 - zwf) * ztauo(iplon,jk)
525 0 : zomco(iplon,jk) = (zomco(iplon,jk) - zwf) / (1.0_r8 - zwf)
526 0 : zgco (iplon,jk) = (zgco(iplon,jk) - zf) / (1.0_r8 - zf)
527 : enddo
528 : endif
529 :
530 : ! End of layer loop
531 : enddo
532 :
533 : ! Clear sky reflectivities
534 : call reftra_sw (klev,ncol, &
535 : lrtchkclr, zgcc, prmu0, ztauc, zomcc, &
536 4300800 : zrefc, zrefdc, ztrac, ztradc)
537 :
538 : ! Total sky reflectivities
539 : call reftra_sw (klev, ncol, &
540 : lrtchkcld, zgco, prmu0, ztauo, zomco, &
541 4300800 : zrefo, zrefdo, ztrao, ztrado)
542 :
543 146227200 : do jk=1,klev
544 1163796480 : do iplon=1,ncol
545 : ! Combine clear and cloudy contributions for total sky
546 1021870080 : ikl = klev+1-jk
547 1021870080 : zclear(iplon) = 1.0_r8 - pcldfmc(iplon,ikl,iw(iplon))
548 1021870080 : zcloud(iplon) = pcldfmc(iplon,ikl,iw(iplon))
549 :
550 1021870080 : zref(iplon,jk) = zclear(iplon)*zrefc(iplon,jk) + zcloud(iplon)*zrefo(iplon,jk)
551 1021870080 : zrefd(iplon,jk)= zclear(iplon)*zrefdc(iplon,jk) + zcloud(iplon)*zrefdo(iplon,jk)
552 1021870080 : ztra(iplon,jk) = zclear(iplon)*ztrac(iplon,jk) + zcloud(iplon)*ztrao(iplon,jk)
553 1021870080 : ztrad(iplon,jk)= zclear(iplon)*ztradc(iplon,jk) + zcloud(iplon)*ztrado(iplon,jk)
554 :
555 : ! Direct beam transmittance
556 :
557 : ! Clear
558 : ! zdbtmc(iplon) = exp(-ztauc(iplon,jk) / prmu0)
559 :
560 : ! Use exponential lookup table for transmittance, or expansion of
561 : ! exponential for low tau
562 1163796480 : ze1(iplon) = ztauc(iplon,jk) / prmu0(iplon)
563 : enddo
564 1168097280 : do iplon=1,ncol
565 1021870080 : if (ze1(iplon) .le. od_lo) then
566 555632885 : zdbtmc(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
567 : else
568 466237195 : tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
569 466237195 : itind(iplon) = tblint * tblind(iplon) + 0.5_r8
570 466237195 : zdbtmc(iplon) = exp_tbl(itind(iplon))
571 : endif
572 :
573 1021870080 : zdbtc(iplon,jk) = zdbtmc(iplon)
574 1021870080 : ztdbtc(iplon,jk+1) = zdbtc(iplon,jk)*ztdbtc(iplon,jk)
575 :
576 : ! Clear + Cloud
577 : ! zdbtmo = exp(-ztauo(jk) / prmu0)
578 :
579 : ! Use exponential lookup table for transmittance, or expansion of
580 : ! exponential for low tau
581 1021870080 : ze1(iplon) = ztauo(iplon,jk) / prmu0(iplon)
582 1021870080 : if (ze1(iplon) .le. od_lo) then
583 511977374 : zdbtmo(iplon) = 1._r8 - ze1(iplon) + 0.5_r8 * ze1(iplon) * ze1(iplon)
584 : else
585 509892706 : tblind(iplon) = ze1(iplon) / (bpade + ze1(iplon))
586 509892706 : itind(iplon) = tblint * tblind(iplon) + 0.5_r8
587 509892706 : zdbtmo(iplon) = exp_tbl(itind(iplon))
588 : endif
589 :
590 1021870080 : zdbt(iplon,jk) = zclear(iplon)*zdbtmc(iplon) + zcloud(iplon)*zdbtmo(iplon)
591 1163796480 : ztdbt(iplon,jk+1) = zdbt(iplon,jk)*ztdbt(iplon,jk)
592 :
593 : enddo
594 : enddo
595 : ! Vertical quadrature for clear-sky fluxes
596 :
597 : call vrtqdr_sw(ncol,klev, iw, &
598 : zrefc, zrefdc, ztrac, ztradc, &
599 : zdbtc, zrdndc, zrupc, zrupdc, ztdbtc, &
600 4300800 : zcd, zcu)
601 :
602 : ! Vertical quadrature for cloudy fluxes
603 :
604 : call vrtqdr_sw(ncol,klev, iw, &
605 : zref, zrefd, ztra, ztrad, &
606 : zdbt, zrdnd, zrup, zrupd, ztdbt, &
607 4300800 : zfd, zfu)
608 :
609 : ! Upwelling and downwelling fluxes at levels
610 : ! Two-stream calculations go from top to bottom;
611 : ! layer indexing is reversed to go bottom to top for output arrays
612 :
613 151065600 : do jk=1,klev+1
614 1199063040 : do iplon=1,ncol
615 1052835840 : ikl=klev+2-jk
616 :
617 : ! Accumulate spectral fluxes over bands - inactive
618 : ! zbbfu(ikl) = zbbfu(ikl) + zincflx(iw)*zfu(jk,iw)
619 : ! zbbfd(ikl) = zbbfd(ikl) + zincflx(iw)*zfd(jk,iw)
620 : ! zbbcu(ikl) = zbbcu(ikl) + zincflx(iw)*zcu(jk,iw)
621 : ! zbbcd(ikl) = zbbcd(ikl) + zincflx(iw)*zcd(jk,iw)
622 : ! zbbfddir(ikl) = zbbfddir(ikl) + zincflx(iw)*ztdbt_nodel(jk)
623 : ! zbbcddir(ikl) = zbbcddir(ikl) + zincflx(iw)*ztdbtc_nodel(jk)
624 :
625 1052835840 : pbbfsu(ibm,iplon,ikl) = pbbfsu(ibm,iplon,ikl) + zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
626 1052835840 : pbbfsd(ibm,iplon,ikl) = pbbfsd(ibm,iplon,ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
627 :
628 : ! Accumulate spectral fluxes over whole spectrum
629 1052835840 : pbbfu(iplon,ikl) = pbbfu(iplon,ikl) + zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
630 1052835840 : pbbfd(iplon,ikl) = pbbfd(iplon,ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
631 1052835840 : pbbcu(iplon,ikl) = pbbcu(iplon,ikl) + zincflx(iplon,iw(iplon))*zcu(iplon,jk,iw(iplon))
632 1199063040 : pbbcd(iplon,ikl) = pbbcd(iplon,ikl) + zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
633 : enddo
634 146227200 : if (idelm .eq. 0) then
635 0 : do iplon=1,ncol
636 0 : pbbfddir(iplon,ikl) = pbbfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
637 0 : pbbcddir(iplon,ikl) = pbbcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
638 : enddo
639 146227200 : elseif (idelm .eq. 1) then
640 1199063040 : do iplon=1,ncol
641 1052835840 : pbbfddir(iplon,ikl) = pbbfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
642 1199063040 : pbbcddir(iplon,ikl) = pbbcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
643 : enddo
644 : endif
645 :
646 : ! Accumulate direct fluxes for UV/visible bands
647 150528000 : if (ibm >= 10 .and. ibm <= 13) then
648 278353920 : do iplon=1,ncol
649 244408320 : puvcd(iplon,ikl) = puvcd(iplon,ikl) + zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
650 278353920 : puvfd(iplon,ikl) = puvfd(iplon,ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
651 : enddo
652 33945600 : if (idelm .eq. 0) then
653 0 : do iplon=1,ncol
654 0 : puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
655 0 : puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
656 : enddo
657 33945600 : elseif (idelm .eq. 1) then
658 278353920 : do iplon=1,ncol
659 244408320 : puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
660 278353920 : puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
661 : enddo
662 : endif
663 : ! band 9 is half-NearIR and half-Visible
664 112281600 : else if (ibm == 9) then
665 85647360 : do iplon=1,ncol
666 75202560 : puvcd(iplon,ikl) = puvcd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
667 75202560 : puvfd(iplon,ikl) = puvfd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
668 75202560 : pnicd(iplon,ikl) = pnicd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
669 85647360 : pnifd(iplon,ikl) = pnifd(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
670 : enddo
671 10444800 : if (idelm .eq. 0) then
672 0 : do iplon=1,ncol
673 0 : puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
674 0 : puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
675 0 : pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
676 0 : pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
677 : enddo
678 10444800 : elseif (idelm .eq. 1) then
679 85647360 : do iplon=1,ncol
680 75202560 : puvfddir(iplon,ikl) = puvfddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
681 75202560 : puvcddir(iplon,ikl) = puvcddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
682 75202560 : pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
683 85647360 : pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
684 : enddo
685 : endif
686 85647360 : do iplon=1,ncol
687 75202560 : pnicu(iplon,ikl) = pnicu(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zcu(iplon,jk,iw(iplon))
688 85647360 : pnifu(iplon,ikl) = pnifu(iplon,ikl) + 0.5_r8*zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
689 : ! Accumulate direct fluxes for near-IR bands
690 : enddo
691 101836800 : else if (ibm == 14 .or. ibm <= 8) then
692 835061760 : do iplon=1,ncol
693 733224960 : pnicd(iplon,ikl) = pnicd(iplon,ikl) + zincflx(iplon,iw(iplon))*zcd(iplon,jk,iw(iplon))
694 835061760 : pnifd(iplon,ikl) = pnifd(iplon,ikl) + zincflx(iplon,iw(iplon))*zfd(iplon,jk,iw(iplon))
695 : enddo
696 101836800 : if (idelm .eq. 0) then
697 0 : do iplon=1,ncol
698 0 : pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt_nodel(iplon,jk)
699 0 : pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc_nodel(iplon,jk)
700 : enddo
701 101836800 : elseif (idelm .eq. 1) then
702 835061760 : do iplon=1,ncol
703 733224960 : pnifddir(iplon,ikl) = pnifddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbt(iplon,jk)
704 835061760 : pnicddir(iplon,ikl) = pnicddir(iplon,ikl) + zincflx(iplon,iw(iplon))*ztdbtc(iplon,jk)
705 : enddo
706 : endif
707 : ! Added for net near-IR flux diagnostic
708 835061760 : do iplon=1,ncol
709 733224960 : pnicu(iplon,ikl) = pnicu(iplon,ikl) + zincflx(iplon,iw(iplon))*zcu(iplon,jk,iw(iplon))
710 835061760 : pnifu(iplon,ikl) = pnifu(iplon,ikl) + zincflx(iplon,iw(iplon))*zfu(iplon,jk,iw(iplon))
711 : enddo
712 : endif
713 :
714 :
715 : ! End loop on jg, g-point interval
716 : enddo
717 :
718 : ! End loop on jb, spectral band
719 : enddo
720 : end do
721 :
722 38400 : end subroutine spcvmc_sw
723 :
724 : end module rrtmg_sw_spcvmc
725 :
726 :
|
