!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE RADIATION                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE radiation(nx,ny,nz,rbufsz,                                   & 2,3
           ptprt,pprt,qv,qc,qr,qi,qs,qh,                                &
           ptbar,pbar,ppi,rhostr,                                       &
           x,y,z,zp,j3inv,                                              &
           soiltyp, tsfc, wetsfc,snowdpth,                              &
           radfrc, radsw, rnflx,                                        &
           rsirbm,rsirdf,rsuvbm,rsuvdf,                                 &
           cosz, cosss,                                                 &
           fdirir,fdifir,fdirpar,fdifpar,                               &
           tem1, tem2, tem3, tem4, tem5,                                &
           tem6, tem7, tem8, tem9, tem10,                               &
           tem11,tem12,tem13,tem14,tem15,tem16,                         &
           radbuf)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine is the main driver of ARPS radiation package. To
!  control the driver, the parameters and variables should be set:
!
!  Input control variables (in namelist &radiation, arps.input)
!
!  radopt    Option to switch on/off radiation physics
!       = 0, No radiation physics;
!       = 1, Simplified surface radiation physics;
!       = 2, Atmospheric radiation transfer parameterization.
!
!  Notes:    When sfcphy is chosen to 3 or 4, radopt=0 will be adjusted
!         to 1 in order to compute the surface energy balance for
!         soil model. However radopt=2 will not be changed.
!
!  radstgr   Option for radiation computing at staggering points
!       = 0, No staggering; Radiation calculation on x-y plane is at
!            all points;
!       = 1, staggering; Radiation calculation on x-y plane is at
!            (even,even) and (odd,odd) points. The values at
!            (even,odd) (odd,even) points are averaged from the
!            surrounding four points. For example for nx=ny=9, the
!            directly calculation are performed at the "x" points,
!            then calculate radiation variables at "o" by averaging
!            from their surrounding "x" points. This scheme can
!            reduce ALMOST HALF of radiation calculation.
!
!
!                      j
!
!                    9 | x o x o x o x o x
!                    8 | o x o x o x o x o
!                    7 | x o x o x o x o x
!                    6 | o x o x o x o x o
!                    5 | x o x o x o x o x
!                    4 | o x o x o x o x o
!                    3 | x o x o x o x o x
!                    2 | o x o x o x o x o
!                    1 | x o x o x o x o x
!                      +-------------------  i
!                        1 2 3 4 5 6 7 8 9
!
!            On boundary, the zero-gradient is assumed.
!
!  rlwopt    Option to choose the longwave schemes.
!       = 0, high = .false. in code, transmission functions are
!            computed using the k-distribution method with linear
!            pressure scaling.  cooling rates are not calculated
!            accurately for pressures less than 20 mb. The
!            computation is faster with high=.false. than with
!            high=.true.
!       = 1, high = .true. in code, transmission functions in the
!            co2, o3 in the co2, o3, and the three water vapor bands
!            with strong absorption are computed using table look-up.
!            cooling rates are computed accurately from the surface
!            up to 0.01 mb.
!
!  dtrad     Time interval (seconds) to update the radiation forcing
!
!  raddiag   Option to dump radiation variables to a file in GrADS
!            format for diagnostic review. The frequency is
!            controled by dtrad. (Effective when radopt=2)
!       = 0, no such dump
!       = 1, dump to a file with a name like 'runname.radout'
!            and its control file has a name like 'runname.radctl'
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  03/11/1996
!
!  MODIFICATION HISTORY:
!
!  4/13/98 (D.Weber and V. Wong)
!  Added precipitable water term in the radopt=1 air emissivity
!  coefficient.
!
!  11/18/98 (Keith Brewster)
!  Changed pibar to ppi (full pi).
!
!  12/8/1998 (Donghai Wang and Vince Wong)
!  Added a new 2-D permanent array, snowcvr(nx,ny), for snow cover.
!  We just used a simple scheme to consider the snow cover process.
!
!  2000/01/10 (Gene Bassett)
!  Snow cover (0 or 1) changed to a fractional value (0 to 1)
!  determined from snow depth.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz
  INTEGER :: rbufsz
