source: palm/trunk/LIB/rrtmg/rrtmg_sw_rad.nomcica.f90 @ 3775

Last change on this file since 3775 was 3272, checked in by suehring, 6 years ago

Merge with branch radition - improved representation of diffuse radiation

File size: 57.3 KB
Line 
1!     path:      $Source$
2!     author:    $Author: miacono $
3!     revision:  $Revision: 23308 $
4!     created:   $Date: 2013-12-27 17:23:51 -0500 (Fri, 27 Dec 2013) $
5!
6
7       module rrtmg_sw_rad
8
9!  --------------------------------------------------------------------------
10! |                                                                          |
11! |  Copyright 2002-2009, Atmospheric & Environmental Research, Inc. (AER).  |
12! |  This software may be used, copied, or redistributed as long as it is    |
13! |  not sold and this copyright notice is reproduced on each copy made.     |
14! |  This model is provided as is without any express or implied warranties. |
15! |                       (http://www.rtweb.aer.com/)                        |
16! |                                                                          |
17!  --------------------------------------------------------------------------
18!
19! ****************************************************************************
20! *                                                                          *
21! *                             RRTMG_SW                                     *
22! *                                                                          *
23! *                                                                          *
24! *                                                                          *
25! *                 a rapid radiative transfer model                         *
26! *                  for the solar spectral region                           *
27! *           for application to general circulation models                  *
28! *                                                                          *
29! *                                                                          *
30! *           Atmospheric and Environmental Research, Inc.                   *
31! *                       131 Hartwell Avenue                                *
32! *                       Lexington, MA 02421                                *
33! *                                                                          *
34! *                                                                          *
35! *                          Eli J. Mlawer                                   *
36! *                       Jennifer S. Delamere                               *
37! *                        Michael J. Iacono                                 *
38! *                        Shepard A. Clough                                 *
39! *                                                                          *
40! *                                                                          *
41! *                                                                          *
42! *                                                                          *
43! *                                                                          *
44! *                                                                          *
45! *                      email:  emlawer@aer.com                             *
46! *                      email:  jdelamer@aer.com                            *
47! *                      email:  miacono@aer.com                             *
48! *                                                                          *
49! *       The authors wish to acknowledge the contributions of the           *
50! *       following people:  Steven J. Taubman, Patrick D. Brown,            *
51! *       Ronald E. Farren, Luke Chen, Robert Bergstrom.                     *
52! *                                                                          *
53! ****************************************************************************
54
55! --------- Modules ---------
56      use parkind, only : im => kind_im, rb => kind_rb
57      use rrsw_vsn
58      use rrtmg_sw_cldprop, only: cldprop_sw
59! *** Move the required call to rrtmg_sw_ini below and the following
60! use association to GCM initialization area ***
61!      use rrtmg_sw_init, only: rrtmg_sw_ini
62      use rrtmg_sw_setcoef, only: setcoef_sw
63      use rrtmg_sw_spcvrt, only: spcvrt_sw
64
65      implicit none
66
67! public interfaces/functions/subroutines
68      public :: rrtmg_sw, inatm_sw, earth_sun
69
70!------------------------------------------------------------------
71      contains
72!------------------------------------------------------------------
73
74!------------------------------------------------------------------
75! Public subroutines
76!------------------------------------------------------------------
77
78      subroutine rrtmg_sw &
79            (ncol    ,nlay    ,icld    ,iaer    , &
80             play    ,plev    ,tlay    ,tlev    , &
81             tsfc    ,h2ovmr  ,o3vmr   ,co2vmr  , &
82             ch4vmr  ,n2ovmr  ,o2vmr   ,asdir   , &
83             asdif   ,aldir   ,aldif   ,coszen  , &
84             adjes   ,dyofyr  ,scon    ,inflgsw , &
85             iceflgsw,liqflgsw,cldfr   ,taucld  , &
86             ssacld  ,asmcld  ,fsfcld  ,cicewp  , &
87             cliqwp  ,reice   ,reliq   ,tauaer  , &
88             ssaaer  ,asmaer  ,ecaer   ,swuflx  , &
89             swdflx  ,swhr    ,swuflxc ,swdflxc , &
90             swhrc   ,dirdflux,difdflux)
91
92! ------- Description -------
93
94! This program is the driver for RRTMG_SW, the AER SW radiation model for
95!  application to GCMs, that has been adapted from RRTM_SW for improved
96!  efficiency and to provide fractional cloudiness and cloud overlap
97!  capability using McICA.
98!
99! Note: The call to RRTMG_SW_INI should be moved to the GCM initialization
100!  area, since this has to be called only once.
101!
102! This routine
103!    b) calls INATM_SW to read in the atmospheric profile from GCM;
104!       all layering in RRTMG is ordered from surface to toa.
105!    c) calls CLDPROP_SW to set cloud optical depth based on input
106!       cloud properties
107!    d) calls SETCOEF_SW to calculate various quantities needed for
108!       the radiative transfer algorithm
109!    e) calls SPCVRT to call the two-stream model that in turn
110!       calls TAUMOL to calculate gaseous optical depths for each
111!       of the 16 spectral bands and to perform the radiative transfer;
112!    f) passes the calculated fluxes and cooling rates back to GCM
113!
114! Two modes of operation are possible:
115!     The mode is chosen by using either rrtmg_sw.nomcica.f90 (to not use
116!     McICA) or rrtmg_sw.f90 (to use McICA) to interface with a GCM.
117!
118!    1) Standard, single forward model calculation (imca = 0); this is
119!       valid only for clear sky or fully overcast clouds
120!    2) Monte Carlo Independent Column Approximation (McICA, Pincus et al.,
121!       JC, 2003) method is applied to the forward model calculation (imca = 1)
122!       This method is valid for clear sky and full or partial cloud conditions.
123!
124! Two methods of cloud property input are possible:
125!     Cloud properties can be input in one of two ways (controlled by input
126!     flags inflag, iceflag and liqflag; see text file rrtmg_sw_instructions
127!     and subroutine rrtmg_sw_cldprop.f90 for further details):
128!
129!    1) Input cloud fraction, cloud optical depth, single scattering albedo
130!       and asymmetry parameter directly (inflgsw = 0)
131!    2) Input cloud fraction and cloud physical properties: ice fracion,
132!       ice and liquid particle sizes (inflgsw = 1 or 2); 
133!       cloud optical properties are calculated by cldprop or cldprmc based
134!       on input settings of iceflgsw and liqflgsw
135!
136! Two methods of aerosol property input are possible:
137!     Aerosol properties can be input in one of two ways (controlled by input
138!     flag iaer, see text file rrtmg_sw_instructions for further details):
139!
140!    1) Input aerosol optical depth, single scattering albedo and asymmetry
141!       parameter directly by layer and spectral band (iaer=10)
142!    2) Input aerosol optical depth and 0.55 micron directly by layer and use
143!       one or more of six ECMWF aerosol types (iaer=6)
144!
145!
146! ------- Modifications -------
147!
148! This version of RRTMG_SW has been modified from RRTM_SW to use a reduced
149! set of g-point intervals and a two-stream model for application to GCMs.
150!
151!-- Original version (derived from RRTM_SW)
152!     2002: AER. Inc.
153!-- Conversion to F90 formatting; addition of 2-stream radiative transfer
154!     Feb 2003: J.-J. Morcrette, ECMWF
155!-- Additional modifications for GCM application
156!     Aug 2003: M. J. Iacono, AER Inc.
157!-- Total number of g-points reduced from 224 to 112.  Original
158!   set of 224 can be restored by exchanging code in module parrrsw.f90
159!   and in file rrtmg_sw_init.f90.
160!     Apr 2004: M. J. Iacono, AER, Inc.
161!-- Modifications to include output for direct and diffuse
162!   downward fluxes.  There are output as "true" fluxes without
163!   any delta scaling applied.  Code can be commented to exclude
164!   this calculation in source file rrtmg_sw_spcvrt.f90.
165!     Jan 2005: E. J. Mlawer, M. J. Iacono, AER, Inc.
166!-- Reformatted for consistency with rrtmg_lw.
167!     Feb 2007: M. J. Iacono, AER, Inc.
168!-- Modifications to formatting to use assumed-shape arrays.
169!     Aug 2007: M. J. Iacono, AER, Inc.
170!-- Modified to output direct and diffuse fluxes either with or without
171!   delta scaling based on setting of idelm flag.
172!     Dec 2008: M. J. Iacono, AER, Inc.
173!-- Added the output direct and diffuse fluxes to the argument of rrtmg_sw
174!     Aug 2018: M. Salim, Humboldt Uni, Berlin, Germany
175
176! --------- Modules ---------
177
178      use parrrsw, only : nbndsw, ngptsw, naerec, nstr, nmol, mxmol, &
179                          jpband, jpb1, jpb2
180      use rrsw_aer, only : rsrtaua, rsrpiza, rsrasya
181      use rrsw_con, only : heatfac, oneminus, pi
182      use rrsw_wvn, only : wavenum1, wavenum2
183
184! ------- Declarations
185
186! ----- Input -----
187! Note: All volume mixing ratios are in dimensionless units of mole fraction obtained
188! by scaling mass mixing ratio (g/g) with the appropriate molecular weights (g/mol)
189      integer(kind=im), intent(in) :: ncol            ! Number of horizontal columns     
190      integer(kind=im), intent(in) :: nlay            ! Number of model layers
191      integer(kind=im), intent(inout) :: icld         ! Cloud overlap method
192                                                      !    0: Clear only
193                                                      !    1: Random
194                                                      !    2: Maximum/random
195                                                      !    3: Maximum
196      integer(kind=im), intent(inout) :: iaer         ! Aerosol option flag
197                                                      !    0: No aerosol
198                                                      !    6: ECMWF method
199                                                      !    10:Input aerosol optical
200                                                      !       properties
201
202      real(kind=rb), intent(in) :: play(:,:)          ! Layer pressures (hPa, mb)
203                                                      !    Dimensions: (ncol,nlay)
204      real(kind=rb), intent(in) :: plev(:,:)          ! Interface pressures (hPa, mb)
205                                                      !    Dimensions: (ncol,nlay+1)
206      real(kind=rb), intent(in) :: tlay(:,:)          ! Layer temperatures (K)
207                                                      !    Dimensions: (ncol,nlay)
208      real(kind=rb), intent(in) :: tlev(:,:)          ! Interface temperatures (K)
209                                                      !    Dimensions: (ncol,nlay+1)
210      real(kind=rb), intent(in) :: tsfc(:)            ! Surface temperature (K)
211                                                      !    Dimensions: (ncol)
212      real(kind=rb), intent(in) :: h2ovmr(:,:)        ! H2O volume mixing ratio
213                                                      !    Dimensions: (ncol,nlay)
214      real(kind=rb), intent(in) :: o3vmr(:,:)         ! O3 volume mixing ratio
215                                                      !    Dimensions: (ncol,nlay)
216      real(kind=rb), intent(in) :: co2vmr(:,:)        ! CO2 volume mixing ratio
217                                                      !    Dimensions: (ncol,nlay)
218      real(kind=rb), intent(in) :: ch4vmr(:,:)        ! Methane volume mixing ratio
219                                                      !    Dimensions: (ncol,nlay)
220      real(kind=rb), intent(in) :: n2ovmr(:,:)        ! Nitrous oxide volume mixing ratio
221                                                      !    Dimensions: (ncol,nlay)
222      real(kind=rb), intent(in) :: o2vmr(:,:)         ! Oxygen volume mixing ratio
223                                                      !    Dimensions: (ncol,nlay)
224      real(kind=rb), intent(in) :: asdir(:)           ! UV/vis surface albedo direct rad
225                                                      !    Dimensions: (ncol)
226      real(kind=rb), intent(in) :: aldir(:)           ! Near-IR surface albedo direct rad
227                                                      !    Dimensions: (ncol)
228      real(kind=rb), intent(in) :: asdif(:)           ! UV/vis surface albedo: diffuse rad
229                                                      !    Dimensions: (ncol)
230      real(kind=rb), intent(in) :: aldif(:)           ! Near-IR surface albedo: diffuse rad
231                                                      !    Dimensions: (ncol)
232
233      integer(kind=im), intent(in) :: dyofyr          ! Day of the year (used to get Earth/Sun
234                                                      !  distance if adjflx not provided)
235      real(kind=rb), intent(in) :: adjes              ! Flux adjustment for Earth/Sun distance
236      real(kind=rb), intent(in) :: coszen(:)          ! Cosine of solar zenith angle
237                                                      !    Dimensions: (ncol)
238      real(kind=rb), intent(in) :: scon               ! Solar constant (W/m2)
239
240      integer(kind=im), intent(in) :: inflgsw         ! Flag for cloud optical properties
241      integer(kind=im), intent(in) :: iceflgsw        ! Flag for ice particle specification
242      integer(kind=im), intent(in) :: liqflgsw        ! Flag for liquid droplet specification
243
244      real(kind=rb), intent(in) :: cldfr(:,:)         ! Cloud fraction
245                                                      !    Dimensions: (ncol,nlay)
246      real(kind=rb), intent(in) :: taucld(:,:,:)      ! In-cloud optical depth
247                                                      !    Dimensions: (nbndsw,ncol,nlay)
248      real(kind=rb), intent(in) :: ssacld(:,:,:)      ! In-cloud single scattering albedo
249                                                      !    Dimensions: (nbndsw,ncol,nlay)
250      real(kind=rb), intent(in) :: asmcld(:,:,:)      ! In-cloud asymmetry parameter
251                                                      !    Dimensions: (nbndsw,ncol,nlay)
252      real(kind=rb), intent(in) :: fsfcld(:,:,:)      ! In-cloud forward scattering fraction
253                                                      !    Dimensions: (nbndsw,ncol,nlay)
254      real(kind=rb), intent(in) :: cicewp(:,:)        ! In-cloud ice water path (g/m2)
255                                                      !    Dimensions: (ncol,nlay)
256      real(kind=rb), intent(in) :: cliqwp(:,:)        ! In-cloud liquid water path (g/m2)
257                                                      !    Dimensions: (ncol,nlay)
258      real(kind=rb), intent(in) :: reice(:,:)         ! Cloud ice effective radius (microns)
259                                                      !    Dimensions: (ncol,nlay)
260                                                      ! specific definition of reice depends on setting of iceflgsw:
261                                                      ! iceflgsw = 0: (inactive)
262                                                      !             
263                                                      ! iceflgsw = 1: ice effective radius, r_ec, (Ebert and Curry, 1992),
264                                                      !               r_ec range is limited to 13.0 to 130.0 microns
265                                                      ! iceflgsw = 2: ice effective radius, r_k, (Key, Streamer Ref. Manual, 1996)
266                                                      !               r_k range is limited to 5.0 to 131.0 microns
267                                                      ! iceflgsw = 3: generalized effective size, dge, (Fu, 1996),
268                                                      !               dge range is limited to 5.0 to 140.0 microns
269                                                      !               [dge = 1.0315 * r_ec]
270      real(kind=rb), intent(in) :: reliq(:,:)         ! Cloud water drop effective radius (microns)
271                                                      !    Dimensions: (ncol,nlay)
272      real(kind=rb), intent(in) :: tauaer(:,:,:)      ! Aerosol optical depth (iaer=10 only)
273                                                      !    Dimensions: (ncol,nlay,nbndsw)
274                                                      ! (non-delta scaled)     
275      real(kind=rb), intent(in) :: ssaaer(:,:,:)      ! Aerosol single scattering albedo (iaer=10 only)
276                                                      !    Dimensions: (ncol,nlay,nbndsw)
277                                                      ! (non-delta scaled)     
278      real(kind=rb), intent(in) :: asmaer(:,:,:)      ! Aerosol asymmetry parameter (iaer=10 only)
279                                                      !    Dimensions: (ncol,nlay,nbndsw)
280                                                      ! (non-delta scaled)     
281      real(kind=rb), intent(in) :: ecaer(:,:,:)       ! Aerosol optical depth at 0.55 micron (iaer=6 only)
282                                                      !    Dimensions: (ncol,nlay,naerec)
283                                                      ! (non-delta scaled)     
284
285! ----- Output -----
286
287      real(kind=rb), intent(out) :: swuflx(:,:)       ! Total sky shortwave upward flux (W/m2)
288                                                      !    Dimensions: (ncol,nlay+1)
289      real(kind=rb), intent(out) :: swdflx(:,:)       ! Total sky shortwave downward flux (W/m2)
290                                                      !    Dimensions: (ncol,nlay+1)
291      real(kind=rb), intent(out) :: swhr(:,:)         ! Total sky shortwave radiative heating rate (K/d)
292                                                      !    Dimensions: (ncol,nlay)
293      real(kind=rb), intent(out) :: swuflxc(:,:)      ! Clear sky shortwave upward flux (W/m2)
294                                                      !    Dimensions: (ncol,nlay+1)
295      real(kind=rb), intent(out) :: swdflxc(:,:)      ! Clear sky shortwave downward flux (W/m2)
296                                                      !    Dimensions: (ncol,nlay+1)
297      real(kind=rb), intent(out) :: swhrc(:,:)        ! Clear sky shortwave radiative heating rate (K/d)
298                                                      !    Dimensions: (ncol,nlay)
299
300! ----- Local -----
301
302! Control
303      integer(kind=im) :: nlayers             ! total number of layers
304      integer(kind=im) :: istart              ! beginning band of calculation
305      integer(kind=im) :: iend                ! ending band of calculation
306      integer(kind=im) :: icpr                ! cldprop/cldprmc use flag
307      integer(kind=im) :: iout                ! output option flag
308      integer(kind=im) :: idelm               ! delta-m scaling flag
309                                              ! [0 = direct and diffuse fluxes are unscaled]
310                                              ! [1 = direct and diffuse fluxes are scaled]
311                                              ! (total downward fluxes are always delta scaled)
312      integer(kind=im) :: isccos              ! instrumental cosine response flag (inactive)
313      integer(kind=im) :: iplon               ! column loop index
314      integer(kind=im) :: i                   ! layer loop index                       ! jk
315      integer(kind=im) :: ib                  ! band loop index                        ! jsw
316      integer(kind=im) :: ia, ig              ! indices
317      integer(kind=im) :: k                   ! layer loop index
318      integer(kind=im) :: ims                 ! value for changing mcica permute seed
319      integer(kind=im) :: imca                ! flag for mcica [0=off, 1=on]
320
321      real(kind=rb) :: zepsec, zepzen         ! epsilon
322      real(kind=rb) :: zdpgcp                 ! flux to heating conversion ratio
323
324! Atmosphere
325      real(kind=rb) :: pavel(nlay+1)          ! layer pressures (mb)
326      real(kind=rb) :: tavel(nlay+1)          ! layer temperatures (K)
327      real(kind=rb) :: pz(0:nlay+1)           ! level (interface) pressures (hPa, mb)
328      real(kind=rb) :: tz(0:nlay+1)           ! level (interface) temperatures (K)
329      real(kind=rb) :: tbound                 ! surface temperature (K)
330      real(kind=rb) :: pdp(nlay+1)            ! layer pressure thickness (hPa, mb)
331      real(kind=rb) :: coldry(nlay+1)         ! dry air column amount
332      real(kind=rb) :: wkl(mxmol,nlay+1)      ! molecular amounts (mol/cm-2)
333
334!      real(kind=rb) :: earth_sun             ! function for Earth/Sun distance factor
335      real(kind=rb) :: cossza                 ! Cosine of solar zenith angle
336      real(kind=rb) :: adjflux(jpband)        ! adjustment for current Earth/Sun distance
337      real(kind=rb) :: solvar(jpband)         ! solar constant scaling factor from rrtmg_sw
338                                              !  default value of 1368.22 Wm-2 at 1 AU
339      real(kind=rb) :: albdir(nbndsw)         ! surface albedo, direct          ! zalbp
340      real(kind=rb) :: albdif(nbndsw)         ! surface albedo, diffuse         ! zalbd
341
342      real(kind=rb) :: taua(nlay+1,nbndsw)    ! Aerosol optical depth
343      real(kind=rb) :: ssaa(nlay+1,nbndsw)    ! Aerosol single scattering albedo
344      real(kind=rb) :: asma(nlay+1,nbndsw)    ! Aerosol asymmetry parameter
345
346! Atmosphere - setcoef
347      integer(kind=im) :: laytrop             ! tropopause layer index
348      integer(kind=im) :: layswtch            ! tropopause layer index
349      integer(kind=im) :: laylow              ! tropopause layer index
350      integer(kind=im) :: jp(nlay+1)          !
351      integer(kind=im) :: jt(nlay+1)          !
352      integer(kind=im) :: jt1(nlay+1)         !
353
354      real(kind=rb) :: colh2o(nlay+1)         ! column amount (h2o)
355      real(kind=rb) :: colco2(nlay+1)         ! column amount (co2)
356      real(kind=rb) :: colo3(nlay+1)          ! column amount (o3)
357      real(kind=rb) :: coln2o(nlay+1)         ! column amount (n2o)
358      real(kind=rb) :: colch4(nlay+1)         ! column amount (ch4)
359      real(kind=rb) :: colo2(nlay+1)          ! column amount (o2)
360      real(kind=rb) :: colmol(nlay+1)         ! column amount
361      real(kind=rb) :: co2mult(nlay+1)        ! column amount
362
363      integer(kind=im) :: indself(nlay+1)
364      integer(kind=im) :: indfor(nlay+1)
365      real(kind=rb) :: selffac(nlay+1)
366      real(kind=rb) :: selffrac(nlay+1)
367      real(kind=rb) :: forfac(nlay+1)
368      real(kind=rb) :: forfrac(nlay+1)
369
370      real(kind=rb) :: &                      !
371                         fac00(nlay+1), fac01(nlay+1), &
372                         fac10(nlay+1), fac11(nlay+1) 
373
374! Atmosphere/clouds - cldprop
375      integer(kind=im) :: ncbands             ! number of cloud spectral bands
376      integer(kind=im) :: inflag              ! flag for cloud property method
377      integer(kind=im) :: iceflag             ! flag for ice cloud properties
378      integer(kind=im) :: liqflag             ! flag for liquid cloud properties
379
380      real(kind=rb) :: cldfrac(nlay+1)        ! layer cloud fraction
381      real(kind=rb) :: tauc(nbndsw,nlay+1)    ! in-cloud optical depth (non-delta scaled)
382      real(kind=rb) :: ssac(nbndsw,nlay+1)    ! in-cloud single scattering albedo (non-delta scaled)
383      real(kind=rb) :: asmc(nbndsw,nlay+1)    ! in-cloud asymmetry parameter (non-delta scaled)
384      real(kind=rb) :: fsfc(nbndsw,nlay+1)    ! in-cloud forward scattering fraction (non-delta scaled)
385      real(kind=rb) :: ciwp(nlay+1)           ! in-cloud ice water path
386      real(kind=rb) :: clwp(nlay+1)           ! in-cloud liquid water path
387      real(kind=rb) :: rel(nlay+1)            ! cloud liquid particle effective radius (microns)
388      real(kind=rb) :: rei(nlay+1)            ! cloud ice particle effective size (microns)
389
390      real(kind=rb) :: taucloud(nlay+1,jpband)  ! in-cloud optical depth
391      real(kind=rb) :: taucldorig(nlay+1,jpband)! in-cloud optical depth (non-delta scaled)
392      real(kind=rb) :: ssacloud(nlay+1,jpband)  ! in-cloud single scattering albedo
393      real(kind=rb) :: asmcloud(nlay+1,jpband)  ! in-cloud asymmetry parameter
394
395! Atmosphere/clouds/aerosol - spcvrt,spcvmc
396      real(kind=rb) :: ztauc(nlay+1,nbndsw)     ! cloud optical depth
397      real(kind=rb) :: ztaucorig(nlay+1,nbndsw) ! unscaled cloud optical depth
398      real(kind=rb) :: zasyc(nlay+1,nbndsw)     ! cloud asymmetry parameter
399                                                !  (first moment of phase function)
400      real(kind=rb) :: zomgc(nlay+1,nbndsw)     ! cloud single scattering albedo
401      real(kind=rb) :: ztaua(nlay+1,nbndsw)     ! total aerosol optical depth
402      real(kind=rb) :: zasya(nlay+1,nbndsw)     ! total aerosol asymmetry parameter
403      real(kind=rb) :: zomga(nlay+1,nbndsw)     ! total aerosol single scattering albedo
404
405      real(kind=rb) :: zbbfu(nlay+2)          ! temporary upward shortwave flux (w/m2)
406      real(kind=rb) :: zbbfd(nlay+2)          ! temporary downward shortwave flux (w/m2)
407      real(kind=rb) :: zbbcu(nlay+2)          ! temporary clear sky upward shortwave flux (w/m2)
408      real(kind=rb) :: zbbcd(nlay+2)          ! temporary clear sky downward shortwave flux (w/m2)
409      real(kind=rb) :: zbbfddir(nlay+2)       ! temporary downward direct shortwave flux (w/m2)
410      real(kind=rb) :: zbbcddir(nlay+2)       ! temporary clear sky downward direct shortwave flux (w/m2)
411      real(kind=rb) :: zuvfd(nlay+2)          ! temporary UV downward shortwave flux (w/m2)
412      real(kind=rb) :: zuvcd(nlay+2)          ! temporary clear sky UV downward shortwave flux (w/m2)
413      real(kind=rb) :: zuvfddir(nlay+2)       ! temporary UV downward direct shortwave flux (w/m2)
414      real(kind=rb) :: zuvcddir(nlay+2)       ! temporary clear sky UV downward direct shortwave flux (w/m2)
415      real(kind=rb) :: znifd(nlay+2)          ! temporary near-IR downward shortwave flux (w/m2)
416      real(kind=rb) :: znicd(nlay+2)          ! temporary clear sky near-IR downward shortwave flux (w/m2)
417      real(kind=rb) :: znifddir(nlay+2)       ! temporary near-IR downward direct shortwave flux (w/m2)
418      real(kind=rb) :: znicddir(nlay+2)       ! temporary clear sky near-IR downward direct shortwave flux (w/m2)
419
420! Optional output fields
421      real(kind=rb) :: swnflx(nlay+2)         ! Total sky shortwave net flux (W/m2)
422      real(kind=rb) :: swnflxc(nlay+2)        ! Clear sky shortwave net flux (W/m2)
423      real(kind=rb), intent(out) :: dirdflux(nlay+2)       ! Direct downward shortwave surface flux
424      real(kind=rb), intent(out) :: difdflux(nlay+2)       ! Diffuse downward shortwave surface flux
425      real(kind=rb) :: uvdflx(nlay+2)         ! Total sky downward shortwave flux, UV/vis 
426      real(kind=rb) :: nidflx(nlay+2)         ! Total sky downward shortwave flux, near-IR
427      real(kind=rb) :: dirdnuv(nlay+2)        ! Direct downward shortwave flux, UV/vis
428      real(kind=rb) :: difdnuv(nlay+2)        ! Diffuse downward shortwave flux, UV/vis
429      real(kind=rb) :: dirdnir(nlay+2)        ! Direct downward shortwave flux, near-IR
430      real(kind=rb) :: difdnir(nlay+2)        ! Diffuse downward shortwave flux, near-IR
431
432! Output - inactive
433!      real(kind=rb) :: zuvfu(nlay+2)         ! temporary upward UV shortwave flux (w/m2)
434!      real(kind=rb) :: zuvfd(nlay+2)         ! temporary downward UV shortwave flux (w/m2)
435!      real(kind=rb) :: zuvcu(nlay+2)         ! temporary clear sky upward UV shortwave flux (w/m2)
436!      real(kind=rb) :: zuvcd(nlay+2)         ! temporary clear sky downward UV shortwave flux (w/m2)
437!      real(kind=rb) :: zvsfu(nlay+2)         ! temporary upward visible shortwave flux (w/m2)
438!      real(kind=rb) :: zvsfd(nlay+2)         ! temporary downward visible shortwave flux (w/m2)
439!      real(kind=rb) :: zvscu(nlay+2)         ! temporary clear sky upward visible shortwave flux (w/m2)
440!      real(kind=rb) :: zvscd(nlay+2)         ! temporary clear sky downward visible shortwave flux (w/m2)
441!      real(kind=rb) :: znifu(nlay+2)         ! temporary upward near-IR shortwave flux (w/m2)
442!      real(kind=rb) :: znifd(nlay+2)         ! temporary downward near-IR shortwave flux (w/m2)
443!      real(kind=rb) :: znicu(nlay+2)         ! temporary clear sky upward near-IR shortwave flux (w/m2)
444!      real(kind=rb) :: znicd(nlay+2)         ! temporary clear sky downward near-IR shortwave flux (w/m2)
445
446
447! Initializations
448
449      zepsec = 1.e-06_rb
450      zepzen = 1.e-10_rb
451      oneminus = 1.0_rb - zepsec
452      pi = 2._rb * asin(1._rb)
453
454      istart = jpb1
455      iend = jpb2
456      iout = 0
457      icpr = 0
458
459! In a GCM with or without McICA, set nlon to the longitude dimension
460!
461! Set imca to select calculation type:
462!  imca = 0, use standard forward model calculation (clear and overcast only)
463!  imca = 1, use McICA for Monte Carlo treatment of sub-grid cloud variability
464!            (clear, overcast or partial cloud conditions)
465
466! *** This version does not use McICA (imca = 0) ***
467
468! Set icld to select of clear or cloud calculation and cloud
469! overlap method (read by subroutine readprof from input file INPUT_RRTM): 
470! Without McICA, SW calculation is limited to clear or fully overcast conditions.
471! icld = 0, clear only
472! icld = 1, with clouds using random cloud overlap (McICA only)
473! icld = 2, with clouds using maximum/random cloud overlap (McICA only)
474! icld = 3, with clouds using maximum cloud overlap (McICA only)
475      if (icld.lt.0.or.icld.gt.3) icld = 2
476
477! Set iaer to select aerosol option
478! iaer = 0, no aerosols
479! iaer = 6, use six ECMWF aerosol types
480!           input aerosol optical depth at 0.55 microns for each aerosol type (ecaer)
481! iaer = 10, input total aerosol optical depth, single scattering albedo
482!            and asymmetry parameter (tauaer, ssaaer, asmaer) directly
483      if (iaer.ne.0.and.iaer.ne.6.and.iaer.ne.10) iaer = 0
484
485! Set idelm to select between delta-M scaled or unscaled output direct and diffuse fluxes
486! NOTE: total downward fluxes are always delta scaled
487! idelm = 0, output direct and diffuse flux components are not delta scaled
488!            (direct flux does not include forward scattering peak)
489! idelm = 1, output direct and diffuse flux components are delta scaled (default)
490!            (direct flux includes part or most of forward scattering peak)
491      idelm = 1
492
493! Call model and data initialization, compute lookup tables, perform
494! reduction of g-points from 224 to 112 for input absorption
495! coefficient data and other arrays.
496!
497! In a GCM this call should be placed in the model initialization
498! area, since this has to be called only once. 
499!      call rrtmg_sw_ini(cpdair)
500
501! This is the main longitude/column loop in RRTMG.
502! Modify to loop over all columns (nlon) or over daylight columns
503
504      do iplon = 1, ncol
505
506! Prepare atmosphere profile from GCM for use in RRTMG, and define
507! other input parameters
508
509         call inatm_sw (iplon, nlay, icld, iaer, &
510              play, plev, tlay, tlev, tsfc, h2ovmr, &
511              o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr, &
512              adjes, dyofyr, scon, inflgsw, iceflgsw, liqflgsw, &
513              cldfr, taucld, ssacld, asmcld, fsfcld, cicewp, cliqwp, &
514              reice, reliq, tauaer, ssaaer, asmaer, &
515              nlayers, pavel, pz, pdp, tavel, tz, tbound, coldry, wkl, &
516              adjflux, solvar, inflag, iceflag, liqflag, cldfrac, tauc, &
517              ssac, asmc, fsfc, ciwp, clwp, rei, rel, taua, ssaa, asma)
518
519!  For cloudy atmosphere, use cldprop to set cloud optical properties based on
520!  input cloud physical properties.  Select method based on choices described
521!  in cldprop.  Cloud fraction, water path, liquid droplet and ice particle
522!  effective radius must be passed in cldprop.  Cloud fraction and cloud
523!  optical properties are transferred to rrtmg_sw arrays in cldprop. 
524
525!  Without McICA, SW calculation is limited to clear or fully overcast conditions.
526!  Stop model if partial cloudiness is present. 
527
528         do i = 1, nlayers
529            if (cldfrac(i).gt.zepsec .and. cldfrac(i).lt.oneminus) then
530               stop 'PARTIAL CLOUD NOT ALLOWED'
531            endif
532         enddo
533         call cldprop_sw(nlayers, inflag, iceflag, liqflag, cldfrac, &
534                         tauc, ssac, asmc, fsfc, ciwp, clwp, rei, rel, &
535                         taucldorig, taucloud, ssacloud, asmcloud)
536         icpr = 1
537
538! Calculate coefficients for the temperature and pressure dependence of the
539! molecular absorption coefficients by interpolating data from stored
540! reference atmospheres.
541
542         call setcoef_sw(nlayers, pavel, tavel, pz, tz, tbound, coldry, wkl, &
543                         laytrop, layswtch, laylow, jp, jt, jt1, &
544                         co2mult, colch4, colco2, colh2o, colmol, coln2o, &
545                         colo2, colo3, fac00, fac01, fac10, fac11, &
546                         selffac, selffrac, indself, forfac, forfrac, indfor)
547
548
549! Cosine of the solar zenith angle
550!  Prevent using value of zero; ideally, SW model is not called from host model when sun
551!  is below horizon
552
553         cossza = coszen(iplon)
554         if (cossza .lt. zepzen) cossza = zepzen
555
556
557! Transfer albedo, cloud and aerosol properties into arrays for 2-stream radiative transfer
558
559! Surface albedo
560!  Near-IR bands 16-24 and 29 (1-9 and 14), 820-16000 cm-1, 0.625-12.195 microns
561         do ib=1,9
562            albdir(ib) = aldir(iplon)
563            albdif(ib) = aldif(iplon)
564         enddo
565         albdir(nbndsw) = aldir(iplon)
566         albdif(nbndsw) = aldif(iplon)
567!  UV/visible bands 25-28 (10-13), 16000-50000 cm-1, 0.200-0.625 micron
568         do ib=10,13
569            albdir(ib) = asdir(iplon)
570            albdif(ib) = asdif(iplon)
571         enddo
572
573
574! Clouds
575         if (icld.eq.0) then
576
577            ztauc(:,:) = 0._rb
578            ztaucorig(:,:) = 0._rb
579            zasyc(:,:) = 0._rb
580            zomgc(:,:) = 1._rb
581
582         elseif (icld.ge.1) then
583            do i=1,nlayers
584               do ib=1,nbndsw
585                  ztauc(i,ib) = taucloud(i,jpb1-1+ib)
586                  ztaucorig(i,ib) = taucldorig(i,jpb1-1+ib)
587                  zasyc(i,ib) = asmcloud(i,jpb1-1+ib)
588                  zomgc(i,ib) = ssacloud(i,jpb1-1+ib)
589               enddo
590            enddo
591
592         endif   
593
594! Aerosol
595! IAER = 0: no aerosols
596         if (iaer.eq.0) then
597
598            ztaua(:,:) = 0._rb
599            zasya(:,:) = 0._rb
600            zomga(:,:) = 1._rb
601
602! IAER = 6: Use ECMWF six aerosol types. See rrsw_aer.f90 for details.
603! Input aerosol optical thickness at 0.55 micron for each aerosol type (ecaer),
604! or set manually here for each aerosol and layer.
605         elseif (iaer.eq.6) then
606
607!            do i = 1, nlayers
608!               do ia = 1, naerec
609!                  ecaer(iplon,i,ia) = 1.0e-15_rb
610!               enddo
611!            enddo
612
613            do i = 1, nlayers
614               do ib = 1, nbndsw
615                  ztaua(i,ib) = 0._rb
616                  zasya(i,ib) = 0._rb
617                  zomga(i,ib) = 0._rb
618                  do ia = 1, naerec
619                     ztaua(i,ib) = ztaua(i,ib) + rsrtaua(ib,ia) * ecaer(iplon,i,ia)
620                     zomga(i,ib) = zomga(i,ib) + rsrtaua(ib,ia) * ecaer(iplon,i,ia) * &
621                                   rsrpiza(ib,ia)
622                     zasya(i,ib) = zasya(i,ib) + rsrtaua(ib,ia) * ecaer(iplon,i,ia) * &
623                                   rsrpiza(ib,ia) * rsrasya(ib,ia)
624                  enddo
625                  if (ztaua(i,ib) == 0._rb) then
626                     ztaua(i,ib) = 0._rb
627                     zasya(i,ib) = 0._rb
628                     zomga(i,ib) = 1._rb
629                  else
630                     if (zomga(i,ib) /= 0._rb) then
631                        zasya(i,ib) = zasya(i,ib) / zomga(i,ib)
632                     endif
633                     if (ztaua(i,ib) /= 0._rb) then
634                        zomga(i,ib) = zomga(i,ib) / ztaua(i,ib)
635                     endif
636                  endif
637               enddo
638            enddo
639
640! IAER=10: Direct specification of aerosol optical properties from GCM
641         elseif (iaer.eq.10) then
642
643            do i = 1 ,nlayers
644               do ib = 1 ,nbndsw
645                  ztaua(i,ib) = taua(i,ib)
646                  zasya(i,ib) = asma(i,ib)
647                  zomga(i,ib) = ssaa(i,ib)
648               enddo
649            enddo
650
651         endif
652
653
654! Call the 2-stream radiation transfer model
655
656         do i=1,nlayers+1
657            zbbcu(i) = 0._rb
658            zbbcd(i) = 0._rb
659            zbbfu(i) = 0._rb
660            zbbfd(i) = 0._rb
661            zbbcddir(i) = 0._rb
662            zbbfddir(i) = 0._rb
663            zuvcd(i) = 0._rb
664            zuvfd(i) = 0._rb
665            zuvcddir(i) = 0._rb
666            zuvfddir(i) = 0._rb
667            znicd(i) = 0._rb
668            znifd(i) = 0._rb
669            znicddir(i) = 0._rb
670            znifddir(i) = 0._rb
671         enddo
672
673         call spcvrt_sw &
674             (nlayers, istart, iend, icpr, idelm, iout, &
675              pavel, tavel, pz, tz, tbound, albdif, albdir, &
676              cldfrac, ztauc, zasyc, zomgc, ztaucorig, &
677              ztaua, zasya, zomga, cossza, coldry, wkl, adjflux, &       
678              laytrop, layswtch, laylow, jp, jt, jt1, &
679              co2mult, colch4, colco2, colh2o, colmol, coln2o, colo2, colo3, &
680              fac00, fac01, fac10, fac11, &
681              selffac, selffrac, indself, forfac, forfrac, indfor, &
682              zbbfd, zbbfu, zbbcd, zbbcu, zuvfd, zuvcd, znifd, znicd, &
683              zbbfddir, zbbcddir, zuvfddir, zuvcddir, znifddir, znicddir)
684
685! Transfer up and down, clear and total sky fluxes to output arrays.
686! Vertical indexing goes from bottom to top; reverse here for GCM if necessary.
687
688         do i = 1, nlayers+1
689            swuflxc(iplon,i) = zbbcu(i)
690            swdflxc(iplon,i) = zbbcd(i)
691            swuflx(iplon,i) = zbbfu(i)
692            swdflx(iplon,i) = zbbfd(i)
693            uvdflx(i) = zuvfd(i)
694            nidflx(i) = znifd(i)
695!  Direct/diffuse fluxes
696            dirdflux(i) = zbbfddir(i)
697            difdflux(i) = swdflx(iplon,i) - dirdflux(i)
698!  UV/visible direct/diffuse fluxes
699            dirdnuv(i) = zuvfddir(i)
700            difdnuv(i) = zuvfd(i) - dirdnuv(i)
701!  Near-IR direct/diffuse fluxes
702            dirdnir(i) = znifddir(i)
703            difdnir(i) = znifd(i) - dirdnir(i)
704         enddo
705
706!  Total and clear sky net fluxes
707         do i = 1, nlayers+1
708            swnflxc(i) = swdflxc(iplon,i) - swuflxc(iplon,i)
709            swnflx(i) = swdflx(iplon,i) - swuflx(iplon,i)
710         enddo
711
712!  Total and clear sky heating rates
713         do i = 1, nlayers
714            zdpgcp = heatfac / pdp(i)
715            swhrc(iplon,i) = (swnflxc(i+1) - swnflxc(i)) * zdpgcp
716            swhr(iplon,i) = (swnflx(i+1) - swnflx(i)) * zdpgcp
717         enddo
718         swhrc(iplon,nlayers) = 0._rb
719         swhr(iplon,nlayers) = 0._rb
720
721! End longitude loop
722      enddo
723
724      end subroutine rrtmg_sw
725
726!*************************************************************************
727      real(kind=rb) function earth_sun(idn)
728!*************************************************************************
729!
730!  Purpose: Function to calculate the correction factor of Earth's orbit
731!  for current day of the year
732
733!  idn        : Day of the year
734!  earth_sun  : square of the ratio of mean to actual Earth-Sun distance
735
736! ------- Modules -------
737
738      use rrsw_con, only : pi
739
740      integer(kind=im), intent(in) :: idn
741
742      real(kind=rb) :: gamma
743
744      gamma = 2._rb*pi*(idn-1)/365._rb
745
746! Use Iqbal's equation 1.2.1
747
748      earth_sun = 1.000110_rb + .034221_rb * cos(gamma) + .001289_rb * sin(gamma) + &
749                   .000719_rb * cos(2._rb*gamma) + .000077_rb * sin(2._rb*gamma)
750
751      end function earth_sun
752
753!***************************************************************************
754      subroutine inatm_sw (iplon, nlay, icld, iaer, &
755            play, plev, tlay, tlev, tsfc, h2ovmr, &
756            o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr, &
757            adjes, dyofyr, scon, inflgsw, iceflgsw, liqflgsw, &
758            cldfr, taucld, ssacld, asmcld, fsfcld, cicewp, cliqwp, &
759            reice, reliq, tauaer, ssaaer, asmaer, &
760            nlayers, pavel, pz, pdp, tavel, tz, tbound, coldry, wkl, &
761            adjflux, solvar, inflag, iceflag, liqflag, cldfrac, tauc, &
762            ssac, asmc, fsfc, ciwp, clwp, rei, rel, taua, ssaa, asma)
763!***************************************************************************
764!
765!  Input atmospheric profile from GCM, and prepare it for use in RRTMG_SW.
766!  Set other RRTMG_SW input parameters. 
767!
768!***************************************************************************
769
770! --------- Modules ----------
771
772      use parrrsw, only : nbndsw, ngptsw, nstr, nmol, mxmol, &
773                          jpband, jpb1, jpb2, rrsw_scon
774      use rrsw_con, only : fluxfac, heatfac, oneminus, pi, grav, avogad
775      use rrsw_wvn, only : ng, nspa, nspb, wavenum1, wavenum2, delwave
776
777! ------- Declarations -------
778
779! ----- Input -----
780! Note: All volume mixing ratios are in dimensionless units of mole fraction obtained
781! by scaling mass mixing ratio (g/g) with the appropriate molecular weights (g/mol)
782      integer(kind=im), intent(in) :: iplon           ! column loop index
783      integer(kind=im), intent(in) :: nlay            ! number of model layers
784      integer(kind=im), intent(in) :: icld            ! clear/cloud flag
785      integer(kind=im), intent(in) :: iaer            ! aerosol option flag
786
787      real(kind=rb), intent(in) :: play(:,:)          ! Layer pressures (hPa, mb)
788                                                      ! Dimensions: (ncol,nlay)
789      real(kind=rb), intent(in) :: plev(:,:)          ! Interface pressures (hPa, mb)
790                                                      ! Dimensions: (ncol,nlay+1)
791      real(kind=rb), intent(in) :: tlay(:,:)          ! Layer temperatures (K)
792                                                      ! Dimensions: (ncol,nlay)
793      real(kind=rb), intent(in) :: tlev(:,:)          ! Interface temperatures (K)
794                                                      ! Dimensions: (ncol,nlay+1)
795      real(kind=rb), intent(in) :: tsfc(:)            ! Surface temperature (K)
796                                                      ! Dimensions: (ncol)
797      real(kind=rb), intent(in) :: h2ovmr(:,:)        ! H2O volume mixing ratio
798                                                      ! Dimensions: (ncol,nlay)
799      real(kind=rb), intent(in) :: o3vmr(:,:)         ! O3 volume mixing ratio
800                                                      ! Dimensions: (ncol,nlay)
801      real(kind=rb), intent(in) :: co2vmr(:,:)        ! CO2 volume mixing ratio
802                                                      ! Dimensions: (ncol,nlay)
803      real(kind=rb), intent(in) :: ch4vmr(:,:)        ! Methane volume mixing ratio
804                                                      ! Dimensions: (ncol,nlay)
805      real(kind=rb), intent(in) :: n2ovmr(:,:)        ! Nitrous oxide volume mixing ratio
806                                                      ! Dimensions: (ncol,nlay)
807      real(kind=rb), intent(in) :: o2vmr(:,:)         ! Oxygen volume mixing ratio
808                                                      ! Dimensions: (ncol,nlay)
809
810      integer(kind=im), intent(in) :: dyofyr          ! Day of the year (used to get Earth/Sun
811                                                      !  distance if adjflx not provided)
812      real(kind=rb), intent(in) :: adjes              ! Flux adjustment for Earth/Sun distance
813      real(kind=rb), intent(in) :: scon               ! Solar constant (W/m2)
814
815      integer(kind=im), intent(in) :: inflgsw         ! Flag for cloud optical properties
816      integer(kind=im), intent(in) :: iceflgsw        ! Flag for ice particle specification
817      integer(kind=im), intent(in) :: liqflgsw        ! Flag for liquid droplet specification
818
819      real(kind=rb), intent(in) :: cldfr(:,:)         ! Cloud fraction
820                                                      ! Dimensions: (ncol,nlay)
821      real(kind=rb), intent(in) :: taucld(:,:,:)      ! In-cloud optical depth (optional)
822                                                      ! Dimensions: (nbndsw,ncol,nlay)
823      real(kind=rb), intent(in) :: ssacld(:,:,:)      ! In-cloud single scattering albedo
824                                                      ! Dimensions: (nbndsw,ncol,nlay)
825      real(kind=rb), intent(in) :: asmcld(:,:,:)      ! In-cloud asymmetry parameter
826                                                      ! Dimensions: (nbndsw,ncol,nlay)
827      real(kind=rb), intent(in) :: fsfcld(:,:,:)      ! In-cloud forward scattering fraction
828                                                      ! Dimensions: (nbndsw,ncol,nlay)
829      real(kind=rb), intent(in) :: cicewp(:,:)        ! In-cloud ice water path (g/m2)
830                                                      ! Dimensions: (ncol,nlay)
831      real(kind=rb), intent(in) :: cliqwp(:,:)        ! In-cloud liquid water path (g/m2)
832                                                      ! Dimensions: (ncol,nlay)
833      real(kind=rb), intent(in) :: reice(:,:)         ! Cloud ice effective size (microns)
834                                                      ! Dimensions: (ncol,nlay)
835      real(kind=rb), intent(in) :: reliq(:,:)         ! Cloud water drop effective radius (microns)
836                                                      ! Dimensions: (ncol,nlay)
837
838      real(kind=rb), intent(in) :: tauaer(:,:,:)      ! Aerosol optical depth
839                                                      ! Dimensions: (ncol,nlay,nbndsw)
840      real(kind=rb), intent(in) :: ssaaer(:,:,:)      ! Aerosol single scattering albedo
841                                                      ! Dimensions: (ncol,nlay,nbndsw)
842      real(kind=rb), intent(in) :: asmaer(:,:,:)      ! Aerosol asymmetry parameter
843                                                      ! Dimensions: (ncol,nlay,nbndsw)
844
845! Atmosphere
846      integer(kind=im), intent(out) :: nlayers        ! number of layers
847
848      real(kind=rb), intent(out) :: pavel(:)          ! layer pressures (mb)
849                                                      ! Dimensions: (nlay)
850      real(kind=rb), intent(out) :: tavel(:)          ! layer temperatures (K)
851                                                      ! Dimensions: (nlay)
852      real(kind=rb), intent(out) :: pz(0:)            ! level (interface) pressures (hPa, mb)
853                                                      ! Dimensions: (0:nlay)
854      real(kind=rb), intent(out) :: tz(0:)            ! level (interface) temperatures (K)
855                                                      ! Dimensions: (0:nlay)
856      real(kind=rb), intent(out) :: tbound            ! surface temperature (K)
857      real(kind=rb), intent(out) :: pdp(:)            ! layer pressure thickness (hPa, mb)
858                                                      ! Dimensions: (nlay)
859      real(kind=rb), intent(out) :: coldry(:)         ! dry air column density (mol/cm2)
860                                                      ! Dimensions: (nlay)
861      real(kind=rb), intent(out) :: wkl(:,:)          ! molecular amounts (mol/cm-2)
862                                                      ! Dimensions: (mxmol,nlay)
863
864      real(kind=rb), intent(out) :: adjflux(:)        ! adjustment for current Earth/Sun distance
865                                                      ! Dimensions: (jpband)
866      real(kind=rb), intent(out) :: solvar(:)         ! solar constant scaling factor from rrtmg_sw
867                                                      ! Dimensions: (jpband)
868                                                      !  default value of 1368.22 Wm-2 at 1 AU
869      real(kind=rb), intent(out) :: taua(:,:)         ! Aerosol optical depth
870                                                      ! Dimensions: (nlay,nbndsw)
871      real(kind=rb), intent(out) :: ssaa(:,:)         ! Aerosol single scattering albedo
872                                                      ! Dimensions: (nlay,nbndsw)
873      real(kind=rb), intent(out) :: asma(:,:)         ! Aerosol asymmetry parameter
874                                                      ! Dimensions: (nlay,nbndsw)
875
876! Atmosphere/clouds - cldprop
877      integer(kind=im), intent(out) :: inflag         ! flag for cloud property method
878      integer(kind=im), intent(out) :: iceflag        ! flag for ice cloud properties
879      integer(kind=im), intent(out) :: liqflag        ! flag for liquid cloud properties
880
881      real(kind=rb), intent(out) :: cldfrac(:)        ! layer cloud fraction
882                                                      ! Dimensions: (nlay)
883      real(kind=rb), intent(out) :: tauc(:,:)         ! in-cloud optical depth (non-delta scaled)
884                                                      ! Dimensions: (nbndsw,nlay)
885      real(kind=rb), intent(out) :: ssac(:,:)         ! in-cloud single scattering albedo (non-delta-scaled)
886                                                      ! Dimensions: (nbndsw,nlay)
887      real(kind=rb), intent(out) :: asmc(:,:)         ! in-cloud asymmetry parameter (non-delta scaled)
888                                                      ! Dimensions: (nbndsw,nlay)
889      real(kind=rb), intent(out) :: fsfc(:,:)         ! in-cloud forward scattering fraction (non-delta scaled)
890                                                      ! Dimensions: (nbndsw,nlay)
891      real(kind=rb), intent(out) :: ciwp(:)           ! in-cloud ice water path
892                                                      ! Dimensions: (nlay)
893      real(kind=rb), intent(out) :: clwp(:)           ! in-cloud liquid water path
894                                                      ! Dimensions: (nlay)
895      real(kind=rb), intent(out) :: rel(:)            ! cloud liquid particle effective radius (microns)
896                                                      ! Dimensions: (nlay)
897      real(kind=rb), intent(out) :: rei(:)            ! cloud ice particle effective size (microns)
898                                                      ! Dimensions: (nlay)
899
900! ----- Local -----
901      real(kind=rb), parameter :: amd = 28.9660_rb    ! Effective molecular weight of dry air (g/mol)
902      real(kind=rb), parameter :: amw = 18.0160_rb    ! Molecular weight of water vapor (g/mol)
903!      real(kind=rb), parameter :: amc = 44.0098_rb   ! Molecular weight of carbon dioxide (g/mol)
904!      real(kind=rb), parameter :: amo = 47.9998_rb   ! Molecular weight of ozone (g/mol)
905!      real(kind=rb), parameter :: amo2 = 31.9999_rb  ! Molecular weight of oxygen (g/mol)
906!      real(kind=rb), parameter :: amch4 = 16.0430_rb ! Molecular weight of methane (g/mol)
907!      real(kind=rb), parameter :: amn2o = 44.0128_rb ! Molecular weight of nitrous oxide (g/mol)
908
909! Set molecular weight ratios (for converting mmr to vmr)
910!  e.g. h2ovmr = h2ommr * amdw)
911      real(kind=rb), parameter :: amdw = 1.607793_rb  ! Molecular weight of dry air / water vapor
912      real(kind=rb), parameter :: amdc = 0.658114_rb  ! Molecular weight of dry air / carbon dioxide
913      real(kind=rb), parameter :: amdo = 0.603428_rb  ! Molecular weight of dry air / ozone
914      real(kind=rb), parameter :: amdm = 1.805423_rb  ! Molecular weight of dry air / methane
915      real(kind=rb), parameter :: amdn = 0.658090_rb  ! Molecular weight of dry air / nitrous oxide
916      real(kind=rb), parameter :: amdo2 = 0.905140_rb ! Molecular weight of dry air / oxygen
917
918      real(kind=rb), parameter :: sbc = 5.67e-08_rb   ! Stefan-Boltzmann constant (W/m2K4)
919
920      integer(kind=im) :: isp, l, ix, n, imol, ib       ! Loop indices
921      real(kind=rb) :: amm, summol                      !
922      real(kind=rb) :: adjflx                           ! flux adjustment for Earth/Sun distance
923!      real(kind=rb) :: earth_sun                        ! function for Earth/Sun distance adjustment
924
925! Add one to nlayers here to include extra model layer at top of atmosphere
926      nlayers = nlay
927
928!  Initialize all molecular amounts to zero here, then pass input amounts
929!  into RRTM array WKL below.
930
931      wkl(:,:) = 0.0_rb
932      cldfrac(:) = 0.0_rb
933      tauc(:,:) = 0.0_rb
934      ssac(:,:) = 1.0_rb
935      asmc(:,:) = 0.0_rb
936      fsfc(:,:) = 0.0_rb
937      ciwp(:) = 0.0_rb
938      clwp(:) = 0.0_rb
939      rei(:) = 0.0_rb
940      rel(:) = 0.0_rb
941      taua(:,:) = 0.0_rb
942      ssaa(:,:) = 1.0_rb
943      asma(:,:) = 0.0_rb
944 
945! Set flux adjustment for current Earth/Sun distance (two options).
946! 1) Use Earth/Sun distance flux adjustment provided by GCM (input as adjes);
947      adjflx = adjes
948!
949! 2) Calculate Earth/Sun distance from DYOFYR, the cumulative day of the year.
950!    (Set adjflx to 1. to use constant Earth/Sun distance of 1 AU).
951      if (dyofyr .gt. 0) then
952         adjflx = earth_sun(dyofyr)
953      endif
954
955! Set incoming solar flux adjustment to include adjustment for
956! current Earth/Sun distance (ADJFLX) and scaling of default internal
957! solar constant (rrsw_scon = 1368.22 Wm-2) by band (SOLVAR).  SOLVAR can be set
958! to a single scaling factor as needed, or to a different value in each
959! band, which may be necessary for paleoclimate simulations.
960!
961      do ib = jpb1,jpb2
962!         solvar(ib) = 1._rb
963         solvar(ib) = scon / rrsw_scon
964         adjflux(ib) = adjflx * solvar(ib)
965      enddo
966
967!  Set surface temperature.
968      tbound = tsfc(iplon)
969
970!  Install input GCM arrays into RRTMG_SW arrays for pressure, temperature,
971!  and molecular amounts. 
972!  Pressures are input in mb, or are converted to mb here.
973!  Molecular amounts are input in volume mixing ratio, or are converted from
974!  mass mixing ratio (or specific humidity for h2o) to volume mixing ratio
975!  here. These are then converted to molecular amount (molec/cm2) below. 
976!  The dry air column COLDRY (in molec/cm2) is calculated from the level
977!  pressures, pz (in mb), based on the hydrostatic equation and includes a
978!  correction to account for h2o in the layer.  The molecular weight of moist
979!  air (amm) is calculated for each layer. 
980!  Note: In RRTMG, layer indexing goes from bottom to top, and coding below
981!  assumes GCM input fields are also bottom to top. Input layer indexing
982!  from GCM fields should be reversed here if necessary.
983
984      pz(0) = plev(iplon,1)
985      tz(0) = tlev(iplon,1)
986      do l = 1, nlayers
987         pavel(l) = play(iplon,l)
988         tavel(l) = tlay(iplon,l)
989         pz(l) = plev(iplon,l+1)
990         tz(l) = tlev(iplon,l+1)
991         pdp(l) = pz(l-1) - pz(l)
992! For h2o input in vmr:
993         wkl(1,l) = h2ovmr(iplon,l)
994! For h2o input in mmr:
995!         wkl(1,l) = h2o(iplon,l)*amdw
996! For h2o input in specific humidity;
997!         wkl(1,l) = (h2o(iplon,l)/(1._rb - h2o(iplon,l)))*amdw
998         wkl(2,l) = co2vmr(iplon,l)
999         wkl(3,l) = o3vmr(iplon,l)
1000         wkl(4,l) = n2ovmr(iplon,l)
1001         wkl(6,l) = ch4vmr(iplon,l)
1002         wkl(7,l) = o2vmr(iplon,l)
1003         amm = (1._rb - wkl(1,l)) * amd + wkl(1,l) * amw           
1004         coldry(l) = (pz(l-1)-pz(l)) * 1.e3_rb * avogad / &
1005                     (1.e2_rb * grav * amm * (1._rb + wkl(1,l)))
1006      enddo
1007
1008! The following section can be used to set values for an additional layer (from
1009! the GCM top level to 1.e-4 mb) for improved calculation of TOA fluxes.
1010! Temperature and molecular amounts in the extra model layer are set to
1011! their values in the top GCM model layer, though these can be modified
1012! here if necessary.
1013! If this feature is utilized, increase nlayers by one above, limit the two
1014! loops above to (nlayers-1), and set the top most (nlayers) layer values here.
1015
1016!      pavel(nlayers) = 0.5_rb * pz(nlayers-1)
1017!      tavel(nlayers) = tavel(nlayers-1)
1018!      pz(nlayers) = 1.e-4_rb
1019!      tz(nlayers-1) = 0.5_rb * (tavel(nlayers)+tavel(nlayers-1))
1020!      tz(nlayers) = tz(nlayers-1)
1021!      pdp(nlayers) = pz(nlayers-1) - pz(nlayers)
1022!      wkl(1,nlayers) = wkl(1,nlayers-1)
1023!      wkl(2,nlayers) = wkl(2,nlayers-1)
1024!      wkl(3,nlayers) = wkl(3,nlayers-1)
1025!      wkl(4,nlayers) = wkl(4,nlayers-1)
1026!      wkl(6,nlayers) = wkl(6,nlayers-1)
1027!      wkl(7,nlayers) = wkl(7,nlayers-1)
1028!      amm = (1._rb - wkl(1,nlayers-1)) * amd + wkl(1,nlayers-1) * amw
1029!      coldry(nlayers) = (pz(nlayers-1)) * 1.e3_rb * avogad / &
1030!                        (1.e2_rb * grav * amm * (1._rb + wkl(1,nlayers-1)))
1031
1032! At this point all molecular amounts in wkl are in volume mixing ratio;
1033! convert to molec/cm2 based on coldry for use in rrtm. 
1034
1035      do l = 1, nlayers
1036         do imol = 1, nmol
1037            wkl(imol,l) = coldry(l) * wkl(imol,l)
1038         enddo
1039      enddo
1040
1041! Transfer aerosol optical properties to RRTM variables;
1042! modify to reverse layer indexing here if necessary.
1043
1044      if (iaer .ge. 1) then
1045         do l = 1, nlayers
1046            do ib = 1, nbndsw
1047               taua(l,ib) = tauaer(iplon,l,ib)
1048               ssaa(l,ib) = ssaaer(iplon,l,ib)
1049               asma(l,ib) = asmaer(iplon,l,ib)
1050            enddo
1051         enddo
1052      endif
1053
1054! Transfer cloud fraction and cloud optical properties to RRTM variables;
1055! modify to reverse layer indexing here if necessary.
1056
1057      if (icld .ge. 1) then
1058         inflag = inflgsw
1059         iceflag = iceflgsw
1060         liqflag = liqflgsw
1061
1062! Move incoming GCM cloud arrays to RRTMG cloud arrays.
1063! For GCM input, incoming reice is defined based on selected ice parameterization (inflglw)
1064
1065         do l = 1, nlayers
1066            cldfrac(l) = cldfr(iplon,l)
1067            ciwp(l) = cicewp(iplon,l)
1068            clwp(l) = cliqwp(iplon,l)
1069            rei(l) = reice(iplon,l)
1070            rel(l) = reliq(iplon,l)
1071            do n = 1,nbndsw
1072               tauc(n,l) = taucld(n,iplon,l)
1073               ssac(n,l) = ssacld(n,iplon,l)
1074               asmc(n,l) = asmcld(n,iplon,l)
1075               fsfc(n,l) = fsfcld(n,iplon,l)
1076            enddo
1077         enddo
1078
1079! If an extra layer is being used in RRTMG, set all cloud properties to zero in the extra layer.
1080
1081!         cldfrac(nlayers) = 0.0_rb
1082!         tauc(:,nlayers) = 0.0_rb
1083!         ssac(:,nlayers) = 1.0_rb
1084!         asmc(:,nlayers) = 0.0_rb
1085!         fsfc(:,nlayers) = 0.0_rb
1086!         ciwp(nlayers) = 0.0_rb
1087!         clwp(nlayers) = 0.0_rb
1088!         rei(nlayers) = 0.0_rb
1089!         rel(nlayers) = 0.0_rb
1090     
1091      endif
1092
1093      end subroutine inatm_sw
1094
1095      end module rrtmg_sw_rad
1096
1097
Note: See TracBrowser for help on using the repository browser.