source: palm/trunk/LIB/rrtmg/rrtmg_sw_spcvrt.f90 @ 1585

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

Added support for RRTMG radiation code

File size: 27.3 KB
Line 
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_spcvrt
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 : hvrspv, hnamspv
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 spcvrt_sw &
35            (nlayers, istart, iend, icpr, idelm, iout, &
36             pavel, tavel, pz, tz, tbound, palbd, palbp, &
37             pclfr, ptauc, pasyc, pomgc, ptaucorig, &
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.
49!
50! Interface:  *spcvrt_sw* is called from *rrtmg_sw.F90* or rrtmg_sw.1col.F90*
51!
52! Method:
53!    Adapted from two-stream model of H. Barker;
54!    Two-stream model options (selected with kmodts in rrtmg_sw_reftra.F90):
55!        1: Eddington, 2: PIFM, Zdunkowski et al., 3: discret ordinates
56!
57! Modifications:
58!
59! Original: H. Barker
60! Revision: Merge with RRTMG_SW: J.-J.Morcrette, ECMWF, Feb 2003
61! Revision: Add adjustment for Earth/Sun distance : MJIacono, AER, Oct 2003
62! Revision: Bug fix for use of PALBP and PALBD: MJIacono, AER, Nov 2003
63! Revision: Bug fix to apply delta scaling to clear sky: AER, Dec 2004
64! Revision: Code modified so that delta scaling is not done in cloudy profiles
65!           if routine cldprop is used; delta scaling can be applied by swithcing
66!           code below if cldprop is not used to get cloud properties.
67!           AER, Jan 2005
68! Revision: Uniform formatting for RRTMG: MJIacono, AER, Jul 2006
69! Revision: Use exponential lookup table for transmittance: MJIacono, AER,
70!           Aug 2007
71!
72! ------------------------------------------------------------------
73
74! ------- Declarations ------
75
76! -------- Input -------
77
78      integer(kind=im), intent(in) :: nlayers
79      integer(kind=im), intent(in) :: istart
80      integer(kind=im), intent(in) :: iend
81      integer(kind=im), intent(in) :: icpr
82      integer(kind=im), intent(in) :: idelm   ! delta-m scaling flag
83                                              ! [0 = direct and diffuse fluxes are unscaled]
84                                              ! [1 = direct and diffuse fluxes are scaled]
85      integer(kind=im), intent(in) :: iout
86      integer(kind=im), intent(in) :: laytrop
87      integer(kind=im), intent(in) :: layswtch
88      integer(kind=im), intent(in) :: laylow
89
90      integer(kind=im), intent(in) :: indfor(:)
91                                                               !   Dimensions: (nlayers)
92      integer(kind=im), intent(in) :: indself(:)
93                                                               !   Dimensions: (nlayers)
94      integer(kind=im), intent(in) :: jp(:)
95                                                               !   Dimensions: (nlayers)
96      integer(kind=im), intent(in) :: jt(:)
97                                                               !   Dimensions: (nlayers)
98      integer(kind=im), intent(in) :: jt1(:)
99                                                               !   Dimensions: (nlayers)
100
101      real(kind=rb), intent(in) :: pavel(:)                    ! layer pressure (hPa, mb)
102                                                               !   Dimensions: (nlayers)
103      real(kind=rb), intent(in) :: tavel(:)                    ! layer temperature (K)
104                                                               !   Dimensions: (nlayers)
105      real(kind=rb), intent(in) :: pz(0:)                      ! level (interface) pressure (hPa, mb)
106                                                               !   Dimensions: (0:nlayers)
107      real(kind=rb), intent(in) :: tz(0:)                      ! level temperatures (hPa, mb)
108                                                               !   Dimensions: (0:nlayers)
109      real(kind=rb), intent(in) :: tbound                      ! surface temperature (K)
110      real(kind=rb), intent(in) :: wkl(:,:)                    ! molecular amounts (mol/cm2)
111                                                               !   Dimensions: (mxmol,nlayers)
112      real(kind=rb), intent(in) :: coldry(:)                   ! dry air column density (mol/cm2)
113                                                               !   Dimensions: (nlayers)
114      real(kind=rb), intent(in) :: colmol(:)
115                                                               !   Dimensions: (nlayers)
116      real(kind=rb), intent(in) :: adjflux(:)                  ! Earth/Sun distance adjustment
117                                                               !   Dimensions: (jpband)
118
119      real(kind=rb), intent(in) :: palbd(:)                    ! surface albedo (diffuse)
120                                                               !   Dimensions: (nbndsw)
121      real(kind=rb), intent(in) :: palbp(:)                    ! surface albedo (direct)
122                                                               !   Dimensions: (nbndsw)
123      real(kind=rb), intent(in) :: prmu0                       ! cosine of solar zenith angle
124      real(kind=rb), intent(in) :: pclfr(:)                    ! cloud fraction
125                                                               !   Dimensions: (nlayers)
126      real(kind=rb), intent(in) :: ptauc(:,:)                  ! cloud optical depth
127                                                               !   Dimensions: (nlayers,nbndsw)
128      real(kind=rb), intent(in) :: pasyc(:,:)                  ! cloud asymmetry parameter
129                                                               !   Dimensions: (nlayers,nbndsw)
130      real(kind=rb), intent(in) :: pomgc(:,:)                  ! cloud single scattering albedo
131                                                               !   Dimensions: (nlayers,nbndsw)
132      real(kind=rb), intent(in) :: ptaucorig(:,:)              ! cloud optical depth, non-delta scaled
133                                                               !   Dimensions: (nlayers,nbndsw)
134      real(kind=rb), intent(in) :: ptaua(:,:)                  ! aerosol optical depth
135                                                               !   Dimensions: (nlayers,nbndsw)
136      real(kind=rb), intent(in) :: pasya(:,:)                  ! aerosol asymmetry parameter
137                                                               !   Dimensions: (nlayers,nbndsw)
138      real(kind=rb), intent(in) :: pomga(:,:)                  ! aerosol single scattering albedo
139                                                               !   Dimensions: (nlayers,nbndsw)
140
141      real(kind=rb), intent(in) :: colh2o(:)
142                                                               !   Dimensions: (nlayers)
143      real(kind=rb), intent(in) :: colco2(:)
144                                                               !   Dimensions: (nlayers)
145      real(kind=rb), intent(in) :: colch4(:)
146                                                               !   Dimensions: (nlayers)
147      real(kind=rb), intent(in) :: co2mult(:)
148                                                               !   Dimensions: (nlayers)
149      real(kind=rb), intent(in) :: colo3(:)
150                                                               !   Dimensions: (nlayers)
151      real(kind=rb), intent(in) :: colo2(:)
152                                                               !   Dimensions: (nlayers)
153      real(kind=rb), intent(in) :: coln2o(:)
154                                                               !   Dimensions: (nlayers)
155
156      real(kind=rb), intent(in) :: forfac(:)
157                                                               !   Dimensions: (nlayers)
158      real(kind=rb), intent(in) :: forfrac(:)
159                                                               !   Dimensions: (nlayers)
160      real(kind=rb), intent(in) :: selffac(:)
161                                                               !   Dimensions: (nlayers)
162      real(kind=rb), intent(in) :: selffrac(:)
163                                                               !   Dimensions: (nlayers)
164      real(kind=rb), intent(in) :: fac00(:)
165                                                               !   Dimensions: (nlayers)
166      real(kind=rb), intent(in) :: fac01(:)
167                                                               !   Dimensions: (nlayers)
168      real(kind=rb), intent(in) :: fac10(:)
169                                                               !   Dimensions: (nlayers)
170      real(kind=rb), intent(in) :: fac11(:)
171                                                               !   Dimensions: (nlayers)
172
173! ------- Output -------
174                                                               !   All Dimensions: (nlayers+1)
175      real(kind=rb), intent(out) :: pbbcd(:)
176      real(kind=rb), intent(out) :: pbbcu(:)
177      real(kind=rb), intent(out) :: pbbfd(:)
178      real(kind=rb), intent(out) :: pbbfu(:)
179      real(kind=rb), intent(out) :: pbbfddir(:)
180      real(kind=rb), intent(out) :: pbbcddir(:)
181
182      real(kind=rb), intent(out) :: puvcd(:)
183      real(kind=rb), intent(out) :: puvfd(:)
184      real(kind=rb), intent(out) :: puvcddir(:)
185      real(kind=rb), intent(out) :: puvfddir(:)
186
187      real(kind=rb), intent(out) :: pnicd(:)
188      real(kind=rb), intent(out) :: pnifd(:)
189      real(kind=rb), intent(out) :: pnicddir(:)
190      real(kind=rb), intent(out) :: pnifddir(:)
191
192! Output - inactive                                            !   All Dimensions: (nlayers+1)
193!      real(kind=rb), intent(out) :: puvcu(:)
194!      real(kind=rb), intent(out) :: puvfu(:)
195!      real(kind=rb), intent(out) :: pnicu(:)
196!      real(kind=rb), intent(out) :: pnifu(:)
197!      real(kind=rb), intent(out) :: pvscd(:)
198!      real(kind=rb), intent(out) :: pvscu(:)
199!      real(kind=rb), intent(out) :: pvsfd(:)
200!      real(kind=rb), intent(out) :: pvsfu(:)
201
202
203! ------- Local -------
204
205      logical :: lrtchkclr(nlayers),lrtchkcld(nlayers)
206
207      integer(kind=im)  :: klev
208      integer(kind=im) :: ib1, ib2, ibm, igt, ikl, ikp, ikx
209      integer(kind=im) :: iw, jb, jg, jl, jk
210!      integer(kind=im), parameter :: nuv = ??
211!      integer(kind=im), parameter :: nvs = ??
212      integer(kind=im) :: itind
213
214      real(kind=rb) :: tblind, ze1
215      real(kind=rb) :: zclear, zcloud
216      real(kind=rb) :: zdbt(nlayers+1), zdbt_nodel(nlayers+1)
217      real(kind=rb) :: zgc(nlayers), zgcc(nlayers), zgco(nlayers)
218      real(kind=rb) :: zomc(nlayers), zomcc(nlayers), zomco(nlayers)
219      real(kind=rb) :: zrdnd(nlayers+1), zrdndc(nlayers+1)
220      real(kind=rb) :: zref(nlayers+1), zrefc(nlayers+1), zrefo(nlayers+1)
221      real(kind=rb) :: zrefd(nlayers+1), zrefdc(nlayers+1), zrefdo(nlayers+1)
222      real(kind=rb) :: zrup(nlayers+1), zrupd(nlayers+1)
223      real(kind=rb) :: zrupc(nlayers+1), zrupdc(nlayers+1)
224      real(kind=rb) :: zs1(nlayers+1)
225      real(kind=rb) :: ztauc(nlayers), ztauo(nlayers)
226      real(kind=rb) :: ztdn(nlayers+1), ztdnd(nlayers+1), ztdbt(nlayers+1)
227      real(kind=rb) :: ztoc(nlayers), ztor(nlayers)
228      real(kind=rb) :: ztra(nlayers+1), ztrac(nlayers+1), ztrao(nlayers+1)
229      real(kind=rb) :: ztrad(nlayers+1), ztradc(nlayers+1), ztrado(nlayers+1)
230      real(kind=rb) :: zdbtc(nlayers+1), ztdbtc(nlayers+1)
231      real(kind=rb) :: zincflx(ngptsw), zdbtc_nodel(nlayers+1) 
232      real(kind=rb) :: ztdbt_nodel(nlayers+1), ztdbtc_nodel(nlayers+1)
233
234      real(kind=rb) :: zdbtmc, zdbtmo, zf, zgw, zreflect
235      real(kind=rb) :: zwf, tauorig, repclc
236!     real(kind=rb) :: zincflux                                   ! inactive
237
238! Arrays from rrtmg_sw_taumoln routines
239
240!      real(kind=rb) :: ztaug(nlayers,16), ztaur(nlayers,16)
241!      real(kind=rb) :: zsflxzen(16)
242      real(kind=rb) :: ztaug(nlayers,ngptsw), ztaur(nlayers,ngptsw)
243      real(kind=rb) :: zsflxzen(ngptsw)
244
245! Arrays from rrtmg_sw_vrtqdr routine
246
247      real(kind=rb) :: zcd(nlayers+1,ngptsw), zcu(nlayers+1,ngptsw)
248      real(kind=rb) :: zfd(nlayers+1,ngptsw), zfu(nlayers+1,ngptsw)
249
250! Inactive arrays
251!     real(kind=rb) :: zbbcd(nlayers+1), zbbcu(nlayers+1)
252!     real(kind=rb) :: zbbfd(nlayers+1), zbbfu(nlayers+1)
253!     real(kind=rb) :: zbbfddir(nlayers+1), zbbcddir(nlayers+1)
254
255! ------------------------------------------------------------------
256
257! Initializations
258
259      ib1 = istart
260      ib2 = iend
261      klev = nlayers
262      iw = 0
263      repclc = 1.e-12_rb
264!      zincflux = 0.0_rb
265
266      do jk=1,klev+1
267         pbbcd(jk)=0._rb
268         pbbcu(jk)=0._rb
269         pbbfd(jk)=0._rb
270         pbbfu(jk)=0._rb
271         pbbcddir(jk)=0._rb
272         pbbfddir(jk)=0._rb
273         puvcd(jk)=0._rb
274         puvfd(jk)=0._rb
275         puvcddir(jk)=0._rb
276         puvfddir(jk)=0._rb
277         pnicd(jk)=0._rb
278         pnifd(jk)=0._rb
279         pnicddir(jk)=0._rb
280         pnifddir(jk)=0._rb
281      enddo
282
283
284! Calculate the optical depths for gaseous absorption and Rayleigh scattering
285
286      call taumol_sw(klev, &
287                     colh2o, colco2, colch4, colo2, colo3, colmol, &
288                     laytrop, jp, jt, jt1, &
289                     fac00, fac01, fac10, fac11, &
290                     selffac, selffrac, indself, forfac, forfrac, indfor, &
291                     zsflxzen, ztaug, ztaur)
292
293
294! Top of shortwave spectral band loop, jb = 16 -> 29; ibm = 1 -> 14
295
296      do jb = ib1, ib2
297         ibm = jb-15
298         igt = ngc(ibm)
299
300! Reinitialize g-point counter for each band if output for each band is requested.
301         if (iout.gt.0.and.ibm.ge.2) iw = ngs(ibm-1)
302
303!        do jk=1,klev+1
304!           zbbcd(jk)=0.0_rb
305!           zbbcu(jk)=0.0_rb
306!           zbbfd(jk)=0.0_rb
307!           zbbfu(jk)=0.0_rb
308!        enddo
309
310! Top of g-point interval loop within each band (iw is cumulative counter)
311         do jg = 1,igt
312            iw = iw+1
313
314! Apply adjustments for correct Earth/Sun distance and zenith angle to incoming solar flux
315            zincflx(iw) = adjflux(jb) * zsflxzen(iw) * prmu0
316!             zincflux = zincflux + adjflux(jb) * zsflxzen(iw) * prmu0           ! inactive
317
318! Compute layer reflectances and transmittances for direct and diffuse sources,
319! first clear then cloudy
320
321! zrefc(jk)  direct albedo for clear
322! zrefo(jk)  direct albedo for cloud
323! zrefdc(jk) diffuse albedo for clear
324! zrefdo(jk) diffuse albedo for cloud
325! ztrac(jk)  direct transmittance for clear
326! ztrao(jk)  direct transmittance for cloudy
327! ztradc(jk) diffuse transmittance for clear
328! ztrado(jk) diffuse transmittance for cloudy
329
330! zref(jk)   direct reflectance
331! zrefd(jk)  diffuse reflectance
332! ztra(jk)   direct transmittance
333! ztrad(jk)  diffuse transmittance
334!
335! zdbtc(jk)  clear direct beam transmittance
336! zdbto(jk)  cloudy direct beam transmittance
337! zdbt(jk)   layer mean direct beam transmittance
338! ztdbt(jk)  total direct beam transmittance at levels
339
340! Clear-sky   
341!   TOA direct beam   
342            ztdbtc(1)=1.0_rb
343            ztdbtc_nodel(1)=1.0_rb
344!   Surface values
345            zdbtc(klev+1) =0.0_rb
346            ztrac(klev+1) =0.0_rb
347            ztradc(klev+1)=0.0_rb
348            zrefc(klev+1) =palbp(ibm)
349            zrefdc(klev+1)=palbd(ibm)
350            zrupc(klev+1) =palbp(ibm)
351            zrupdc(klev+1)=palbd(ibm)
352           
353! Cloudy-sky   
354!   Surface values
355            ztrao(klev+1) =0.0_rb
356            ztrado(klev+1)=0.0_rb
357            zrefo(klev+1) =palbp(ibm)
358            zrefdo(klev+1)=palbd(ibm)
359           
360! Total sky   
361!   TOA direct beam   
362            ztdbt(1)=1.0_rb
363            ztdbt_nodel(1)=1.0_rb
364!   Surface values
365            zdbt(klev+1) =0.0_rb
366            ztra(klev+1) =0.0_rb
367            ztrad(klev+1)=0.0_rb
368            zref(klev+1) =palbp(ibm)
369            zrefd(klev+1)=palbd(ibm)
370            zrup(klev+1) =palbp(ibm)
371            zrupd(klev+1)=palbd(ibm)
372   
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)=(pclfr(ikl) > 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) + ptauc(ikl,ibm)
397!               zomco(jk) = ptauc(ikl,ibm) * pomgc(ikl,ibm) + ztaur(ikl,iw)
398!               zgco (jk) = (ptauc(ikl,ibm) * pomgc(ikl,ibm) * pasyc(ikl,ibm) + &
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 - pclfr(ikl)
413                  zcloud = pclfr(ikl)
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) + ptaucorig(ikl,ibm)
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) + ptauc(ikl,ibm)
463                  zomco(jk) = ztauc(jk) * zomcc(jk) + ptauc(ikl,ibm) * pomgc(ikl,ibm) 
464                  zgco (jk) = (ptauc(ikl,ibm) * pomgc(ikl,ibm) * pasyc(ikl,ibm) + &
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) + ptauc(ikl,ibm)
472                  zomco(jk) = ptaua(ikl,ibm) * pomga(ikl,ibm) + ptauc(ikl,ibm) * pomgc(ikl,ibm) + &
473                              ztaur(ikl,iw) * 1.0_rb
474                  zgco (jk) = (ptauc(ikl,ibm) * pomgc(ikl,ibm) * pasyc(ikl,ibm) + &
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
491! Clear sky reflectivities
492            call reftra_sw (klev, &
493                            lrtchkclr, zgcc, prmu0, ztauc, zomcc, &
494                            zrefc, zrefdc, ztrac, ztradc)
495
496! Total sky reflectivities     
497            call reftra_sw (klev, &
498                            lrtchkcld, zgco, prmu0, ztauo, zomco, &
499                            zrefo, zrefdo, ztrao, ztrado)
500
501
502            do jk=1,klev
503
504! Combine clear and cloudy contributions for total sky
505               ikl = klev+1-jk
506               zclear = 1.0_rb - pclfr(ikl)
507               zcloud = pclfr(ikl)
508
509               zref(jk) = zclear*zrefc(jk) + zcloud*zrefo(jk)
510               zrefd(jk)= zclear*zrefdc(jk) + zcloud*zrefdo(jk)
511               ztra(jk) = zclear*ztrac(jk) + zcloud*ztrao(jk)
512               ztrad(jk)= zclear*ztradc(jk) + zcloud*ztrado(jk)
513
514! Direct beam transmittance       
515
516! Clear
517!                zdbtmc = exp(-ztauc(jk) / prmu0)
518
519! Use exponential lookup table for transmittance, or expansion of
520! exponential for low tau
521               ze1 = ztauc(jk) / prmu0
522               if (ze1 .le. od_lo) then
523                  zdbtmc = 1._rb - ze1 + 0.5_rb * ze1 * ze1
524               else
525                  tblind = ze1 / (bpade + ze1)
526                  itind = tblint * tblind + 0.5_rb
527                  zdbtmc = exp_tbl(itind)
528               endif
529
530               zdbtc(jk) = zdbtmc
531               ztdbtc(jk+1) = zdbtc(jk)*ztdbtc(jk)
532
533! Clear + Cloud
534!                zdbtmo = exp(-ztauo(jk) / prmu0)
535
536! Use exponential lookup table for transmittance, or expansion of
537! exponential for low tau
538               ze1 = ztauo(jk) / prmu0
539               if (ze1 .le. od_lo) then
540                  zdbtmo = 1._rb - ze1 + 0.5_rb * ze1 * ze1
541               else
542                  tblind = ze1 / (bpade + ze1)
543                  itind = tblint * tblind + 0.5_rb
544                  zdbtmo = exp_tbl(itind)
545               endif
546
547               zdbt(jk) = zclear*zdbtmc + zcloud*zdbtmo
548               ztdbt(jk+1) = zdbt(jk)*ztdbt(jk)
549       
550            enddo           
551                 
552! Vertical quadrature for clear-sky fluxes
553
554            call vrtqdr_sw (klev, iw, &
555                            zrefc, zrefdc, ztrac, ztradc, &
556                            zdbtc, zrdndc, zrupc, zrupdc, ztdbtc, &
557                            zcd, zcu)
558     
559! Vertical quadrature for cloudy fluxes
560
561            call vrtqdr_sw (klev, iw, &
562                            zref, zrefd, ztra, ztrad, &
563                            zdbt, zrdnd, zrup, zrupd, ztdbt, &
564                            zfd, zfu)
565
566! Upwelling and downwelling fluxes at levels
567!   Two-stream calculations go from top to bottom;
568!   layer indexing is reversed to go bottom to top for output arrays
569
570            do jk=1,klev+1
571               ikl=klev+2-jk
572
573! Accumulate spectral fluxes over bands - inactive
574!               zbbfu(ikl) = zbbfu(ikl) + zincflx(iw)*zfu(jk,iw) 
575!               zbbfd(ikl) = zbbfd(ikl) + zincflx(iw)*zfd(jk,iw)
576!               zbbcu(ikl) = zbbcu(ikl) + zincflx(iw)*zcu(jk,iw)
577!               zbbcd(ikl) = zbbcd(ikl) + zincflx(iw)*zcd(jk,iw)
578!               zbbfddir(ikl) = zbbfddir(ikl) + zincflx(iw)*ztdbt_nodel(jk)
579!               zbbcddir(ikl) = zbbcddir(ikl) + zincflx(iw)*ztdbtc_nodel(jk)
580
581! Accumulate spectral fluxes over whole spectrum 
582               pbbfu(ikl) = pbbfu(ikl) + zincflx(iw)*zfu(jk,iw)
583               pbbfd(ikl) = pbbfd(ikl) + zincflx(iw)*zfd(jk,iw)
584               pbbcu(ikl) = pbbcu(ikl) + zincflx(iw)*zcu(jk,iw)
585               pbbcd(ikl) = pbbcd(ikl) + zincflx(iw)*zcd(jk,iw)
586               if (idelm .eq. 0) then
587                  pbbfddir(ikl) = pbbfddir(ikl) + zincflx(iw)*ztdbt_nodel(jk)
588                  pbbcddir(ikl) = pbbcddir(ikl) + zincflx(iw)*ztdbtc_nodel(jk)
589               elseif (idelm .eq. 1) then
590                  pbbfddir(ikl) = pbbfddir(ikl) + zincflx(iw)*ztdbt(jk)
591                  pbbcddir(ikl) = pbbcddir(ikl) + zincflx(iw)*ztdbtc(jk)
592               endif
593
594! Accumulate direct fluxes for UV/visible bands
595               if (ibm >= 10 .and. ibm <= 13) then
596                  puvcd(ikl) = puvcd(ikl) + zincflx(iw)*zcd(jk,iw)
597                  puvfd(ikl) = puvfd(ikl) + zincflx(iw)*zfd(jk,iw)
598                  if (idelm .eq. 0) then
599                     puvfddir(ikl) = puvfddir(ikl) + zincflx(iw)*ztdbt_nodel(jk)
600                     puvcddir(ikl) = puvcddir(ikl) + zincflx(iw)*ztdbtc_nodel(jk)
601                  elseif (idelm .eq. 1) then
602                     puvfddir(ikl) = puvfddir(ikl) + zincflx(iw)*ztdbt(jk)
603                     puvcddir(ikl) = puvcddir(ikl) + zincflx(iw)*ztdbtc(jk)
604                  endif
605! Accumulate direct fluxes for near-IR bands
606               else if (ibm == 14 .or. ibm <= 9) then 
607                  pnicd(ikl) = pnicd(ikl) + zincflx(iw)*zcd(jk,iw)
608                  pnifd(ikl) = pnifd(ikl) + zincflx(iw)*zfd(jk,iw)
609                  if (idelm .eq. 0) then
610                     pnifddir(ikl) = pnifddir(ikl) + zincflx(iw)*ztdbt_nodel(jk)
611                     pnicddir(ikl) = pnicddir(ikl) + zincflx(iw)*ztdbtc_nodel(jk)
612                  elseif (idelm .eq. 1) then
613                     pnifddir(ikl) = pnifddir(ikl) + zincflx(iw)*ztdbt(jk)
614                     pnicddir(ikl) = pnicddir(ikl) + zincflx(iw)*ztdbtc(jk)
615                  endif
616               endif
617
618            enddo
619
620! End loop on jg, g-point interval
621         enddo             
622
623! End loop on jb, spectral band
624      enddo                   
625
626      end subroutine spcvrt_sw
627
628      end module rrtmg_sw_spcvrt
629
630
Note: See TracBrowser for help on using the repository browser.