!
!-----------------------------------------------------------------------
!
!  Define ARPS variables
!
!-----------------------------------------------------------------------
!
  REAL :: ptprt (nx,ny,nz)
  REAL :: pprt  (nx,ny,nz)
  REAL :: qv    (nx,ny,nz)
  REAL :: qc    (nx,ny,nz)
  REAL :: qr    (nx,ny,nz)
  REAL :: qi    (nx,ny,nz)
  REAL :: qs    (nx,ny,nz)
  REAL :: qh    (nx,ny,nz)

  REAL :: ptbar (nx,ny,nz)
  REAL :: pbar  (nx,ny,nz)
  REAL :: ppi   (nx,ny,nz)
  REAL :: rhostr(nx,ny,nz)

  REAL :: x      (nx)
  REAL :: y      (ny)
  REAL :: z      (nz)

  REAL :: zp    (nx,ny,nz)  ! The physical height coordinate defined at
                            ! w-point of staggered grid.
  REAL :: j3inv (nx,ny,nz)

  INTEGER :: soiltyp(nx,ny) ! Soil type at each point
  REAL :: tsfc   (nx,ny)
  REAL :: wetsfc (nx,ny)    ! Surface soil moisture in the top 1 cm layer
  REAL :: snowdpth(nx,ny)   ! Snow depth (m)

  REAL :: radfrc(nx,ny,nz)  ! Radiation forcing (K/s)
  REAL :: radsw  (nx,ny)    ! Solar radiation down to the surface
  REAL :: rnflx  (nx,ny)    ! Net radiation flux absorbed by surface

  REAL :: rsirbm(nx,ny)     ! Solar IR surface albedo for beam radiation
  REAL :: rsirdf(nx,ny)     ! Solar IR surface albedo for diffuse radiation
  REAL :: rsuvbm(nx,ny)     ! Solar UV surface albedo for beam radiation
  REAL :: rsuvdf(nx,ny)     ! Solar UV surface albedo for diffuse radiation

  REAL :: cosz  (nx,ny)     ! Cosine of zenith
  REAL :: cosss (nx,ny)     ! Cosine of angle between sun light and
                            ! surface terrain slope

  REAL :: fdirir (nx,ny)    ! all-sky direct downward IR flux
                            ! (0.7-10 micron) at the surface
  REAL :: fdifir (nx,ny)    ! all-sky diffuse downward IR flux
                            ! at the surface
  REAL :: fdirpar(nx,ny)    ! all-sky direct downward par flux
                            ! (0.4-0.7 micron) at the surface
  REAL :: fdifpar(nx,ny)    ! all-sky diffuse downward par flux
                            ! at the surface
  REAL :: radbuf(rbufsz)    ! temporary arrays used for radiation
                            ! transfer computing
!
!-----------------------------------------------------------------------
!
!  Temporary arrays which have the vertical coordinate inversed, that
!  is, k=1 is for top while k=nz is at the surface.
!
!-----------------------------------------------------------------------
!
  REAL :: tem1 (nx,ny,nz)    ! pinv, slpmag, slpdir
  REAL :: tem2 (nx,ny,nz)    ! tinv
  REAL :: tem3 (nx,ny,nz)    ! qvinv
  REAL :: tem4 (nx,ny,nz)    ! o3a
  REAL :: tem5 (nx,ny,nz)    ! ccld
  REAL :: tem6 (nx,ny,nz)    ! flxir
  REAL :: tem7 (nx,ny,nz)    ! flcir
  REAL :: tem8 (nx,ny,nz)    ! flxuv
  REAL :: tem9 (nx,ny,nz)    ! flcuv
  REAL :: tem10(nx,ny,nz)    ! dfdts
  REAL :: tem11(nx,ny,nz)    ! tauir
  REAL :: tem12(nx,ny,nz)    ! taual
  REAL :: tem13(nx,ny,nz)    ! tauswi
  REAL :: tem14(nx,ny,nz)    ! tauswl
  REAL :: tem15(nx,ny,nz)    ! reffi
  REAL :: tem16(nx,ny,nz)    ! reffl
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'radcst.inc'
  INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
!  Local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k

  REAL :: albedo, albedoz

  REAL :: tema,temb
  REAL :: rad2deg

  REAL :: frac_snowcover
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF ( radopt == 1 .AND. ( sfcphy /= 3 .AND. sfcphy /= 4 ) ) THEN
    RETURN
  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate solar zenith angle. Array radsw is used as a2dr2
