source: palm/trunk/LIB/rrtmg/rrtmg_sw_spcvmc.f90 @ 3999

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

Added support for RRTMG radiation code

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