!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITIAL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!EMK BMJ
SUBROUTINE initial(mptr,nx,ny,nz,nxndg,nyndg,nzndg,nstyps,exbcbufsz, & 3,9
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
sdteb,sdtwb,sdtnb,sdtsb, &
ubar,vbar,ptbar,pbar,ptbari,pbari, &
rhostr,rhostri,qvbar,ppi,csndsq, &
x,y,z,zp,hterain, mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
sinlat, kmh,kmv,rprntl, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,ptsfc,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
raing,rainc,prcrate,exbcbuf,bcrlx,radfrc,radsw, &
rnflx,usflx,vsflx,ptsflx,qvsflx, &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr, &
temxy1,tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11,tem12,tem13, &
tem14,tem15,tem16,tem17,tem18,tem19, &
tem20,tem21,tem22,tem23,tem24,tem25,tem26, &
tem1_0,tem2_0,tem3_0)
!EMK BMJ
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Initialize the model parameters and variables, including base state
! variables, dependent variables and grid structure.
!
! This is the main driver for model initializations.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 11/5/92.
!
! MODIFICATION HISTORY:
!
! 5/02/92 (M. Xue)
! Added full documentation.
!
! 5/03/92 (M. Xue)
! Further documentation.
!
! 9/14/1992 (M. Xue)
! Different surface drag coefficients defined for momentum, heat and
! moisture.
! Three options included for the Coriolis force calculations.
!
! 2/12/94 (Yuhe Liu)
! Surface data and variables added for surface energy budget model.
!
! 6/10/94 (M. Xue &AS)
! Added call to initpwr.
!
! 02/07/1995 (Yuhe Liu)
! Added a new 2-D permanent array, veg(nx,ny), to the argument list
!
! 08/30/1995 (Yuhe Liu)
! Moved the initialization of external boundary arrays from the
! main program to this subroutine.
!
! 10/31/95 (D. Weber)
! Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
! radiation condition.
!
! 1/22/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 07/22/97 (D. Weber)
! Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version
! of the upper w-p radiation condition (fftopt=2).
!
! 08/01/97 (Zonghui Huo)
! Added Kain-fritsch cumulus parameterization scheme.
!
! 11/05/97 (D. Weber)
! Added temxy5 array for use in the bottom boundary condition for
! pertrubation pressure (hydrostatic).
!
! 11/06/97 (D. Weber)
! Added three additional levels to the mapfct array. The three
! levels (4,5,6) represent the inverse of the first three in order.
! The inverse map factors are computed to improve efficiency.
!
! 4/15/1998 (Donghai Wang)
! Added the source terms to the right hand terms of the qc,qr,qi,qs
! equations due to the K-F cumulus parameterization.
!
! 4/15/1998 (Donghai Wang)
! Added the running average vertical velocity (array w0avg)
! for the K-F cumulus parameterization scheme.
!
! 9/15/1998 (D. Weber)
! Added ptbari, pbari (inverse) for use in optimizing the code.
!
! 8/31/1998 (K. Brewster)
! Added call to ININUDGE to initialize nudging terms.
!
! 11/18/1998 (K. Brewster)
! Changed pibar to ppi (full pi) and moved initialization.
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
! 07/10/2001 (K. Brewster)
! Added increment arrays to argument list and to call to ininudge.
!
! 13 March 2002 (Eric Kemp)
! Added arrays for WRF BMJ cumulus scheme.
!
! April 2002 (Fanyou Kong)
! Added WRF new Kain-Fritsch (April 2002 version: KF_ETA) scheme
! initialization (lookup table)
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! mptr Grid identifier.
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nxndg Number of x grid points for nudging (1 or nx)
! nyndg Number of y grid points for nudging (1 or ny)
! nzndg Number of z grid points for nudging (1 or nz)
!
! OUTPUT:
!
! u x-component of velocity at all time levels (m/s).
! v y-component of velocity at all time levels (m/s).
! w z-component of velocity at all time levels (m/s).
! wcont Contravariant vertical velocity (m/s)
! ptprt Perturbation potential temperature at all time levels (K)
! pprt Perturbation pressure at all time levels (Pascal)
! qv Water vapor specific humidity at all time levels (kg/kg)
! qc Cloud water mixing ratio at all time levels (kg/kg)
! qr Rainwater mixing ratio at all time levels (kg/kg)
! qi Cloud ice mixing ratio at all time levels (kg/kg)
! qs Snow mixing ratio at all time levels (kg/kg)
! qh Hail mixing ratio at all time levels (kg/kg)
! tke Turbulent Kinetic Energy ((m/s)**2)
!
! udteb Time tendency of u field at east boundary (m/s**2)
! udtwb Time tendency of u field at west boundary (m/s**2)
! udtnb Time tendency of u field at north boundary (m/s**2)
! udtsb Time tendency of u field at south boundary (m/s**2)
!
! vdteb Time tendency of v field at east boundary (m/s**2)
! vdtwb Time tendency of v field at west boundary (m/s**2)
! vdtnb Time tendency of v field at north boundary (m/s**2)
! vdtsb Time tendency of v field at south boundary (m/s**2)
!
! wdteb Time tendency of w field at east boundary (m/s**2)
! wdtwb Time tendency of w field at west boundary (m/s**2)
! wdtnb Time tendency of w field at north boundary (m/s**2)
! wdtsb Time tendency of w field at south boundary (m/s**2)
!
! pdteb Time tendency of pprt field at east boundary (PASCAL/s)
! pdtwb Time tendency of pprt field at west boundary (PASCAL/s)
! pdtnb Time tendency of pprt field at north boundary (PASCAL/s)
! pdtsb Time tendency of pprt field at south boundary (PASCAL/s)
!
! sdteb Time tendency of a scalar field at east boundary (m/s**2)
! sdtwb Time tendency of a scalar field at west boundary (m/s**2)
! sdtnb Time tendency of a scalar field at north boundary (m/s**2)
! sdtsb Time tendency of a scalar field at south boundary (m/s**2)
!
! ubar Base state x-velocity component (m/s)
! vbar Base state y-velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! ptbari Inverse Base state potential temperature (K)
! pbari Inverse Base state pressure (Pascal)
! rhostr Base state density rhobar times j3 (kg/m**3)
! rhostri Inverse base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg).
! ppi Exner function
! csndsq Sound wave speed squared.
!
! x x-coordinate of grid points in computational space (m)
! y y-coordinate of grid points in computational space (m)
! z z-coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space (m)
! hterain Terrain height (m)
! mapfct Map factors at scalar, u and v points
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
! aj3x Avgx of the coordinate transformation Jacobian d(zp)/dz
! aj3y Avgy of the coordinate transformation Jacobian d(zp)/dz
! aj3z Avgz of the coordinate transformation Jacobian d(zp)/dz
! j3inv Inverse of the coordinate transformation j3
!
! trigs1 Array containing pre-computed trig function for fftopt=1.
! trigs2 Array containing pre-computed trig function for fftopt=1.
! ifax1 Array containing the factors of nx for fftopt=1.
! ifax2 Array containing the factors of ny for fftopt=1.
!
! vwork1 2-D work array for fftopt=2 option.
! vwork2 2-D work array for fftopt=2 option.
! wsave1 Work array for fftopt=2 option.
! wsave2 Work array for fftopt=2 option.
!
! sinlat Sin of latitude at each grid point
!
! kmh Horizontal turb. mixing coef. for momentum ( m**2/s )
! kmv Vertical turb. mixing coef. for momentum ( m**2/s )
! rprntl Reciprocal of Prandtl number
!
! soiltyp Soil type
! stypfrct Soil type fraction
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsfc Skin temperature at the ground or ocean surface (K)
! tsoil Deep soil temperature (K) (in deep 1 m layer)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
! ptsfc Ground surface potential temperature (K)
! qvsfc Effective S.H. at sfc.
!
! ptcumsrc Source term in pt-equation due to cumulus parameterization
! qcumsrc Source term in water equations due to cumulus parameterization
!
! nca Counter for CAPE release in the Kain-Fritsch scheme
! kfraincv K-F convective precipitation rate
! cldefi BMJ cloud efficiency
! xland BMJ land/sea mask
! bmjraincv BMJ convective precipitation amount
!
! radfrc Radiation forcing (K/s)
! radsw Solar radiation reaching the surface
! rnflx Net absorbed radiation by the surface
!
! raing Grid supersaturation rain
! rainc Cumulus convective rain
!
! usflx Surface flux of u-momentum (kg/(m*s**2))
! vsflx Surface flux of v-momentum (kg/(m*s**2))
! ptsflx Surface heat flux (K*kg/(m**2 * s ))
! qvsflx Surface moisture flux of (kg/(m**2 * s))
!
! temxy1 Temporary work array
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
! tem4 Temporary work array.
! tem5 Temporary work array.
! tem6 Temporary work array.
! tem7 Temporary work array.
! tem8 Temporary work array.
! tem9 Temporary work array.
! tem10 Temporary work array.
! tem11 Temporary work array.
! tem12 Temporary work array.
! tem13 Temporary work array.
! tem14 Temporary work array.
! tem15 Temporary work array.
! tem16 Temporary work array.
! tem17 Temporary work array.
! tem18 Temporary work array.
! tem19 Temporary work array.
! tem20 Temporary work array.
! tem21 Temporary work array.
! tem22 Temporary work array.
! tem23 Temporary work array.
! tem24 Temporary work array.
! tem25 Temporary work array.
! tem26 Temporary work array.
!
! tem1_0 Temporary work array.
! tem2_0 Temporary work array.
! tem3_0 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: mptr ! Grid identifier
INCLUDE 'timelvls.inc'
INTEGER :: nx,ny,nz ! The number of grid points in 3 directions
INTEGER :: nxndg,nyndg,nzndg ! The number of grid points in 3 directions
REAL :: u (nx,ny,nz,nt) ! Total u-velocity (m/s).
REAL :: v (nx,ny,nz,nt) ! Total v-velocity (m/s).
REAL :: w (nx,ny,nz,nt) ! Total w-velocity (m/s).
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature
! from that of base state atmosphere (Kelvin).
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure from that
! of base state atmosphere (Pascal).
REAL :: qv (nx,ny,nz,nt) ! Water vapor specific humidity (kg/kg).
REAL :: qc (nx,ny,nz,nt) ! Cloud water mixing ratio (kg/kg).
REAL :: qr (nx,ny,nz,nt) ! Rain water mixing ratio (kg/kg).
REAL :: qi (nx,ny,nz,nt) ! Cloud ice mixing ratio (kg/kg).
REAL :: qs (nx,ny,nz,nt) ! Snow mixing ratio (kg/kg).
REAL :: qh (nx,ny,nz,nt) ! Hail mixing ratio (kg/kg).
REAL :: tke (nx,ny,nz,nt) ! Turbulent Kinetic Energy ((m/s)**2)
REAL :: udteb (ny,nz) ! T-tendency of u at e-boundary (m/s**2)
REAL :: udtwb (ny,nz) ! T-tendency of u at w-boundary (m/s**2)
REAL :: udtnb (nx,nz) ! T-tendency of u at n-boundary (m/s**2)
REAL :: udtsb (nx,nz) ! T-tendency of u at s-boundary (m/s**2)
REAL :: vdteb (ny,nz) ! T-tendency of v at e-boundary (m/s**2)
REAL :: vdtwb (ny,nz) ! T-tendency of v at w-boundary (m/s**2)
REAL :: vdtnb (nx,nz) ! T-tendency of v at n-boundary (m/s**2)
REAL :: vdtsb (nx,nz) ! T-tendency of v at s-boundary (m/s**2)
REAL :: wdteb (ny,nz) ! T-tendency of w at e-boundary (m/s**2)
REAL :: wdtwb (ny,nz) ! T-tendency of w at w-boundary (m/s**2)
REAL :: wdtnb (nx,nz) ! T-tendency of w at n-boundary (m/s**2)
REAL :: wdtsb (nx,nz) ! T-tendency of w at s-boundary (m/s**2)
REAL :: pdteb (ny,nz) ! T-tendency of pprt at e-boundary (PASCAL/s)
REAL :: pdtwb (ny,nz) ! T-tendency of pprt at w-boundary (PASCAL/s)
REAL :: pdtnb (nx,nz) ! T-tendency of pprt at n-boundary (PASCAL/s)
REAL :: pdtsb (nx,nz) ! T-tendency of pprt at s-boundary (PASCAL/s)
REAL :: sdteb (ny,nz) ! T-tendency of w at e-boundary (m/s**2)
REAL :: sdtwb (ny,nz) ! T-tendency of w at w-boundary (m/s**2)
REAL :: sdtnb (nx,nz) ! T-tendency of w at n-boundary (m/s**2)
REAL :: sdtsb (nx,nz) ! T-tendency of w at s-boundary (m/s**2)
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s).
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s).
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: ptbari(nx,ny,nz) ! Inverse Base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse Base state pressure (Pascal).
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: rhostri(nx,ny,nz) ! Inverse base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity
! (kg/kg).
REAL :: ppi (nx,ny,nz) ! Exner function.
REAL :: csndsq(nx,ny,nz) ! Sound wave speed squared.
REAL :: x (nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y (ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z (nz) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point on the staggered grid.
REAL :: hterain(nx,ny) ! Terrain height (m).
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian -d(zp)/dx.
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian -d(zp)/dy.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian d(zp)/dz.
REAL :: aj3x (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL :: aj3y (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL :: aj3z (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL :: j3inv (nx,ny,nz) ! Inverse of J3
REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
INTEGER :: ifax1(13) ! Array containing the factors of nx for
! fftopt=1.
INTEGER :: ifax2(13) ! Array containing the factors of ny for
! fftopt=1.
REAL :: vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2.
REAL :: vwork2 (ny,nx+1) ! 2-D work array for fftopt=2.
REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2.
REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2.
REAL :: sinlat(nx,ny) ! Sin of latitude at each grid point
REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: rprntl(nx,ny,nz) ! Reciprocal of Prandtl number
INTEGER :: nstyps ! Number of soil type
INTEGER :: soiltyp(nx,ny,nstyps) ! Soil types at grids
REAL :: stypfrct(nx,ny,nstyps) ! Fraction of soil types
INTEGER :: vegtyp (nx,ny) ! Vegetation type
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: roufns (nx,ny) ! Surface roughness
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: tsfc (nx,ny,0:nstyps) ! Temperature at ground (K)
! (in top 1 cm layer)
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
! (in deep 1 m layer)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture
! in the top 1 cm layer
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture
! in the deep 1 m layer
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: qvsfc (nx,ny,0:nstyps) ! Effective qv at sfc.
REAL :: ptsfc (nx,ny) ! Ground surface potential
! temperature (K)
REAL :: ptcumsrc(nx,ny,nz) ! Source term in pt-equation due
! to cumulus parameterization
REAL :: qcumsrc(nx,ny,nz,5) ! Source term in water equations due
! to cumulus parameterization:
! qcumsrc(1,1,1,1) for qv equation
! qcumsrc(1,1,1,2) for qc equation
! qcumsrc(1,1,1,3) for qr equation
! qcumsrc(1,1,1,4) for qi equation
! qcumsrc(1,1,1,5) for qs equation
REAL :: w0avg(nx,ny,nz) ! a closing running average vertical
! velocity in 10min for K-F scheme
REAL :: kfraincv(nx,ny) ! K-F convective rainfall (cm)
INTEGER :: nca(nx,ny) ! K-F counter for CAPE release
!EMK BMJ
REAL,INTENT(INOUT) :: cldefi(nx,ny) ! BMJ cloud efficiency
REAL,INTENT(INOUT) :: xland(nx,ny) ! BMJ land mask
! (1.0 = land, 2.0 = sea)
REAL,INTENT(INOUT) :: bmjraincv(nx,ny) ! BMJ convective rainfall (cm)
!EMK END
REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s)
REAL :: radsw(nx,ny) ! Solar radiation reacing the surface
REAL :: rnflx(nx,ny) ! Net absorbed radiation by the surface
REAL :: raing(nx,ny) ! Grid supersaturation rain
REAL :: rainc(nx,ny) ! Cumulus convective rain
REAL :: prcrate(nx,ny,4) ! precipitation rates (kg/(m**2*s))
! prcrate(1,1,1) = total precipitation rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2))
REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2))
REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m*s**2))
REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s))
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx(nx,ny) ! EXBC relaxation coefficients
REAL :: uincr(nxndg,nyndg,nzndg) ! Analysis increment for u
REAL :: vincr(nxndg,nyndg,nzndg) ! Analysis increment for v
REAL :: wincr(nxndg,nyndg,nzndg) ! Analysis increment for w
REAL :: pincr(nxndg,nyndg,nzndg) ! Analysis increment for p
REAL :: ptincr(nxndg,nyndg,nzndg) ! Analysis increment for pt
REAL :: qvincr(nxndg,nyndg,nzndg) ! Analysis increment for qv
REAL :: qcincr(nxndg,nyndg,nzndg) ! Analysis increment for qc
REAL :: qrincr(nxndg,nyndg,nzndg) ! Analysis increment for qr
REAL :: qiincr(nxndg,nyndg,nzndg) ! Analysis increment for qi
REAL :: qsincr(nxndg,nyndg,nzndg) ! Analysis increment for qs
REAL :: qhincr(nxndg,nyndg,nzndg) ! Analysis increment for qh
REAL :: temxy1(nx,ny) ! Temporary work array
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
REAL :: tem2 (nx,ny,nz) ! Temporary work array.
REAL :: tem3 (nx,ny,nz) ! Temporary work array.
REAL :: tem4 (nx,ny,nz) ! Temporary work array.
REAL :: tem5 (nx,ny,nz) ! Temporary work array.
REAL :: tem6 (nx,ny,nz) ! Temporary work array.
REAL :: tem7 (nx,ny,nz) ! Temporary work array.
REAL :: tem8 (nx,ny,nz) ! Temporary work array.
REAL :: tem9 (nx,ny,nz) ! Temporary work array.
REAL :: tem10 (nx,ny,nz) ! Temporary work array.
REAL :: tem11 (nx,ny,nz) ! Temporary work array.
REAL :: tem12 (nx,ny,nz) ! Temporary work array.
REAL :: tem13 (nx,ny,nz) ! Temporary work array.
REAL :: tem14 (nx,ny,nz) ! Temporary work array.
REAL :: tem15 (nx,ny,nz) ! Temporary work array.
REAL :: tem16 (nx,ny,nz) ! Temporary work array.
REAL :: tem17 (nx,ny,nz) ! Temporary work array.
REAL :: tem18 (nx,ny,nz) ! Temporary work array.
REAL :: tem19 (nx,ny,nz) ! Temporary work array.
REAL :: tem20 (nx,ny,nz) ! Temporary work array.
REAL :: tem21 (nx,ny,nz) ! Temporary work array.
REAL :: tem22 (nx,ny,nz) ! Temporary work array.
REAL :: tem23 (nx,ny,nz) ! Temporary work array.
REAL :: tem24 (nx,ny,nz) ! Temporary work array.
REAL :: tem25 (nx,ny,nz) ! Temporary work array.
REAL :: tem26 (nx,ny,nz) ! Temporary work array.
REAL :: tem1_0(0:nx,0:ny,0:nz) ! Temporary work array.
REAL :: tem2_0(0:nx,0:ny,0:nz) ! Temporary work array.
REAL :: tem3_0(0:nx,0:ny,0:nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'indtflg.inc'
INCLUDE 'phycst.inc'
INCLUDE 'exbc.inc'
INCLUDE 'nudging.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j
INTEGER :: ireturn
REAL :: tem
!EMK BMJ
LOGICAL :: restart
INTEGER,PARAMETER :: vegwaterflag = 14
INTEGER,PARAMETER :: xland_waterflag = 2
INTEGER,PARAMETER :: xland_landflag = 1
!EMK END
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
mgrid = mptr
grdin = 0
basin = 0
varin = 0
mstin = 0
rainin= 0
prcin = 0
icein = 0
trbin = 0
sfcin = 0
landin= 0
radin = 0
flxin = 0
!wdt update - init0 no longer necessary (arrays set to 0 after allocation)
!
!-----------------------------------------------------------------------
!
! INITialize model array VARiables.
!
!-----------------------------------------------------------------------
!
!EMK BMJ
CALL initgrdvar
(nx,ny,nz,nt,nstyps,exbcbufsz, &
x,y,z,zp,hterain,mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
udteb, udtwb, vdtnb, vdtsb, &
pdteb,pdtwb ,pdtnb ,pdtsb, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
ubar,vbar,ptbar,pbar,ptbari,pbari, &
rhostr,rhostri,qvbar,ppi,csndsq, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
raing,rainc,prcrate,exbcbuf, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3,tem4,tem5, &
tem6,tem7,tem8,tem9)
DO j=1,ny-1
DO i=1,nx-1
tem = 0.5 * ( pprt(i,j,1,2)+pbar(i,j,1) &
+ pprt(i,j,2,2)+pbar(i,j,2) )
ptsfc(i,j)=tsfc(i,j,0)*(p0/tem)**rddcp
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate the sin of the lattitude of each grid point, to be used
! in the calculation of latitude-dependent Coriolis parameters.
!
!-----------------------------------------------------------------------
!
CALL gtsinlat(nx,ny,x,y, sinlat, tem1,tem2, tem3)
!
!-----------------------------------------------------------------------
!
! Initialize arrays that store the lookup table data.
!
!-----------------------------------------------------------------------
!
CALL initlktb
!
!-----------------------------------------------------------------------
!
! Initialize the external boundary data array.
!
!-----------------------------------------------------------------------
!
IF( lbcopt == 2 .AND. mptr == 1 ) THEN
ireturn = 0
CALL extbdtini(nx,ny,nz, &
u,v,w,ptprt,pprt, &
qv,qc,qr,qi,qs,qh,ptbar,pbar, &
exbcbuf(nu0exb), exbcbuf(nv0exb), &
exbcbuf(nw0exb), exbcbuf(npt0exb), &
exbcbuf(npr0exb), exbcbuf(nqv0exb), &
exbcbuf(nqc0exb), exbcbuf(nqr0exb), &
exbcbuf(nqi0exb), exbcbuf(nqs0exb), &
exbcbuf(nqh0exb), exbcbuf(nudtexb), &
exbcbuf(nvdtexb), exbcbuf(nwdtexb), &
exbcbuf(nptdtexb),exbcbuf(nprdtexb), &
exbcbuf(nqvdtexb),exbcbuf(nqcdtexb), &
exbcbuf(nqrdtexb),exbcbuf(nqidtexb), &
exbcbuf(nqsdtexb),exbcbuf(nqhdtexb), &
bcrlx, &
tem1,tem2,tem3,tem4,tem5,tem6, &
tem7,tem8,tem9,tem10,tem11,ireturn)
IF( ireturn == 1 ) THEN
WRITE (6,'(a/a)') &
'Can not find the external boundary data. Dump the', &
'history file and restart file and then STOP the model.'
CALL arpsstop('arpsstop called from initial with ext boundary file',1)
ELSE IF( ireturn == 2 ) THEN
WRITE (6,'(a/a)') &
'Can not open the external boundary data. Dump the history', &
'file and restart file and then STOP the model.'
CALL arpsstop('arpsstop called from initial with opening ext &
& boundary file ',1)
ELSE IF( ireturn == 3 ) THEN
WRITE (6,'(a/a)') &
'Read errors in the external boundary data file. Dump the', &
'history file and restart file and then STOP the model.'
CALL arpsstop('arpsstop called from initial with reading ext &
& boundary file ',1)
ELSE IF( ireturn /= 0 ) THEN
WRITE (6,'(a/a)') &
'Other errors in getting the external boundary data. Dump the', &
'history file and restart file and then STOP the model.'
CALL arpsstop('arpsstop called from initial with the ext &
& boundary file ',1)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Initialize nudging assimilation variables
!
!-----------------------------------------------------------------------
!
IF( nudgopt > 0 ) &
CALL ininudge(nxndg,nyndg,nzndg, &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr,ireturn)
!-----------------------------------------------------------------------
!
! Initialize KF_ETA arrays and look-up tables.
!
!-----------------------------------------------------------------------
!KFY KF_ETA
IF (cnvctopt == 5) THEN
!...Now initialize KF_ETA look-up tables
IF (initopt == 2) THEN
restart = .TRUE.
ELSE
restart = .FALSE.
END IF
CALL interface_wrf_kfinit(nx,ny,nz,nca,restart)
END IF ! IF (cnvctopt == 5) THEN
!KFY KF_ETA END
!-----------------------------------------------------------------------
!
! Initialize BMJ arrays and look-up tables.
!
!-----------------------------------------------------------------------
!EMK BMJ
IF (cnvctopt == 4) THEN
DO j = 1,ny-1
DO i = 1,nx-1
IF (vegtyp(i,j) == vegwaterflag) THEN
xland(i,j) = xland_waterflag
ELSE
xland(i,j) = xland_landflag
END IF
END DO ! DO i = 1,nx
END DO ! DO j = 1,ny
!...Now initialize BMJ look-up tables
IF (initopt == 2) THEN
restart = .TRUE.
ELSE
restart = .FALSE.
END IF
CALL interface_wrf_bmjinit(nx,ny,nz,cldefi,restart)
END IF ! IF (cnvctopt == 4) THEN
!EMK END
!
!-----------------------------------------------------------------------
!
! End of model initialization.
!
!-----------------------------------------------------------------------
!
RETURN
END SUBROUTINE initial
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INIGRD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE inigrd(nx,ny,nz, & 3,18
x,y,z,zp,hterain,mapfct,j1,j2,j3,zp1d,dzp1d,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Initialize the model grid variables.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/17/1991.
!
! MODIFICATION HISTORY:
!
! 5/02/92 (M. Xue)
! Added full documentation.
!
! 5/03/92 (M. Xue)
! Further documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
! 12/22/94 (Yuhe Liu)
! Changed variable tuning, which was hard wired inside this
! subroutine, to strhtune, which is input from namelist input file.
!
! 1/22/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 11/06/97 (D. Weber)
! Added three additional levels to the mapfct array. The three
! levels (4,5,6) represent the inverse of the first three in order.
! The inverse map factors are computed to improve efficiency.
!
! 2000/04/11 (Gene Bassett)
! Only update the terrain data with a call to BCS2D for ternopt=1.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! OUTPUT:
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space
! (m)
! hterain Terrain height (m)
! mapfct Map factors at scalar, u and v points
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! WORK ARRAY:
!
! zp1d Temporary work array.
! dzp1d Temporary work array.
! tem1 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number of grid points in 3 directions
REAL :: x(nx) ! x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y(ny) ! y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z(nz) ! z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL :: zp(nx,ny,nz) ! Physical height coordinate defined at
! w-point of the staggered grid.
REAL :: j1(nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dx.
REAL :: j2(nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dy.
REAL :: j3(nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/dz.
REAL :: hterain(nx,ny) ! Terrain height.
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: zp1d (nz) ! Temporary array
REAL :: dzp1d(nz) ! Temporary array
REAL :: tem1(nx,ny,nz) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
REAL :: xs, ys, zs, radnd, pi2
REAL :: zflat1,z1,z2
REAL :: zpmax
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'phycst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Define a uniform model grid.
!
!-----------------------------------------------------------------------
!
CALL setgrd ( nx,ny, x,y )
!
!-----------------------------------------------------------------------
!
! Set map factor at scalar, u, and v points
!
!-----------------------------------------------------------------------
!
IF ( mpfctopt /= 0 ) THEN
DO j=1,ny-1
DO i=1,nx-1
xs = 0.5*(x(i)+x(i+1))
ys = 0.5*(y(j)+y(j+1))
CALL xytomf( 1,1,xs,ys,mapfct(i,j,1) )
IF(maptest == 1)THEN ! set up a symmetric field...
mapfct(i,j,1) = 1.0 + ABS((xs-0.5*(x(1)+x(nx))) &
/(x(nx)-x(1)) &
+(ys-0.5*(y(1)+y(ny))) &
/(y(ny)-y(1)))
END IF
mapfct(i,j,4) = 1.0/mapfct(i,j,1)
mapfct(i,j,7) = mapfct(i,j,1)*mapfct(i,j,1)
mapfct(i,j,8) = 0.25*mapfct(i,j,1) ! for use in sovlwpim
! and wcontra...
END DO
END DO
DO j=1,ny-1
DO i=1,nx
ys = 0.5*(y(j)+y(j+1))
CALL xytomf( 1,1,x(i),ys,mapfct(i,j,2) )
IF(maptest == 1)THEN ! set up a symmetric field...
mapfct(i,j,2) = 1.0 + ABS((x(i)-0.5*(x(1)+x(nx))) &
/(x(nx)-x(1)) &
+(ys-0.5*(y(1)+y(ny))) &
/(y(ny)-y(1)))
END IF
mapfct(i,j,5) = 1.0/mapfct(i,j,2)
END DO
END DO
DO j=1,ny
DO i=1,nx-1
xs = 0.5*(x(i)+x(i+1))
CALL xytomf( 1,1,xs,y(j),mapfct(i,j,3) )
IF(maptest == 1)THEN ! set up a symmetric field...
mapfct(i,j,3) = 1.0 + ABS((xs-0.5*(x(1)+x(nx))) &
/(x(nx)-x(1)) &
+(y(j)-0.5*(y(1)+y(ny))) &
/(y(ny)-y(1)))
END IF
mapfct(i,j,6) = 1.0/mapfct(i,j,3)
END DO
END DO
ELSE
DO k=1,7
DO j=1,ny
DO i=1,nx
mapfct(i,j,k) = 1.0
mapfct(i,j,8) = 0.25
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Print map factor at scalar, u, and v points for the purpose of
! debug
!
!-----------------------------------------------------------------------
!
! CALL getunit( mpunit )
! CALL asnctl ('NEWLOCAL', 1, ierr)
! CALL asnfile('mapfactor.data', '-F f77 -N ieee', ierr)
!
! OPEN (mpunit,file='mapfactor.data',form='unformatted')
!
! CALL edgfill(mapfct(1,1,1),1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
! write(mpunit) ((mapfct(i,j,1),i=1,nx),j=1,ny)
!
! CALL edgfill(mapfct(1,1,2),1,nx,1,nx, 1,ny,1,ny-1, 1,1,1,1)
! write(mpunit) ((mapfct(i,j,2),i=1,nx),j=1,ny)
!
! CALL edgfill(mapfct(1,1,3),1,nx,1,nx-1, 1,ny,1,ny, 1,1,1,1)
! write(mpunit) ((mapfct(i,j,3),i=1,nx),j=1,ny)
!
! CLOSE( mpunit )
! CALL retunit( mpunit )
!
!-----------------------------------------------------------------------
!
! End of debug print.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
z(k) = zrefsfc + (k-2) * dz
END DO
!
!-----------------------------------------------------------------------
!
! Specify the terrain
!
!-----------------------------------------------------------------------
!
IF( ternopt == 0 ) THEN ! No terrain, the ground is flat
DO j=1,ny-1
DO i=1,nx-1
hterain(i,j) = zrefsfc
END DO
END DO
ELSE IF( ternopt == 1 ) THEN ! Bell-shaped mountain
!
!-----------------------------------------------------------------------
!
! Define the bell-shaped mountain
!
!-----------------------------------------------------------------------
!
pi2 = 1.5707963267949
DO j=1,ny-1
DO i=1,nx-1
xs = (x(i)+x(i+1))*0.5
ys = (y(j)+y(j+1))*0.5
zs = z(2)
IF( mntwidy < 0.0 .OR. runmod == 2 ) THEN ! 2-d terrain in x-z plane.
radnd = 1.0+((xs-mntctrx)/mntwidx)**2
ELSE IF( mntwidx < 0.0 .OR. runmod == 3 ) THEN ! 2-d terrain in y-z plane.
radnd = 1.0+((ys-mntctry)/mntwidy)**2
ELSE ! 3-d terrain
radnd = 1.0+((xs-mntctrx)/mntwidx)**2 &
+((ys-mntctry)/mntwidy)**2
END IF
hterain(i,j) = zrefsfc + hmount/radnd
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Make sure that the terrain satisfies the specified boundary
! conditions.
!
!-----------------------------------------------------------------------
!
!call test_dump (hterain,"XXXinit_hterain",nx,ny,1,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend1dew(hterain,nx,ny,ebc,wbc,0,mptag,tem1)
CALL mprecv1dew(hterain,nx,ny,ebc,wbc,0,mptag,tem1)
CALL mpsend1dns(hterain,nx,ny,nbc,sbc,0,mptag,tem1)
CALL mprecv1dns(hterain,nx,ny,nbc,sbc,0,mptag,tem1)
END IF
!call test_dump (hterain,"XXXBinit_hterain",nx,ny,1,0,1)
CALL acct_interrupt(bc_acct)
CALL bcs2d(nx,ny,hterain, ebc,wbc,nbc,sbc)
CALL acct_stop_inter
!call test_dump (hterain,"XXXAinit_hterain",nx,ny,1,0,1)
ELSE IF( ternopt == 2 ) THEN ! Read from terrain data base
!
!-----------------------------------------------------------------------
!
! Read in the terrain data.
!
!-----------------------------------------------------------------------
!
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL readtrn( nx,ny,dx,dy, terndta, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
hterain )
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Set up a stretched vertical grid.
!
! For strhopt=1, function y = a+b*x**3 is used to specify dz as a
! function of k.
! For strhopt=2, function y = c + a*tanh(b*x) is used to specify dz
! as a function of k.
!
!-----------------------------------------------------------------------
!
IF ( strhopt == 0 ) THEN
DO k=1,nz
zp1d(k) = z(k)
END DO
ELSE IF ( strhopt == 1 .OR.strhopt == 2 ) THEN
z1 = zrefsfc + MAX(0.0, MIN(dlayer1, z(nz-2)-zrefsfc ))
z2 = z1 + MAX(0.0, MIN(dlayer2, z(nz-1)-z1 ))
IF( dlayer1 >= (nz-3)*dzmin ) THEN
WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)') &
'Can not setup a vertical grid with uniform dz=',dzmin, &
' over the depth of ',dlayer1,' please specify a smaller ', &
'value of dlayer1. Program stopped INIGRD.'
CALL arpsstop('arpsstop called from INIGRD with ther vertical grid ' &
,1)
END IF
CALL strhgrd(nz,strhopt,zrefsfc,z1,z2,z(nz-1), &
dzmin,strhtune, zp1d,dzp1d)
ELSE
WRITE(6,'(1x,a,i3,a/)') &
'Invalid vertical grid stretching option, strhopt was ',strhopt, &
'. Program stopped in INIGRD.'
CALL arpsstop('arpsstop called from INIGRD with stretching ' ,1)
END IF
!
!-----------------------------------------------------------------------
!
! Physical height of computational grid defined as
!
! Zp=(z-zrefsfc)*(Zm-hterain)/(Zm-zrefsfc)+hterain for z=<Zm.
! ZP=z for z>Zm
!
! where Zm the height at which the grid levels becomes flat.
! Hm < Zm =< Ztop, hm is the height of mountain and Ztop the height
! of model top.
!
!-----------------------------------------------------------------------
!
DO k=nz-1,2,-1
IF(zp1d(k) <= zflat) THEN
zflat1 = zp1d(k)
EXIT
END IF
END DO
zflat1=MAX(MIN(
DO k=2,nz-1
IF(zp1d(k) > zflat1) THEN
DO j=1,ny-1
DO i=1,nx-1
zp(i,j,k)=zp1d(k)
END DO
END DO
ELSE
DO j=1,ny-1
DO i=1,nx-1
zp(i,j,k)=(zp1d(k)-zrefsfc)*(zflat1-hterain(i,j)) &
/(zflat1-zrefsfc)+hterain(i,j)
END DO
END DO
END IF
END DO
DO j=1,ny-1
DO i=1,nx-1
zp(i,j,2)=hterain(i,j)
zp(i,j,1)=2.0*zp(i,j,2)-zp(i,j,3)
zp(i,j,nz)=2.0*zp(i,j,nz-1)-zp(i,j,nz-2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate transformation Jacobians J1, J2 and J3.
!
!-----------------------------------------------------------------------
!
CALL jacob(nx,ny,nz,x,y,z,zp,j1,j2,j3)
RETURN
END SUBROUTINE inigrd
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE STRHGRD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE strhgrd(nz,strhopt,z0,z1,z2,ztop,dzmin,strhtune, z,dzk) 3,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To construct a vertically stretched grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/17/94.
!
! MODIFICATION HISTORY:
!
! 05/11/95 (Jinxing Zong and MX)
!
! A bug fix for the case of nonzero zrefsfc, the reference height of
! the surface. Results not affected for zrefsfc=0 (default value).
!
!-----------------------------------------------------------------------
!
!
! INPUT:
!
! nz The vertical dimension of ARPS grid.
!
!
! OUTPUT:
!
! z Array containing the height of veritical grid levels.
! dzk Array containing the grid spacing between vertical levels
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nz
INTEGER :: strhopt
REAL :: z0
REAL :: z1
REAL :: z2
REAL :: ztop
REAL :: dzmin
REAL :: strhtune
REAL :: z (nz)
REAL :: dzk (nz)
REAL :: rnzh,dzm
REAL :: a,b,c,hnew,zkmid,dzu
INTEGER :: nzh,nzl,k
REAL :: dz
IF( (z1-z0) == (nz-3)*dzmin.AND.(z1-z0) == (ztop-z0) ) THEN
dz = (ztop-z0)/(nz-3)
DO k=1,nz-1
dzk(k)= dz
END DO
DO k=1,nz
z(k)=z0 + (k-2) * dz
END DO
WRITE(6,'(/1x,a,f13.3,/a,f13.3)') &
'Layer 1 depth was as deep as the entire domain. i', &
'A uniform vertical grid is assumed with dz=',dz, &
' over the model depth of ',ztop-z0
RETURN
END IF
IF(z1 < z0) z1 = z0
IF(z2 > ztop) z2 = ztop
nzl = (z1-z0)/dzmin
IF( (z1-z0) >= (nz-3)*dzmin ) THEN
WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)') &
'Can not setup a vertical grid with uniform dz=',dzmin, &
' over the depth of ',z1-z0,' please specify a smaller', &
' value of dlayer1 '
CALL arpsstop('arpsstop called from STRHGRD with stretching ' ,1)
END IF
IF( z2 >= ztop ) THEN
dzm = (ztop-z0-nzl*dzmin)/(nz-3-nzl)
! print*, nzl*dzmin + (nz-3-nzl)*dzm
nzh = 0
dzu = 2*dzm - dzmin
ELSE
a = 2*(nz-3-nzl)
b = 2*z0-ztop-z2-(nz-3-3*nzl)*dzmin
c = dzmin*(z2-z0-nzl*dzmin)
dzm = (-b + SQRT(b*b-4*a*c) )/(2*a)
rnzh = (ztop-z2)/(2*dzm-dzmin)
nzh = INT(rnzh)
hnew = nzl*dzmin + nzh*(2*dzm-dzmin) + &
(nz-3-nzl-nzh)*dzm + z0
IF( nzh /= 0 ) THEN
dzu = (2*dzm-dzmin) + (ztop-hnew)/nzh
ELSE IF( nz-3-nzl-nzh /= 0 ) THEN
dzm = dzm + (ztop-hnew)/(nz-3-nzl-nzh)
dzu = (2*dzm-dzmin)
END IF
END IF
DO k=1,nzl+1
dzk(k)=dzmin
END DO
IF( strhopt == 1 ) THEN
a = (dzm-dzmin)
DO k=nzl+2,nz-2-nzh
dzk(k)= dzm+a* &
((2.0*FLOAT(k-nzl-2)/FLOAT(nz-4-nzh-nzl)-1.0) )**3
END DO
ELSE
zkmid=0.5*FLOAT( nz-nzh+nzl)
IF( nzl+2-zkmid == 0.0 ) THEN
b = 0.0
ELSE
b= strhtune* 2.0/(nzl+2-zkmid)
END IF
a=(dzmin-dzm)/TANH( strhtune* 2.0)
DO k=nzl+2,nz-2-nzh
dzk(k)=dzm + a*TANH(b*(FLOAT(k)-zkmid))
END DO
END IF
DO k=nz-2-nzh+1, nz-2
dzk(k)= dzu
END DO
dzk(nz-1) = dzk(nz-2)
dzk(nz ) = dzk(nz-1)
z(2) = z0
DO k=3,nz-1
z(k) = z(k-1)+dzk(k-1)
END DO
z(1) =z(2)-dzk(1)
z(nz-1)=ztop
z(nz)=z(nz-1)+dzk(nz-1)
RETURN
END SUBROUTINE strhgrd
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE JACOB ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE jacob(nx,ny,nz,x,y,z,zp, j1,j2,j3) 6,10
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate transformation Jacobians J1, J2 and J3.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/17/1991.
!
! MODIFICATION HISTORY:
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/2/94 (M. Xue)
! Loop 710 that resets j2 on north and south boundary to
! zero gradient deleted. It shouldn't have been there.
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! OUTPUT:
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space
! (m)
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number of grid points in 3 directions
REAL :: x(nx) ! x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y(ny) ! y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z(nz) ! z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL :: zp(nx,ny,nz) ! the physical height coordinate defined at
! w-point of the staggered grid.
REAL :: j1(nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dx.
REAL :: j2(nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dy.
REAL :: j3(nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/dz.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: mptag
REAL, allocatable :: mp_tem(:)
INTEGER :: istat
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'bndry.inc'
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (mp_opt > 0) THEN
ALLOCATE (mp_tem(MAX(nx,ny)*nz),stat=istat)
IF (istat /= 0) THEN
CALL arpsstop ("arpsstop called from JACOB: ERROR allocating mp_tem" &
,1)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate transformation Jacobian J1 defined as -del(zp)/del(x)
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny-1
DO i=2,nx-1
j1(i,j,k) = -2 * (zp(i,j,k)-zp(i-1,j,k)) / (x(i+1)-x(i-1))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! X - boundary conditions of j1
!
!-----------------------------------------------------------------------
!
!call test_dump (zp,"XXXinit_zp",nx,ny,nz,0,1)
!call test_dump (x,"XXXinit_x",nx,1,1,1,1)
!call test_dump (j1,"XXXinit_j1",nx,ny,nz,1,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(j1,nx,ny,nz,ebc,wbc,1,mptag,mp_tem)
CALL mprecv2dew(j1,nx,ny,nz,ebc,wbc,1,mptag,mp_tem)
END IF
!call test_dump (j1,"XXXBinit_j1",nx,ny,nz,1,1)
CALL acct_interrupt(bc_acct)
IF(wbc == 1) THEN ! Rigid wall boundary condition
DO k=1,nz
DO j=1,ny-1
j1( 1,j,k)=j1( 3 ,j,k)
END DO
END DO
ELSE IF(wbc == 2) THEN ! Periodic boundary condition.
DO k=1,nz
DO j=1,ny-1
j1( 1,j,k)=j1(nx-2,j,k)
END DO
END DO
ELSE IF(wbc /= 0) THEN
DO k=1,nz
DO j=1,ny-1
j1( 1,j,k)=j1( 2 ,j,k)
END DO
END DO
END IF
IF(ebc == 1) THEN ! Rigid wall boundary condition
DO k=1,nz
DO j=1,ny-1
j1(nx,j,k)=j1(nx-2,j,k)
END DO
END DO
ELSE IF(ebc == 2) THEN ! Periodic boundary condition.
DO k=1,nz
DO j=1,ny-1
j1(nx,j,k)=j1( 3 ,j,k)
END DO
END DO
ELSE IF(ebc /= 0) THEN
DO k=1,nz
DO j=1,ny-1
j1(nx,j,k)=j1(nx-1,j,k)
END DO
END DO
END IF
!call test_dump (j1,"XXXCinit_j1",nx,ny,nz,1,1)
DO k=1,nz
DO i=1,nx
j1(i,ny,k) = j1(i,ny-1,k)
END DO
END DO
!call test_dump (j1,"XXXAinit_j1",nx,ny,nz,1,1)
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! Calculate transformation Jacobian J2 defined as -del(zp)/del(y)
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO i=1,nx-1
DO j=2,ny-1
j2(i,j,k) = -2 * (zp(i,j,k)-zp(i,j-1,k)) / (y(j+1)-y(j-1))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Y - boundary conditions of j2
!
!-----------------------------------------------------------------------
!
!call test_dump (j2,"XXXinit_j2",nx,ny,nz,2,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dns(j2,nx,ny,nz,nbc,sbc,2,mptag,mp_tem)
CALL mprecv2dns(j2,nx,ny,nz,nbc,sbc,2,mptag,mp_tem)
END IF
CALL acct_interrupt(bc_acct)
IF(sbc == 1) THEN ! Rigid wall boundary condition
DO k=1,nz
DO i=1,nx-1
j2(i, 1,k)=j2(i, 3 ,k)
END DO
END DO
ELSE IF(sbc == 2) THEN ! Periodic boundary condition.
IF (mp_opt == 0) THEN
DO k=1,nz
DO i=1,nx-1
j2(i, 1,k)=j2(i,ny-2,k)
END DO
END DO
END IF
ELSE IF(sbc /= 0) THEN
DO k=1,nz
DO i=1,nx-1
j2(i, 1,k)=j2(i, 2 ,k)
END DO
END DO
END IF
IF(nbc == 1) THEN ! Rigid wall boundary condition
DO k=1,nz
DO i=1,nx-1
j2(i,ny,k)=j2(i,ny-2,k)
END DO
END DO
ELSE IF(nbc == 2) THEN ! Periodic boundary condition.
IF (mp_opt == 0) THEN
DO k=1,nz
DO i=1,nx-1
j2(i,ny,k)=j2(i, 3 ,k)
END DO
END DO
END IF
ELSE IF(nbc /= 0) THEN
DO k=1,nz
DO i=1,nx-1
j2(i,ny,k)=j2(i,ny-1,k)
END DO
END DO
END IF
DO k=1,nz
DO j=1,ny
j2(nx,j,k) = j2(nx-1,j,k)
END DO
END DO
!call test_dump (j2,"XXXAinit_j2",nx,ny,nz,2,1)
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! Calculate transformation Jacobian J3 defined as del(zp)/del(z)
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
j3(i,j,k) = (zp(i,j,k+1)-zp(i,j,k))/(z(k+1)-z(k))
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
j3(nx,j,k) = j3(nx-1,j,k)
END DO
END DO
DO k=1,nz-1
DO i=1,nx
j3(i,ny,k) = j3(i,ny-1,k)
END DO
END DO
DO j=1,ny
DO i=1,nx
j3(i,j,nz) = j3(i,j,nz-1)
END DO
END DO
IF (mp_opt > 0) THEN
DEALLOCATE (mp_tem,stat=istat)
END IF
RETURN
END SUBROUTINE jacob
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITGRDVAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!EMK BMJ
SUBROUTINE initgrdvar(nx,ny,nz,nts,nstyps,exbcbufsz, & 2,46
x,y,z,zp,hterain,mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
udteb, udtwb, vdtnb, vdtsb, &
pdteb,pdtwb ,pdtnb ,pdtsb, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
ubar,vbar,ptbar,pbar,ptbari,pbari, &
rhostr,rhostri,qvbar,ppi,csndsq, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
raing,rainc,prcrate,exbcbuf, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3,tem4,tem5, &
tem6,tem7,tem8,tem9)
!EMK END
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Initialize the model array variables.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/17/1992.
!
! MODIFICATION HISTORY:
!
! 5/02/92 (M. Xue)
! Added full documentation.
!
! 5/03/92 (M. Xue)
! Further documentation.
!
! 10/7/1992 (M. Xue)
! Added call to subroutine extinit, the option three of
! initialization.
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
! 02/07/1995 (Yuhe Liu)
! Added a new 2-D permanent array, veg(nx,ny), to the argument list
!
! 10/31/95 (D. Weber)
! Added trigs1,trigs2,ifax1,ifax2 for use in the fft code for the
! upper radiation condition.
!
! 1/22/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 07/22/97 (D. Weber)
! Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version
! of the upper w-p radiation condition (fftopt=2).
!
! 08/01/97 (Zonghui Huo)
! Added Kain-fritsch cumulus parameterization scheme.
!
! 11/06/97 (D. Weber)
! Added three additional levels to the mapfct array. The three
! levels (4,5,6) represent the inverse of the first three in order.
! The inverse map factors are computed to improve efficiency.
!
! 12/05/97 (K. Brewster)
! Added argument, nt, so that routines that do not require more
! than one time level can be initialized using less memory.
!
! 4/15/1998 (Donghai Wang)
! Added the source terms to the right hand terms of the qc,qr,qi,qs
! equations due to the K-F cumulus parameterization.
!
! 4/15/1998 (Donghai Wang)
! Added the running average vertical velocity (array w0avg)
! for the K-F cumulus parameterization scheme.
!
! 11/18/98 (Keith Brewster)
! Changed pibar to ppi (full pi).
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
! 11/06/2001 (Yunheng Wang)
! Added mpupdatei calling for ict/icb to solve radiation forcing
! difference for MPI run.
!
! 2002/02/28 (Gene Bassett)
! Replaced mpupdatei for ict/icb to a call to mpmax0.
!
! 13 March 2002 (Eric Kemp)
! Added arrays for WRF BMJ cumulus scheme.
!
! 04/10/2002 (Yunheng Wang)
! Subsituted mpupdatei calls for mpmax0 calls again
! because mpmax0 calls were not correct.
!
!-----------------------------------------------------------------------
!
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nts Number of time levels to be initialized.
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space
! (m)
! hterain The height of terrain (m)
!
! mapfct Map factors at scalar, u and v points
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
! aj3x Avgx of the coordinate transformation Jacobian d(zp)/dz
! aj3y Avgy of the coordinate transformation Jacobian d(zp)/dz
! aj3z Avgz of the coordinate transformation Jacobian d(zp)/dz
! trigs1 Array containing pre-computed trig function for fftopt=1.
! trigs1 Array containing pre-computed trig function for fftopt=1.
! ifax1 Array containing the factors of nx for fftopt=1.
! ifax2 Array containing the factors of ny for fftopt=1.
!
! vwork1 2-D work array for fftopt=2.
! vwork2 2-D work array for fftopt=2.
! wsave1 Work array for fftopt=2.
! wsave2 Work array for fftopt=2.
!
! OUTPUT:
!
! u x component of velocity at all time levels (m/s)
! v y component of velocity at all time levels (m/s)
! w Vertical component of Cartesian velocity
! at all time levels (m/s)
! wcont Contravariant vertical velocity (m/s)
! ptprt Perturbation potential temperature at all time levels
! (K)
! pprt Perturbation pressure at all time levels (Pascal)
! qv Water vapor specific humidity at all time levels
! (kg/kg)
! qc Cloud water mixing ratio at all time levels (kg/kg)
! qr Rainwater mixing ratio at all time levels (kg/kg)
! qi Cloud ice mixing ratio at all time levels (kg/kg)
! qs Snow mixing ratio at all time levels (kg/kg)
! qh Hail mixing ratio at all time levels (kg/kg)
! tke Turbulent Kinetic Energy ((m/s)**2)
!
! udteb Time tendency of u field at east boundary (m/s**2)
! udtwb Time tendency of u field at west boundary (m/s**2)
!
! vdtnb Time tendency of v field at north boundary (m/s**2)
! vdtsb Time tendency of v field at south boundary (m/s**2)
!
! pdteb Time tendency of pprt field at east boundary (PASCAL/s)
! pdtwb Time tendency of pprt field at west boundary (PASCAL/s)
! pdtnb Time tendency of pprt field at north boundary
! (PASCAL/s)
! pdtsb Time tendency of pprt field at south boundary
! (PASCAL/s)
!
! ubar Base state zonal velocity component (m/s)
! vbar Base state meridional velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! ptbari Inverse Base state potential temperature (K)
! pbari Inverse Base state pressure (Pascal)
! rhostr Base state density rhobar times j3 (kg/m**3)
! rhostri Inverse base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg)
! ppi Exner function.
! csndsq Sound wave speed squared.
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsfc Skin temperature at the ground or ocean surface (K)
! tsoil Deep soil temperature (K) (in deep 1 m layer)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
! snowdpth Snow depth (m)
! qvsfc Effective S.H. at sfc.
!
! ptcumsrc Source term in pt-equation due to cumulus parameterization
! qcumsrc Source term in water equations due to cumulus parameterization
! kfraincv K-F convective rainfall (cm)
! nca K-F counter for CAPE release
! cldefi BMJ cloud efficiency
! xland BMJ land/sea mask
! bmjraincv BMJ convective rainfall (cm)
!
! radfrc Radiation forcing (K)
! radsw Solar radiation reaching the surface
! rnflx Net absorbed radiation by the surface
!
! raing Grid supersaturation rain
! rainc Cumulus convective rain
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3 Temporary work array.
! tem4 Temporary work array.
! tem5 Temporary work array.
! tem6 Temporary work array.
! tem7 Temporary work array.
! tem8 Temporary work array.
! tem9 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nts ! Number of time levels to be initialized.
INTEGER :: tpast ! Index of time level for the past time.
INTEGER :: tpresent ! Index of time level for the present
! time.
INTEGER :: tfuture ! Index of time level for the future
! time.
INTEGER :: nx,ny,nz ! The number of grid points in 3
! directions
REAL :: x (nx) ! x-coord. of the physical and compu-
! tational grid. Defined at u-point.
REAL :: y (ny) ! y-coord. of the physical and compu-
! tational grid. Defined at v-point.
REAL :: z (nz) ! z-coord. of the computational grid.
! Defined at w-point on the staggered
! grid.
REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at
! w-point of the staggered grid.
REAL :: hterain(nx,ny) ! The height of the terrain.
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as - d( zp )/d( x ).
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as - d( zp )/d( y ).
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! defined as d( zp )/d( z ).
REAL :: aj3x (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL :: aj3y (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL :: aj3z (nx,ny,nz) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL :: j3inv (nx,ny,nz) ! Inverse of j3
REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
! function for fftopt=1.
INTEGER :: ifax1(13) ! Array containing the factors of nx for
! fftopt=1.
INTEGER :: ifax2(13) ! Array containing the factors of ny for
! fftopt=1.
REAL :: vwork1 (nx+1,ny+1) ! 2-D work array for fftopt=2.
REAL :: vwork2 (ny,nx+1) ! 2-D work array for fftopt=2.
REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2.
REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2.
REAL :: u (nx,ny,nz,nts) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nts) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nts) ! Total w-velocity (m/s)
REAL :: wcont (nx,ny,nz) ! Contravariant vertical velocity (m/s)
REAL :: ptprt (nx,ny,nz,nts) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nts) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nts) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nts) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nts) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nts) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nts) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nts) ! Hail mixing ratio (kg/kg)
REAL :: tke (nx,ny,nz,nts) ! Turbulent Kinetic Energy ((m/s)**2)
REAL :: udteb (ny,nz) ! T-tendency of u at e-boundary (m/s**2)
REAL :: udtwb (ny,nz) ! T-tendency of u at w-boundary (m/s**2)
REAL :: vdtnb (nx,nz) ! T-tendency of v at n-boundary (m/s**2)
REAL :: vdtsb (nx,nz) ! T-tendency of v at s-boundary (m/s**2)
REAL :: pdteb (ny,nz) ! T-tendency of pprt at e-boundary
! (PASCAL/s)
REAL :: pdtwb (ny,nz) ! T-tendency of pprt at w-boundary
! (PASCAL/s)
REAL :: pdtnb (nx,nz) ! T-tendency of pprt at n-boundary
! (PASCAL/s)
REAL :: pdtsb (nx,nz) ! T-tendency of pprt at s-boundary
! (PASCAL/s)
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: ptbari(nx,ny,nz) ! Inverse Base state pot. temperature (K)
REAL :: pbari (nx,ny,nz) ! Inverse Base state pressure (Pascal).
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: rhostri(nx,ny,nz) ! Inv. base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity(kg/kg)
REAL :: ppi (nx,ny,nz) ! Exner function.
REAL :: csndsq(nx,ny,nz) ! Sound wave speed squared.
INTEGER :: nstyps ! Number of soil types
INTEGER :: soiltyp(nx,ny,nstyps) ! Soil types at grid points
REAL :: stypfrct(nx,ny,nstyps) ! Fraction of soil types
INTEGER :: vegtyp (nx,ny) ! Vegetation type
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: roufns (nx,ny) ! Surface roughness
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: tsfc (nx,ny,0:nstyps) ! Temperature at ground (K)
! (in top 1 cm layer)
REAL :: qvsfc (nx,ny,0:nstyps) ! Effective qv at sfc.
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
! (in deep 1 m layer)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture in the top
! 1 cm layer
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture in the deep
! 1 m layer
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: ptcumsrc(nx,ny,nz) ! Source term in pt-equation due
! to cumulus parameterization
REAL :: qcumsrc(nx,ny,nz,5) ! Source term in water equations due
! to cumulus parameterization:
! qcumsrc(1,1,1,1) for qv equation
! qcumsrc(1,1,1,2) for qc equation
! qcumsrc(1,1,1,3) for qr equation
! qcumsrc(1,1,1,4) for qi equation
! qcumsrc(1,1,1,5) for qs equation
REAL :: w0avg(nx,ny,nz) ! a closing running average vertical
! velocity in 10min for K-F scheme
REAL :: kfraincv(nx,ny) ! K-F convective rainfall (cm)
INTEGER :: nca(nx,ny) ! K-F counter for CAPE release
!EMK BMJ
REAL,INTENT(INOUT) :: cldefi(nx,ny) ! BMJ cloud efficiency
REAL,INTENT(INOUT) :: xland(nx,ny) ! BMJ land mask
! (1.0 = land, 2.0 = sea)
REAL,INTENT(INOUT) :: bmjraincv(nx,ny) ! BMJ convective rainfall (cm)
!EMK END
REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s)
REAL :: radsw (nx,ny) ! Solar radiation reaching the surface
REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface
REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2))
REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2))
REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m*s**2))
REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s))
REAL :: raing(nx,ny) ! Grid supersaturation rain
REAL :: rainc(nx,ny) ! Cumulus convective rain
REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precipitation rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: tem1 (nx,ny,nz) ! Temporary work array.
REAL :: tem2 (nx,ny,nz) ! Temporary work array.
REAL :: tem3 (nx,ny,nz) ! Temporary work array.
REAL :: tem4 (nx,ny,nz) ! Temporary work array.
REAL :: tem5 (nx,ny,nz) ! Temporary work array.
REAL :: tem6 (nx,ny,nz) ! Temporary work array.
REAL :: tem7 (nx,ny,nz) ! Temporary work array.
REAL :: tem8 (nx,ny,nz) ! Temporary work array.
REAL :: tem9 (nx,ny,nz) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
REAL :: temp
REAL :: alatpro(2)
REAL :: sclf
REAL :: dxscl ! Model x-direction grid spacing
! normalized by the map scale
! dxscl=dx/sclf
REAL :: dyscl ! Model y-direction grid spacing
! normalized by the map scale
! dyscl=dy/sclf
REAL :: xs,ys, swx,swy, ctrx, ctry
REAL zpmax
INTEGER :: mptag
!wdt update
REAL :: rmin,rmax
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
INCLUDE 'phycst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! If initopt = 1, initialize the model fields using intialization
! routines. Typically from an analytical
! definition of initial perturbations.
!
! If initopt = 2, initialize the model fields from a restart file
! produced by a previous model run.
!
! If initopt = 3, initialize the model fields by reading in an
! external data file.
!
!-----------------------------------------------------------------------
!
IF( nts > 1 ) THEN
tpast=1
tpresent=2
tfuture=3
ELSE
tpast=1
tpresent=1
tfuture=1
END IF
IF( initopt == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! Initialization of model GRiD definition arrays.
!
!-----------------------------------------------------------------------
!
CALL inigrd(nx,ny,nz, &
x,y,z,zp,hterain,mapfct,j1,j2,j3,tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! Define the base state atmospheric variables.
!
!-----------------------------------------------------------------------
!
IF ((inibasopt == 1) .AND. (max_fopen < nprocs)) THEN
CALL wrtcomment("ERROR: for message passing version with "// &
"inibasopt=1, max_fopen (in arps.input)",1)
CALL arpsstop ("arpsstop called from initgrdvar due to nproc_x-y &
& nproc_x*nproc_y (in arps.input).",1)
END IF
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL inibase(nx,ny,nz, ubar,vbar,ptbar,pbar,ptbari,pbari, &
rhostr,rhostri,qvbar, &
x,y,z,zp,j3, tem1,tem2,tem3,tem4)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
!
!-----------------------------------------------------------------------
!
! Initialize time dependent model variables.
!
!-----------------------------------------------------------------------
!
CALL initdvr(nx,ny,nz,nts, &
ubar,vbar,ptbar,pbar,rhostr,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptcumsrc,qcumsrc,raing,rainc,prcrate,tem1)
!
!-----------------------------------------------------------------------
!
! Set the time tendencies of u, v and pprt on the lateral boundaries
! to zero for the initial time step
!
! These tendencies will be used by lateral boundary condition options
! 4 and 5 only.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
udteb(j,k) = 0.0
udtwb(j,k) = 0.0
pdteb(j,k) = 0.0
pdtwb(j,k) = 0.0
END DO
END DO
DO k=1,nz
DO i=1,nx
vdtnb(i,k) = 0.0
vdtsb(i,k) = 0.0
pdtnb(i,k) = 0.0
pdtsb(i,k) = 0.0
END DO
END DO
ELSE IF( initopt == 2) THEN ! restart run
!
!-----------------------------------------------------------------------
!
! Read in restart data from a restart file to initialize fields
! u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh at time tpast and tpresent,
! the base state variables ubar,vbar,ptbar,pbar,rhostr,qvbar,
! and the time tendencies of variables at the lateral boundries.
!
! Fields at tfuture are set to equal to those at tpresent.
!
! This subroutine also sets the value of tstart to the restart
! data time. The value from input file in over-written.
!
!-----------------------------------------------------------------------
!
!EMK BMJ
CALL rstin(nx,ny,nz,nts, nstyps, exbcbufsz, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
udteb, udtwb, vdtnb, vdtsb, &
pdteb, pdtwb, pdtnb, pdtsb, &
ubar,vbar,ptbar,pbar,rhostr,qvbar, &
x,y,z,zp,hterain,mapfct,j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
radfrc,radsw,rnflx, &
raing,rainc,prcrate,exbcbuf,tem1,tem2)
!EMK END
!
!-----------------------------------------------------------------------
!
! Set map projection parameters which were not stored in restart
! data file.
!
!-----------------------------------------------------------------------
!
alatpro(1) = trulat1
alatpro(2) = trulat2
IF( sclfct /= 1.0) THEN
sclf = 1.0/sclfct
dxscl = dx*sclf
dyscl = dy*sclf
ELSE
sclf = 1.0
dxscl = dx
dyscl = dy
END IF
CALL setmapr( mapproj,sclf,alatpro,trulon )
CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
! swx = ctrx - (float(nx-3)/2.) * dxscl
! swy = ctry - (float(ny-3)/2.) * dyscl
swx = ctrx - (FLOAT(nproc_x*(nx-3))/2.) * dxscl
swy = ctry - (FLOAT(nproc_y*(ny-3))/2.) * dyscl
CALL setorig( 1, swx, swy)
! set up the model origin to the coord.
CALL setcornerll(nx,ny, x, y)
ELSE IF( initopt == 3) THEN ! External data input.
!
!-----------------------------------------------------------------------
!
! Read in externally supplied initial fields. These fields include
! u, v, w, ptprt, pprt, qv, qc, qr, qi, qs, and qh at time level
! tpresent, and the base state variables ubar, vbar, ptbar, pbar,
! rhostr,qvbar.
!
! Fields at tpast and tfuture are set to their values at tpresent.
!
!-----------------------------------------------------------------------
!
CALL extinit(nx,ny,nz,nts,nstyps, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
ubar,vbar,ptbar,pbar,rhostr,qvbar, &
x,y,z,zp,hterain,j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
!
!-----------------------------------------------------------------------
!
! Set map projection parameters which were not stored in history
! data file.
!
!-----------------------------------------------------------------------
!
alatpro(1) = trulat1
alatpro(2) = trulat2
IF( sclfct /= 1.0) THEN
sclf = 1.0/sclfct
dxscl = dx*sclf
dyscl = dy*sclf
ELSE
sclf = 1.0
dxscl = dx
dyscl = dy
END IF
CALL setmapr( mapproj,sclf,alatpro,trulon )
CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
! swx = ctrx - (float(nx-3)/2.) * dxscl
! swy = ctry - (float(ny-3)/2.) * dyscl
swx = ctrx - (FLOAT(nproc_x*(nx-3))/2.) * dxscl
swy = ctry - (FLOAT(nproc_y*(ny-3))/2.) * dyscl
CALL setorig( 1, swx, swy)
! set up the model origin to the coord.
IF ( mpfctopt /= 0 ) THEN
DO j=1,ny-1
DO i=1,nx-1
xs = 0.5*(x(i)+x(i+1))
ys = 0.5*(y(j)+y(j+1))
CALL xytomf( 1,1,xs,ys,mapfct(i,j,1) )
mapfct(i,j,4) = 1.0/mapfct(i,j,1)
mapfct(i,j,7) = mapfct(i,j,1)*mapfct(i,j,1)
mapfct(i,j,8) = 0.25*mapfct(i,j,1)
END DO
END DO
DO j=1,ny-1
DO i=1,nx
ys = 0.5*(y(j)+y(j+1))
CALL xytomf( 1,1,x(i),ys,mapfct(i,j,2) )
mapfct(i,j,5) = 1.0/mapfct(i,j,2)
END DO
END DO
DO j=1,ny
DO i=1,nx-1
xs = 0.5*(x(i)+x(i+1))
CALL xytomf( 1,1,xs,y(j),mapfct(i,j,3) )
mapfct(i,j,6) = 1.0/mapfct(i,j,3)
END DO
END DO
ELSE
DO k=1,7
DO j=1,ny
DO i=1,nx
mapfct(i,j,k) = 1.0
mapfct(i,j,8) = 0.25
END DO
END DO
END DO
END IF
CALL setcornerll(nx,ny, x, y)
!
!-----------------------------------------------------------------------
!
! Set the time tendencies of u, v and pprt on the lateral boundaries
! to zero for the initial time step
!
! These tendencies will be used by lateral boundary condition options
! 4 and 5 only.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
udteb(j,k) = 0.0
udtwb(j,k) = 0.0
pdteb(j,k) = 0.0
pdtwb(j,k) = 0.0
END DO
END DO
DO k=1,nz
DO i=1,nx
vdtnb(i,k) = 0.0
vdtsb(i,k) = 0.0
pdtnb(i,k) = 0.0
pdtsb(i,k) = 0.0
END DO
END DO
ELSE
WRITE(6,'(1x,a,i3,a)') &
'Invalid option for intialization, INITOPT was ',initopt, &
' Job stopped.'
CALL arpsstop ("arpsstop called from initgrdvar due to initopt",1)
END IF
!
!-----------------------------------------------------------------------
!
! Define the reversed vertical indeces of height cldh2m and cldm2l
! which represent the levels that separate high, middle, and low
! clouds in the computation of the solar radiation transfer
!
!-----------------------------------------------------------------------
!
ict = nz-2
icb = nz-2
DO k=nz-2,2,-1
tem1(1,1,k) = (zp(1,1,k)-zp(1,1,2))*(zflat-zrefsfc) &
/(zflat-zp(1,1,2))+zrefsfc
END DO
DO k=nz-2,2,-1
IF ( tem1(1,1,k) <= cldh2m) THEN
ict = k
EXIT
END IF
END DO
! for bit-for-bit MP accuracy:
CALL mpupdatei(ict, 1)
DO k=nz-2,2,-1
IF ( tem1(1,1,k) <= cldm2l) THEN
icb = k
EXIT
END IF
END DO
! for bit-for-bit MP accuracy:
CALL mpupdatei(icb, 1)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
j3inv(i,j,k)=1.0/j3(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=2,nx-1
aj3x(i,j,k)=0.5*(j3(i,j,k)+j3(i-1,j,k))
END DO
END DO
END DO
!call test_dump (aj3x,"XXXinit_aj3x",nx,ny,nz,1,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(aj3x,nx,ny,nz,ebc,wbc,1,mptag,tem2)
CALL mprecv2dew(aj3x,nx,ny,nz,ebc,wbc,1,mptag,tem2)
! CALL mpsend2dns(aj3x,nx,ny,nz,1,mptag,tem2)
! CALL mprecv2dns(aj3x,nx,ny,nz,1,mptag,tem2)
END IF
!call test_dump (aj3x,"XXXBinit_aj3x",nx,ny,nz,1,1)
CALL acct_interrupt(bc_acct)
CALL bcsu(nx,ny,nz,1,ny-1,1,nz-1,ebc,wbc,aj3x)
CALL acct_stop_inter
!call test_dump (aj3x,"XXXAinit_aj3x",nx,ny,nz,1,1)
DO k=1,nz-1
DO j=2,ny-1
DO i=1,nx-1
aj3y(i,j,k)=0.5*(j3(i,j,k)+j3(i,j-1,k))
END DO
END DO
END DO
!call test_dump (aj3y,"XXXinit_aj3y",nx,ny,nz,2,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
! CALL mpsend2dew(aj3y,nx,ny,nz,2,mptag,tem2)
! CALL mprecv2dew(aj3y,nx,ny,nz,2,mptag,tem2)
CALL mpsend2dns(aj3y,nx,ny,nz,nbc,sbc,2,mptag,tem2)
CALL mprecv2dns(aj3y,nx,ny,nz,nbc,sbc,2,mptag,tem2)
END IF
CALL acct_interrupt(bc_acct)
CALL bcsv(nx,ny,nz,1,nx-1,1,nz-1,nbc,sbc,aj3y)
CALL acct_stop_inter
!call test_dump (aj3y,"XXXAinit_aj3y",nx,ny,nz,2,1)
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
aj3z(i,j,k)=0.5*(j3(i,j,k)+j3(i,j,k-1))
END DO
END DO
END DO
CALL acct_interrupt(bc_acct)
CALL bcsw(nx,ny,nz,1,nx-1,1,ny-1,tbc,bbc,aj3z)
CALL acct_stop_inter
!call test_dump (aj3z,"XXXAinit_aj3z",nx,ny,nz,3,1)
!
!-----------------------------------------------------------------------
!
! Calculate the trigs1,trigs2,ifax1,ifax2 arrays by calling set99
!
! OR
!
! calculate wsave1,wsave2 by calling vcosi.
!
!-----------------------------------------------------------------------
!
DO i=1,13
ifax1(i) = 0
ifax2(i) = 0
END DO
DO i=1,3*(nx-1)/2+1
trigs1(i) = 0
END DO
DO j=1,3*(ny-1)/2+1
trigs2(j) = 0
END DO
DO j=1,ny+1
DO i=1,nx+1
vwork1(i,j) = 0.0
END DO
END DO
DO i=1,nx+1
DO j=1,ny
vwork2(j,i) = 0.0
END DO
END DO
DO i=1,3*(nx-1)+15
wsave2(i) = 0.0
END DO
DO i=1,3*(ny-1)+15
wsave1(i) = 0.0
END DO
IF ( tbc == 4 ) THEN ! set up the fft work arrays for use in the
! upper radiation boundary condition.
IF ( fftopt == 1 ) THEN ! set up periodic work arrays...
IF ( ny == 4 ) THEN ! set up trigs in x direction only
CALL set99(trigs1,ifax1,nx-1) ! NOTE: nx must be ODD!!!!
! and of special character...
! see fft99f.f for details....
ELSE IF ( nx == 4 ) THEN ! set up trigs in y direction only
CALL set99(trigs2,ifax2,ny-1) ! NOTE: ny must be ODD!!!!
! and of special character...
! see fft99f.f for details....
ELSE ! set up for 2-d transform...
CALL set99(trigs1,ifax1,nx-1) ! NOTE: nx must be ODD!!!!
! and of special character...
! see fft99f.f for details....
CALL set99(trigs2,ifax2,ny-1) ! NOTE: ny must be ODD!!!!
! and of special character...
! see fft99f.f for details....
END IF
ELSE IF ( fftopt == 2 ) THEN ! set up the cos fft arrays...
IF(ny == 4)THEN ! set up function in x direction only...
CALL vcosti(nx-1,wsave2) ! nx should be even.
ELSE IF(nx == 4)THEN ! set up function in y direction only...
CALL vcosti(ny-1,wsave1) ! ny should be even.
ELSE ! set up functions for 2-d transform...
CALL vcosti(ny-1,wsave1) ! ny should be even.
CALL vcosti(nx-1,wsave2) ! nx should be even.
END IF ! end of run type if block...
END IF ! end of fftopt if block.....
END IF
!
!-----------------------------------------------------------------------
!
! Find the lowest model layer (index rayklow) that is entirely or
! partially contained in the specified Rayleigh damping (sponge)
! layers.
!
! The Rayleigh damping is then applied only to layers with
! k greater than or equal to rayklow.
!
!-----------------------------------------------------------------------
!
rayklow = nz-1
DO k=nz-1,2,-1
zpmax = zp(1,1,k)
DO j=1,ny-1
DO i=1,nx-1
zpmax = MAX( zp(i,j,k), zpmax )
END DO
END DO
! for bit-for-bit accuracy with MP version:
rmin = zpmax
call mpmax0(zpmax,rmin)
IF( zpmax < zbrdmp ) THEN
rayklow = MAX(2, k+1)
EXIT
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Calculate Exner function and store in ppi
!
!-----------------------------------------------------------------------
!
CALL setppi(nx,ny,nz,nts,tpresent,pprt,pbar,ppi)
!
!-----------------------------------------------------------------------
!
! Calculate and store the sound wave speed squared in csndsq.
!
!-----------------------------------------------------------------------
!
IF(csopt == 1) THEN ! Original definition of sound speed.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
csndsq(i,j,k)= cpdcv*pbar(i,j,k)*j3(i,j,k)/rhostr(i,j,k)
END DO
END DO
END DO
ELSE IF(csopt == 2) THEN ! Original sound speed multiplied
! by a factor
temp = csfactr**2*cpdcv
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
csndsq(i,j,k)= temp * pbar(i,j,k)*j3(i,j,k)/rhostr(i,j,k)
END DO
END DO
END DO
ELSE ! Specified constant sound speed.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
csndsq(i,j,k)= csound * csound
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Fill the edges of base state arrays that are otherwise undefined.
! This is done for safty reason only.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny
vbar (nx,j,k) = vbar (nx-1,j,k)
END DO
END DO
DO k=1,nz-1
DO i=1,nx
ubar (i,ny,k) = ubar (i,ny-1,k)
END DO
END DO
DO i=1,nx
DO j=1,ny
ubar (i,j,nz) = ubar (i,j,nz-1)
vbar (i,j,nz) = vbar (i,j,nz-1)
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
ptbar (nx,j,k) = ptbar (nx-1,j,k)
pbar (nx,j,k) = pbar (nx-1,j,k)
ppi (nx,j,k) = ppi (nx-1,j,k)
qvbar (nx,j,k) = qvbar (nx-1,j,k)
csndsq(nx,j,k) = csndsq(nx-1,j,k)
rhostr(nx,j,k) = rhostr(nx-1,j,k)
END DO
END DO
DO k=1,nz-1
DO i=1,nx
ptbar (i,ny,k) = ptbar (i,ny-1,k)
pbar (i,ny,k) = pbar (i,ny-1,k)
ppi (i,ny,k) = ppi (i,ny-1,k)
qvbar (i,ny,k) = qvbar (i,ny-1,k)
csndsq(i,ny,k) = csndsq(i,ny-1,k)
rhostr(i,ny,k) = rhostr(i,ny-1,k)
END DO
END DO
DO i=1,nx
DO j=1,ny
ptbar (i,j,nz) = ptbar (i,j,nz-1)
pbar (i,j,nz) = pbar (i,j,nz-1)
ppi (i,j,nz) = ppi (i,j,nz-1)
qvbar (i,j,nz) = qvbar (i,j,nz-1)
csndsq(i,j,nz) = csndsq(i,j,nz-1)
rhostr(i,j,nz) = rhostr(i,j,nz-1)
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ptbari(i,j,k) = 1.0/ptbar(i,j,k)
pbari(i,j,k) = 1.0/pbar(i,j,k)
rhostri(i,j,k) = 1.0/rhostr(i,j,k)
END DO
END DO
END DO
IF ( sfcphy > 0 ) THEN
CALL initsfc(nx,ny,nz,nstyps, &
pbar,pprt(1,1,1,1), &
ptbar,ptprt(1,1,1,1), &
qvbar,qv(1,1,1,1), &
soiltyp,stypfrct, vegtyp, lai,roufns,veg,tem1, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc)
END IF
RETURN
END SUBROUTINE initgrdvar
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITDVR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE initdvr(nx,ny,nz,nts, & 1,9
ubar,vbar,ptbar,pbar,rhostr,qvbar, &
x,y,z,zp,hterain, j1,j2,j3, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptcumsrc,qcumsrc,raing,rainc,prcrate,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Initialize the model time dependent variables.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/1991.
!
! MODIFICATION HISTORY:
!
! 5/25/92 (M. Xue)
! Added full documentation.
!
! 5/03/92 (M. Xue)
! Further documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
! 01/28/95 (G. Bassett)
! Added pt0opt=5 (soup can shaped perturbation to ptprt).
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nts Number of time levels to be initialized.
!
! ubar Base state x-velocity component (m/s)
! vbar Base state y-velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhostr Base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg).
!
! x x-coordinate of grid points in computational space (m)
! y y-coordinate of grid points in computational space (m)
! z z-coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space(m)
! hterain Terrain height (m)
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! OUTPUT:
!
! u x-component of velocity at all time levels (m/s).
! v y-component of velocity at all time levels (m/s).
! w z-component of velocity at all time levels (m/s).
! ptprt Perturbation potential temperature at all time levels
! (K)
! pprt Perturbation pressure at all time levels (Pascal)
! qv Water vapor specific humidity at all time levels
! (kg/kg)
! qc Cloud water mixing ratio at all time levels (kg/kg)
! qr Rainwater mixing ratio at all time levels (kg/kg)
! qi Cloud ice mixing ratio at all time levels (kg/kg)
! qs Snow mixing ratio at all time levels (kg/kg)
! qh Hail mixing ratio at all time levels (kg/kg)
!
! ptcumsrc Source term in pt-equation due to cumulus parameterization
! qcumsrc Source term in water equations due to cumulus parameterization
! raing Grid supersaturation rain
! rainc Cumulus convective rain
! prcrate Precipitation ratesrain
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number of grid points in 3
! directions
INTEGER :: nts ! Number of time levels to be initialized.
INTEGER :: tpast ! Index of time level for the past time.
INTEGER :: tpresent ! Index of time level for the present
! time.
INTEGER :: tfuture ! Index of time level for the future
! time.
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal).
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity(kg/kg)
REAL :: x (nx) ! x-coord. of the physical and compu-
! tational grid. Defined at u-point.
REAL :: y (ny) ! y-coord. of the physical and compu-
! tational grid. Defined at v-point.
REAL :: z (nz) ! z-coord. of the computational grid.
! Defined at w-point on the staggered
! grid.
REAL :: zp (nx,ny,nz) ! Physical height coordinate defined at
! w-point of the staggered grid.
REAL :: hterain(nx,ny) ! Terrain height (m).
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dx.
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dy.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/dz.
REAL :: u (nx,ny,nz,nts) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nts) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nts) ! Total w-velocity (m/s)
REAL :: ptprt (nx,ny,nz,nts) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nts) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nts) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nts) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nts) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nts) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nts) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nts) ! Hail mixing ratio (kg/kg)
REAL :: ptcumsrc(nx,ny,nz) ! Source term in pt-equation due
! to cumulus parameterization
REAL :: qcumsrc(nx,ny,nz,5) ! Source term in water equations due
! to cumulus parameterization:
! qcumsrc(1,1,1,1) for qv equation
! qcumsrc(1,1,1,2) for qc equation
! qcumsrc(1,1,1,3) for qr equation
! qcumsrc(1,1,1,4) for qi equation
! qcumsrc(1,1,1,5) for qs equation
REAL :: raing (nx,ny) ! Grid supersaturation rain
REAL :: rainc (nx,ny) ! Cumulus convective rain
REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL :: tem1(nx,ny,nz) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: xs, ys, zs
REAL :: us, vs, ws, rhobar
REAL :: radnd , pi,pi2,pi4
INTEGER :: i,j,k, n
INTEGER :: iseed,ibgn,iend,jbgn,jend,kbgn,kend
INTEGER :: ebc1,wbc1,nbc1,sbc1
REAL :: amplitud
REAL :: knumx,lnumy,mnumz
REAL :: lnthx,lnthy,lnthz
REAL :: lambda,lambdah,lambda2
INTEGER :: mptag
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'phycst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF( nts > 1 ) THEN
tpast=1
tpresent=2
tfuture=3
ELSE
tpast=1
tpresent=1
tfuture=1
END IF
!
!-----------------------------------------------------------------------
!
! Specify the initial potential temperature perturbation.
!
!-----------------------------------------------------------------------
!
IF( pt0opt == 1 .OR. pt0opt == 6 ) THEN ! Bubble shaped perturbation
!
!-----------------------------------------------------------------------
!
! Define a potential temperature perturbation for a bubble-shaped
! disturbance.
!
!-----------------------------------------------------------------------
!
pi2 = 1.5707963267949
IF( ptpert0 == 0.0) THEN
DO k= 1,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
ptprt(i,j,k,1) = 0.0
END DO
END DO
END DO
ELSE
DO k= 1,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
!wdt multiple bubbles for timing tests:
!xs = (x(i)+x(i+1))*0.5
!ys = (y(j)+y(j+1))*0.5
xs = (x(i)+x(i+1))*0.5-x(1)
ys = (y(j)+y(j+1))*0.5-y(1)
zs = (zp(i,j,k)+zp(i,j,k+1))*0.5
IF( pt0rady < 0.0 .OR. runmod == 2 ) THEN
! 2-d bubble in x-z plane.
radnd = SQRT( ((xs-pt0ctrx)/pt0radx)**2 &
+ ((zs-pt0ctrz)/pt0radz)**2 )
ELSE IF( pt0radx < 0.0 .OR. runmod == 3 ) THEN
! 2-d bubble in y-z plane.
radnd = SQRT( ((ys-pt0ctry)/pt0rady)**2 &
+ ((zs-pt0ctrz)/pt0radz)**2 )
ELSE ! 3-d bubble
radnd = SQRT( ((xs-pt0ctrx)/pt0radx)**2 + &
((ys-pt0ctry)/pt0rady)**2 + &
((zs-pt0ctrz)/pt0radz)**2 )
END IF
IF(radnd >= 1.0) THEN
ptprt(i,j,k,1) = 0.0
ELSE
ptprt(i,j,k,1) = ptpert0*(COS(pi2*radnd )**2)
END IF
END DO
END DO
END DO
IF(pt0opt == 6) THEN ! Perturbation speficied in T'.
DO k= 1,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
ptprt(i,j,k,1) = ptprt(i,j,k,1)*(p0/pbar(i,j,k))**rddcp
END DO
END DO
END DO
END IF
END IF
ELSE IF( pt0opt == 2 .OR. pt0opt == 3 ) THEN
! Random field
!
!-----------------------------------------------------------------------
!
! Define a potential temperature perturbation by a random function.
! This ensures that the horizontal average of the perturbation in a
! specified domain is zero.
!
! Fill an array with zeros.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ptprt(i,j,k,1) = 0.0
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! The following parameters define the portion of domain to be
! initialized with a random potential temperture perturbation.
! Users can modify them to fit their needs.
!
!-----------------------------------------------------------------------
!
ibgn = 1
iend = nx-1
jbgn = 1
jend = ny-1
kbgn = 1
kend = nz-1
iseed = -100
DO k= kbgn,kend
CALL ranary(nx,ny,ibgn,iend,jbgn,jend, &
iseed,ptpert0,ptprt(1,1,k,1))
END DO
IF( pt0opt == 3 ) THEN ! Symmetric random perturbation
!
!-----------------------------------------------------------------------
!
! Create a random perturbation field symmetric about both the x and y
! axes.
!
!-----------------------------------------------------------------------
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx/2
ptprt(i,j,k,1) = ptprt(nx-i,j,k,1)
END DO
END DO
END DO
DO k=1,nz-1
DO i=1,nx-1
DO j=1,ny/2
ptprt(i,j,k,1) = ptprt(i,ny-j,k,1)
END DO
END DO
END DO
IF( nx == ny ) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = (ptprt(i,j,k,1)+ptprt(j,i,k,1))*0.5
END DO
END DO
END DO
DO k=1,nz-1
DO i=1,nx-1
DO j=1,ny-1
ptprt(i,j,k,1) = tem1(i,j,k)
END DO
END DO
END DO
END IF
END IF
ELSE IF( pt0opt == 4 ) THEN ! Bubble given in Skamarock and
! Klemp (1994)
pi2 = 1.5707963267949
DO k=1,nz-1
DO i=1,nx-1
DO j=1,ny-1
xs = (x(i)+x(i+1))*0.5
zs = (zp(i,j,k)+zp(i,j,k+1))*0.5
ptprt(i,j,k,1) = ptpert0 &
*SIN(pi2*2*zs/((nz-3)*dz))/(1+((xs-pt0ctrx)/pt0radx)**2)
END DO
END DO
END DO
ELSE IF( pt0opt == 5 ) THEN ! Soup can shaped perturbation
DO k= 1,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
xs = (x(i)+x(i+1))*0.5
ys = (y(j)+y(j+1))*0.5
zs = (zp(i,j,k)+zp(i,j,k+1))*0.5
IF ( ABS(zs-pt0ctrz)/pt0radz < 1 ) THEN
IF( pt0rady < 0.0 .OR. runmod == 2 ) THEN
! 2-d bubble in x-z plane.
radnd = ABS(xs-pt0ctrx)/pt0radx
ELSE IF( pt0radx < 0.0 .OR. runmod == 3 ) THEN
! 2-d bubble in y-z plane.
radnd = ABS(ys-pt0ctry)/pt0rady
ELSE ! 3-d bubble
radnd = SQRT( ((xs-pt0ctrx)/pt0radx)**2 + &
((ys-pt0ctry)/pt0rady)**2 )
END IF
IF(radnd >= 1.0) THEN
ptprt(i,j,k,1) = 0.0
ELSE
ptprt(i,j,k,1) = ptpert0
END IF
ELSE
ptprt(i,j,k,1) = 0.0
END IF
END DO
END DO
END DO
END IF
ebc1=0
wbc1=0
sbc1=0
nbc1=0
IF( ebc == 1 .OR.ebc == 2 .OR. ebc == 3 ) ebc1=ebc
IF( wbc == 1 .OR.wbc == 2 .OR. wbc == 3 ) wbc1=wbc
IF( sbc == 1 .OR.sbc == 2 .OR. sbc == 3 ) sbc1=sbc
IF( nbc == 1 .OR.nbc == 2 .OR. nbc == 3 ) nbc1=nbc
!call test_dump (ptprt,"XXXinit_ptprt",nx,ny,nz,0,0)
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend2dew(ptprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
CALL mprecv2dew(ptprt,nx,ny,nz,ebc1,wbc1,0,mptag,tem1)
CALL mpsend2dns(ptprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
CALL mprecv2dns(ptprt,nx,ny,nz,nbc1,sbc1,0,mptag,tem1)
END IF
CALL acct_interrupt(bc_acct)
CALL bcsclr(nx,ny,nz,dtbig, &
ptprt(1,1,1,1),ptprt(1,1,1,1), &
ptprt(1,1,1,1),tem1,tem1,tem1,tem1, &
ebc1,wbc1,nbc1,sbc1,tbc,bbc)
CALL acct_stop_inter
!call test_dump (ptprt,"XXXAinit_ptprt",nx,ny,nz,0,1)
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ptprt(i,j,k,n)=ptprt(i,j,k,1)
pprt(i,j,k,n)=0.0
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
u(i,j,k,n)=ubar(i,j,k)
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
v(i,j,k,n)=vbar(i,j,k)
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
w(i,j,k,n)=0.0
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qv(i,j,k,n)= qvbar(i,j,k)
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qc(i,j,k,n)=0.0
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qr(i,j,k,n)= 0.0
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qi(i,j,k,n)= 0.0
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qs(i,j,k,n)= 0.0
END DO
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qh(i,j,k,n)= 0.0
END DO
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ptcumsrc(i,j,k)= 0.0
qcumsrc(i,j,k,1)= 0.0
qcumsrc(i,j,k,2)= 0.0
qcumsrc(i,j,k,3)= 0.0
qcumsrc(i,j,k,4)= 0.0
qcumsrc(i,j,k,5)= 0.0
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
raing(i,j)= 0.0
rainc(i,j)= 0.0
prcrate(i,j,1)= 0.0
prcrate(i,j,2)= 0.0
prcrate(i,j,3)= 0.0
prcrate(i,j,4)= 0.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! The following setup is to overwrite the total u, v, w, pprt,
! ptprt, and qv for Beltrami flow initial conditions.
!
!-----------------------------------------------------------------------
!
IF ( pt0opt == 7 ) THEN
pi = 3.1415926535898
pi2 = 2*pi
pi4 = 4*pi
amplitud = 2.0 ! amplitude A=2 m/s
tmixopt = 1 ! constant viscosity option
tmixcst = 1.0 ! viscosity=1 m**2/s
tmixvert = 1.0 ! compute all mixing terms
wbc = 2 ! reset boundary conditions
ebc = 2 ! to periodical condition
nbc = 2
sbc = 2
tbc = 2
bbc = 2
lnthx = dx*(nx-3) ! length in x
lnthy = dy*(ny-3) ! length in y
lnthz = dz*(nz-3) ! length in z
knumx = pi4/lnthx ! wave number in x-dir
lnumy = pi2/lnthy ! wave number in y-dir
mnumz = pi2/lnthz ! wave number in z-dir
lambdah = knumx*knumx + lnumy*lnumy
lambda2 = lambdah + mnumz*mnumz
lambda = SQRT( lambda2 )
! print *, ' amplitude = ',amplitud
PRINT *, ' lnthx = ',lnthx, &
' lnthy = ',lnthy, &
' lnthz = ',lnthz
PRINT *, ' knumx = ',knumx, &
' lnumy = ',lnumy, &
' mnumz = ',mnumz
PRINT *, ' lambda1 = ',lambdah, &
' lambda2 = ',lambda2, &
' lambda = ',lambda
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
ys = 0.5*(y(j+1)+y(j))
zs = 0.5*(z(k+1)+z(k))
u(i,j,k,1) = -amplitud/lambdah &
*(lambda*lnumy &
*COS(knumx*x(i))*SIN(lnumy*ys)*SIN(mnumz*zs) &
+mnumz*knumx &
*SIN(knumx*x(i))*COS(lnumy*ys)*COS(mnumz*zs))
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
xs = 0.5*(x(i+1)+x(i))
zs = 0.5*(z(k+1)+z(k))
v(i,j,k,1)= amplitud/lambdah &
* (lambda*knumx &
*SIN(knumx*xs)*COS(lnumy*y(j))*SIN(mnumz*zs) &
-mnumz*lnumy &
*COS(knumx*xs)*SIN(lnumy*y(j))*COS(mnumz*zs))
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
xs = 0.5*(x(i+1)+x(i))
ys = 0.5*(y(j+1)+y(j))
w(i,j,k,1)=amplitud &
*COS(knumx*xs)*COS(lnumy*ys)*SIN(mnumz*z(k))
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
zs = 0.5*(z(k+1)+z(k))
us = 0.5*(u(i+1,j,k,1)+u(i,j,k,1))
vs = 0.5*(v(i,j+1,k,1)+v(i,j,k,1))
ws = 0.5*(w(i,j,k+1,1)+w(i,j,k,1))
rhobar = rhostr(i,j,k)/j3(i,j,k)
pprt(i,j,k,1) = p0-rhobar*(0.5*(us*us+vs*vs+ws*ws)+g*zs) &
- pbar(i,j,k)
END DO
END DO
END DO
DO n=1,nts
DO k=1,nz
DO j=1,ny
DO i=1,nx
u (i,j,k,n) = u(i,j,k,1)
v (i,j,k,n) = v(i,j,k,1)
w (i,j,k,n) = w(i,j,k,1)
pprt (i,j,k,n) = pprt(i,j,k,1)
ptprt(i,j,k,n) = 0.0
qv (i,j,k,n) = 0.0
END DO
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE initdvr
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EXTINIT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE extinit(nx,ny,nz,nts, nstyps, & 1,36
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
ubar,vbar,ptbar,pbar,rhostr,qvbar, &
x,y,z,zp,hterain,j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
uprt,vprt,qvprt,kmh,kmv,wbar, &
tem6,tem7,tem8)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read in the initial fields from an externally supplied initial
! data file.
!
! These fields include u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh
! at time level tpresent, and the base state variables
! ubar,vbar,ptbar,pbar,rhostr,qvbar.
!
! Fields at tpast and tfuture are set to their values at tpresent.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/7/1992
!
! MODIFICATION HISTORY:
!
! 3/25/94 (G. Bassett, M. Xue)
! Modified to read in regular binary history dumps.
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
! 01/14/1995 (M. Xue)
! Corrected where jacob was called. It should be called
! before do loop 90.
!
! 03/29/1997 (M. Xue)
! Modification made so that when values of mapproj,trulat1,trulat2,
! trulon,ctrlat,ctrlon in the input data do not match those in
! input file, the values in the data are used instead of those
! in input, as was done before. Warning messages will be printed.
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
! 03/19/2002 (Keith Brewster)
! Corrected print statement related to mis-match in data and user times.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nts Number of time levels to be initialized.
!
! OUTPUT:
!
! u x component of velocity at times tpast and tpresent
! (m/s)
! v y component of velocity at times tpast and tpresent
! (m/s)
! w Vertical component of Cartesian velocity at times
! tpast and tpresent (m/s)
! ptprt Perturbation potential temperature at times tpast and
! tpresent (K)
! pprt Perturbation pressure at times tpast and tpresent
! (Pascal)
!
! qv Water vapor specific humidity at times tpast and
! tpresent (kg/kg)
! qc Cloud water mixing ratio at times tpast and tpresent
! (kg/kg)
! qr Rainwater mixing ratio at times tpast and tpresent
! (kg/kg)
! qi Cloud ice mixing ratio at times tpast and tpresent
! (kg/kg)
! qs Snow mixing ratio at times tpast and tpresent (kg/kg)
! qh Hail mixing ratio at times tpast and tpresent (kg/kg)
! tke Turbulence kinetic energy (m**2/s)
!
! ubar Base state zonal velocity component (m/s)
! vbar Base state meridional velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhostr Base state density rhobar times j3 (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space(m)
! hterain The height of terrain (m)
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsfc Skin temperature at the ground or ocean surface (K)
! tsoil Deep soil temperature (K) (in deep 1 m layer)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
! snowdpth Snow depth (m)
! qvsfc Effective S.H. at sfc.
!
! ptcumsrc Source term in pt-equation due to cumulus parameterization
! qcumsrc Source term in water equations due to cumulus parameterization
! raing Grid supersaturation rain
! rainc Cumulus convective rain
! prcrate Precipitation rates
!
! radfrc Radiation forcing (K/s)
! radsw Solar radiation reaching the surface
! rnflx Net radiation flux absorbed by surface
!
! usflx Surface flux of u-momentum (kg/(m*s**2))
! vsflx Surface flux of v-momentum (kg/(m*s**2))
! ptsflx Surface heat flux (K*kg/(m**2 * s ))
! qvsflx Surface moisture flux of (kg/(m**2 * s))
!
! tstart The time when the time integration starts, which is set
! to the time of the restart data
!
! WORK ARRAYS:
!
! tem6 Temporary work array.
! tem7 Temporary work array.
! tem8 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nts ! Number of time levels to be initialized.
INTEGER :: tpast ! Index of time level for the past time.
INTEGER :: tpresent ! Index of time level for the present
! time.
INTEGER :: tfuture ! Index of time level for the future
! time.
REAL :: u (nx,ny,nz,nts) ! Total u-velocity (m/s)
REAL :: v (nx,ny,nz,nts) ! Total v-velocity (m/s)
REAL :: w (nx,ny,nz,nts) ! Total w-velocity (m/s)
REAL :: ptprt (nx,ny,nz,nts) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nts) ! Perturbation pressure (Pascal)
REAL :: qv (nx,ny,nz,nts) ! Water vapor specific humidity (kg/kg)
REAL :: qc (nx,ny,nz,nts) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz,nts) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz,nts) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz,nts) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz,nts) ! Hail mixing ratio (kg/kg)
REAL :: tke (nx,ny,nz,nts) ! Turbulence kinetic energy (m**2/s)
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: rhostr(nx,ny,nz) ! Base state density rhobar times j3.
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific
! humidity(kg/kg)
REAL :: x(nx) ! x-coord. of the physical and compu-
! tational grid. Defined at u-point.
REAL :: y(ny) ! y-coord. of the physical and compu-
! tational grid. Defined at v-point.
REAL :: z(nz) ! z-coord. of the computational grid.
! Defined at w-point on the staggered
! grid.
REAL :: zp(nx,ny,nz) ! Physical height coordinate defined at
! w-point of the staggered grid.
REAL :: hterain(nx,ny) ! Terrain height (m).
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dx.
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian
! -d(zp)/dy.
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian
! d(zp)/dz.
INTEGER :: nstyps ! Number of soil types
INTEGER :: soiltyp(nx,ny,nstyps) ! Soil types at grid points
REAL :: stypfrct(nx,ny,nstyps) ! Fraction of soil types
INTEGER :: vegtyp (nx,ny) ! Vegetation type
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: roufns (nx,ny) ! Surface roughness
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: tsfc (nx,ny,0:nstyps) ! Temperature at ground (K)
! (in top 1 cm layer)
REAL :: qvsfc (nx,ny,0:nstyps) ! Effective qv at sfc.
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
! (in deep 1 m layer)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture in the top
! 1 cm layer
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture in the deep
! 1 m layer
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: ptcumsrc(nx,ny,nz) ! Source term in pt-equation due
! to cumulus parameterization
REAL :: qcumsrc(nx,ny,nz,5) ! Source term in water equations due
! to cumulus parameterization:
! qcumsrc(1,1,1,1) for qv equation
! qcumsrc(1,1,1,2) for qc equation
! qcumsrc(1,1,1,3) for qr equation
! qcumsrc(1,1,1,4) for qi equation
! qcumsrc(1,1,1,5) for qs equation
REAL :: raing(nx,ny) ! Grid supersaturation rain
REAL :: rainc(nx,ny) ! Cumulus convective rain
REAL :: prcrate(nx,ny,4) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s)
REAL :: radsw (nx,ny) ! Solar radiation reaching the surface
REAL :: rnflx (nx,ny) ! Net radiation flux absorbed by surface
REAL :: usflx (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2))
REAL :: vsflx (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2))
REAL :: ptsflx(nx,ny) ! Surface heat flux (K*kg/(m*s**2))
REAL :: qvsflx(nx,ny) ! Surface moisture flux (kg/(m**2*s))
REAL :: uprt (nx,ny,nz) ! Temporary array
REAL :: vprt (nx,ny,nz) ! Temporary array
REAL :: qvprt(nx,ny,nz) ! Temporary array
REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: wbar (nx,ny,nz) ! Temporary array
REAL :: tem6 (nx,ny,nz) ! Temporary array
REAL :: tem7 (nx,ny,nz) ! Temporary array
REAL :: tem8 (nx,ny,nz) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: lengbf, lenfil
INTEGER :: i, j, k, ireturn, tim
REAL :: amax,amin
CHARACTER (LEN=80) :: runnamesv
INTEGER :: nocmnt_old ! Number of comment lines
CHARACTER (LEN=80 ) :: cmnt_old(50) ! String of comments on this model run
INTEGER :: nchin
REAL :: time_tmp
REAL :: tstopsv,thisdmpsv,mapprojsv,latitudsv,ctrlatsv,ctrlonsv
INTEGER :: monthsv,daysv,yearsv,hoursv,minutesv,secondsv
REAL :: trulat1sv,trulat2sv,trulonsv
REAL :: dxsv,dysv,strhoptsv,dzsv,dzminsv,zrefsfcsv,dlayer1sv,dlayer2sv, &
strhtunesv,zflatsv
REAL :: p0inv,eps,tvbar
INTEGER :: abstsec0,abstsec1
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'phycst.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
p0inv=1./p0
IF( nts > 1 ) THEN
tpast=1
tpresent=2
tfuture=3
ELSE
tpast=1
tpresent=1
tfuture=1
END IF
!
!-----------------------------------------------------------------------
!
! Read in the initial fields from the input data file.
!
!-----------------------------------------------------------------------
!
lengbf = 80
CALL strlnth( inigbf, lengbf )
lenfil = 80
CALL strlnth( inifile, lenfil )
IF (mp_opt > 0) THEN
runnamesv = inigbf
WRITE(inigbf,'(a,a,2i2.2)') &
runnamesv(1:lengbf),'_',loc_x,loc_y
lengbf = lengbf + 5
runnamesv = inifile
WRITE(inifile,'(a,a,2i2.2)') &
runnamesv(1:lenfil),'_',loc_x,loc_y
lenfil = lenfil + 5
END IF
tim = tpresent
runnamesv = runname
tstopsv = tstop
thisdmpsv = thisdmp
mapprojsv = mapproj
trulat1sv = trulat1
trulat2sv = trulat2
trulonsv = trulon
latitudsv = latitud
ctrlatsv = ctrlat
ctrlonsv = ctrlon
monthsv = month
daysv = day
yearsv = year
hoursv = hour
minutesv = minute
secondsv = second
dxsv = dx
dysv = dy
strhoptsv = strhopt
dzsv = dz
dzminsv = dzmin
zrefsfcsv = zrefsfc
dlayer1sv = dlayer1
dlayer2sv = dlayer2
strhtunesv = strhtune
zflatsv = zflat
nocmnt_old = nocmnt
DO i=1,nocmnt_old
cmnt_old(i)=cmnt(i)
END DO
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL dtaread(nx,ny,nz,nstyps, inifmt, nchin ,inigbf,lengbf, &
inifile,lenfil,time_tmp, &
x,y,z,zp,uprt,vprt,w(1,1,1,tim),ptprt(1,1,1,tim), &
pprt(1,1,1,tim),qvprt,qc(1,1,1,tim),qr(1,1,1,tim), &
qi(1,1,1,tim),qs(1,1,1,tim),qh(1,1,1,tim), &
tke(1,1,1,tim),kmh,kmv, &
ubar,vbar,wbar,ptbar,pbar,rhostr,qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc, tsoil, wetsfc, wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem6,tem7,tem8)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
!
!-----------------------------------------------------------------------
!
! To restore the original runname and comments, etc.
! from ARPS input file.
!
!-----------------------------------------------------------------------
eps = 0.0001
IF(mapproj /= mapprojsv .OR.ABS(trulat1sv-trulat1) > eps &
.OR.ABS(trulat2sv-trulat2) > eps &
.OR.ABS(trulonsv -trulon ) > eps) THEN
WRITE(6,'(/2(/5x,a),2(/5x,3(a,f13.3))/)') &
'WARNING: Map projection parameters in the input data do not ', &
'match those specified in input file. Values in data used.', &
'In data, trulat1=',trulat1 ,', trulat2=',trulat2, &
', trulon=',trulon, &
'In input, trulat1=',trulat1sv,', trulat2=',trulat2sv, &
', trulon=',trulonsv
END IF
IF(ABS(ctrlatsv-ctrlat) > eps .OR. ABS(ctrlonsv-ctrlon) > eps ) THEN
WRITE(6,'(/3(/5x,a),2(/5x,2(a,f13.3))/)') &
'WARNING: Central latitude and/or longitude of the grid ', &
'in the input data do not match those specified in input file.', &
'Values in data used.', &
'In data , ctrlat=',ctrlat ,', ctrlon=',ctrlon, &
'In input, ctrlat=',ctrlatsv,', ctrlon=',ctrlonsv
END IF
CALL ctim2abss(year,month,day,hour,minute,second, abstsec0)
CALL ctim2abss(yearsv,monthsv,daysv,hoursv,minutesv,secondsv, &
abstsec1)
abstsec0 = abstsec0 + nint(time_tmp)
abstsec1 = abstsec1 + nint(tstart)
IF ( abstsec0 /= abstsec1 ) THEN
WRITE(6,'(a,2(/a,1x,i4.4,5(a,i2.2),a,f10.3))') &
'WARNING: Data time is inconsistent with user input time.', &
' Data initime:', year,'-',month,'-',day,'.', &
hour,':',minute,':',second, &
' model data time: ', tstart, &
' User initime:', yearsv,'-',monthsv,'-',daysv,'.', &
hoursv,':',minutesv,':',secondsv, &
' model starting time: ', time_tmp
IF ( timeopt == 1 ) THEN
WRITE(6,'(a)') 'Program continues using user specified time'
year = yearsv
month = monthsv
day = daysv
hour = hoursv
minute = minutesv
second = secondsv
ELSE IF ( timeopt == 2 ) THEN
WRITE(6,'(a)') 'Program continues using data time'
tstart = time_tmp
ELSE
WRITE(6,'(a)') 'Program stopped in subroutine EXTINIT'
CALL arpsstop ("arpsstop called from extinit due to timeopt",1)
END IF
ELSE
year = yearsv
month = monthsv
day = daysv
hour = hoursv
minute = minutesv
second = secondsv
IF (myproc == 0) &
WRITE(6,'(a,i4.4,5(a,i2.2),3x,a,f10.3,a)') &
'Use specified initial time: ', &
year,'-',month,'-',day,'.',hour,':',minute,':',second, &
'and model starting time: ', tstart, 'seconds'
END IF
runname = runnamesv
tstop = tstopsv
thisdmp = thisdmpsv
dx = dxsv
dy = dysv
strhopt = strhoptsv
dz = dzsv
dzmin = dzminsv
zrefsfc = zrefsfcsv
dlayer1 = dlayer1sv
dlayer2 = dlayer2sv
strhtune = strhtunesv
zflat = zflatsv
nocmnt = nocmnt_old
DO i=1,nocmnt
cmnt(i)=cmnt_old(i)
END DO
IF( ireturn /= 0) THEN
WRITE(6,'(3(/1x,a))') &
'Error occured when reading history data ', &
inifile(1:lenfil)//' or '//inigbf(1:lengbf), &
'Program stopped in EXTINIT.'
CALL arpsstop ("arpsstop called from extinit in reading file",1)
END IF
!
!-----------------------------------------------------------------------
!
! Calculate rhostr & jacobians, set hterain.
!
!-----------------------------------------------------------------------
!
CALL jacob(nx,ny,nz,x,y,z,zp,j1,j2,j3)
DO k= 1,nz-1
DO j= 1,ny-1
DO i= 1,nx-1
tvbar=(ptbar(i,j,k)*((pbar(i,j,k)*p0inv)**rddcp))* &
((1.0+rvdrd*qvbar(i,j,k))/(1.0+qvbar(i,j,k)))
rhostr(i,j,k)= ABS(j3(i,j,k))*pbar(i,j,k)/(rd*tvbar)
!
! rhostr(i,j,k)= abs(j3(i,j,k))*pbar(i,j,k)/
! : ( Rd * ptbar(i,j,k)*(pbar(i,j,k)*p0inv)**rddcp )
!
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
rhostr(i,j, 1)=rhostr(i,j,2)
rhostr(i,j,nz-1)=rhostr(i,j,nz-2)
END DO
END DO
DO i=1,nx
DO j=1,ny
hterain(i,j) = zp(i,j,2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Put the 3-d variables into their 4-d arrays.
!
!-----------------------------------------------------------------------
!
tim = tpresent
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
u(i,j,k,tim)=uprt(i,j,k)+ubar(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
v(i,j,k,tim)=vprt(i,j,k)+vbar(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
qv(i,j,k,tim)=qvprt(i,j,k)+qvbar(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set the future values of variables to their current values.
! This is done primarily for safety reasons.
!
!-----------------------------------------------------------------------
!
IF(nts > 1) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
u (i,j,k,tfuture) = u (i,j,k,tpresent)
v (i,j,k,tfuture) = v (i,j,k,tpresent)
w (i,j,k,tfuture) = w (i,j,k,tpresent)
ptprt(i,j,k,tfuture) = ptprt(i,j,k,tpresent)
pprt (i,j,k,tfuture) = pprt (i,j,k,tpresent)
qv (i,j,k,tfuture) = qv (i,j,k,tpresent)
qc (i,j,k,tfuture) = qc (i,j,k,tpresent)
qr (i,j,k,tfuture) = qr (i,j,k,tpresent)
qi (i,j,k,tfuture) = qi (i,j,k,tpresent)
qs (i,j,k,tfuture) = qs (i,j,k,tpresent)
qh (i,j,k,tfuture) = qh (i,j,k,tpresent)
tke (i,j,k,tfuture) = tke (i,j,k,tpresent)
u (i,j,k,tpast ) = u (i,j,k,tpresent)
v (i,j,k,tpast ) = v (i,j,k,tpresent)
w (i,j,k,tpast ) = w (i,j,k,tpresent)
ptprt(i,j,k,tpast ) = ptprt(i,j,k,tpresent)
pprt (i,j,k,tpast ) = pprt (i,j,k,tpresent)
qv (i,j,k,tpast ) = qv (i,j,k,tpresent)
qc (i,j,k,tpast ) = qc (i,j,k,tpresent)
qr (i,j,k,tpast ) = qr (i,j,k,tpresent)
qi (i,j,k,tpast ) = qi (i,j,k,tpresent)
qs (i,j,k,tpast ) = qs (i,j,k,tpresent)
qh (i,j,k,tpast ) = qh (i,j,k,tpresent)
tke (i,j,k,tpast ) = tke (i,j,k,tpresent)
END DO
END DO
END DO
END IF
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
ptcumsrc(i,j,k)=0.0
qcumsrc(i,j,k,1)=0.0
qcumsrc(i,j,k,2)=0.0
qcumsrc(i,j,k,3)=0.0
qcumsrc(i,j,k,4)=0.0
qcumsrc(i,j,k,5)=0.0
END DO
END DO
END DO
!-----------------------------------------------------------------------
!
! Print out the max/min of initial arrays read in.
!
!-----------------------------------------------------------------------
!
tim = tpresent
IF (myproc ==0) &
WRITE(6,'(1x,a/)') 'Min. and max. of the data arrays read in:'
CALL a3dmax0(x,1,nx,1,nx,1,1,1,1, 1,1,1,1, amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'xmin = ', amin,', xmax =',amax
CALL a3dmax0(y,1,ny,1,ny,1,1,1,1, 1,1,1,1, amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'ymin = ', amin,', ymax =',amax
CALL a3dmax0(z,1,nz,1,nz,1,1,1,1, 1,1,1,1, amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'zmin = ', amin,', zmax =',amax
CALL a3dmax0(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'zpmin = ', amin,', zpmax =',amax
CALL a3dmax0(j3,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'j3min = ', amin,', j3max =',amax
CALL a3dmax0(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,', ubarmax =',amax
CALL a3dmax0(vbar,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,', vbarmax =',amax
CALL a3dmax0(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,', ptbarmax=',amax
CALL a3dmax0(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,', pbarmax =',amax
CALL a3dmax0(rhostr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'rhostrmin=', amin,', rhostrmax=',amax
CALL a3dmax0(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin,', qvbarmax=',amax
CALL a3dmax0( amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'umin = ', amin,', umax =',amax
CALL a3dmax0( amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'vmin = ', amin,', vmax =',amax
CALL a3dmax0( amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'wmin = ', amin,', wmax =',amax
CALL a3dmax0(ptprt(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,', ptprtmax=',amax
CALL a3dmax0(pprt(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,', pprtmax =',amax
CALL a3dmax0(qv(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qvmin = ', amin,', qvmax =',amax
CALL a3dmax0(qc(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qcmin = ', amin,', qcmax =',amax
CALL a3dmax0(qr(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qrmin = ', amin,', qrmax =',amax
CALL a3dmax0(qi(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qimin = ', amin,', qimax =',amax
CALL a3dmax0(qs(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qsmin = ', amin,', qsmax =',amax
CALL a3dmax0(qh(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qhmin = ', amin,', qhmax =',amax
CALL a3dmax0(tke(1,1,1,tim),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'tkemin = ', amin,', tkemax =',amax
! IF ( sfcphy.gt.0 ) THEN
! CALL a3dmax0(tsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
! write(6,'(1x,2(a,e13.6))') 'tsfcmin = ',amin,', tsfcmax =',amax
! CALL a3dmax0(tsoil,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
! write(6,'(1x,2(a,e13.6))') 'tsoilmin= ',amin,', tsoilmax =',amax
! CALL a3dmax0(wetsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
! write(6,'(1x,2(a,e13.6))') 'wetsmin = ',amin,', wetsmax =',amax
! CALL a3dmax0(wetdp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
! write(6,'(1x,2(a,e13.6))') 'wetdmin = ',amin,', wetdmax =',amax
! CALL a3dmax0(wetcanp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
! write(6,'(1x,2(a,e13.6))') 'wetcmin = ',amin,', wetcmax =',amax
! ENDIF
CALL a3dmax0(ptcumsrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'ptcummin= ', amin,', ptcummax=',amax
CALL a3dmax0(qcumsrc(1,1,1,1),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qvcummin= ', amin,', qvcummax=',amax
CALL a3dmax0(qcumsrc(1,1,1,2),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qccummin= ', amin,', qccummax=',amax
CALL a3dmax0(qcumsrc(1,1,1,3),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qrcummin= ', amin,', qrcummax=',amax
CALL a3dmax0(qcumsrc(1,1,1,4),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qicummin= ', amin,', qicummax=',amax
CALL a3dmax0(qcumsrc(1,1,1,5),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
IF (myproc ==0) &
WRITE(6,'(1x,2(a,e13.6))') 'qscummin= ', amin,', qscummax=',amax
RETURN
END SUBROUTINE extinit
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITSFC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE initsfc(nx,ny,nz,nstyps, & 1,79
pbar,pprt,ptbar,ptprt,qvbar,qv, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
ndvi, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Initialize the surface characteristics data and soil model
! variables according to option parameters sfcdat and soilinit.
!
! The surface and soil data files are sequential binary files.
!
! Their structures should be:
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 02/17/94
!
! MODIFICATION:
!
! 02/07/1995 (Yuhe Liu)
! Added a new 2-D permanent array, veg(nx,ny), to the soil model and
! at the same time delete the table data array veg(14).
!
! 01/29/1995 (V. Wong and X. Song)
! Add a flag wtrexist in "initsfc".
!
! 12/04/1997 (Yuhe Liu)
! Set the soil variables to their default value read in from input
! file if they are not included in soilinit file.
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
! 2000/01/07 (Gene Bassett)
! Modified the behavior for the case where soil data read in from file
! and also present in history file. If soilinit=2, initopt=3 and
! sfcin=1, for soil variables that are not in soildata file the values
! in the history file are used.
!
!-----------------------------------------------------------------------
!
! INPUT and OUTPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the z-direction (sfc/top)
!
! pbar Base state pressure
! pprt Preturbation pressure
!
! ptbar Base state potential temperature
! ptprt Preturbation potential temperature
!
! qvbar Base state water vapor mixing ratio
! qv Water vapor mixing ratio
!
! soiltyp Soil type at the horizontal grid points
! stypfrct Soil type fraction
! vegtyp Vegetation type at the horizontal grid points
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsfc Temperature at ground (K) (in top 1 cm layer)
! tsoil Deep soil temperature (K) (in deep 1 m layer)
! wetsfc Surface soil moisture in the top 1 cm layer
! wetdp Deep soil moisture in the deep 1 m layer
! wetcanp Canopy water amount
! qvsfc Effective S.H. at sfc.
!
! TEMPORARY WORKING ARRAY
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations. (Local Variables)
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL :: pi
PARAMETER ( pi = 3.141592654)
INTEGER :: nx, ny, nz
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: ptbar (nx,ny,nx) ! Base state potential temperature (K)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: qvbar (nx,ny,nz) ! Base state specific humidity (kg/kg)
REAL :: qv (nx,ny,nz) ! Total specific humidity (kg/kg)
INTEGER :: nstyps ! Number of soil types
INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type at each point
REAL :: stypfrct(nx,ny,nstyps) ! Fraction of soil types
INTEGER :: vegtyp (nx,ny) ! Vegetation type at each point
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: roufns (nx,ny) ! Surface roughness
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: ndvi (nx,ny) ! NDVI
REAL :: tsfc (nx,ny,0:nstyps) ! Ground surface temperature (K)
REAL :: tsoil (nx,ny,0:nstyps) ! Deep soil temperature (K)
REAL :: wetsfc (nx,ny,0:nstyps) ! Ground surface soil moisture
REAL :: wetdp (nx,ny,0:nstyps) ! Deep soil moisture
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy moisture
REAL :: snowdpth(nx,ny) ! Snow depth (m)
REAL :: qvsfc (nx,ny,0:nstyps) ! Effective surface specific
! humidity (kg/kg)
!
!-----------------------------------------------------------------------
!
! Define local variables
!
!-----------------------------------------------------------------------
!
REAL :: psfc ! Surface pressure
REAL :: rhgs ! The relative humidity at ground surface
REAL :: qvsatts ! Saturated specific humidity at surface
REAL :: wrmax ! Maximum value for canopy moisture, wetcanp
REAL :: pterm ! Temporary variable, real positive term flag
REAL :: p0inv ! Inverse of p0 (1000 mb)
REAL :: tema ! Temporary variable
REAL :: temb ! Temporary variable
!
!-----------------------------------------------------------------------
!
! Global constants and parameters, most of them specify the
! model run options.
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'indtflg.inc'
INCLUDE 'phycst.inc'
INCLUDE 'soilcst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'mp.inc' ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,is,imid,jmid
INTEGER :: tsfcin
INTEGER :: tsoilin
INTEGER :: wsfcin
INTEGER :: wdpin
INTEGER :: wcanpin
INTEGER :: snowdin
INTEGER :: mptag, astat
INTEGER, allocatable :: mpitmp(:)
REAL, allocatable :: mprtmp(:)
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF (mp_opt > 0 .AND. (sfcdat == 1 .OR. soilinit == 1)) THEN
ALLOCATE(mpitmp(MAX(nx,ny)),stat=astat)
IF (astat /= 0) THEN
WRITE (6,*) &
"INITSFC: ERROR allocating temporary array mpitmp."
CALL arpsstop ("arpsstop called from initsfc allocating mpitmp",1)
END IF
ALLOCATE(mprtmp(MAX(nx,ny)),stat=astat)
IF (astat /= 0) THEN
WRITE (6,*) &
"INITSFC: ERROR allocating temporary array mprtmp."
CALL arpsstop ("arpsstop called from initsfc allocating mprtmp",1)
END IF
END IF
p0inv = 1.0/p0
IF ( initopt == 2 ) THEN
WRITE (6, '(a,a/a/a/a/a)') &
'Model initialization option was set to restart. ', &
'All the surface and soil arrays should have been read', &
'in from the restart file. No more initialization is ', &
'done in subroutine INITSFC.'
wtrexist=0
DO j = 1, ny
DO i = 1, nx
DO is = 1,nstyp
IF (soiltyp(i,j,is) == 13) THEN
wtrexist=1
RETURN
END IF
END DO
END DO
END DO
RETURN
END IF
IF ( sfcdat == 1 ) THEN
nstyp = 1
DO j=1, ny-1
DO i=1, nx-1
soiltyp(i,j,1) = styp
vegtyp (i,j) = vtyp
lai (i,j) = lai0
roufns (i,j) = roufns0
veg (i,j) = veg0
END DO
END DO
ELSE IF (sfcdat == 2 .OR. (sfcdat == 3.AND.landin /= 1) ) THEN
!
!-----------------------------------------------------------------------
!
! Read the surface property data from file sfcdtfl (passed
! in globcst.inc).
!
!-----------------------------------------------------------------------
!
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL readsfcdt( nx,ny,nstyps,sfcdtfl,dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi )
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
ELSE IF(sfcdat == 3 .AND. landin == 1) THEN
WRITE(6,'(1x,a/,a/)') &
'Surface property data in the initialization file was used.', &
'No more initialization performed.'
ELSE
WRITE(6,'(1x,a,i3,a/)') &
'Invalid surface data input option. sfcdat =',sfcdat, &
'. Program stopped in INITSFC.'
CALL arpsstop ("arpsstop called from initsfc with sfcdat option",1)
END IF
print *, ' nstyps = ',nstyps
imid=(nx/2)+1
jmid=(ny/2)+1
DO is=1,nstyps
print *, ' Sample soil ( ',imid,jmid,') = ', &
soiltyp(imid,jmid,is),stypfrct(imid,jmid,is)
END DO
IF ( nstyp == 1 ) THEN
DO j=1,ny
DO i=1,nx
stypfrct(i,j,1) = 1.0
END DO
END DO
END IF
IF ( soilinit == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! All surface variables are specified uniformly in horizontal from
! the input namelist.
!
!-----------------------------------------------------------------------
DO j=1, ny-1
DO i=1, nx-1
psfc = .5 * ( (pbar(i,j,1)+pprt(i,j,1)) &
+ (pbar(i,j,2)+pprt(i,j,2)) )
IF( soiltyp(i,j,1) == 13) THEN ! For water
tsfc(i,j,0) = ptswtr0*(psfc*p0inv)**rddcp
ELSE ! For land
tsfc(i,j,0) = ptslnd0*(psfc*p0inv)**rddcp
END IF
tsoil (i,j,0) = tsoil0
wetsfc (i,j,0) = wetsfc0
wetdp (i,j,0) = wetdp0
wetcanp(i,j,0) = wetcanp0
snowdpth(i,j) = snowdpth0
END DO
END DO
DO is=1,nstyp
DO j=1,ny-1
DO i=1,nx-1
tsfc (i,j,is) = tsfc (i,j,0)
tsoil (i,j,is) = tsoil (i,j,0)
wetsfc (i,j,is) = wetsfc (i,j,0)
wetdp (i,j,is) = wetdp (i,j,0)
wetcanp(i,j,is) = wetcanp(i,j,0)
END DO
END DO
END DO
ELSE IF ((soilinit == 2).OR. &
(soilinit == 3.AND.sfcin /= 1) ) THEN
!
!-----------------------------------------------------------------------
!
! Read the soil variables from file soilinfl. soilinfl is
! passed in globcst.inc.
!
!-----------------------------------------------------------------------
!
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL readsoil(nx,ny,nstyps,soilinfl,dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,soiltyp)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
IF (sfcin /= 1) THEN
IF ( tsfcin /= 1 ) THEN
DO is=1,nstyp
DO j=1, ny-1
DO i=1, nx-1
psfc = .5 * ( (pbar(i,j,1)+pprt(i,j,1)) &
+ (pbar(i,j,2)+pprt(i,j,2)) )
IF( soiltyp(i,j,1) == 13) THEN ! For water
tsfc(i,j,is) = ptswtr0*(psfc*p0inv)**rddcp
ELSE ! For land
tsfc(i,j,is) = ptslnd0*(psfc*p0inv)**rddcp
END IF
END DO
END DO
END DO
END IF
IF ( tsoilin /= 1 ) THEN
DO is=1,nstyp
DO j=1, ny-1
DO i=1, nx-1
tsoil(i,j,is) = tsoil0
END DO
END DO
END DO
END IF
IF ( wsfcin /= 1 ) THEN
DO is=1,nstyp
DO j=1, ny-1
DO i=1, nx-1
wetsfc(i,j,is) = wetsfc0
END DO
END DO
END DO
END IF
IF ( wdpin /= 1 ) THEN
DO is=1,nstyp
DO j=1, ny-1
DO i=1, nx-1
wetdp(i,j,is) = wetdp0
END DO
END DO
END DO
END IF
IF ( wcanpin /= 1 ) THEN
DO is=1,nstyp
DO j=1, ny-1
DO i=1, nx-1
wetcanp(i,j,is) = wetcanp0
END DO
END DO
END DO
END IF
IF ( snowdin /= 1 ) THEN
DO j=1, ny-1
DO i=1, nx-1
snowdpth(i,j) = snowdpth0
END DO
END DO
END IF
END IF
ELSE IF(soilinit == 3.AND.sfcin == 1) THEN
WRITE(6,'(1x,a/,a/)') &
'Surface variable in the initialization file was used.', &
'No more initialization performed.'
IF ( nstyp == 1 ) THEN
DO j=1,ny-1
DO i=1,nx-1
tsfc (i,j,1) = tsfc (i,j,0)
tsoil (i,j,1) = tsoil (i,j,0)
wetsfc (i,j,1) = wetsfc (i,j,0)
wetdp (i,j,1) = wetdp (i,j,0)
wetcanp(i,j,1) = wetcanp(i,j,0)
END DO
END DO
END IF
ELSE IF(soilinit == 4) THEN
p0inv=1./p0
DO j=1,ny-1
DO i=1,nx-1
tema=0.5*((ptbar(i,j,2)+ptprt(i,j,2)) &
+(ptbar(i,j,1)+ptprt(i,j,1)))
psfc=0.5*((pbar(i,j,2)+pprt(i,j,2)) &
+(pbar(i,j,1)+pprt(i,j,1)))
temb = tema*(psfc*p0inv)**rddcp
tsfc (i,j,0) = temb + tsprt
tsoil (i,j,0) = temb + t2prt
wetsfc (i,j,0) = 0.0
wetdp (i,j,0) = 0.0
wetcanp(i,j,0) = wetcanp0
snowdpth(i,j) = snowdpth0
END DO
END DO
DO is=1,nstyp
DO j=1,ny-1
DO i=1,nx-1
tsfc (i,j,is) = tsfc (i,j,0)
tsoil (i,j,is) = tsoil(i,j,0)
wetsfc (i,j,is) = wgrat*wsat(soiltyp(i,j,is))
wetdp (i,j,is) = w2rat*wsat(soiltyp(i,j,is))
wetcanp(i,j,is) = wetcanp(i,j,0)
wetsfc(i,j,0 ) = wetsfc(i,j,0)+stypfrct(i,j,is)*wetsfc(i,j,is)
wetdp (i,j,0 ) = wetdp(i,j,0)+stypfrct(i,j,is)*wetdp(i,j,is)
END DO
END DO
END DO
ELSE
WRITE(6,'(1x,a,i3,a/)') &
'Invalid surface variable input option. soilinit =',soilinit, &
'. Program stopped in INITSFC.'
CALL arpsstop ("arpsstop called from initsfc with soilinit option",1)
END IF
IF ( moist /= 0 ) THEN
DO is=1,nstyp
DO j = 1, ny-1
DO i = 1, nx-1
wrmax = 0.2*1.0E-3*veg(i,j)*lai(i,j)
wetsfc (i,j,is) = MAX( wetsfc(i,j,is), 0.0 )
wetsfc (i,j,is) = MIN( wetsfc(i,j,is), wsat(soiltyp(i,j,is)) )
wetdp (i,j,is) = MAX( wetdp(i,j,is), 0.0 )
wetdp (i,j,is) = MIN( wetdp(i,j,is), wsat(soiltyp(i,j,is)) )
wetcanp(i,j,is) = MAX( wetcanp(i,j,is), 0.0 )
wetcanp(i,j,is) = MIN( wetcanp(i,j,is), wrmax )
END DO
END DO
END DO
DO j = 1, ny-1
DO i = 1, nx-1
wetsfc (i,j,0) = MAX( wetsfc(i,j,0), 0.0 )
wetdp (i,j,0) = MAX( wetdp(i,j,0), 0.0 )
wetcanp(i,j,0) = MAX( wetcanp(i,j,0), 0.0 )
END DO
END DO
ELSE
DO is=0,nstyp
DO j = 1, ny-1
DO i = 1, nx-1
wetsfc(i,j,is) = 0.0
wetdp(i,j,is) = 0.0
wetcanp(i,j,is) = 0.0
END DO
END DO
END DO
END IF
DO is=1,nstyp
DO j = 1, ny-1
DO i = 1, nx-1
psfc = .5 * ( (pbar(i,j,1)+pprt(i,j,1)) &
+ (pbar(i,j,2)+pprt(i,j,2)) )
qvsatts = f_qvsat( psfc, tsfc(i,j,is) )
IF ( soiltyp(i,j,is) == 12 .OR. soiltyp(i,j,is) == 13) THEN
! Ice and water
rhgs = 1.0
ELSE IF ( soiltyp(i,j,is) > 0 ) THEN
pterm = .5+SIGN(.5,wetsfc(i,j,is)-1.1*wfc(soiltyp(i,j,is)))
rhgs = pterm + (1.-pterm) &
* 0.25 * ( 1.-COS(wetsfc(i,j,is)*pi &
/(1.1*wfc(soiltyp(i,j,is)))) )**2
ELSE
rhgs = 0.0
END IF
qvsfc(i,j,is) = rhgs * qvsatts
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Ensure soil values are consistent with each other.
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
tsfc (i,j,0) = 0.0
tsoil (i,j,0) = 0.0
wetsfc (i,j,0) = 0.0
wetdp (i,j,0) = 0.0
wetcanp(i,j,0) = 0.0
END DO
END DO
DO is = 1,nstyp
DO j=1,ny-1
DO i=1,nx-1
tsfc (i,j,0) = tsfc (i,j,0) &
+ tsfc (i,j,is) * stypfrct(i,j,is)
tsoil (i,j,0) = tsoil (i,j,0) &
+ tsoil (i,j,is) * stypfrct(i,j,is)
wetsfc (i,j,0) = wetsfc (i,j,0) &
+ wetsfc (i,j,is) * stypfrct(i,j,is)
wetdp (i,j,0) = wetdp (i,j,0) &
+ wetdp (i,j,is) * stypfrct(i,j,is)
wetcanp(i,j,0) = wetcanp(i,j,0) &
+ wetcanp(i,j,is) * stypfrct(i,j,is)
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
IF ((vegtyp(i,j) == 14) .OR. (tsfc(i,j,0) > 273.16)) THEN
snowdpth(i,j) = 0. ! Make sure that snow cover is consistent
END IF ! with vegetation type and soil temperature.
END DO
END DO
DO j = 1, ny-1
DO i = 1, nx-1
qvsfc(i,j,0) = 0.0
END DO
END DO
DO is=1,nstyp
DO j = 1, ny-1
DO i = 1, nx-1
qvsfc(i,j,0) = qvsfc(i,j,0) + qvsfc(i,j,is)*stypfrct(i,j,is)
END DO
END DO
END DO
IF ( sfcdat == 1 ) THEN
DO is=1,nstyp
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend1diew(soiltyp(1,1,is),nx,ny,ebc,wbc,0,mptag,mpitmp)
CALL mprecv1diew(soiltyp(1,1,is),nx,ny,ebc,wbc,0,mptag,mpitmp)
CALL mpsend1dins(soiltyp(1,1,is),nx,ny,nbc,sbc,0,mptag,mpitmp)
CALL mprecv1dins(soiltyp(1,1,is),nx,ny,nbc,sbc,0,mptag,mpitmp)
CALL mpsend1dew(stypfrct(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(stypfrct(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(stypfrct(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(stypfrct(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
END IF
CALL acct_interrupt(bc_acct)
CALL bcis2d(nx,ny, soiltyp (1,1,is), ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, stypfrct(1,1,is), ebc,wbc,nbc,sbc)
CALL acct_stop_inter
END DO
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend1diew(vegtyp,nx,ny,ebc,wbc,0,mptag,mpitmp)
CALL mprecv1diew(vegtyp,nx,ny,ebc,wbc,0,mptag,mpitmp)
CALL mpsend1dins(vegtyp,nx,ny,nbc,sbc,0,mptag,mpitmp)
CALL mprecv1dins(vegtyp,nx,ny,nbc,sbc,0,mptag,mpitmp)
CALL mpsend1dew(lai,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(lai,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(lai,nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(lai,nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(roufns,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(roufns,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(roufns,nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(roufns,nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(veg,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(veg,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(veg,nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(veg,nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(snowdpth,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(snowdpth,nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(snowdpth,nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(snowdpth,nx,ny,nbc,sbc,0,mptag,mprtmp)
END IF
CALL acct_interrupt(bc_acct)
CALL bcis2d(nx,ny, vegtyp, ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, lai, ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, roufns, ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, veg, ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, snowdpth,ebc,wbc,nbc,sbc)
CALL acct_stop_inter
END IF
!call test_dump (roufns,"XXXAinit_roufns",nx,ny,1,0,1)
IF ( soilinit == 1 ) THEN
DO is=0,nstyp
IF (mp_opt > 0) THEN
CALL acct_interrupt(mp_acct)
CALL mpsend1dew(tsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(tsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(tsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(tsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(tsoil(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(tsoil(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(tsoil(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(tsoil(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(wetsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(wetsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(wetsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(wetsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(wetdp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(wetdp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(wetdp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(wetdp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(wetcanp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(wetcanp(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(wetcanp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(wetcanp(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mpsend1dew(qvsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mprecv1dew(qvsfc(1,1,is),nx,ny,ebc,wbc,0,mptag,mprtmp)
CALL mpsend1dns(qvsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
CALL mprecv1dns(qvsfc(1,1,is),nx,ny,nbc,sbc,0,mptag,mprtmp)
END IF
CALL acct_interrupt(bc_acct)
CALL bcs2d (nx,ny, tsfc (1,1,is), ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, tsoil (1,1,is), ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, wetsfc (1,1,is), ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, wetdp (1,1,is), ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, wetcanp(1,1,is), ebc,wbc,nbc,sbc)
CALL bcs2d (nx,ny, qvsfc (1,1,is), ebc,wbc,nbc,sbc)
CALL acct_stop_inter
END DO
END IF
IF (mp_opt > 0 .AND. (sfcdat == 1 .OR. soilinit == 1)) THEN
DEALLOCATE (mpitmp,stat=astat)
DEALLOCATE (mprtmp,stat=astat)
END IF
wtrexist=0
DO j = 1, ny
DO i = 1, nx
DO is = 1,nstyp
IF (soiltyp(i,j,is) == 13) THEN
wtrexist=1
RETURN
END IF
END DO
END DO
END DO
RETURN
WRITE (6,'(/a,i2/a)') &
' Read error in surface data file ' &
//soilinfl//' with the I/O unit ',sfcunit, &
'The model will STOP in subroutine INITSFC.'
CALL arpsstop ("arpsstop called from initsfc reading sfc data file",1)
WRITE (6,'(/a,i2/a)') &
' Read error in surface initial data file ' &
//soilinfl//' with the I/O unit ',sfcunit, &
'The model will STOP in subroutine INITSFC.'
CALL arpsstop ("arpsstop called from initsfc reading sfc data file-2",1)
END SUBROUTINE initsfc
!wdt update - init0 no longer necessary (arrays set to 0 after allocation)
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FLZERO ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE flzero(a,n) 103
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Fill in an entire array with zeros.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/1991.
!
! MODIFICATION HISTORY:
!
! 5/02/92 (M. Xue)
! Added full documentation.
!
! 5/03/92 (M. Xue)
! Further documentation.
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! n The dimension of 1-D array a.
!
! OUTPUT:
!
! a 1-D array filled with zeros.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, n
REAL :: a(n) ! 1-D array filled with zeros
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO i=1,n
a(i)=0.0
END DO
RETURN
END SUBROUTINE flzero
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INITLKTB ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE initlktb 2
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Initializes arrays used for lookup table functions.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Adwait Sathye
! 06/02/94
!
! MODIFICATION HISTORY:
!
! 2/17/97 (J. Zong)
! Corrected definitions of the power functions in loop 100.
!
! 2/24/97 (Jinxing Zong, Ming Xue and Yuhe Liu)
! Defined five pwr arrays for lookup tables to replace fractional
! power calculations in Tao ice microphysics.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
! pwr arrays filled with evenly spaced data points
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INCLUDE FILES
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! LOCAL VARIABLES
!
!-----------------------------------------------------------------------
!
INTEGER :: i
REAL :: temper
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO i = 0, 1001
pwr2046(i) = ( 0.05 * i / 1000.0 ) ** 0.2046
pwr525(i) = ( 0.05 * i / 1000.0 ) ** 0.525
pwr875(i) = ( 0.05 * i / 1000.0 ) ** 0.875
pwr1364(i) = ( 0.05 * i / 1000000.0 ) ** 0.1364
END DO
!-----------------------------------------------------------------------
!
! Build a look up table for the Latent heat of vaporation, calculated
! in the temperature range of (-100,50) degree Celcius
! according to formula
!
! Lv = lathv * (273.15/tbar)**(0.167+3.67E-4*tbar)
!
! Lv for t(K) is given by
!
! index = max( min(151, NINT(t-172.15)), 1)
! Lv = latheatv(index)
!
! where lathv = 2.500780e6 is set in phycst.inc, and t is the
!
!-----------------------------------------------------------------------
DO i=-100,50
temper = i+273.15
latheatv(i+101)=lathv * (273.15/temper)**(0.167+3.67E-4*temper)
END DO
!-----------------------------------------------------------------------
!
! Build lookup tables for T**0.81 and (rhobar/mu/psi**2)**one6th,
! where mu and psi stand for dynamic viscosity of air and diffusivity
! of water vapor in air, respectively. The temerature ranges from
! -100 to 50 degree Celcius, pressure from 0 to 1500 hPa, and air
! density from 0 to 1.5 kg/m**3 in the calculations of the tables.
! With the above ranges of T, p and air density, (rhobar/mu/psi**2)
! is between 0.0 and 3000.0.
!
!-----------------------------------------------------------------------
DO i = 0, 151
pwr81(i) = ( 173.15 + i ) ** 0.81
END DO
DO i = 0, 10001
pwr1666(i) = ( 3000.0 * i / 10000.0 ) ** 0.1666667
END DO
!-----------------------------------------------------------------------
!
! Lookup tables for (rhobar*q)**a, where q can be qr, qs or qh.
! rhobar is in g/cm***3. rhobar*q ranges from 0 to 50.e-6. The
! increment of the tables is 50.e-9.
!
!-----------------------------------------------------------------------
DO i = 0, 10001
pwr2 (i) = ( 0.05 * i / 10000000.0 ) ** 0.2
pwr0625 (i) = ( 0.05 * i / 10000000.0 ) ** 0.0625
pwr15625(i) = ( 0.05 * i / 10000000.0 ) ** 0.15625
END DO
RETURN
END SUBROUTINE initlktb
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GTSINLAT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE gtsinlat(nx,ny, x,y, sinlat, xs, ys, temxy) 1,1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculate the sin of the lattitude of each grid point, to be used
! in the calculation of latitude-dependent Coriolis parameters.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Ming Xue
! 3/21/95.
!
! 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)
! x x-coordinate of grid points in computational space (m)
! y y-coordinate of grid points in computational space (m)
!
! OUTPUT:
!
! sinlat Sin of latitude at each grid point
!
! WORK ARRAYS:
!
! xs x-coordinate at scalar points.
! ys y-coordinate at scalar points.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny ! Dimensions of model grid.
REAL :: x (nx) ! The x-coord. of the u points.
REAL :: y (ny) ! The y-coord. of the v points.
REAL :: sinlat(nx,ny) ! Sin of latitude at each grid point
REAL :: xs (nx) ! x-coord. of scalar points.
REAL :: ys (ny) ! y-coord. of scalar points.
REAL :: temxy (nx,ny) ! Work array.
REAL :: r2d
INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO i=1,nx-1
xs(i) = (x(i)+x(i+1))*0.5
END DO
xs(nx) = -xs(nx-2)+2.0*xs(nx-1)
DO j=1,ny-1
ys(j) = (y(j)+y(j+1))*0.5
END DO
ys(ny) = -ys(ny-2)+2.0*ys(ny-1)
!
CALL xytoll(nx,ny,xs,ys,sinlat,temxy)
r2d = ATAN(1.0)*4.0/180.0
DO j=1,ny
DO i=1,nx
sinlat(i,j) = SIN( r2d * sinlat(i,j) )
END DO
END DO
RETURN
END SUBROUTINE gtsinlat
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SETPPI ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE setppi(nx,ny,nz,nt,tlvl,pprt,pbar,ppi) 2
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nt,tlvl
REAL :: pprt(nx,ny,nz,nt)
REAL :: pbar(nx,ny,nz)
REAL :: ppi (nx,ny,nz)
!
INCLUDE 'phycst.inc'
!
INTEGER :: i,j,k
REAL :: p0inv
!
p0inv=1./p0
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ppi(i,j,k)=((pprt(i,j,k,tlvl)+pbar(i,j,k))*p0inv)**rddcp
END DO
END DO
END DO
!
RETURN
END SUBROUTINE setppi