!  temporarily.
!
!    zp(*,*,2) = surface terrain
!    cosz    = cosine of zenith
!    cosss   = cosine of angle between sun light and terrain slope
!
!    radsw = a2dr2, square ratio of average distance to the time
!                   dependent distance from the earth to the sun
!
!  Array tem1 and tem2 are used as temporary 2-D arrays
!
!    tem1(*,*,1) = rjday
!    tem1(*,*,2) = tloc
!    tem1(*,*,3) = latscl
!    tem1(*,*,4) = lonscl
!
!    tem2(*,*,1) = slpmag
!    tem2(*,*,2) = slpdir
!    tem2(*,*,3) = tem1(*,*)
!    tem2(*,*,4) = tem2(*,*)
!
!-----------------------------------------------------------------------
!
  CALL zenangl( nx,ny, x,y, zp(1,1,2), cosz, cosss, radsw,              &
                tem1(1,1,1),tem1(1,1,2),tem1(1,1,3),tem1(1,1,4),        &
                tem2(1,1,1),tem2(1,1,2),tem2(1,1,3),tem2(1,1,4) )
!
!-----------------------------------------------------------------------
!
!  Calculate surface albedo which is dependent on solar zenith angle
!  and soil moisture. Set the albedo for different types of solar
!  flux to be same.
!
!    rsirbm   Solar IR surface albedo for beam radiation
!    rsirdf   Solar IR surface albedo for diffuse radiation
!    rsuvbm   Solar UV surface albedo for beam radiation
!    rsuvdf   Solar UV surface albedo for diffuse radiation
!
!-----------------------------------------------------------------------
!
  rad2deg = 180.0/3.141592654

  DO j=1,ny-1
    DO i=1,nx-1

      albedoz = 0.01 * ( EXP( 0.003286         & ! zenith dependent albedo
          * SQRT( ( ACOS(cosz(i,j))*rad2deg ) ** 3 ) ) - 1.0 )

      IF ( sfcphy == 0 ) THEN             ! soil type not defined
        tema = 0
      ELSE
        tema = wetsfc(i,j)/wsat(soiltyp(i,j))
      END IF

      frac_snowcover = MIN(snowdpth(i,j)/snowdepth_crit, 1.0)

      IF ( tema > 0.5 ) THEN
        albedo = albedoz + (1.-frac_snowcover)*0.14                     &
                         + frac_snowcover*snow_albedo
      ELSE
        albedo = albedoz + (1.-frac_snowcover)*(0.31 - 0.34 * tema)     &
                         + frac_snowcover*snow_albedo
      END IF

      rsirbm(i,j) = albedo
      rsirdf(i,j) = albedo
      rsuvbm(i,j) = albedo
      rsuvdf(i,j) = albedo

    END DO
  END DO

  IF ( radopt == 1 ) THEN

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k) = (ptbar(i,j,k)+ptprt(i,j,k))*ppi(i,j,k)
          tem3(i,j,k) = pbar(i,j,k) + pprt(i,j,k)
        END DO
      END DO
    END DO

!
!-----------------------------------------------------------------------
!
!  Calculate solar radiation at the surface using a simplified
!  scheme.
!
!  tem1(*,*,1) = prcpln, precipitation path length
!
!  tem2        = temperature
!  tem3        = total pressure
!
!  fdirir      = fraction of solar radiation reaching the surface
!
!  Note: added precipitable water contribution to the emissa term
!        for radopt=1.
!
!-----------------------------------------------------------------------
!
    CALL solrsfc( nx,ny,nz,                                             &
                  tem2, tem3, qv, cosz, tem1(1,1,1), fdirir )

    DO j=1,ny-1
      DO i=1,nx-1
        tema = 0.5 * ( tem2(i,j,1) + tem2(i,j,2) ) ! surface air temp.

!
!  new code to include the water part of the downward emissa
!  coefficient....
!
        temb = 0.0

        DO k = 2,nz-2
          temb = temb+rhostr(i,j,k)*qv(i,j,k)*(zp(i,j,k+1)-zp(i,j,k))
        END DO

        temb = 0.17*ALOG10(temb*0.1)
        temb = MIN(1.0,temb+emissa)

!  end of new code for emissa..

        radsw(i,j) = solarc * radsw(i,j) * cosss(i,j) * fdirir(i,j)

        rnflx(i,j) = radsw(i,j) * (1.0-rsirbm(i,j))                     &
                   + temb   * sbcst * tema**4                           &
                   - emissg * sbcst * tsfc(i,j)**4
      END DO
    END DO

  ELSE IF ( radopt == 2 ) THEN
