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

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

Added support for RRTMG radiation code

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