!
!##################################################################
!##################################################################
!###### ######
!###### 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