!
!-----------------------------------------------------------------------
!
!  Compute the atmospheric radiation forcing
!
!  Temporary arrays are used for
!
!    tem1  = pinv,  Pressure in mb at scalar points
!    tem2  = tinv,  Temperature
!    tem3  = qvinv, Water vapor mixing ratio (g/g)
!    tem4  = o3a,   Ozone (o3) mixing ratio (g/g)
!    tem5  = ccld,  Cloud coverage (fraction)
!    tem6  = flxir, all-sky net downward LW flux
!    tem7  = flcir, clear-sky net downward LW flux
!    tem8  = flxuv, all-sky solar flux (downward minus upward)
!    tem9  = flcuv, clear-sky solar flux (downward minus upward)
!    tem10 = dfdts, Sensitivity of net downward flux to surface temperature
!    tem11 = tauir,  Cloud optical depth for LW IR
!    tem12 = taual,  Aerosol optical thickness
!    tem13 = tauswi, Cloud optical depth for solar IR for ice particles
!    tem14 = tauswl, Cloud optical depth for solar IR for liquid particles
!    tem15 = reffi,  Effective cloud-particle size for ice particles
!    tem16 = reffl,  Effective cloud-particle size for liquid particles
!
!    cosss = cosine of angle between sun light and terrain slope
!
!    radsw = a2dr2, input,  square ratio of average distance to the
!                           time dependent distance from the earth to
!                           the sun
!          = radsw, output,    solar radiation reaching the surface
!    rnflx = st4,   temporary, longwave upward flux at surface
!          = rnflx, output,    net radiation flux absorbed by surface
!
!-----------------------------------------------------------------------
!
    CALL radtrns(nx,ny,nz, rbufsz,                                      &
                 ptprt,pprt,qv,qc,qr,qi,qs,qh,                          &
                 ptbar,pbar,ppi,rhostr, tsfc,                           &
                 x,y,z,zp, j3inv,                                       &
                 radfrc, radsw,rnflx, cosss,                            &
                 rsirbm,rsirdf,rsuvbm,rsuvdf, cosz,                     &
                 fdirir,fdifir,fdirpar,fdifpar, rnflx,                  &
                 tem1, tem2, tem3, tem4, tem5,                          &
                 tem6, tem7, tem8, tem9, tem10,                         &
                 tem11,tem12,tem13,tem14,tem15,tem16,                   &
                 radbuf)

  END IF

  RETURN

END SUBROUTINE radiation
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SOLRSFC                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE solrsfc(nx,ny,nz,                                            & 1
           temp, pres, qv, cosz, prcpln, trsw )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the atmosphere transmittance for solar radiation
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  03/25/1996
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz
!
!-----------------------------------------------------------------------
!
!  Define ARPS variables
!
!-----------------------------------------------------------------------
!
  REAL :: temp  (nx,ny,nz)     ! temperature
  REAL :: pres  (nx,ny,nz)     ! total pressure
  REAL :: qv    (nx,ny,nz)     ! Mixing ratio

  REAL :: cosz  (nx,ny)

  REAL :: prcpln(nx,ny)        ! Precipitation path length

  REAL :: trsw  (nx,ny)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'soilcst.inc'
!
!-----------------------------------------------------------------------
!
!  Local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k

  REAL :: dirf
  REAL :: trrg
  REAL :: trwv
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Calculate the precipitation path length by a vertical integral if
!  moist > 0.
!
!-----------------------------------------------------------------------
!
  DO j = 1, ny-1
    DO i = 1, nx-1
      prcpln(i,j) = 0.0
    END DO
  END DO

  IF ( moist /= 0 ) THEN
    DO j = 1, ny-1
      DO i = 1, nx-1

        prcpln(i,j) = 0.0

        DO k = 2, nz-2
          prcpln(i,j) = prcpln(i,j)                                     &
                      + 0.5*(qv(i,j,k)+qv(i,j,k+1))                     &
                      * 0.5*(pres(i,j,k)+pres(i,j,k+1))/101300.0        &
                      * SQRT(2.*273.16/(temp(i,j,k)+temp(i,j,k+1)))     &
                      * (pres(i,j,k)-pres(i,j,k+1))
        END DO

        prcpln(i,j) = prcpln(i,j) / g
      END DO
    END DO
  END IF

  DO j = 1, ny-1
    DO i = 1, nx-1
