chem/module_subs_tuv.F

   1       MODULE tuv_subs
   2
   3       use params_mod, only : dp
   4
   5       IMPLICIT none
   6
   7       private
   8       public :: tuv_radfld, sundis, calc_zenith
   9
  10       integer :: nstr = 1        ! stream count
  11
  12       CONTAINS
  13
  14       SUBROUTINE tuv_radfld( nlambda_start, cld_od_opt, cldfrac, nlyr, nwave, &
  15                              zenith, z, albedo, &
  16                              aircol, o2col, o3col, so2col, no2col, &
  17                              tauaer300, tauaer400, tauaer600, tauaer999, &
  18                              waer300, waer400, waer600, waer999, &
  19                              gaer300, gaer400, gaer600, gaer999, &
  20                              dtaer, omaer, gaer, dtcld, omcld, gcld, &
  21                              has_aer_ra_feedback, &
  22                              qll, dobsi, o3_xs, no2_xs, o2_xs, &
  23                              so2_xs, wmin, wc, tlev, srb_o2_xs, radfld, efld, &
  24                              e_dir, e_dn, e_up, &
  25 #ifdef SW_DEBUG
  26                              dir_fld, dwn_fld, up_fld, dt_cld, tuv_diags )
  27 #else
  28                              dir_fld, dwn_fld, up_fld, dt_cld )
  29 #endif
  30 !-----------------------------------------------------------------------------
  31 !     ... calculate the radiation field
  32 !-----------------------------------------------------------------------------
  33
  34       use srb,       only : la_srb
  35       use rad_trans, only : rtlink
  36
  37 !-----------------------------------------------------------------------------
  38 !     ... dummy arguments
  39 !-----------------------------------------------------------------------------
  40       integer, intent(in)  :: nlambda_start
  41       integer, intent(in)  :: nlyr
  42       integer, intent(in)  :: nwave
  43       integer, intent(in)  :: cld_od_opt
  44       real, intent(in)  :: zenith
  45       real, intent(in)  :: dobsi
  46       real, intent(in)  :: wmin
  47       real, intent(in)  :: z(:)
  48       real, intent(in)  :: albedo(:)
  49       real, intent(in)  :: aircol(:)
  50       real, intent(in)  :: o2col(:)
  51       real, intent(in)  :: o3col(:)
  52       real, intent(in)  :: so2col(:)
  53       real, intent(in)  :: no2col(:)
  54       real, intent(in)  :: tauaer300(:)
  55       real, intent(in)  :: tauaer400(:)
  56       real, intent(in)  :: tauaer600(:)
  57       real, intent(in)  :: tauaer999(:)
  58       real, intent(in)  :: waer300(:)
  59       real, intent(in)  :: waer400(:)
  60       real, intent(in)  :: waer600(:)
  61       real, intent(in)  :: waer999(:)
  62       real, intent(in)  :: gaer300(:)
  63       real, intent(in)  :: gaer400(:)
  64       real, intent(in)  :: gaer600(:)
  65       real, intent(in)  :: gaer999(:)
  66       real, intent(in)  :: qll(:)
  67       real, intent(in)  :: wc(:)
  68       real, intent(in)  :: tlev(:)
  69       real, intent(in)  :: cldfrac(:)
  70       real, intent(in)  :: o2_xs(:)
  71       real, intent(in)  :: so2_xs(:)
  72       real, intent(in)  :: o3_xs(:,:)
  73       real, intent(in)  :: no2_xs(:,:)
  74       real, intent(out) :: srb_o2_xs(:,:)
  75       real, intent(out) :: radfld(:,:)
  76       real, intent(out) :: efld(:,:)
  77       real, intent(inout)  :: dir_fld(:,:), dwn_fld(:,:), up_fld(:,:)
  78       real, intent(inout)  :: e_dir(:,:), e_dn(:,:), e_up(:,:)
  79       real, intent(inout)  :: dt_cld(:)
  80       real, intent(inout)  :: dtaer(:,:), omaer(:,:), gaer(:,:)
  81       real, intent(inout)  :: dtcld(:,:), omcld(:,:), gcld(:,:)
  82       logical, intent(in)  :: has_aer_ra_feedback
  83 #ifdef SW_DEBUG
  84       logical, intent(in) :: tuv_diags
  85 #endif
  86
  87 !-----------------------------------------------------------------------------
  88 !     ... local variables
  89 !-----------------------------------------------------------------------------
  90       integer :: k, n
  91       integer :: wn
  92       integer :: n_radlev, n_radlevp1
  93       integer :: nid(0:nlyr)
  94       real :: dtrl(nlyr,nwave)
  95       real :: dto2(nlyr,nwave)
  96       real :: dto3(nlyr,nwave)
  97       real :: dtso2(nlyr,nwave)
  98       real :: dtno2(nlyr,nwave)
  99 !     real :: dtcld(nlyr,nwave)
 100 !     real :: dtaer(nlyr,nwave)
 101       real :: dtsnw(nlyr,nwave)
 102
 103 !     real :: omcld(nlyr,nwave)
 104 !     real :: gcld(nlyr,nwave)
 105 !     real :: omaer(nlyr,nwave)
 106 !     real :: gaer(nlyr,nwave)
 107       real :: omsnw(nlyr,nwave)
 108       real :: gsnw(nlyr,nwave)
 109
 110       real :: edir(nlyr+1)
 111       real :: edn(nlyr+1)
 112       real :: eup(nlyr+1)
 113       real :: fdir(nlyr+1)
 114       real :: fdn(nlyr+1)
 115       real :: fup(nlyr+1)
 116
 117       real :: vcol(nlyr)
 118       real :: scol(nlyr)
 119
 120       real :: dsdh(0:nlyr,nlyr)
 121
 122       n_radlev = size( radfld,dim=2 )
 123       n_radlevp1 = n_radlev + 1
 124
 125       do wn = 1,nwave
 126         omcld(:,wn) = 0.
 127         omaer(:,wn) = 0.
 128         omsnw(:,wn) = 0.
 129         gcld(:,wn)  = 0.
 130         gaer(:,wn)  = 0.
 131         gsnw(:,wn)  = 0.
 132         dtcld(:,wn) = 0.
 133         dtaer(:,wn) = 0.
 134         dtsnw(:,wn) = 0.
 135       end do
 136
 137       call odrl( wc, aircol, dtrl )
 138       call seto2( o2col, o2_xs, dto2 )
 139       call odo3( o3col, o3_xs, dto3, dobsi )
 140       call setso2( so2col, so2_xs, dtso2 )
 141       call setno2( no2col, no2_xs, dtno2 )
 142 !-------------------------------------------------------------
 143 ! aerosol optical depths
 144 !-------------------------------------------------------------
 145       if( has_aer_ra_feedback ) then
 146         call setaer( nlambda_start, wc, tauaer300, tauaer400, &
 147                      tauaer600, tauaer999, waer300, &
 148                      waer400, waer600, waer999,     &
 149                      gaer300, gaer400, gaer600,     &
 150                      gaer999, dtaer, omaer, gaer )
 151       endif
 152 !-------------------------------------------------------------
 153 ! cloud optical depths (cloud water units = g/m3)
 154 !-------------------------------------------------------------
 155       call setcld( nlambda_start, cld_od_opt, z, qll, cldfrac, &
 156                    dtcld, omcld, gcld )
 157       dt_cld(:n_radlev) = dtcld(2:n_radlevp1,1)
 158
 159       call sphers( nlyr, z, zenith, dsdh, nid )
 160
 161 #ifdef SW_DEBUG
 162       if( tuv_diags ) then
 163         open(unit=33,file='WRF-TUV.dbg.out')
 164         write(33,*) 'tuv_radfld: tuv_diags'
 165         write(33,'(''nlyr = '',i4)') nlyr
 166         write(33,'(''dsdh(1,1) = '',1p,g15.7)') dsdh(1,1)
 167         write(33,*) 'dsdh(nlyr,:)'
 168         do n = 1,nlyr,5
 169           write(33,'(1p,5g15.7)') dsdh(nlyr,n:min(n+4,nlyr))
 170         end do
 171         close(33)
 172       endif
 173 #endif
 174
 175       call airmas( nlyr, dsdh, nid, aircol, vcol, scol )
 176       call la_srb( nlyr, z, tlev, wmin, &
 177                    vcol, scol, o2_xs, dto2, srb_o2_xs )
 178
 179       do wn = nlambda_start,nwave
 180         call rtlink( &
 181            nstr, nlyr+1, nlyr, nwave, &
 182            wn, albedo(wn), zenith, &
 183            dsdh, nid, &
 184            dtrl,  &
 185            dto3,  &
 186            dto2, &
 187            dtso2, &
 188            dtno2,  &
 189            dtcld, omcld, gcld, &
 190            dtaer, omaer, gaer, &
 191            dtsnw, omsnw, gsnw, &
 192 #ifdef SW_DEBUG
 193            edir, edn, eup, fdir, fdn, fup, tuv_diags )
 194 #else
 195            edir, edn, eup, fdir, fdn, fup )
 196 #endif
 197 !       radfld(wn,1:nlyr-1)  = fdir(2:nlyr) + fdn(2:nlyr) + fup(2:nlyr)
 198 !       efld(1:nlyr-1,wn)    = edir(2:nlyr) + edn(2:nlyr) + eup(2:nlyr)
 199 !       dir_fld(1:nlyr-1,wn) = fdir(2:nlyr)
 200 !       dwn_fld(1:nlyr-1,wn) = fdn(2:nlyr)
 201 !       up_fld(1:nlyr-1,wn)  = fup(2:nlyr)
 202 !       e_dir(1:nlyr-1,wn)   = edir(2:nlyr)
 203 !       e_dn(1:nlyr-1,wn)    = edn(2:nlyr)
 204 !       e_up(1:nlyr-1,wn)    = eup(2:nlyr)
 205         radfld(wn,1:n_radlev) = fdir(2:n_radlevp1) + fdn(2:n_radlevp1) + fup(2:n_radlevp1)
 206         efld(1:n_radlev,wn)    = edir(2:n_radlevp1) + edn(2:n_radlevp1) + eup(2:n_radlevp1)
 207         dir_fld(1:n_radlev,wn) = fdir(2:n_radlevp1)
 208         dwn_fld(1:n_radlev,wn) = fdn(2:n_radlevp1)
 209         up_fld(1:n_radlev,wn)  = fup(2:n_radlevp1)
 210         e_dir(1:n_radlev,wn)   = edir(2:n_radlevp1)
 211         e_dn(1:n_radlev,wn)    = edn(2:n_radlevp1)
 212         e_up(1:n_radlev,wn)    = eup(2:n_radlevp1)
 213       end do
 214
 215       END SUBROUTINE tuv_radfld
 216
 217       SUBROUTINE odrl( wc, aircol, dtrl )
 218 !-----------------------------------------------------------------------------*
 219 !=  PURPOSE:                                                                 =*
 220 !=  Compute Rayleigh optical depths as a function of altitude and wavelength =*
 221 !-----------------------------------------------------------------------------*
 222 !=  PARAMETERS:                                                              =*
 223 !=  C       - REAL, number of air molecules per cm^2 at each specified    (O)=*
 224 !=            altitude layer                                                 =*
 225 !=  DTRL    - REAL, Rayleigh optical depth at each specified altitude     (O)=*
 226 !=            and each specified wavelength                                  =*
 227 !-----------------------------------------------------------------------------*
 228
 229 !-----------------------------------------------------------------------------*
 230 !     ...dummy arguments
 231 !-----------------------------------------------------------------------------*
 232       REAL,    intent(in)  :: aircol(:)
 233       REAL,    intent(in)  :: wc(:)
 234       REAL,    intent(out) :: dtrl(:,:)
 235
 236 !-----------------------------------------------------------------------------*
 237 !     ...local variables
 238 !-----------------------------------------------------------------------------*
 239       INTEGER :: nwave, nlyr
 240       INTEGER :: wn
 241       REAL    :: srayl, wmicrn, xx
 242
 243       nwave = size( wc )
 244       nlyr  = size( aircol )
 245 !-----------------------------------------------------------------------------*
 246 ! compute Rayleigh cross sections and depths:
 247 !-----------------------------------------------------------------------------*
 248       DO wn = 1,nwave
 249 !-----------------------------------------------------------------------------*
 250 ! Rayleigh scattering cross section from WMO 1985 (originally from
 251 ! Nicolet, M., On the molecular scattering in the terrestrial atmosphere:
 252 ! An empirical formula for its calculation in the homoshpere, Planet.
 253 ! Space Sci., 32, 1467-1468, 1984.
 254 !-----------------------------------------------------------------------------*
 255         wmicrn =  wc(wn)*1.E-3
 256         IF( wmicrn <= 0.55 ) THEN
 257           xx = 3.6772 + 0.389*wmicrn + 0.09426/wmicrn
 258         ELSE
 259           xx = 4.04
 260         ENDIF
 261         srayl = 4.02e-28/(wmicrn)**xx
 262 !-----------------------------------------------------------------------------*
 263 ! alternate (older) expression from
 264 ! Frohlich and Shaw, Appl.Opt. v.11, p.1773 (1980).
 265 !-----------------------------------------------------------------------------*
 266         dtrl(:nlyr,wn) = aircol(:nlyr)*srayl
 267       END DO
 268
 269       END SUBROUTINE odrl
 270
 271       SUBROUTINE odo3( o3col, o3xs, dto3, dobsi )
 272 !-----------------------------------------------------------------------------
 273 !=  NAME:  Optical Depths of O3
 274 !=  PURPOSE:
 275 !=  Compute ozone optical depths as a function of altitude and wavelength
 276 !-----------------------------------------------------------------------------
 277 !=  PARAMETERS:
 278 !=  O3XS   - REAL, molecular absoprtion cross section (cm^2) of O3 at     (I)
 279 !=           each specified wavelength and altitude
 280 !=  C      - REAL, ozone vertical column increments, molec cm-2, for each (I)
 281 !=           layer
 282 !=  DTO3   - REAL, optical depth due to ozone absorption at each          (O)
 283 !=           specified altitude at each specified wavelength
 284 !-----------------------------------------------------------------------------
 285
 286 !-----------------------------------------------------------------------------
 287 !     ... dummy arguments
 288 !-----------------------------------------------------------------------------
 289       REAL, intent(in)    :: dobsi
 290       REAL, intent(in)    :: o3col(:)
 291       REAL, intent(in)    :: o3xs(:,:)
 292       REAL, intent(inout) :: dto3(:,:)
 293
 294       INTEGER :: nlyr, nwave
 295       INTEGER :: wn
 296       REAL    :: dob_at_grnd, scale_fac
 297
 298       nwave = size(o3xs,dim=1)
 299       nlyr  = size(o3col)
 300
 301       if( dobsi == 0. ) then
 302 !-----------------------------------------------------------------------------
 303 !  no scaling
 304 !-----------------------------------------------------------------------------
 305         DO wn = 1,nwave
 306           dto3(:nlyr,wn) = o3col(:nlyr) * o3xs(wn,:nlyr)
 307         END DO
 308       else
 309 !-----------------------------------------------------------------------------
 310 !  scale model o3 column to dobsi
 311 !-----------------------------------------------------------------------------
 312         dob_at_grnd = sum( o3col(:nlyr) )/2.687e16
 313         scale_fac   = dobsi/dob_at_grnd
 314         DO wn = 1,nwave
 315           dto3(:nlyr,wn) = scale_fac * o3col(:nlyr) * o3xs(wn,:nlyr)
 316         END DO
 317       endif
 318
 319       END SUBROUTINE odo3
 320
 321       SUBROUTINE seto2( o2col, o2xs1, dto2 )
 322 !-----------------------------------------------------------------------------
 323 !=  PURPOSE:
 324 !=  Set up an altitude profile of air molecules.  Subroutine includes a
 325 !=  shape-conserving scaling method that allows scaling of the entire
 326 !=  profile to a given sea-level pressure.
 327 !-----------------------------------------------------------------------------
 328 !=  PARAMETERS:
 329 !=  CZ      - REAL, number of air molecules per cm^2 at each specified    (O)
 330 !=            altitude layer
 331 !-----------------------------------------------------------------------------
 332
 333 !-----------------------------------------------------------------------------
 334 !     ... dummy arguments
 335 !-----------------------------------------------------------------------------
 336       REAL, intent(in)  :: o2col(:)
 337       REAL, intent(in)  :: o2xs1(:)
 338       REAL, intent(out) :: dto2(:,:)
 339
 340 !-----------------------------------------------------------------------------
 341 !     ... local variables
 342 !-----------------------------------------------------------------------------
 343       INTEGER :: nlyr, nwave
 344       INTEGER :: wn
 345
 346       nwave = size(o2xs1)
 347       nlyr  = size(o2col)
 348
 349 !-----------------------------------------------------------------------------
 350 !  Assumes that O2 = 20.95 % of air density.  If desire different O2
 351 !  profile (e.g. for upper atmosphere) then can load it here.
 352 !-----------------------------------------------------------------------------
 353       DO wn = 1,nwave
 354         dto2(:nlyr,wn) = o2col(:nlyr) * o2xs1(wn)
 355       ENDDO
 356
 357       END SUBROUTINE seto2
 358
 359       SUBROUTINE setso2( colso2, so2_xs, dtso2 )
 360 !-----------------------------------------------------------------------------
 361 !=  PURPOSE:
 362 !=  Set up an altitude profile of SO2 molecules, and corresponding absorption
 363 !=  optical depths.  Subroutine includes a shape-conserving scaling method
 364 !=  that allows scaling of the entire profile to a given overhead SO2
 365 !=  column amount.
 366 !-----------------------------------------------------------------------------
 367 !=  PARAMETERS:
 368 !=  SO2_XS - REAL, molecular absoprtion cross section (cm^2) of O2 at     (I)
 369 !=           each specified wavelength
 370 !=  DTSO2  - REAL, optical depth due to SO2 absorption at each            (O)
 371 !=           specified altitude at each specified wavelength
 372 !-----------------------------------------------------------------------------
 373
 374 !-----------------------------------------------------------------------------
 375 !     ... dummy arguments
 376 !-----------------------------------------------------------------------------
 377       REAL,    intent(in)  :: colso2(:)
 378       REAL,    intent(in)  :: so2_xs(:)
 379       REAL,    intent(out) :: dtso2(:,:)
 380
 381 !-----------------------------------------------------------------------------
 382 !     ... local variables
 383 !-----------------------------------------------------------------------------
 384       integer :: nwave, nlyr
 385       integer :: wn
 386
 387       nwave = size( so2_xs )
 388       nlyr  = size( colso2 )
 389
 390       DO wn = 1,nwave
 391         dtso2(:nlyr,wn) = colso2(:nlyr)*so2_xs(wn)
 392       END DO
 393
 394       END SUBROUTINE setso2
 395
 396       SUBROUTINE setno2( colno2, no2_xs, dtno2 )
 397 !-----------------------------------------------------------------------------
 398 !=  NAME:  Optical Depths of no2
 399 !=  PURPOSE:
 400 !=  Compute no2 optical depths as a function of altitude and wavelength
 401 !-----------------------------------------------------------------------------
 402 !=  PARAMETERS:
 403 !=  NO2_XS - REAL, molecular absoprtion cross section (cm^2) of no2 at    (I)
 404 !=           each specified wavelength and altitude
 405 !=  COLNO2 - REAL, no2 vertical column increments, molec cm-2, for each   (I)
 406 !=           layer
 407 !=  DTNO2  - REAL, optical depth due to no2 absorption at each            (O)
 408 !=           specified altitude at each specified wavelength
 409 !-----------------------------------------------------------------------------
 410
 411 !-----------------------------------------------------------------------------
 412 !     ... dummy arguments
 413 !-----------------------------------------------------------------------------
 414       REAL, intent(in)    :: colno2(:)
 415       REAL, intent(in)    :: no2_xs(:,:)
 416       REAL, intent(inout) :: dtno2(:,:)
 417
 418       INTEGER :: nlyr, nwave
 419       INTEGER :: wn
 420
 421       nwave = size(no2_xs,dim=1)
 422       nlyr  = size(colno2)
 423
 424       DO wn = 1,nwave
 425         dtno2(:nlyr,wn) = colno2(:nlyr) * no2_xs(wn,:nlyr)
 426       END DO
 427
 428       END SUBROUTINE setno2
 429
 430       subroutine setaer( nlambda_start, wc, tauaer300, tauaer400, &
 431                          tauaer600, tauaer999,               &
 432                          waer300, waer400, waer600, waer999, &
 433                          gaer300, gaer400, gaer600, gaer999, &
 434                          dtaer, omaer, gaer )
 435 !----------------------------------------------------------------------
 436 ! The routine is based on aerosol treatment in module_ra_rrtmg_sw.F
 437 ! INPUT:
 438 ! nzlev: number of specified altitude levels in the working grid
 439 ! z: specified altitude working grid
 440 ! Aerosol optical properties at 300, 400, 600 and 999 nm.
 441 !   tauaer300, tauaer400, tauaer600, tauaer999: Layer AODs
 442 !   waer300, waer400, waer600, waer999: Layer SSAs
 443 !   gaer300, gaer400, gaer600, gaer999: Layer asymmetry parameters
 444
 445 ! OUTPUT:
 446 ! dtaer: Layer AOD at FTUV wavelengths
 447 ! omaer: Layer SSA at FTUV wavelengths
 448 ! gaer : Layer asymmetry parameters at FTUV wavelengths
 449 !------------------------------------------------------------------------
 450
 451 !-----------------------------------------------------------------------------
 452 ! Dummy arguments
 453 !-----------------------------------------------------------------------------
 454       integer, intent(in) :: nlambda_start
 455       real, intent(in)  :: wc(:)
 456       real, intent(in)  :: tauaer300(:), tauaer400(:),    &
 457                            tauaer600(:), tauaer999(:)
 458       real, intent(in)  :: waer300(:), waer400(:),        &
 459                            waer600(:), waer999(:)
 460       real, intent(in)  :: gaer300(:), gaer400(:),        &
 461                            gaer600(:), gaer999(:)
 462       real, intent(out) :: dtaer(:,:), omaer(:,:), gaer(:,:)
 463
 464 !-----------------------------------------------------------------------------
 465 ! Local Variables
 466 !-----------------------------------------------------------------------------
 467       real, parameter :: thresh = 1.e-9
 468       integer     :: k, wn, nlyr, nwave
 469       real        :: ang, slope, wfac
 470
 471       nlyr =  size(dtaer,dim=1)
 472       nwave = size(dtaer,dim=2)
 473
 474 wave_loop: &
 475       do wn = nlambda_start,nwave
 476         wfac = wc(wn)*1.e-3 - .6
 477         do k = 1,nlyr-1
 478 !-----------------------------------------------------------------------------
 479 ! use angstrom exponent to calculate aerosol optical depth; wc is in nm.
 480 !-----------------------------------------------------------------------------
 481           if( tauaer300(k) > thresh .and. tauaer999(k) > thresh ) then
 482             ang = log(tauaer300(k)/tauaer999(k))/log(0.999/0.3)
 483             dtaer(k,wn) = tauaer400(k)*(0.4/(wc(wn)*1.e-3))**ang
 484 !-----------------------------------------------------------------------------
 485 ! ssa - use linear interpolation/extrapolation
 486 !-----------------------------------------------------------------------------
 487             slope = 5.*(waer600(k) - waer400(k))
 488             omaer(k,wn) = slope*wfac + waer600(k)
 489             omaer(k,wn) = max( .4,min( 1.,omaer(k,wn) ) )
 490 !-----------------------------------------------------------------------------
 491 ! asymmetry parameter - use linear interpolation/extrapolation
 492 !-----------------------------------------------------------------------------
 493             slope = 5.*(gaer600(k) - gaer400(k))
 494             gaer(k,wn) = slope*wfac + gaer600(k)
 495             gaer(k,wn) = max( .5,min( 1.,gaer(k,wn) ) )
 496           endif
 497         end do
 498       end do wave_loop
 499
 500       end subroutine setaer
 501
 502       subroutine setcld( nlambda_start, cld_od_opt, z, xlwc, cldfrac, &
 503                          dtcld, omcld, gcld )
 504 !-----------------------------------------------------------------------------
 505 != PURPOSE:
 506 != Set up cloud optical depth, single albedo and g
 507 !-----------------------------------------------------------------------------
 508 != PARAMETERS:
 509 != PARAMETERS:
 510 != NZ - INTEGER, number of specified altitude levels in the working (I)
 511 != grid
 512 != Z - real(dp), specified altitude working grid (km) (I)
 513 != XLWC Cloud water content g/M3 (I)
 514 !=
 515 != dtcld - cloud optical depth
 516 != omcld - cloud droplet single albedo
 517 != gcld  - g
 518 !-----------------------------------------------------------------------------
 519 !
 520 ! VERTICAL DOMAIN is from bottom(1) to TOP (TOP=nz)
 521 ! CCM from top(1) to bottom(nz)
 522 !-----------------------------------------------------------------------------
 523
 524 !-----------------------------------------------------------------------------
 525 ! ... Dummy arguments
 526 !-----------------------------------------------------------------------------
 527       integer, intent(in) :: nlambda_start
 528       integer, intent(in) :: cld_od_opt
 529       real, intent(in)  :: z(:)
 530       real, intent(in)  :: xlwc(:)
 531       real, intent(in)  :: cldfrac(:)
 532       real, intent(inout) :: dtcld(:,:)
 533       real, intent(inout) :: omcld(:,:)
 534       real, intent(inout) :: gcld(:,:)
 535 !-----------------------------------------------------------------------------
 536 ! ... Local variables
 537 !-----------------------------------------------------------------------------
 538       real, parameter :: km2m = 1.e3        ! kilometer to meter
 539       real, parameter :: wden = 1.e6        ! g/m3 (1 m3 water = 1e6 g water)
 540       real, parameter :: re = 10.0 * 1.e-6  ! assuming cloud drop radius = 10 um to M
 541       real, parameter :: fac = 1./(wden*re)
 542
 543       integer  :: astat
 544       integer  :: wn
 545       integer  :: nlyr, nwave
 546       real, allocatable :: wrk(:), layer_cldfrac(:)
 547
 548       nlyr  = size(dtcld,dim=1)
 549       nwave = size(dtcld,dim=2)
 550
 551       allocate( wrk(nlyr),layer_cldfrac(nlyr),stat=astat )
 552       if( astat /= 0 ) then
 553         call wrf_error_fatal( 'setcld: failed to allocate wrk' )
 554       endif
 555
 556 !-----------------------------------------------------------------------------
 557 ! ... calculate optical depth
 558 !-----------------------------------------------------------------------------
 559       wrk(1:nlyr-1) = (z(2:nlyr) - z(1:nlyr-1))*km2m   !  (km -> m)
 560       wrk(1:nlyr-1) = 1.5 * .5*(xlwc(1:nlyr-1) + xlwc(2:nlyr))*wrk(1:nlyr-1)*fac
 561       wrk(1:nlyr-1) = max( wrk(1:nlyr-1),0. )
 562       if( cld_od_opt == 2 ) then
 563         layer_cldfrac(1:nlyr-1) = .5*(cldfrac(1:nlyr-1) + cldfrac(2:nlyr))
 564         wrk(1:nlyr-1) = wrk(1:nlyr-1)*layer_cldfrac(1:nlyr-1)*sqrt( layer_cldfrac(1:nlyr-1) )
 565       endif
 566 !----------------------------------------------------
 567 ! ....calculate cloud optical depth T
 568 ! following Liao et al. JGR, 104, 23697, 1999
 569 !----------------------------------------------------
 570       if( any( wrk(1:nlyr-1) > 0. ) ) then
 571         do wn = nlambda_start,nwave
 572           dtcld(1:nlyr-1,wn) = wrk(1:nlyr-1)
 573           omcld(1:nlyr-1,wn) = .9999
 574           gcld (1:nlyr-1,wn) = .85
 575         end do
 576       endif
 577
 578       if( allocated( wrk ) ) then
 579         deallocate( wrk )
 580       endif
 581       if( allocated( layer_cldfrac ) ) then
 582         deallocate( layer_cldfrac )
 583       endif
 584
 585       end subroutine setcld
 586
 587       REAL FUNCTION sundis( julday )
 588 !-----------------------------------------------------------------------------
 589 ! purpose:
 590 ! calculate earth-sun distance variation for a given date. based on
 591 ! fourier coefficients originally from: spencer, j.w., 1971, fourier
 592 ! series representation of the position of the sun, search, 2:172
 593 !-----------------------------------------------------------------------------
 594 ! parameters:
 595 ! idate - integer, specification of the date, from yymmdd (i)
 596 ! esrm2 - real(dp), variation of the earth-sun distance (o)
 597 ! esrm2 = (average e/s dist)^2 / (e/s dist on day idate)^2
 598 !-----------------------------------------------------------------------------
 599
 600       implicit none
 601
 602 !-----------------------------------------------------------------------------
 603 ! ... dummy arguments
 604 !-----------------------------------------------------------------------------
 605       integer, intent(in) :: julday
 606
 607 !-----------------------------------------------------------------------------
 608 ! ... local variables
 609 !-----------------------------------------------------------------------------
 610       real(dp), parameter :: pi = 3.1415926_dp
 611
 612       real(dp) :: dayn, thet0
 613       real(dp) :: sinth, costh, sin2th, cos2th
 614
 615 !-----------------------------------------------------------------------------
 616 ! ... parse date to find day number (julian day)
 617 !-----------------------------------------------------------------------------
 618       dayn = real(julday - 1,kind=dp) + .5_dp
 619 !-----------------------------------------------------------------------------
 620 ! ... define angular day number and compute esrm2:
 621 !-----------------------------------------------------------------------------
 622       thet0 = 2._dp*pi*dayn/365._dp
 623 !-----------------------------------------------------------------------------
 624 ! ... calculate sin(2*thet0), cos(2*thet0) from
 625 ! addition theorems for trig functions for better
 626 ! performance; the computation of sin2th, cos2th
 627 ! is about 5-6 times faster than the evaluation
 628 ! of the intrinsic functions sin and cos
 629 !-----------------------------------------------------------------------------
 630       sinth  = sin( thet0 )
 631       costh  = cos( thet0 )
 632       sin2th = 2._dp*sinth*costh
 633       cos2th = costh*costh - sinth*sinth
 634       sundis  = real( 1.000110_dp + .034221_dp*costh + .001280_dp*sinth &
 635                                   + .000719_dp*cos2th + .000077_dp*sin2th )
 636
 637       END FUNCTION sundis
 638
 639       subroutine calc_zenith( lat, long, julday, gmt, zenith, &
 640                               its, ite, jts, jte, &
 641                               ims, ime, jms, jme )
 642 !-------------------------------------------------------------------
 643 ! this subroutine calculates solar zenith and azimuth angles for a
 644 ! part  time and location.  must specify:
 645 ! input:
 646 ! lat - latitude in decimal degrees
 647 ! long - longitude in decimal degrees
 648 ! gmt  - greenwich mean time - decimal military eg.
 649 ! 22.75 = 45 min after ten pm gmt
 650 ! output:
 651 ! zenith
 652 ! azimuth
 653 ! .. Scalar Arguments ..
 654 !-------------------------------------------------------------------
 655         integer, intent(in)  :: julday
 656         integer, intent(in)  :: its,ite
 657         integer, intent(in)  :: jts,jte
 658         integer, intent(in)  :: ims,ime
 659         integer, intent(in)  :: jms,jme
 660         real(dp), intent(in) :: gmt
 661         real, intent(in)     :: lat(ims:ime,jms:jme)
 662         real, intent(in)     :: long(ims:ime,jms:jme)
 663         real, intent(out)    :: zenith(ims:ime,jms:jme)
 664
 665 !-------------------------------------------------------------------
 666 ! .. Local variables
 667 !-------------------------------------------------------------------
 668       real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp
 669       real(dp), parameter :: r2d = 1.0_dp/d2r
 670
 671       integer  :: i, j
 672       real(dp) :: caz, csz, cw, d, ec, epsi, eqt, eyt, feqt, feqt1, &
 673           feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi, &
 674           pi, ra, raz, rdecl, reqt, rlt, rml, rra, ssw, sw, tab, w, wr, &
 675           yt, zpt, zr
 676
 677       d = real(julday,dp) + gmt/24.0_dp
 678 !-------------------------------------------------------------------
 679 ! calc geom mean longitude
 680 !-------------------------------------------------------------------
 681       ml  = 279.2801988_dp + d*(.9856473354_dp + 2.267E-13_dp*d)
 682       rml = ml*d2r
 683 !-------------------------------------------------------------------
 684 ! calc equation of time in sec
 685 ! w = mean long of perigee
 686 ! e = eccentricity
 687 ! epsi = mean obliquity of ecliptic
 688 !-------------------------------------------------------------------
 689       w = 282.4932328_dp + d*(4.70684E-5_dp + 3.39E-13_dp*d)
 690       wr = w*d2r
 691       ec   = 1.6720041E-2_dp - d*(1.1444E-9_dp + 9.4E-17_dp*d)
 692       epsi = 23.44266511_dp - d*(3.5626E-7_dp + 1.23E-15_dp*d)
 693       pepsi = epsi*d2r
 694       yt = (tan(pepsi/2.0_dp))**2
 695       cw = cos(wr)
 696       sw = sin(wr)
 697       ssw = sin(2.0_dp*wr)
 698       eyt = 2._dp*ec*yt
 699       feqt1 = -sin(rml)*cw*(eyt + 2._dp*ec)
 700       feqt2 = cos(rml)*sw*(2._dp*ec - eyt)
 701       feqt3 = sin(2._dp*rml)*(yt - (5._dp*ec**2/4._dp)*(cw**2 - sw**2))
 702       feqt4 = cos(2._dp*rml)*(5._dp*ec**2*ssw/4._dp)
 703       feqt5 = sin(3._dp*rml)*(eyt*cw)
 704       feqt6 = -cos(3._dp*rml)*(eyt*sw)
 705       feqt7 = -sin(4._dp*rml)*(.5_dp*yt**2)
 706       feqt  = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7
 707       eqt   = feqt*13751.0_dp
 708
 709 !-------------------------------------------------------------------
 710 ! convert eq of time from sec to deg
 711 !-------------------------------------------------------------------
 712       reqt = eqt/240._dp
 713 !-------------------------------------------------------------------
 714 ! calc right ascension in rads
 715 !-------------------------------------------------------------------
 716       ra  = ml - reqt
 717       rra = ra*d2r
 718 !-------------------------------------------------------------------
 719 ! calc declination in rads, deg
 720 !-------------------------------------------------------------------
 721       tab = 0.43360_dp*sin(rra)
 722       rdecl = atan(tab)
 723       do j = jts,jte
 724         do i = its,ite
 725 !-------------------------------------------------------------------
 726 ! calc local hour angle
 727 !-------------------------------------------------------------------
 728           lbgmt = 12.0_dp - eqt/3600._dp + real(long(i,j),dp)*24._dp/360._dp
 729           lzgmt = 15.0_dp*(gmt - lbgmt)
 730           zpt   = lzgmt*d2r
 731           rlt   = real(lat(i,j),dp)*d2r
 732           csz   = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
 733           csz   = min( 1._dp,csz )
 734           zr    = acos(csz)
 735           zenith(i,j) = real( zr/d2r,4 )
 736         end do
 737       end do
 738
 739       end subroutine calc_zenith
 740
 741       SUBROUTINE sphers( nlyr, z, zen, dsdh, nid )
 742 !-----------------------------------------------------------------------------
 743 !=  PURPOSE:
 744 !=  Calculate slant path over vertical depth ds/dh in spherical geometry.
 745 !=  Calculation is based on:  A.Dahlback, and K.Stamnes, A new spheric model
 746 !=  for computing the radiation field available for photolysis and heating
 747 !=  at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B)
 748 !-----------------------------------------------------------------------------
 749 !=  PARAMETERS:
 750 !=  NZ      - INTEGER, number of specified altitude levels in the working (I)
 751 !=            grid
 752 !=  Z       - REAL, specified altitude working grid (km)                  (I)
 753 !=  ZEN     - REAL, solar zenith angle (degrees)                          (I)
 754 !=  DSDH    - REAL, slant path of direct beam through each layer crossed  (O)
 755 !=            when travelling from the top of the atmosphere to layer i;
 756 !=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1
 757 !=  NID     - INTEGER, number of layers crossed by the direct beam when   (O)
 758 !=            travelling from the top of the atmosphere to layer i;
 759 !=            NID(i), i = 0..NZ-1
 760 !-----------------------------------------------------------------------------
 761 !  This program is free software;  you can redistribute it and/or modify
 762 !  it under the terms of the GNU General Public License as published by the
 763 != Free Software Foundation;  either version 2 of the license, or (at your
 764 != option) any later version.
 765 != The TUV package is distributed in the hope that it will be useful, but
 766 != WITHOUT ANY WARRANTY;  without even the implied warranty of MERCHANTIBI-
 767 != LITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
 768 != License for more details.
 769 != To obtain a copy of the GNU General Public License, write to:
 770 != Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 771 !-----------------------------------------------------------------------------
 772 != To contact the authors, please mail to:
 773 != Sasha Madronich, NCAR/ACD, P.O.Box 3000, Boulder, CO, 80307-3000, USA
 774 != send email to:  sasha@ucar.edu
 775 !-----------------------------------------------------------------------------
 776
 777       use params_mod, only : pi, radius
 778
 779       IMPLICIT NONE
 780
 781 !-----------------------------------------------------------------------------
 782 !     ... dummy arguments
 783 !-----------------------------------------------------------------------------
 784       INTEGER, intent(in) :: nlyr
 785       REAL,    intent(in) :: zen
 786       REAL,    intent(in) :: z(:)
 787
 788       INTEGER, intent(inout) :: nid(0:nlyr)
 789       REAL,    intent(inout) :: dsdh(0:nlyr,nlyr)
 790
 791
 792 !-----------------------------------------------------------------------------
 793 !     ... local variables
 794 !-----------------------------------------------------------------------------
 795       REAL, parameter ::  dr = pi/180.
 796
 797       INTEGER :: j, jm1, k
 798       INTEGER :: id
 799       REAL    :: re
 800       REAL    :: zd(0:nlyr)
 801       REAL    :: ds_dh(1:nlyr)
 802       REAL(8) :: zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
 803
 804       zenrad = REAL( zen*dr,8 )
 805 !-----------------------------------------------------------------------------
 806 ! include the elevation above sea level to the radius of the earth:
 807 !-----------------------------------------------------------------------------
 808       re = radius + z(1)
 809
 810 !-----------------------------------------------------------------------------
 811 ! from bottom-up to top-down
 812 ! note zd is the elevation above earth surface:
 813 !-----------------------------------------------------------------------------
 814       zd(0:nlyr) = z(nlyr+1:1:-1) - z(1)
 815
 816 !-----------------------------------------------------------------------------
 817 ! initialize nid
 818 !-----------------------------------------------------------------------------
 819       nid(0:nlyr) = 0
 820
 821 !-----------------------------------------------------------------------------
 822 ! calculate ds/dh of every layer
 823 !-----------------------------------------------------------------------------
 824 layer_loop : &
 825       DO k = 0, nlyr
 826         ds_dh(:) = 0.
 827         rpsinz = real(re + zd(k),8) * SIN(zenrad)
 828 !       IF( zen > 90.0 .AND. rpsinz < real(re,8) ) THEN
 829         IF( zen <= 90.0 .or. rpsinz >= real(re,8) ) THEN
 830 !-----------------------------------------------------------------------------
 831 ! Find index of layer in which the screening height lies
 832 !-----------------------------------------------------------------------------
 833           id = k
 834           IF( zen > 90.0 ) THEN
 835             DO j = 1,nlyr
 836               IF( rpsinz < real(zd(j-1) + re,8) .AND. &
 837                   rpsinz >= real(zd(j) + re,8) ) then
 838                 id = j
 839               ENDIF
 840             END DO
 841           END IF
 842
 843           DO j = 1, id
 844             jm1 = j - 1
 845 !           IF( j == id .AND. id == k .AND. zen > 90.0 ) then
 846             IF( j /= id .or. k /= id .or. zen <= 90.0 ) then
 847               sm = 1.0_8
 848             ELSE
 849               sm = -1.0_8
 850             ENDIF
 851             rj   = real(re + zd(jm1),8)
 852             rjp1 = real(re + zd(j),8)
 853             dhj  = zd(jm1) - zd(j)
 854
 855             ga = max( rj*rj - rpsinz*rpsinz,0.0_8 )
 856             gb = max( rjp1*rjp1 - rpsinz*rpsinz,0.0_8 )
 857
 858             IF( id > k .AND. j == id ) THEN
 859               dsj = SQRT( ga )
 860             ELSE
 861               dsj = SQRT( ga ) - sm*SQRT( gb )
 862             END IF
 863             ds_dh(j) = real( dsj/dhj,4 )
 864           END DO
 865           nid(k) = id
 866         ELSE
 867           nid(k) = -1
 868         ENDIF
 869         dsdh(k,:) = ds_dh(:)
 870
 871       END DO layer_loop
 872
 873       END SUBROUTINE sphers
 874
 875       SUBROUTINE airmas( nlyr, dsdh, nid, cz, vcol, scol )
 876 !-----------------------------------------------------------------------------
 877 !=  PURPOSE:
 878 !=  Calculate vertical and slant air columns, in spherical geometry, as a
 879 !=  function of altitude.
 880 !-----------------------------------------------------------------------------
 881 !=  PARAMETERS:
 882 !=  NZ      - INTEGER, number of specified altitude levels in the working (I)
 883 !=            grid
 884 !=  DSDH    - REAL, slant path of direct beam through each layer crossed  (O)
 885 !=            when travelling from the top of the atmosphere to layer i;
 886 !=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1
 887 !=  NID     - INTEGER, number of layers crossed by the direct beam when   (O)
 888 !=            travelling from the top of the atmosphere to layer i;
 889 !=            NID(i), i = 0..NZ-1
 890 !=  VCOL    - REAL, output, vertical air column, molec cm-2, above level iz
 891 !=  SCOL    - REAL, output, slant air column in direction of sun, above iz
 892 !=            also in molec cm-2
 893 !-----------------------------------------------------------------------------
 894
 895       use params_mod, only : largest
 896
 897       IMPLICIT NONE
 898
 899 !-----------------------------------------------------------------------------
 900 !     ... dummy arguments
 901 !-----------------------------------------------------------------------------
 902       INTEGER, intent(in) :: nlyr
 903       INTEGER, intent(in) :: nid(0:nlyr)
 904       REAL,    intent(in) :: dsdh(0:nlyr,nlyr)
 905       REAL,    intent(in) :: cz(nlyr)
 906
 907       REAL, intent(inout) :: vcol(nlyr), scol(nlyr)
 908
 909 !-----------------------------------------------------------------------------
 910 !     ... local variables
 911 !-----------------------------------------------------------------------------
 912       INTEGER :: lyr, j, nlev, nlevi
 913       REAL    :: sum, vsum
 914
 915 !-----------------------------------------------------------------------------
 916 ! calculate vertical and slant column from each level: work downward
 917 !-----------------------------------------------------------------------------
 918       nlev = nlyr + 1
 919       vsum = 0.
 920       DO lyr = 1, nlyr
 921         nlevi = nlev - lyr
 922         vsum = vsum + cz(nlevi)
 923         vcol(nlevi) = vsum
 924         sum = 0.
 925         IF( nid(lyr) < 0 ) THEN
 926           sum = largest
 927         ELSE
 928 !-----------------------------------------------------------------------------
 929 ! single pass layers:
 930 !-----------------------------------------------------------------------------
 931           DO j = 1, MIN(nid(lyr), lyr)
 932             sum = sum + cz(nlev-j)*dsdh(lyr,j)
 933           END DO
 934 !-----------------------------------------------------------------------------
 935 ! double pass layers:
 936 !-----------------------------------------------------------------------------
 937            DO j = MIN(nid(lyr),lyr)+1, nid(lyr)
 938              sum = sum + 2.*cz(nlev-j)*dsdh(lyr,j)
 939            END DO
 940         ENDIF
 941         scol(nlevi) = sum
 942       END DO
 943
 944       END SUBROUTINE airmas
 945
 946       END MODULE tuv_subs