!
!-----------------------------------------------------------------------
!
!  Calculate the direction fractor of transmisivity
!
!-----------------------------------------------------------------------
!
      dirf = 35.0 / SQRT( 1224.0 * cosz(i,j)**2 + 1.0 )
!
!-----------------------------------------------------------------------
!
!  Calculate Rayleigh sccattering and absorption transmisivity, trrg.
!
!-----------------------------------------------------------------------
!
      trrg = 1.021 - 0.084 * SQRT( dirf                                 &
                           * ( 949. * 0.5*(pres(i,j,2)+pres(i,j,1))     &
                           * 1.e-8 + .051 ) )
                                ! Rayleigh scatt. and abs. transmission
! AB, Eq. 2, psfc in Pascal
!
!-----------------------------------------------------------------------
!
!  Calculate the transmisivity of water vapor, trwv. The
!  precipitation path length used in the calculation has been
!  calculated and stored in prcpln(i,j).
!
!  For cloud cover, use constant dirfc = 5/3 as the direction
!  fractor. (How to determine if cloudy or not?)
!
!  For clear sky, use dirf as the direction fractor, instead of dirfc
!
!-----------------------------------------------------------------------
!
      trwv = 1.0 - 2.9 * prcpln(i,j) * dirf                             &
           / ( 5.925 * prcpln(i,j) * dirf                               &
           + ( 1. + 141.5 * prcpln(i,j) * dirf ) ** 0.634 )
!
!-----------------------------------------------------------------------
!
!  Calculate the fraction of solar radiation reaching the surface.
!
!-----------------------------------------------------------------------
!
      trsw(i,j) = trrg * trwv

    END DO
  END DO

  RETURN
END SUBROUTINE solrsfc
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ZENANGL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE zenangl(nx,ny, x,y, hterain, cosz, cosss, a2dr2,             & 1,2
           rjday,tloc, latscl,lonscl, slpmag,slpdir,                    &
           tem1,tem2)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate cosine of solar zenith angle.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu and Vince Wong
!  11/16/93
!
!  MODIFICATION HISTORY:
!  2/2/99  Vince Wong and Jik Leong
!  This modification calculates the solar declination angle and
!  equation of time using a method found on page C24 of the
!  1996 Astronomical Almanac.
!  The mothod is good to 0.01 degrees in the sky over the
!  period 1950 to 2050.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!
!    x        X coordinates at scalar points
!    y        Y coordinates at scalar points
!    hterain  Surface terrain
!
!  OUTPUT:
!
!    cosz     Cosine of zenith
!    cosss    Cosine of angle between sun light and terrain slope
!    a2dr2    Square ratio of average distance to the time
!             dependent distance from the earth to the sun
!
!  WORK ARRAY:
!
!    rjday    Julian day at each grid point
!    tloc     Local time at each grid point
!    latscl   Latitudes  at scalar points
!    lonscl   Longitudes at scalar points
!    slpmag   Surface terrain slope magnitude
!    slpdir   Surface terrain slope direction
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny

  REAL :: x(nx)
  REAL :: y(ny)
  REAL :: hterain(nx,ny)

  REAL :: cosz(nx,ny),   cosss(nx,ny), a2dr2(nx,ny)
  REAL :: rjday(nx,ny),  tloc(nx,ny)
  REAL :: latscl(nx,ny), lonscl(nx,ny)
  REAL :: slpmag(nx,ny), slpdir(nx,ny)

  REAL :: tem1(nx,ny), tem2(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Include file:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j

  REAL :: xs, ys

  REAL :: hour0, yrday
  REAL :: deg2rad, pi, pi2

  REAL :: etau, shrangl, sdeclin
  REAL :: azimuth, sinz
  REAL :: dpsi, sinpsi, cospsi

  REAL :: anncyc

  REAL :: hr, days2k, lsun, gsun, obliq, lambda, xsun, ysun
  REAL :: asun, alpha, rad2deg

  LOGICAL :: firstcall        ! First call flag of this subroutine

  SAVE firstcall, hour0, pi, pi2, deg2rad, yrday, rad2deg
  DATA firstcall/.true./
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (firstcall) THEN
    pi  = 3.14159265358979
    pi2 = 2.0 * pi
    deg2rad = pi/180.0
    rad2deg = 1./deg2rad

    hour0 = FLOAT(hour)                                                 &
          + FLOAT(minute)/60.0                                          &
          + FLOAT(second)/3600.0

    IF ( MOD(year, 4) == 0 ) THEN
      yrday = 366.
    ELSE
      yrday = 365.
    END IF

    firstcall = .false.
  END IF

  CALL sfcslp( nx,ny, hterain, slpmag,slpdir, tem1,tem2 )

  IF ( mapproj == 0 ) THEN
    DO j=1,ny-1
      DO i=1,nx-1
        latscl(i,j) = ctrlat
        lonscl(i,j) = ctrlon
      END DO
    END DO
  ELSE
    DO j=1,ny-1
      ys = 0.5*(y(j)+y(j+1))
      DO i=1,nx-1
        xs = 0.5*(x(i)+x(i+1))
        CALL xytoll(1,1, xs,ys, latscl(i,j), lonscl(i,j))
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate the local time at each grid point. The
!  following formula is based on that the input time is the GMT
!  time at the reference grid point of the center.
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1

      latscl(i,j) = deg2rad * latscl(i,j)        ! lat: -90 to 90

      tloc(i,j) = hour0 + (curtim-dtbig)/3600.0                         &
                + lonscl(i,j)/15.0

      rjday(i,j) = jday + INT( tloc(i,j)/24.0 )

      tloc(i,j) = MOD( tloc(i,j), 24.0 )

      IF ( tloc(i,j) < 0. ) THEN
        tloc(i,j) = tloc(i,j) + 24.0            ! Local time
        rjday(i,j) = MOD( rjday(i,j)-1, yrday ) ! Julian day at each pts
      END IF

    END DO
  END DO

  DO j=1,ny-1
    DO i=1,nx-1
      anncyc = pi2 * ( rjday(i,j) - 1.0 ) / yrday

      a2dr2(i,j) = 1.000110                                             &
                 + 0.034221 * COS(anncyc)                               &
                 + 0.00128  * SIN(anncyc)                               &
                 + 0.000719 * COS(2.*anncyc)                            &
                 + 0.000077 * SIN(2.*anncyc)  ! PX, Eq. 17

      hr = hour + minute / 60.

!  days before (-ve) or after (+ve) 1/1/2000
      days2k = 367 * year - 7 * ( year + ( month + 9 ) / 12 ) / 4       &
              + 275 * month / 9 + day - 730531.5 + hr / 24.

      lsun = 280.461 + 0.9856474 * days2k     ! Mean Longitude of the Sun
      950     IF ( lsun < 0 ) THEN
        lsun = lsun + 360.
        GO TO 950
      ELSE IF ( lsun > 360 ) THEN
        lsun = lsun - 360.
        GO TO 950
      END IF

      gsun = 357.528 + 0.9856003 * days2k     ! Mean anomaly of the Sun
      960     IF ( gsun < 0 ) THEN
        gsun = gsun + 360.
        GO TO 960
      ELSE IF ( gsun > 360 ) THEN
        gsun = gsun - 360.
        GO TO 960
      END IF

      lambda = lsun + 1.915 * SIN(gsun*deg2rad)  & ! Ecliptic longitude
          + 0.02 * SIN(2*gsun*deg2rad)
      970     IF ( lambda < 0 ) THEN
        lambda = lambda + 360.
        GO TO 970
      ELSE IF ( lambda > 360 ) THEN
        lambda = lambda - 360.
        GO TO 970
      END IF

      obliq = 23.439 - 0.0000004 * days2k     ! Obliquity of the ecliptic

      xsun = COS(lambda*deg2rad)
      ysun = COS(obliq*deg2rad) * SIN(lambda*deg2rad)
      asun = ATAN(ysun/xsun)*rad2deg
      IF ( xsun < 0. ) THEN
        alpha = asun + 180   ! Right Ascension (RA)
      ELSE IF ( ( ysun < 0. ) .AND. ( xsun > 0. ) ) THEN
        alpha = asun + 360
      ELSE
        alpha = asun
      END IF

      etau = ( lsun - alpha ) * 4. / 60.      ! Equation of time in hour

!    etau = 0.158 * sin( pi*(rjday(i,j)+10.)/91.25 ) ! Equation of time
!    :       + 0.125 * sin( pi*rjday(i,j)/182.5 )       ! Wong, Eq. 8

      shrangl = 15.0 * deg2rad                            & ! Hour angle
          * ( tloc(i,j) + etau - 12.0)            ! Wong, Eq. 7

!    sdeclin = 23.5 * deg2rad
!    :          * cos( 2.0*pi*(rjday(i,j)-173.)/yrday ) ! Wong, Eq. 6
      sdeclin = ASIN(SIN(obliq*deg2rad)*SIN(lambda*deg2rad))
                                                ! Declination (in radian)

      cosz(i,j) = COS(latscl(i,j)) * COS(sdeclin) * COS(shrangl)        &
                + SIN(latscl(i,j)) * SIN(sdeclin)

!    print *, cos(latscl(i,j)),cos(sdeclin),cos(shrangl)
!    print *, sin(latscl(i,j)),sin(sdeclin)
!    print *,sdeclin,shrangl

      sinz = SIN ( ACOS(cosz(i,j)) )
!
!-----------------------------------------------------------------------
!
!  Consider the effects of the terrain slope on the solar radiation.
!  The slope magnitude and direction has been computed by subroutine
!  SFCSLP and passed in by slpmag and slpdir.
!
!-----------------------------------------------------------------------
!
      sinpsi = COS(sdeclin) * SIN(shrangl) / sinz

      cospsi = ( cosz(i,j) * SIN( latscl(i,j) ) - SIN( sdeclin ) )      &
             / ( sinz * COS( latscl(i,j) ) )

      azimuth = ATAN2( sinpsi, cospsi )

      dpsi = azimuth - slpdir(i,j)

      cosss(i,j) = COS( slpmag(i,j) ) * cosz(i,j)                       &
                 + SIN( slpmag(i,j) ) * sinz * COS( dpsi )

      cosz (i,j) = MAX( cosz (i,j), 0.0 )
      cosss(i,j) = MAX( cosss(i,j), 0.0 )

    END DO
  END DO

  RETURN
END SUBROUTINE zenangl
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE SFCSLP                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE sfcslp( nx,ny, hterain, slpmag,slpdir, tem1,tem2 ) 1,7
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This program is to compute the terrain slope vector in terms of
!  magnitude and direction.
!
!  Magnitude slpmag:
!
!                                   1
!    cos(slpmag) = -----------------------------------,   0 <= slpmag <= pi/2
!                  sqrt( 1 + (dz/dx)**2 + (dz/dy)**2 )
!
!  Direction slpdir:
!
!                  dz/dx
!    tan(slpdir) = -----,                              - pi <= slpdir <= pi
!                  dz/dy
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu and Vince Wong
!  2/16/94
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction (east/west)
!  ny       Number of grid points in the y-direction (north/south)
!
!  hterain  The heights of terrain
!
!  OUTPUT:
!
!  slpmag   Magnitude of terrain surface slope vector
!  slpdir   Direction of terrain surface slope vector from north
!           clockwise
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx,ny

  REAL :: hterain(nx,ny)    ! The heights of terrain

  REAL :: slpmag (nx,ny)    ! Terrain surface slope vector magnitude
  REAL :: slpdir (nx,ny)    ! Terrain surface slope vector direction
                            ! from north clockwise

  REAL :: tem1(nx,ny)       ! 2-D temporary array
  REAL :: tem2(nx,ny)       ! 2-D temporary array
!
!-----------------------------------------------------------------------
!
!  Local declarations
!
!-----------------------------------------------------------------------
!
  REAL :: pi
  PARAMETER ( pi = 3.141592654)

  REAL :: twdxinv, twdyinv, dzsds2

  INTEGER :: i, j

  INTEGER :: mptag1,mptag2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'bndry.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF ( ternopt == 0 ) THEN

    DO j = 1, ny-1
      DO i = 1, nx-1
        slpmag(i,j) = 0.
        slpdir(i,j) = 0.
      END DO
    END DO

  ELSE
!
!-----------------------------------------------------------------------
!
!                                           dz       dz
!  Use the central difference to calculate ---- and ----.
!                                           dx       dy
!
!-----------------------------------------------------------------------
!
    twdxinv = 1.0 / (2.0*dx)

    DO j = 1, ny-1
      DO i = 2, nx-2

        tem1(i,j) = ( hterain(i+1,j) - hterain(i-1,j) ) * twdxinv

      END DO
    END DO

    twdyinv = 1.0 / (2.0*dy)

    DO j = 2, ny-2
      DO i = 1, nx-1

        tem2(i,j) = ( hterain(i,j+1) - hterain(i,j-1) ) * twdyinv

      END DO
    END DO

!call test_dump (tem1,"XXXradfrc_tem1",nx,ny,1,0,0)
!call test_dump (tem2,"XXXradfrc_tem2",nx,ny,1,0,0)
    IF (mp_opt > 0) THEN
      CALL acct_interrupt(mp_acct)
      CALL mpsend1dew(tem1,nx,ny,ebc,wbc,0,mptag1,slpmag)  ! use slpmag as temp
      CALL mpsend1dns(tem2,nx,ny,nbc,sbc,0,mptag2,slpmag)  ! use slpmag as temp
      CALL mprecv1dew(tem1,nx,ny,ebc,wbc,0,mptag1,slpmag)  ! use slpmag as temp
      CALL mprecv1dns(tem2,nx,ny,nbc,sbc,0,mptag2,slpmag)  ! use slpmag as temp
    END IF

    CALL acct_interrupt(bc_acct)

    IF ( wbc == 1 ) THEN

      DO j = 1,ny-1
        tem1(1,j) = - tem1(2,j)
      END DO

    ELSE IF ( wbc == 2 ) THEN

      IF (mp_opt == 0) THEN
        DO j = 1,ny-1
          tem1(1,j) = tem1(nx-2,j)
        END DO
      END IF

    ELSE IF ( wbc /= 0 ) THEN

      DO j = 1, ny-1
        tem1(1,j) = ( hterain(2,j) - hterain(1,j) ) / dx
      END DO

    END IF

    IF ( ebc == 1 ) THEN

      DO j = 1, ny-1
        tem1(nx-1,j) = - tem1(nx-2,j)
      END DO

    ELSE IF ( ebc == 2 ) THEN

      IF (mp_opt == 0) THEN
        DO j = 1,ny-1
          tem1(nx-1,j) = tem1(2,j)
        END DO
      END IF

    ELSE IF ( ebc /= 0 ) THEN

      DO j = 1, ny-1

        tem1(nx-1,j) = ( hterain(nx-1,j) - hterain(nx-2,j) ) / dx

      END DO

    END IF

    IF ( sbc == 1 ) THEN

      DO i = 1, nx-1
        tem2(i,1) = - tem2(i,2)
      END DO

    ELSE IF ( sbc == 2 ) THEN

      IF (mp_opt == 0) THEN
        DO i = 1, nx-1
          tem2(i,1) = tem2(i,ny-2)
        END DO
      END IF

    ELSE IF ( sbc /= 0 ) THEN

      DO i = 1, nx-1
        tem2(i,1) = ( hterain(i,2) - hterain(i,1) ) / dy
      END DO

    END IF

    IF ( nbc == 1 ) THEN

      DO i = 1, nx-1
        tem2(i,ny-1) = - tem2(i,ny-2)
      END DO

    ELSE IF ( nbc == 2 ) THEN

      IF (mp_opt == 0) THEN
        DO i = 1, nx-1
          tem2(i,ny-1) = tem2(i,2)
        END DO
      END IF

    ELSE IF ( nbc /= 0 ) THEN

      DO i = 1, nx-1
        tem2(i,ny-1) = ( hterain(i,ny-1) - hterain(i,ny-2) ) / dy
      END DO

    END IF
!call test_dump (tem1,"XXXAradfrc_tem1",nx,ny,1,0,1)
!call test_dump (tem2,"XXXAradfrc_tem2",nx,ny,1,0,1)

    CALL acct_stop_inter

    DO j = 1, ny-1
      DO i = 1, nx-1

        dzsds2 = tem1(i,j)**2 + tem2(i,j)**2

        slpmag(i,j) = ACOS( 1. / SQRT( 1. + dzsds2 ) )

        IF ( dzsds2 == 0. ) THEN
          slpdir(i,j) = 0.
        ELSE
          slpdir(i,j) = ATAN2( tem1(i,j), tem2(i,j) )
        END IF

      END DO
    END DO

  END IF

  RETURN
END SUBROUTINE sfcslp