!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CORDINTG ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!EMK BMJ
SUBROUTINE cordintg(mptr,frstep, nx,ny,nz,nxndg,nyndg,nzndg, & 3,69
rbufsz,nstyps,exbcbufsz, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth, &
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, 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,raing,rainc,prcrate,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
radfrc,radsw,rnflx,rad2d,radbuf,exbcbuf,bcrlx, &
usflx,vsflx,ptsflx,qvsflx, &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr, &
phydro,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 END
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the forward integration of all model time-dependent
! variables.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/01/1991.
!
! MODIFICATION HISTORY:
!
! 4/26/92 (K. Droegemeier)
! Added full documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 01/24/94 (Yuhe Liu)
! Added array tsfc, qvsfc, tsoil, wetsfc, wetdp, wetcanp, and some of
! temporary arrays to the arguments of subroutine SFCFLX, for
! prediction of the surface temperature and specific humidity
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 02/07/1995 (Yuhe Liu)
! Added a new 2-D array, veg(nx,ny), to the argument list
!
! 7/6/1995 (M. Xue)
! Pressure detrending is applied to the base grid (as opposed to
! nested grids) only.
!
! 10/31/95 (D. Weber)
! Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
! radiation condition.
!
! 01/25/96 (Donghai Wang & Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 01/31/1996 (V. Wong and X. Song)
! Added a parameter, qpfgfrq, that controls the frequency of calling
! the subroutine QPFGRID.
!
! 08/14/1996 (Yuhe Liu)
! Moved the definition of radiation buffer to the ARPS main program
! and passed it through argument list.
!
! 07/22/97 (D. Weber)
! Added wsave1,wsave2,vwork1,vwork2 for use in the upper w-p
! radiation condition (fftopt=2).
!
! 08/01/97 (Zonghui Huo)
! Added Kain-fritsch cumulus parameterization scheme.
!
! 10/21/97 (Donghai Wang)
! Using total density (rho) in the calculation of the pressure
! gradient force terms, and added the second order terms
! in the linerized buoyancy terms.
!
! 11/05/97 (D. Weber)
! Added phydro array for use in the bottom boundary condition for
! perturbation 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.
!
! 2/11/1998 (M. Xue)
! Corrections made to ensure no moist process is activatived when
! moist=0.
!
! 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/18/98 (D. Weber)
! Added precomputed averages of j3 in the x,y, and z directions
! to improve code efficiency.
!
! 8/31/1998 (K. Brewster)
! Added call to NUDGEALL for nudging to observation increments.
!
! 11/18/98 (Keith Brewster)
! Changed pibar to ppi (full pi).
!
! 12/09/1998 (Donghai Wang)
! Added the snow cover.
!
! 12/05/2001 (Yunheng Wang)
! Replaced sfcflx call with sfcphysics call which was a structure
! change by M. Xue.
!
! 07/10/2001 (K. Brewster)
! Added increment arrays to argument list and to call to nudgeall
!
! 07/23/2001 (K. Brewster)
! Added mptr to call to tinteg, needed for printing of diagnostic noise
! parameter.
!
! 13 March 2002 (Eric Kemp)
! Added arrays for WRF BMJ cumulus scheme.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! mptr Grid identifier.
! frstep Flag to determine if this is the initial time step.
! 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
!
! 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)
! wcont Contravariant vertical velocity (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 Turbulent Kinetic Energy ((m/s)**2)
! pbldpth Planetary boundary layer depth (m)
!
! 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)
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! 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.
!
! 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)
! 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 Jacobian d(zp)/dz
!
! 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.
! vwork2 2-D work array for fftopt=2.
! wsave1 Work array for fftopt=2.
! wsave2 Work array for fftopt=2.
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsfc Ground surface temperature (K)
! tsoil Deep soil temperature (K) (in deep 1 m layer)
! qvsfc Effective qv at sfc.
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
!
! ptsfc Ground surface potential temperature (K)
! qvsfc Effective S.H. at sfc.
!
! OUTPUT:
! u Updated (by dtbig) x velocity component at times tpast
! and tpresent (m/s).
! v Updated (by dtbig) y velocity component at times tpast
! and tpresent (m/s).
! w Updated (by dtbig) z velocity component at times tpast
! and tpresent (m/s).
! ptprt Updated (by dtbig) perturbation potential temperature
! at times tpast and tpresent (K).
! pprt Updated (by dtbig) perturbation pressure at times tpast
! and tpresent (Pascal).
! qv Updated (by dtbig) water vapor specific humidity at times
! tpast and tpresent (kg/kg).
! qc Updated (by dtbig) cloud water mixing ratio at times
! tpast and tpresent (kg/kg).
! qr Updated (by dtbig) rainwater mixing ratio at times tpast
! and tpresent (kg/kg).
! qi Updated (by dtbig) cloud ice mixing ratio at times tpast
! and tpresent (kg/kg).
! qs Updated (by dtbig) snow mixing ratio at times tpast and
! tpresent (kg/kg).
! qh Updated (by dtbig) hail mixing ratio at times tpast and
! tpresent (kg/kg).
!
! 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)
!
! 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
!
! tsfc Ground surface temperature (K)
! qvsfc Effective qv at sfc.
! 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
! raing Grid supersaturation rain
! rainc Cumulus convective rain
! 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/s)
! radsw Solar radiation reaching the ground surface (W/m**2)
! rnflx Net radiation flux at the ground surface (W/m**2)
! rad2d 2-D arrays for radiation calculation
! radbuf Buffer array to carry temporary working arrays for
! radiation calculation
!
! WORK ARRAYS:
!
! 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)
!
! 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)
!
! ptdteb Time tendency of ptprt field at east boundary (K/s)
! ptdtwb Time tendency of ptprt field at west boundary (K/s)
! ptdtnb Time tendency of ptprt field at north boundary(K/s)
! ptdtsb Time tendency of ptprt field at south boundary(K/s)
!
! sdteb Time tendency of a scalar at e-boundary
! sdtwb Time tendency of a scalar at w-boundary
! sdtnb Time tendency of a scalar at n-boundary
! sdtsb Time tendency of a scalar at s-boundary
!
! 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))
!
! 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 ! Force explicit declarations
INTEGER :: mptr ! Grid identifier.
INCLUDE 'timelvls.inc'
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
INTEGER :: nxndg,nyndg,nzndg ! Number of grid points in 3 directions
INTEGER :: frstep ! Indicator of the first time step call
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 (K)
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 :: pbldpth(nx,ny,nt) ! Planetary boundary layer depth (m)
REAL :: udteb (ny,nz) ! Time tendency of u at e-boundary (m/s**2)
REAL :: udtwb (ny,nz) ! Time tendency of u at w-boundary (m/s**2)
REAL :: udtnb (nx,nz) ! Time tendency of u at n-boundary (m/s**2)
REAL :: udtsb (nx,nz) ! Time tendency of u at s-boundary (m/s**2)
REAL :: vdteb (ny,nz) ! Time tendency of v at e-boundary (m/s**2)
REAL :: vdtwb (ny,nz) ! Time tendency of v at w-boundary (m/s**2)
REAL :: vdtnb (nx,nz) ! Time tendency of v at n-boundary (m/s**2)
REAL :: vdtsb (nx,nz) ! Time tendency of v at s-boundary (m/s**2)
REAL :: wdteb (ny,nz) ! Time tendency of w at e-boundary (m/s**2)
REAL :: wdtwb (ny,nz) ! Time tendency of w at w-boundary (m/s**2)
REAL :: wdtnb (nx,nz) ! Time tendency of w at n-boundary (m/s**2)
REAL :: wdtsb (nx,nz) ! Time tendency of w at s-boundary (m/s**2)
REAL :: pdteb (ny,nz) ! Time tendency of pprt at e-boundary (Pascal/s)
REAL :: pdtwb (ny,nz) ! Time tendency of pprt at w-boundary (Pascal/s)
REAL :: pdtnb (nx,nz) ! Time tendency of pprt at n-boundary (Pascal/s)
REAL :: pdtsb (nx,nz) ! Time tendency of pprt at s-boundary (Pascal/s)
REAL :: phydro(nx,ny) ! Big time step forcing for computing
! hydrostatic pprt at k=1.
REAL :: sdteb (ny,nz) ! Time tendency of a scalar at e-boundary
REAL :: sdtwb (ny,nz) ! Time tendency of a scalar at w-boundary
REAL :: sdtnb (nx,nz) ! Time tendency of a scalar at n-boundary
REAL :: sdtsb (nx,nz) ! Time tendency of a scalar at s-boundary
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 of staggered grid.
REAL :: mapfct(nx,ny,8) ! Map factors at scalar, u and v points
REAL :: j1 (nx,ny,nz) ! Coordinate transformation Jacobian -d(zp)/d(x)
REAL :: j2 (nx,ny,nz) ! Coordinate transformation Jacobian -d(zp)/d(y)
REAL :: j3 (nx,ny,nz) ! Coordinate transformation Jacobian 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) ! Coordinate transformation Jacobian d(zp)/d(z)
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 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 :: 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 :: ptsfc (nx,ny) ! Potential temperature at the ground level (K)
REAL :: qvsfc (nx,ny,0:nstyps) ! Effective qv at the surface (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 :: 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(IN) :: 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 :: raing (nx,ny) ! Gridscale rainfall (mm)
REAL :: rainc (nx,ny) ! Cumulus rainfall (mm)
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
!
!-----------------------------------------------------------------------
!
! Arrays for radiation
!
!-----------------------------------------------------------------------
!
INCLUDE 'radcst.inc' ! Radiation constants and parameters
INTEGER :: rbufsz ! Radiation buffer size
REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s)
REAL :: radsw (nx,ny) ! Solar radiation reaching the surface
REAL :: rnflx (nx,ny) ! Net absorbed radiation flux at surface
REAL :: rad2d(nx,ny,nrad2d)
! Buffur array to carry the variables calculated or used in
! radiation calculation. The last index defines the variables:
! 1 = nrsirbm, Solar IR surface albedo for beam
! 2 = nrsirdf, Solar IR surface albedo for diffuse
! 3 = nrsuvbm, Solar UV surface albedo for beam
! 4 = nrsuvdf, Solar UV surface albedo for diffuse
! 5 = ncosz, Cosine of zenith
! 6 = ncosss, Cosine of angle between sun light and slope
! 7 = nfdirir, all-sky direct downward IR flux (0.7-10 micron)
! at the surface
! 8 = nfdifir, all-sky diffuse downward IR flux
! at the surface
! 9 = nfdirpar, all-sky direct downward and par flux
! (0.4-0.7 micron) at the surface
! 10 = nfdifpar, all-sky diffuse downward and par flux
! at thediffuse downward and par flux
! at the surface
REAL :: radbuf(rbufsz) ! Buffer array storing temporary working
! arrays for radiation computing
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx(nx,ny) ! EXBC relaxation coefficients
!
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 :: 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
!
!-----------------------------------------------------------------------
!
! Temporary arrays
!
!-----------------------------------------------------------------------
!
REAL :: cdm (nx,ny) ! Drag coefficient for surface momentum flux
REAL :: cdh (nx,ny) ! Drag coefficient for surface heat flux
REAL :: cdq (nx,ny) ! Drag coefficient for surface moisture flux
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' ! Global constants that control model execution
INCLUDE 'bndry.inc' ! Boundary condition control parameters
INCLUDE 'phycst.inc' ! Physical constants
INCLUDE 'exbc.inc' ! EXBC constants and parameters
INCLUDE 'nudging.inc' ! variables for nudging assimilation
INCLUDE 'mp.inc' ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
! Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: nocalls
SAVE nocalls
DATA nocalls / 0 /
INTEGER :: i,j,k
INTEGER :: ireturn ! Return status indicator
REAL :: dtbig1
REAL :: pisum
INTEGER :: tstrtlvl
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
mgrid = mptr
!
!-----------------------------------------------------------------------
!
! Calculate the radiation forcing and surface fluxes if desired.
!
!-----------------------------------------------------------------------
!
CALL set_acct
(rad_acct)
IF ( radopt /= 0 .AND. MOD(nstep-1, nradstp) == 0 ) THEN
IF (myproc == 0) &
WRITE(6,'(a,i8,a,f10.2,a)') &
' Compute radiation at time step,', nstep, &
', model time=',curtim,' (s)'
CALL radiation(nx,ny,nz,rbufsz, &
ptprt(1,1,1,1),pprt(1,1,1,1), &
qv(1,1,1,1),qc(1,1,1,1),qr(1,1,1,1), &
qi(1,1,1,1),qs(1,1,1,1),qh(1,1,1,1), &
ptbar,pbar,ppi,rhostr, &
x,y,z,zp, j3inv, &
soiltyp(1,1,1),tsfc(1,1,0),wetsfc(1,1,0), &
snowdpth,radfrc,radsw,rnflx, &
rad2d(1,1,nrsirbm), rad2d(1,1,nrsirdf), &
rad2d(1,1,nrsuvbm), rad2d(1,1,nrsuvdf), &
rad2d(1,1,ncosz), rad2d(1,1,ncosss), &
rad2d(1,1,nfdirir), rad2d(1,1,nfdifir), &
rad2d(1,1,nfdirpar),rad2d(1,1,nfdifpar), &
tem1, tem2, tem3, tem4, tem5, &
tem6, tem7, tem8, tem9, tem10, &
tem11,tem12,tem13,tem14,tem15,tem16, &
radbuf)
IF( lvldbg >= 2 ) THEN
CALL checkshx(radfrc, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'radfrc', tem1)
CALL checkshy(radfrc, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'radfrc', tem1)
CALL checkshx(radsw,nx,ny,1,1,nx-1,1,ny-1,1,1, &
'radsw', tem1)
CALL checkshy(radsw,nx,ny,1,1,nx-1,1,ny-1,1,1, &
'radsw', tem1)
CALL checkshx(rnflx,nx,ny,1,1,nx-1,1,ny-1,1,1, &
'rnflx', tem1)
CALL checkshy(rnflx,nx,ny,1,1,nx-1,1,ny-1,1,1, &
'rnflx', tem1)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate the surface fluxes of momentum, heat, and water
! vapor for use in the turbulent mixing subroutine
!
!-----------------------------------------------------------------------
!
CALL set_acct(misc_acct)
IF(sadvopt /= 4) THEN ! Leapfrop scheme
tstrtlvl = tpast
ELSE ! Forward-based scheme
tstrtlvl = tpresent
END IF
CALL sfcphysics(nx,ny,nz,nstyps, &
u(1,1,1,tstrtlvl),v(1,1,1,tstrtlvl), &
w(1,1,1,tstrtlvl),ptprt(1,1,1,tstrtlvl), &
pprt(1,1,1,tstrtlvl),qv(1,1,1,tstrtlvl), &
qr(1,1,1,tstrtlvl), &
rhostr,ptbar,pbar,qvbar, &
x,y,z, zp, j1,j2,j3,j3inv, prcrate(1,1,1), &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
radsw,rnflx, &
cdm,cdh,cdq,usflx,vsflx,ptsflx,qvsflx,pbldpth )
!call test_dump(usflx,"XXXAsfcflx_usflx",nx,ny,1,1,1)
!call test_dump(vsflx,"XXXAsfcflx_vsflx",nx,ny,1,2,1)
!call test_dump(ptsflx,"XXXAsfcflx_ptsflx",nx,ny,1,0,1)
!call test_dump(qvsflx,"XXXAsfcflx_qvsflx",nx,ny,1,0,1)
!call test_dump(pbldpth,"XXXAsfcflx_pbldpth",nx,ny,3,0,1)
!
CALL set_acct(bc_acct)
IF( lbcopt == 2 .AND. mgrid == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! Read in variables from an external boundary file, and
! calculate the linear time tendencies of the external data.
!
!-----------------------------------------------------------------------
!
CALL extbdt(nx,ny,nz, ptbar,pbar, ireturn, &
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), &
tem1,tem2,tem3,tem4,tem5,tem6, &
tem7,tem8,tem9,tem10,tem11)
IF( ireturn == 0 ) THEN
GO TO 112
ELSE 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.'
GO TO 111
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.'
GO TO 111
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.'
GO TO 111
ELSE
WRITE (6,'(a/a)') &
'Other errors in getting the external boundary data. Dump the', &
'history file and restart file and then STOP the model.'
GO TO 111
END IF
111 CONTINUE
CALL set_acct(output_acct)
!
!-----------------------------------------------------------------------
!
! Write out external data arrays and stop the program if data read
! failed.
!
!-----------------------------------------------------------------------
!
IF (mp_opt > 0) THEN
CALL arpsstop("arpstop called from CORDINTG ",1)
END IF
CALL abortdmp(mptr,nx,ny,nz,nstyps, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
ubar,vbar,ptbar,pbar,rhostr,qvbar,kmh,kmv, &
x,y,z,zp,zp(1,1,2),mapfct, j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3, tem4)
! 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
!EMK BMJ
CALL rstout(nx,ny,nz,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,zp(1,1,2),mapfct, &
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)
!EMK END
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
CALL arpsstop("arpsstop called from CORDINTG run unstable, aborted",1)
112 CONTINUE
END IF
!
!-----------------------------------------------------------------------
!
! Adjust potential temperature and Qv fields with grid condensation
! or warm rain microphysics
!
!-----------------------------------------------------------------------
!
CALL set_acct(qpfgrd_acct)
IF ( (moist /= 0) .AND. (cnvctopt == 1) .AND. &
(MOD(curtim+0.001,qpfgfrq) < (0.5*dtbig)) ) THEN
CALL qpfgrid(nx,ny,nz,prcrate(1,1,2), &
pprt(1,1,1,2),ptprt(1,1,1,2),qv(1,1,1,2),pbar,ptbar, &
rhostr,zp,j3,j3inv,raing, tem1(1,1,1),tem1(1,1,2), &
tem1(1,1,3),tem1(1,1,4))
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL acct_interrupt(bc_acct)
CALL exbcpt( nx,ny,nz, curtim, ptprt(1,1,1,2), &
exbcbuf(npt0exb),exbcbuf(nptdtexb) )
CALL exbcq( nx,ny,nz, 1,curtim, qv(1,1,1,2), &
exbcbuf(nqv0exb),exbcbuf(nqvdtexb) )
CALL acct_stop_inter
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Integrate the conservation equations of momentum, heat, pressure
! and water variables forward by one timestep
!
!-----------------------------------------------------------------------
!
CALL set_acct(tinteg_acct)
!EMK BMJ
CALL tinteg(mptr, frstep, nx,ny,nz,exbcbufsz, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth, &
udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
phydro,sdteb,sdtwb,sdtnb,sdtsb, &
ubar,vbar,ptbar,pbar,ptbari,pbari, &
rhostr,rhostri,qvbar,ppi,csndsq, &
x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2,sinlat, &
ptsfc,qvsfc(1,1,0),prcrate, radfrc, &
usflx,vsflx,ptsflx,qvsflx,kmh,kmv,rprntl, &
ptcumsrc,qcumsrc,raing,rainc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
exbcbuf, bcrlx, &
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 END
!call test_dump (pprt(1,1,1,2),"XXXAtinteg_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,3),"XXXAtinteg_pprt3",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpast),"XXXAtinteg_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAtinteg_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAtinteg_qr3",nx,ny,nz,0,1)
CALL set_acct(misc_acct)
IF( pdetrnd == 1 .AND. mptr == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! Detrend pressure. This should be used only when unwanted trend
! develops in the pressure field. The precedure forces the
! domain averaged Exner function to zero.
!
! Pressure detrending is sometimes necessary when open lateral
! boundary condition is used. Pressure detrending is applied
! to the base grid (as opposed to nested grids) only.
!
!-----------------------------------------------------------------------
!
pisum = 0.0
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
pisum = pisum + ppi(i,j,k)/(pbar(i,j,k)+pprt(i,j,k,3)) &
* pprt(i,j,k,3)
END DO
END DO
END DO
IF (mp_opt > 0) THEN
CALL mptotal(pisum)
pisum = pisum/(nproc_x*nproc_y)
END IF
pisum = pisum/((nx-3)*(ny-3)*(nz-3))
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
pprt(i,j,k,3) = pprt(i,j,k,3) &
- pisum*(pbar(i,j,k)+pprt(i,j,k,3))/ppi(i,j,k)
END DO
END DO
END DO
!call test_dump (pprt(1,1,1,3),"XXXApdtrnmd_pprt3",nx,ny,nz,0,1)
END IF
!
!-----------------------------------------------------------------------
!
! Call the microphysics routines which adjust the thermodynamic
! and water variables.
!
! Upon exiting this subroutine, supersaturation is removed.
!
!-----------------------------------------------------------------------
!
IF( moist /= 0 ) THEN
IF( mphyopt == 1 .OR. (mphyopt == 0.AND.moist == 1) ) THEN
IF ( cnvctopt /= 1 ) THEN
!-----------------------------------------------------------------------
!
! Kessler warm rain microphysics or saturation adjustment only.
!
!-----------------------------------------------------------------------
IF(frstep == 1) THEN
dtbig1 = dtbig/2
ELSE
dtbig1 = dtbig
END IF
CALL set_acct(warmrain_acct)
!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_qr3",nx,ny,nz,0,1)
CALL microph(nx,ny,nz,dtbig1, &
ptprt,pprt,qv,qc,qr,qi,qs,qh, &
raing,prcrate(1,1,4), &
rhostr,pbar,ptbar,ppi,j3,j3inv, &
tem1,tem2,tem3,tem4,tem5,tem6)
!call test_dump (qr(1,1,1,tpast),"XXXAmicroph_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAmicroph_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAmicroph_qr3",nx,ny,nz,0,1)
END IF
!
!-----------------------------------------------------------------------
!
! The microphysics may change the boundary values of both the
! potential temperature and specific humidity. For the choice of
! external boundary conditions, we must adjust them to the external
! BC.
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL acct_interrupt(bc_acct)
CALL exbcpt(nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture), &
exbcbuf(npt0exb), exbcbuf(nptdtexb))
CALL exbcq( nx,ny,nz, 1,curtim+dtbig, qv(1,1,1,tfuture), &
exbcbuf(nqv0exb),exbcbuf(nqvdtexb))
CALL exbcq( nx,ny,nz, 2,curtim+dtbig, qc(1,1,1,tfuture), &
exbcbuf(nqc0exb),exbcbuf(nqcdtexb) )
CALL exbcq( nx,ny,nz, 3,curtim+dtbig, qr(1,1,1,tfuture), &
exbcbuf(nqr0exb),exbcbuf(nqrdtexb) )
CALL acct_stop_inter
END IF
ELSE IF( mphyopt == 2 ) THEN ! Ice microphysics.
IF(frstep == 1) THEN
dtbig1 = dtbig/2
ELSE
dtbig1 = dtbig
END IF
!
!-----------------------------------------------------------------------
!
! Lin ice microphysics scheme with modifications by W.G. Tao.
! (Lin, Y.-L., R. D. Farley, and H. D. Orville, 1983:
! Bulk parameterization of the snow field in a cloud model.
! J. Climate Appl. Meteor., 22, 1065-1092.
! Tao, W.-K., and J. Simpson, 1993: Goddard cumulus ensemble model.
! Part I: Model description. Terres. Atmos. Ocean Sci., 4, 35-72.)
!
!-----------------------------------------------------------------------
!
CALL set_acct(linice_acct)
!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_ice_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_ice_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_ice_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXBmicroph_ice_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXBmicroph_ice_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXBmicroph_ice_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
CALL microph_ice(nx,ny,nz,dtbig1, &
ptprt,pprt,qv,qc,qr,qi,qs,qh, &
raing,prcrate(1,1,4), &
rhostr,pbar,ptbar,qvbar,ppi,j3,j3inv, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11,tem12,tem13,tem14,tem15)
!call test_dump (qr(1,1,1,tpast),"XXXAmicroph_ice_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAmicroph_ice_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAmicroph_ice_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXAmicroph_ice_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAmicroph_ice_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAmicroph_ice_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!
!-----------------------------------------------------------------------
!
! The microphysics may change the boundary values of both the
! potential temperature and specific humidity. For the choice of
! external boundary conditions, we must adjust them to the external
! BC.
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
CALL acct_interrupt(bc_acct)
CALL exbcpt( nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture), &
exbcbuf(npt0exb),exbcbuf(nptdtexb) )
CALL exbcq( nx,ny,nz, 1,curtim+dtbig, qv(1,1,1,tfuture), &
exbcbuf(nqv0exb),exbcbuf(nqvdtexb) )
CALL exbcq( nx,ny,nz, 2,curtim+dtbig, qc(1,1,1,tfuture), &
exbcbuf(nqc0exb),exbcbuf(nqcdtexb) )
CALL exbcq( nx,ny,nz, 3,curtim+dtbig, qr(1,1,1,tfuture), &
exbcbuf(nqr0exb),exbcbuf(nqrdtexb) )
CALL exbcq( nx,ny,nz, 4,curtim+dtbig, qi(1,1,1,tfuture), &
exbcbuf(nqi0exb),exbcbuf(nqidtexb) )
CALL exbcq( nx,ny,nz, 5,curtim+dtbig, qs(1,1,1,tfuture), &
exbcbuf(nqs0exb),exbcbuf(nqsdtexb) )
CALL exbcq( nx,ny,nz, 6,curtim+dtbig, qh(1,1,1,tfuture), &
exbcbuf(nqh0exb),exbcbuf(nqhdtexb) )
CALL acct_stop_inter
END IF
!call test_dump (qr(1,1,1,tpast),"XXXAexbcq_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAexbcq_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAexbcq_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXAexbcq_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAexbcq_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAexbcq_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
ELSE IF ( mphyopt == 3 ) THEN
IF(frstep == 1) THEN
dtbig1 = dtbig/2
ELSE
dtbig1 = dtbig
END IF
!
!-----------------------------------------------------------------------
!
! Schultz ice microphysics scheme (Schultz, P., 1995: An explicit
! cloud physics parameterization for operational numerical
! weather prediction. Mon. Wea. Rev., 123, 3331-3343).
!
!-----------------------------------------------------------------------
!
CALL set_acct(nemice_acct)
!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_nem_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_nem_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_nem_qr3",nx,ny,nz,0,1)
CALL microph_nem(nx,ny,nz,dtbig1,ptprt,pprt,qv,qc, &
qr,qi,qs,qh,raing,prcrate(1,1,4),rhostr, &
pbar,ptbar,ppi,j3,j3inv, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11,tem12,tem13, &
tem14,tem15,tem16)
!call test_dump2(qr(1,1,1,tpast),"XXXAmicroph_nem_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAmicroph_nem_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAmicroph_nem_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!
!-----------------------------------------------------------------------
!
! The microphysics may change the boundary values of both the
! potential temperature and specific humidity. For the choice of
! external boundary conditions, we must adjust them to the external
! BC.
!
!-----------------------------------------------------------------------
!
IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_nembc_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_nembc_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_nembc_qr3",nx,ny,nz,0,1)
CALL acct_interrupt(bc_acct)
CALL exbcpt( nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture), &
exbcbuf(npt0exb),exbcbuf(nptdtexb) )
CALL exbcq( nx,ny,nz, 1,curtim+dtbig, qv(1,1,1,tfuture), &
exbcbuf(nqv0exb),exbcbuf(nqvdtexb) )
CALL exbcq( nx,ny,nz, 2,curtim+dtbig, qc(1,1,1,tfuture), &
exbcbuf(nqc0exb),exbcbuf(nqcdtexb) )
CALL exbcq( nx,ny,nz, 3,curtim+dtbig, qr(1,1,1,tfuture), &
exbcbuf(nqr0exb),exbcbuf(nqrdtexb) )
CALL exbcq( nx,ny,nz, 4,curtim+dtbig, qi(1,1,1,tfuture), &
exbcbuf(nqi0exb),exbcbuf(nqidtexb) )
CALL exbcq( nx,ny,nz, 5,curtim+dtbig, qs(1,1,1,tfuture), &
exbcbuf(nqs0exb),exbcbuf(nqsdtexb) )
CALL exbcq( nx,ny,nz, 6,curtim+dtbig, qh(1,1,1,tfuture), &
exbcbuf(nqh0exb),exbcbuf(nqhdtexb) )
CALL acct_stop_inter
END IF
!call test_dump (qr(1,1,1,tpast),"XXXAmicroph_nembc_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAmicroph_nembc_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAmicroph_nembc_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXAmicroph_nembc_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAmicroph_nembc_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAmicroph_nembc_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
ELSE ! Invalid option
WRITE(6,'(1x,a,/1x,a,i3)') &
'Invalid option for microphysics parameterization.', &
'Job stopped here. MPHYOPT was ',mphyopt
CALL arpsstop("arpsstop called from CORDINTG improper mphyopt",1)
END IF
DO j=1,ny-1
DO i=1,nx-1
prcrate(i,j,1)=prcrate(i,j,2)+prcrate(i,j,3)+prcrate(i,j,4)
END DO
END DO
END IF ! moist.ne.0
!
!-----------------------------------------------------------------------
!
! Apply nudging
!
!-----------------------------------------------------------------------
!
CALL set_acct(misc_acct)
IF(nudgopt > 0 .AND. curtim < ndstop .AND. curtim >= ndstart .AND. &
MOD(nstep, nudgstp) == 0 ) THEN
WRITE(6,'(a,f10.1)') &
' Nudging forecast at ',(curtim+dtbig)
CALL nudgeall(nx,ny,nz, nxndg,nyndg,nzndg, &
u(1,1,1,tpresent),v(1,1,1,tpresent), &
w(1,1,1,tpresent),pprt(1,1,1,tpresent),ptprt(1,1,1,tpresent), &
qv(1,1,1,tpresent),qc(1,1,1,tpresent),qr(1,1,1,tpresent), &
qi(1,1,1,tpresent),qs(1,1,1,tpresent),qh(1,1,1,tpresent), &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr)
CALL nudgeall(nx,ny,nz, nxndg,nyndg,nzndg, &
u(1,1,1,tfuture),v(1,1,1,tfuture), &
w(1,1,1,tfuture),pprt(1,1,1,tfuture),ptprt(1,1,1,tfuture), &
qv(1,1,1,tfuture),qc(1,1,1,tfuture),qr(1,1,1,tfuture), &
qi(1,1,1,tfuture),qs(1,1,1,tfuture),qh(1,1,1,tfuture), &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr)
END IF
!
!-----------------------------------------------------------------------
!
! Apply the Asselin time filter to all variables
!
!-----------------------------------------------------------------------
!
IF( flteps /= 0 .AND. frstep /= 1 ) THEN
CALL tfilt(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke)
!call test_dump2(qr(1,1,1,tpast),"XXXAtfilt_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAtfilt_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAtfilt_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
END IF
!
!-----------------------------------------------------------------------
!
!
! Call the Klemp-Lilly (1978)/Durran (1983) version of the radiation
!
!
!-----------------------------------------------------------------------
!
IF(rbcopt == 2.OR.rbcopt == 3.OR.rbcopt == 4)THEN
CALL set_acct(bc_acct)
IF(frstep == 1) THEN
dtbig1 = dtbig/2
ELSE
dtbig1 = dtbig
END IF
CALL bdtu(nx,ny,nz,dtbig1, u,ubar,udteb,udtwb)
CALL bdtv(nx,ny,nz,dtbig1, v,vbar,vdtnb,vdtsb)
END IF
!
!-----------------------------------------------------------------------
!
! To prepare for the next timestep integration, swap time levels
! for all time-dependent variables.
!
! The swap is as follows:
!
! Fields at time tpresent are assigned to tpast,
! Fields at time tfuture are assigned to tpresent,
! Fields at time tfuture are not changed.
!
!-----------------------------------------------------------------------
!
CALL set_acct(misc_acct)
!call test_dump (u(1,1,1,tpast),"XXXBflip_u1",nx,ny,nz,0,1)
!call test_dump (u(1,1,1,tpresent),"XXXBflip_u2",nx,ny,nz,0,1)
!call test_dump (u(1,1,1,tfuture),"XXXBflip_u3",nx,ny,nz,0,1)
!call test_dump (v(1,1,1,tpast),"XXXBflip_v1",nx,ny,nz,0,1)
!call test_dump (v(1,1,1,tpresent),"XXXBflip_v2",nx,ny,nz,0,1)
!call test_dump (v(1,1,1,tfuture),"XXXBflip_v3",nx,ny,nz,0,1)
!call test_dump (w(1,1,1,tpast),"XXXBflip_w1",nx,ny,nz,0,1)
!call test_dump (w(1,1,1,tpresent),"XXXBflip_w2",nx,ny,nz,0,1)
!call test_dump (w(1,1,1,tfuture),"XXXBflip_w3",nx,ny,nz,0,1)
!call test_dump (ptprt(1,1,1,tpast),"XXXBflip_ptprt1",nx,ny,nz,0,1)
!call test_dump (ptprt(1,1,1,tpresent),"XXXBflip_ptprt2",nx,ny,nz,0,1)
!call test_dump (ptprt(1,1,1,tfuture),"XXXBflip_ptprt3",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,tpast),"XXXBflip_pprt1",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,tpresent),"XXXBflip_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,tfuture),"XXXBflip_pprt3",nx,ny,nz,0,1)
!call test_dump (qv(1,1,1,tpast),"XXXBflip_qv1",nx,ny,nz,0,1)
!call test_dump (qv(1,1,1,tpresent),"XXXBflip_qv2",nx,ny,nz,0,1)
!call test_dump (qv(1,1,1,tfuture),"XXXBflip_qv3",nx,ny,nz,0,1)
!call test_dump (qc(1,1,1,tpast),"XXXBflip_qc1",nx,ny,nz,0,1)
!call test_dump (qc(1,1,1,tpresent),"XXXBflip_qc2",nx,ny,nz,0,1)
!call test_dump (qc(1,1,1,tfuture),"XXXBflip_qc3",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpast),"XXXBflip_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBflip_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBflip_qr3",nx,ny,nz,0,1)
!call test_dump (qi(1,1,1,tpast),"XXXBflip_qi1",nx,ny,nz,0,1)
!call test_dump (qi(1,1,1,tpresent),"XXXBflip_qi2",nx,ny,nz,0,1)
!call test_dump (qi(1,1,1,tfuture),"XXXBflip_qi3",nx,ny,nz,0,1)
!call test_dump (qs(1,1,1,tpast),"XXXBflip_qs1",nx,ny,nz,0,1)
!call test_dump (qs(1,1,1,tpresent),"XXXBflip_qs2",nx,ny,nz,0,1)
!call test_dump (qs(1,1,1,tfuture),"XXXBflip_qs3",nx,ny,nz,0,1)
!call test_dump (qh(1,1,1,tpast),"XXXBflip_qh1",nx,ny,nz,0,1)
!call test_dump (qh(1,1,1,tpresent),"XXXBflip_qh2",nx,ny,nz,0,1)
!call test_dump (qh(1,1,1,tfuture),"XXXBflip_qh3",nx,ny,nz,0,1)
!call test_dump (tke(1,1,1,tpast),"XXXBflip_tke1",nx,ny,nz,0,0)
!call test_dump (tke(1,1,1,tpresent),"XXXBflip_tke2",nx,ny,nz,0,0)
!call test_dump (tke(1,1,1,tfuture),"XXXBflip_tke3",nx,ny,nz,0,0)
CALL tflip(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke)
!
!-----------------------------------------------------------------------
!
! Print max./min. for the nested grid case.
!
!-----------------------------------------------------------------------
! IF( mgrid.eq.1 .and. nestgrd.eq.1 ) THEN
curtim = curtim + dtbig
! CALL maxmin(mptr,nx,ny,nz,tlevel,rhobar,
! : u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv,
! : x,y,z,zp,mapfct,
! : tsfc,tsoil,wetsfc,wetdp,wetcanp,
! : tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
! Check the stability of the integration. If the model solution
! appears to be exceeding the linear stability condition, then
! perform a data dump for post-mortem.
!
!-----------------------------------------------------------------------
!
CALL chkstab(mgrid,nx,ny,nz,nstyps, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
ubar,vbar,ptbar,pbar,rhostr,qvbar,kmh,kmv, &
x,y,z,zp,zp(1,1,2),mapfct,j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3, tem4)
curtim = curtim - dtbig
! ENDIF
RETURN
END SUBROUTINE cordintg
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TINTEG ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!EMK BMJ
SUBROUTINE tinteg(mptr, frstep, nx,ny,nz,exbcbufsz, & 1,66
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth, &
udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
phydro,sdteb,sdtwb,sdtnb,sdtsb, &
ubar,vbar,ptbar,pbar,ptbari,pbari, &
rhostr,rhostri,qvbar,ppi,csndsq, &
x,y,z,zp, mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2,sinlat, &
ptsfc,qvsfc,prcrate, radfrc, &
usflx,vsflx,ptsflx,qvsflx,kmh,kmv,rprntl, &
ptcumsrc,qcumsrc,raing,rainc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
exbcbuf,bcrlx, &
rhofct,defsq,lenscl, &
uforce,vforce,wforce,pforce,ptforce, &
tem1,tem2,tem3,tem4,tem5,tem6, &
tem7,tem8,tem9,tem10,tem11,tem12,tem13,tem14,tem15, &
tem16,tem17,tem18, &
tem1_0,tem2_0,tem3_0)
!EMK END
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Orchestrate the time integration of the dynamics of the basic governing
! equations for a single time step. Note that the time tendencies of the
! non-conservative processes (e.g., microphysics, radiation, etc.)
! are not included in this routine but are applied in their
! corresponding routines.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/21/92.
!
! MODIFICATION HISTORY:
!
! 5/03/92 (K. Droegemeier)
! Added full documentation.
!
! 5/20/92 (M. Xue)
! Reformatted certain documentations.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 5/26/93 (M. Xue)
! Added w into the solvpt argument list.
!
! 9/7/94 (M.Xue)
! Call to ADVP modified.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 8/20/95 (M. Xue)
! Bug fix with the use of ptcumsrc.
!
! 10/31/95 (D. Weber)
! Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
! radiation condition.
!
! 01/23/1996 (Donghai Wang, Ming Xue and Yuhe Liu)
! Added the map projection factor to ARPS governing equations.
!
! 03/11/96 (M. Xue)
! Modified the call to SOLVTKE, so that all three time levels of
! ptprt, pprt etc are passed into the routine.
!
! 03/19/96 (Yuhe Liu)
! Added radiation forcing
!
! 07/10/1997 (Fanyou Kong - CMRP)
! Fixed a bug in 'solvtke' with corrected temporary arga
!
! 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.
!
! 10/21/97 (Donghai Wang)
! Using total density (rho) in the calculation of the pressure
! gradient force terms, and added the second order terms
! in the linerized buoyancy terms.
!
! 11/05/97 (D. Weber)
! Added phydro array for use in the bottom boundary condition for
! perturbation 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/18/98 (D. Weber)
! Added precomputed averages of j3 in the x, y, and z directions
! to improve the code efficiency.
!
! 07/23/2001 (K. Brewster)
! Added mptr to argument list, needed for printing of diagnostic
! noise parameter. Also added mptr to calls to smlstep, solvuw,
! solvwpex, and solvwpim.
!
! 13 March 2002 (Eric Kemp)
! Added arrays for WRF BMJ cumulus scheme.
!
! April 2002 (Fanyou Kong)
! Added cnvctopt=5 option for the new WRF K-F (KF_ETA) scheme
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! frstep Initial time step flag.
! 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
!
! 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)
! wcont Contravariant vertical velocity (m/s)
! computational coordinates (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)
!
! 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)
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! 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)
! cnsd Sound wave speed.
!
! 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)
! 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.
! vwork2 2-D work array for fftopt=2.
! wsave1 Work array for fftopt=2.
! wsave2 Work array for fftopt=2.
!
! sinlat Sin of latitude at each grid point
!
! 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 (kg/(m**2 * s))
!
! 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)
!
! 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
!
! radfrc Radiation forcing (K/s)
!
! OUTPUT:
!
! u x component of velocity at time tfuture (m/s)
! v y component of velocity at time tfuture (m/s)
! w Vertical component of Cartesian velocity at tfuture (m/s)
! ptprt Perturbation potential temperature at time tfuture (K)
! pprt Perturbation pressure at time tfuture (Pascal)
! qv Water vapor specific humidity at time tfuture (kg/kg)
! qc Cloud water mixing ratio at time tfuture (kg/kg)
! qr Rainwater mixing ratio at time tfuture (kg/kg)
! qi Cloud ice mixing ratio at time tfuture (kg/kg)
! qs Snow mixing ratio at time tfuture (kg/kg)
! qh Hail mixing ratio at time tfuture (kg/kg)
!
! 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)
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! tke Turbulent Kinetic Energy ((m/s)**2)
! pbldpth Planetary boundary layer depth (m)
!
! 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
!
! WORK ARRAYS:
!
! defsq Deformation squared (1/s**2)
! lenscl Turbulent mixing length scale (m)
! uforce Acoustically inactive forcing terms in u-Eq. ((kg/(m*s)**2)
! vforce Acoustically inactive forcing terms in v-Eq. ((kg/(m*s)**2)
! wforce Acoustically inactive forcing terms in w-Eq. ((kg/(m*s)**2)
! pforce Acoustically inactive forcing terms in p-eq. (Pascal/s)
! ptforce Gravity wave inactive forcing terms in pt-eq. (K*kg/(m**3*s))
!
! rhofct rho-factor: rhobar/rho
!
! 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)
!
! 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)
!
! ptdteb Time tendency of ptprt field at east boundary (K/s)
! ptdtwb Time tendency of ptprt field at west boundary (K/s)
! ptdtnb Time tendency of ptprt field at north boundary(K/s)
! ptdtsb Time tendency of ptprt field at south boundary(K/s)
!
! sdteb Time tendency of a scalar at e-boundary
! sdtwb Time tendency of a scalar at w-boundary
! sdtnb Time tendency of a scalar at n-boundary
! sdtsb Time tendency of a scalar at s-boundary
!
! 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
!
! tem1_0 Temporary work array.
! tem2_0 Temporary work array.
! tem3_0 Temporary work array.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'timelvls.inc'
INTEGER :: mptr
INTEGER :: nx, ny, nz ! 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 (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: tke (nx,ny,nz,nt) ! Turbulent Kinetic Energy ((m/s)**2)
REAL :: pbldpth(nx,ny,nt) ! Planetary boundary layer depth (m)
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 :: udteb (ny,nz) ! Time tendency of u at e-boundary (m/s**2)
REAL :: udtwb (ny,nz) ! Time tendency of u at w-boundary (m/s**2)
REAL :: udtnb (nx,nz) ! Time tendency of u at n-boundary (m/s**2)
REAL :: udtsb (nx,nz) ! Time tendency of u at s-boundary (m/s**2)
REAL :: vdteb (ny,nz) ! Time tendency of v at e-boundary (m/s**2)
REAL :: vdtwb (ny,nz) ! Time tendency of v at w-boundary (m/s**2)
REAL :: vdtnb (nx,nz) ! Time tendency of v at n-boundary (m/s**2)
REAL :: vdtsb (nx,nz) ! Time tendency of v at s-boundary (m/s**2)
REAL :: wdteb (ny,nz) ! Time tendency of w at e-boundary (m/s**2)
REAL :: wdtwb (ny,nz) ! Time tendency of w at w-boundary (m/s**2)
REAL :: wdtnb (nx,nz) ! Time tendency of w at n-boundary (m/s**2)
REAL :: wdtsb (nx,nz) ! Time tendency of w at s-boundary (m/s**2)
REAL :: pdteb (ny,nz) ! Time tendency of pprt at e-boundary (Pascal/s)
REAL :: pdtwb (ny,nz) ! Time tendency of pprt at w-boundary (Pascal/s)
REAL :: pdtnb (nx,nz) ! Time tendency of pprt at n-boundary (Pascal/s)
REAL :: pdtsb (nx,nz) ! Time tendency of pprt at s-boundary (Pascal/s)
REAL :: phydro(nx,ny) ! Pressure at k=1 computed using pert.
! hydrostatic relation.
REAL :: sdteb (ny,nz) ! Time tendency of a variable at e-boundary
REAL :: sdtwb (ny,nz) ! Time tendency of a variable at w-boundary
REAL :: sdtnb (nx,nz) ! Time tendency of a variable at n-boundary
REAL :: sdtsb (nx,nz) ! Time tendency of a variable at s-boundary
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 of the staggered grid.
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) ! Coordinate transformation Jacobian defined as
! d( zp )/d( z ).
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 :: ptsfc (nx,ny) ! Potential temperature at the ground level (K)
REAL :: qvsfc (nx,ny) ! Effective qv at the surface (kg/kg)
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
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 :: 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
REAL :: defsq (nx,ny,nz) ! Deformation squared (1/s**2)
REAL :: lenscl(nx,ny,nz) ! Turbulent mixing length scale (m)
REAL :: uforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in u-momentum equation (kg/(m*s)**2)
REAL :: vforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in v-momentum equation (kg/(m*s)**2)
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
REAL :: pforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in pressure equation (Pascal/s)
REAL :: ptforce(nx,ny,nz) ! Gravity wave inactive forcing terms
! in pressure equation (Pascal/s)
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(IN) :: 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 :: raing(nx,ny) ! Grid supersaturation rain
REAL :: rainc(nx,ny) ! Cumulus convective rain
REAL :: radfrc(nx,ny,nz) ! Radiation forcing (K/s)
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: bcrlx (nx,ny) ! EXBC relaxation coefficients
REAL :: rhofct(nx,ny,nz) ! rho-factor: rhobar/rho
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 :: 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.
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
REAL :: dtbig1 ! Local value of big time step size
REAL :: dtsml1 ! Local value of small time step size
INTEGER :: frstep ! Flag for the initial time step
INTEGER :: i,j,k
INTEGER :: qflag ! Indicator for the water/ice type
! when calling SOLVQ.
INTEGER :: ntst
REAL :: tst
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! To initialize the three-time level scheme, we assume for
! the first step that the values of all variables at time
! past equal the values at time present. We then perform a
! forward-in-time integration for the first step only and
! then, using the same algorithm, we perform a centered-in-time
! integration for all subsequent steps. The size of the
! timestep is adjusted accordingly, such that the leapfrog step
! size is twice that of the first forward-in-time step.
!
!-----------------------------------------------------------------------
!
IF(frstep == 1) THEN ! frstep=1 on the first time step
! only - indicating a forward-in-
! time integration.
dtbig1 = dtbig/2
dtsml1 = dtsml/2
ELSE
dtbig1 = dtbig
dtsml1 = dtsml
END IF
!
!-----------------------------------------------------------------------
!
! Calculate wcont at time tpresent. wcont will be used in FRCUVW
! and FRCP. rhostr averaged to u, v, w points. They are stored
! in tem1, tem2 and tem3.
!
!-----------------------------------------------------------------------
!
CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
!call test_dump (w(1,1,1,tpresent),"XXX1Bwc_w1",nx,ny,nz,1,0)
CALL wcontra(nx,ny,nz, &
u(1,1,1,tpresent),v(1,1,1,tpresent), &
w(1,1,1,tpresent),mapfct,j1,j2,j3,aj3z, &
rhostr,tem1,tem2,tem3,wcont,tem4,tem5)
!
!-----------------------------------------------------------------------
!
! Compute the acoustically inactive terms in the momentum and
! pressure equations that are held fixed during the small time
! step computations. This includes advection, buoyancy, mixing
! (both physical and computational), and the Coriolis terms.
! These forcing terms are accumulated into arrays for each
! of the momentum equations, e.g., uforce for the u-equation,
! vforce for the v-equation, wforce for the w-equation and
! pforce for the p-equation.
!
! Note: Arrays, pforce and ptforce, are used as work arrays in
! subroutine frcuvw.
!
!-----------------------------------------------------------------------
!
!call test_dump2(rhofct,"XXXBfrcuvw_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
CALL frcuvw(nx,ny,nz,exbcbufsz,dtbig1, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth, &
ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,qvbar, &
usflx,vsflx, x,y,z,zp, mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, sinlat, ptsfc, &
uforce,vforce,wforce,kmh,kmv,rprntl,lenscl,defsq, &
exbcbuf, bcrlx,rhofct,phydro, &
pforce,ptforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11,tem12,tem13,tem14)
!call test_dump (uforce,"XXX2uforce",nx,ny,nz,1,0)
!call test_dump2(rhofct,"XXXAfrcuvw_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
IF (tmixopt == 4) THEN
!
!-----------------------------------------------------------------------
!
! Integrate the TKE equation.
! NOTE: Arrays lenscl and defsq should NOT be changed between
! calls to frcuvw and solvtke.
! Note: ptforce is used as a work array in the following call.
!
!-----------------------------------------------------------------------
!
CALL solvtke(nx,ny,nz,dtbig1,u,v,wcont,ptprt,pprt, &
qv,qc,qr,qi,qs,qh,tke, &
ubar,vbar,ptbar,pbar,ptbari,rhostr,rhostri,qvbar, &
x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,j3inv, &
kmh,kmv,rprntl,lenscl,defsq, &
ptsflx,qvsflx, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem9,tem10,tem11,ptforce, & ! ptforce used as a tem array
tem1_0,tem2_0,tem3_0)
END IF
CALL frcp(nx,ny,nz,exbcbufsz, dtbig1, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptbar,pbar,rhostr,qvbar,mapfct,j1,j2,j3,aj3x,aj3y,aj3z, &
pforce, &
exbcbuf, bcrlx, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
!call test_dump (pprt(1,1,1,2),"XXXAfrcp_pprt2",nx,ny,nz,0,0)
!call test_dump (pprt(1,1,1,3),"XXXAfrcp_pprt3",nx,ny,nz,0,0)
!call test_dump (pforce,"XXXAfrcp_pforce",nx,ny,nz,0,0)
!
!-----------------------------------------------------------------------
!
! Calculate gravity wave or acoustic wave inactive terms in the
! potential temperature equation.
!
! The force terms are stored in ptforce.
!
!-----------------------------------------------------------------------
!
CALL frcpt(nx,ny,nz, exbcbufsz, dtbig1,ptprt,u,v,w,wcont, &
ptbar,rhostr,rhostri,kmh,kmv,rprntl, &
usflx,vsflx,ptsflx,pbldpth, &
x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,j3inv,ptsfc, &
ptforce, &
exbcbuf, bcrlx, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11, &
tem1_0,tem2_0,tem3_0,tem12)
CALL set_acct(tinteg_acct)
IF ( radopt == 2 ) THEN
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
ptforce(i,j,k) = ptforce(i,j,k) &
+ rhostr(i,j,k)*radfrc(i,j,k)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Calculate source and sink terms in temperature (ptprt) and
! moisture (qv) equations that are due to subgrid scale cumulus
! convection.
!
!-----------------------------------------------------------------------
!
!EMK BMJ
! IF (cnvctopt == 1 .OR. cnvctopt == 2 .OR. cnvctopt == 3) THEN
IF (cnvctopt == 1 .OR. cnvctopt == 2 .OR. cnvctopt == 3 .OR. &
cnvctopt == 4 .OR. cnvctopt == 5) THEN
!END EMK
IF (cnvctopt == 3 .OR. cnvctopt == 5) THEN
CALL set_acct(kfcum_acct)
!EMK BMJ
ELSE IF (cnvctopt == 4 ) THEN
CALL set_acct(bmjcum_acct)
!EMK END
ELSE
CALL set_acct(kuocum_acct)
ENDIF
!-----------------------------------------------------------------------
!
! Calculate w0avg, is close to a running mean vertical velocity,
! tst is the number of time steps in 10min for K-F scheme
!
!-----------------------------------------------------------------------
!
IF (cnvctopt == 3 .OR. cnvctopt == 5) THEN
ntst=nint( 600.0/dtbig )
tst=FLOAT(ntst)
IF( (curtim-tstart) <= 300.0 .AND. initopt /= 2 ) THEN
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
w0avg(i,j,k)= (w(i,j,k,1)+w(i,j,k,2))*0.5
END DO
END DO
END DO
ELSE IF ( (curtim-tstart) > 300.0 .OR. initopt == 2) THEN
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
w0avg(i,j,k)=(w0avg(i,j,k)*(tst-1.)+w(i,j,k,2))/tst
END DO
END DO
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate vertical velocity for KUO scheme, or running-average
! vertical velocity (m/s) for Kain Fritsch scheme
!
!-----------------------------------------------------------------------
!
! IF( cnvctopt == 3 .AND. MOD(curtim+0.001,confrq) <= (0.5*dtbig) &
! .AND. ((curtim-tstart) > 300.0 &
IF( (cnvctopt == 3 .OR. cnvctopt == 5) .AND. MOD(curtim+0.001,confrq) &
<= (0.5*dtbig) .AND. ((curtim-tstart) > 300.0 &
.OR. initopt == 2)) THEN
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
! tem1(i,j,k) = (w(i,j,k,1)+w(i,j,k,2))*0.5
tem1(i,j,k) = w0avg(i,j,k)
END DO
END DO
END DO
ELSE IF( MOD(curtim+0.001,confrq) <= (0.5*dtbig).OR. &
nstep == 1 ) THEN
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = wcont(i,j,k)*j3(i,j,k)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Call cumulus parameterization schemes
! Make sure to reset kfraincv, ptcumsrc and qcumsrc to 0 once nca<=0
!
!-----------------------------------------------------------------------
!
IF( (moist /= 0) .AND. ((curtim-tstart) > 300.0 .OR. initopt == 2) .AND. &
(MOD(curtim+0.001,confrq) <= (0.5*dtbig))) THEN
DO j=1,ny-1
DO i=1,nx-1
IF (nca(i,j) <= 0) THEN
kfraincv(i,j) = 0.0
prcrate(i,j,3) = 0.0
END IF
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
IF (nca(i,j) <= 0) THEN
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 IF
END DO
END DO
END DO
!EMK BMJ
CALL qpfcums(nx,ny,nz,prcrate(1,1,3),qvsflx, &
u(1,1,1,2),v(1,1,1,2),tem1, &
pprt(1,1,1,2),ptprt(1,1,1,2),qv(1,1,1,2), &
pbar,ptbar,qvbar,rhostr,zp,j3, &
ptcumsrc,qcumsrc,rainc,nca,kfraincv, &
cldefi,xland,bmjraincv, &
tem2,tem3,tem4,tem5,tem6, &
tem7,tem8,tem9,tem10,tem11)
!EMK END
END IF
!
! Accumulate rainc and update nca. Note: raincv is in cm.
!
DO j=1,ny-1
DO i=1,nx-1
rainc(i,j) = rainc(i,j) + 10.0*kfraincv(i,j) + 10.0*bmjraincv(i,j)
IF ( nca(i,j) >= 1 ) nca(i,j) = nca(i,j) - 1
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
ptforce(i,j,k) = ptforce(i,j,k) &
+ rhostr(i,j,k)*ptcumsrc(i,j,k)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! End of cumulus parameterization
!
!-----------------------------------------------------------------------
!
CALL set_acct(misc_acct)
IF( lvldbg >= 2 ) THEN
CALL checkuhx(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'ufrcx', tem1)
CALL checkuhy(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2, &
'ufrcy', tem1)
CALL checkvhx(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vfrcx', tem1)
CALL checkvhy(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2, &
'vfrcy', tem1)
CALL checkwhx(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wfrcx', tem1)
CALL checkwhy(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wfrcy', tem1)
CALL checkshx(pforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'pfrcx', tem1)
CALL checkshy(pforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'pfrcy', tem1)
CALL checkshx(ptforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'ptfrcx', tem1)
CALL checkshy(ptforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2, &
'ptfrcy', tem1)
END IF
!
!-----------------------------------------------------------------------
!
! Integrate the momentum and pressure equations (i.e., the
! acoustically-active equations) in time using a mode-splitting
! approach. The momentum components are advanced in time using
! a forward scheme (relative to the pressure gradient force terms)
! and the pressure is advanced using a backward scheme (relative to
! the divergence term, which is evaluated using the newly updated
! velocities hence the backward scheme). These small time steps
! bring the momentum and pressure from time tpast to tfuture.
! During these steps, the acoustically-inactive forcing terms (e.g.,
! uforce, pforce, etc.) are held fixed at time tpresent (i.e. they
! are evaluated only once, at time tpresent, for each large time
! step integration), hence the time integration is leapfrog-
! in-time relative to these forcing terms.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! When this solver is called for a nested grid (not the base grid),
! the information stored in u, v, w and pprt is used to set the
! time tendency of these variables at the lateral boundaries.
!
!-----------------------------------------------------------------------
!
IF( mgrid /= 1 .AND. nestgrd == 1 ) THEN
CALL nestbdt(nx,ny,nz,u, 1,nx,1,ny-1,1,nz-1,dtbig, &
udteb,udtwb,udtnb,udtsb)
CALL nestbdt(nx,ny,nz,v, 1,nx-1,1,ny,1,nz-1,dtbig, &
vdteb,vdtwb,vdtnb,vdtsb)
CALL nestbdt(nx,ny,nz,w, 1,nx-1,1,ny-1,1,nz,dtbig, &
wdteb,wdtwb,wdtnb,wdtsb)
CALL nestbdt(nx,ny,nz,pprt, 1,nx-1,1,ny-1,1,nz-1,dtbig, &
pdteb,pdtwb,pdtnb,pdtsb)
CALL nestbdt(nx,ny,nz,ptprt,1,nx-1,1,ny-1,1,nz-1,dtbig, &
sdteb,sdtwb,sdtnb,sdtsb)
!
!-----------------------------------------------------------------------
!
! For the first forward time step, make sure that the fields at
! tpast equal those at tpresent. This can only be done after
! the time tendencies on the boundaries are calculated and stored
! in the time tendency arrays.
!
!-----------------------------------------------------------------------
!
IF( frstep == 1 ) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
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)
pprt (i,j,k,tpast)=pprt (i,j,k,tpresent)
ptprt(i,j,k,tpast)=ptprt(i,j,k,tpresent)
END DO
END DO
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Note here defsq is used as a work array by ACOUST.
!
!-----------------------------------------------------------------------
!
CALL set_acct(smlstp_acct)
CALL smlstep(mptr, nx,ny,nz, exbcbufsz, dtbig1,dtsml1, &
u,v,w,wcont,pprt,ptprt, &
udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
phydro,sdteb,sdtwb,sdtnb,sdtsb, &
ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,rhostri, &
csndsq, &
uforce,vforce,wforce,pforce,ptforce, &
x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
exbcbuf,rhofct, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,defsq, &
tem10,tem11,tem12,tem13,tem14,tem15,tem16,tem17, &
tem18,tem1_0,tem2_0,tem3_0)
!call test_dump (pprt(1,1,1,2),"XXXAsmlstep_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,3),"XXXAsmlstep_pprt3",nx,ny,nz,0,1)
CALL set_acct(misc_acct)
IF( lvldbg >= 1 ) THEN
CALL checkuhx( nx,ny,nz,1,nx,1,ny-1,1,nz-1, 'uxbig', tem1)
CALL checkuhy( nx,ny,nz,1,nx,1,ny-1,1,nz-1, 'uybig', tem1)
CALL checkvhx( nx,ny,nz,1,nx-1,1,ny,1,nz-1, 'vxbig', tem1)
CALL checkvhy( nx,ny,nz,1,nx-1,1,ny,1,nz-1, 'vybig', tem1)
CALL checkwhx( nx,ny,nz,1,nx-1,1,ny-1,1,nz, 'wxbig', tem1)
CALL checkwhy( nx,ny,nz,1,nx-1,1,ny-1,1,nz, 'wybig', tem1)
CALL checkshx(pprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pxbig', tem1)
CALL checkshy(pprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pybig', tem1)
CALL checkshx(ptprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptxbig', tem1)
CALL checkshy(ptprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptybig', tem1)
END IF
!
!-----------------------------------------------------------------------
!
! Since wcont was reset in SMLSTP. Now need to re-calculate its
! value at tpresent. Which will be used in SOLVQV and SOLVQ.
!
!-----------------------------------------------------------------------
!
CALL set_acct(tinteg_acct)
CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
!call test_dump (w(1,1,1,tpresent),"XXX2Bwc_w1",nx,ny,nz,1,0)
CALL wcontra(nx,ny,nz, &
u(1,1,1,tpresent),v(1,1,1,tpresent), &
w(1,1,1,tpresent),mapfct,j1,j2,j3,aj3z, &
rhostr,tem1,tem2,tem3,wcont,tem4,tem5)
!
!-----------------------------------------------------------------------
!
! Since pressure was reset in SMLSTP. Now need to re-calculate its
! value at tfuture which will be used in microphysics.
!
!-----------------------------------------------------------------------
!
CALL setppi(nx,ny,nz,nt,tfuture,pprt,pbar,ppi)
!call test_dump (pprt(1,1,1,2),"XXXAsetppi_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,3),"XXXAsetppi_pprt3",nx,ny,nz,0,1)
IF( lvldbg >= 1 ) THEN
CALL acct_interrupt(misc_acct)
CALL checkwhx(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wcxbig', tem1)
CALL checkwhy(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wcybig', tem1)
CALL acct_stop_inter
END IF
!
!-----------------------------------------------------------------------
!
! Integrate the liquid water substance continuity equations
! forward one timestep.
!
! Sources and sinks due to phase changes and radiative processes
! are handled separately.
!
! If the simulation is for dry dynamics, then we skip over
! the moisture computations to save computer resources.
! The flag for this choice is "moist", i.e.,
!
! Moist = 0 for dry run
! Moist = 1 for moist run
!
!-----------------------------------------------------------------------
!
IF( moist == 1) THEN ! Determine if this run is dry or moist.
!
!-----------------------------------------------------------------------
!
! Water vapor equation. Array uforce and vforce are used as work
! arrays.
! Before we call the solve routine for the water varaibles, we
! need to compute ustr,vstr,wstr and store them into tem9,10,and 11
! These arrays cannot be altered and will be used in the
! following six moisture solve calls.
!
!
!-----------------------------------------------------------------------
!
IF (sadvopt == 1 .OR. sadvopt == 2 .OR. sadvopt == 3 ) THEN
! 2nd or 4th order advection
CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
tem9(i,j,k)=u(i,j,k,2)*tem1(i,j,k)
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
tem10(i,j,k)=v(i,j,k,2)*tem2(i,j,k)
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
tem11(i,j,k)=wcont(i,j,k)*tem3(i,j,k)
END DO
END DO
END DO
ELSE IF( sadvopt == 4.OR.sadvopt == 5) THEN ! FCT advection
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem3_0(i,j,k)=rhostr(i,j,k)
END DO
END DO
END DO
CALL acct_interrupt(bc_acct)
CALL extndsbc(tem3_0,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)
CALL acct_stop_inter
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx
tem9(i,j,k)=u(i,j,k,2)*(tem3_0(i-1,j,k)+tem3_0(i,j,k)) &
*mapfct(i,j,5)*0.5
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx-1
tem10(i,j,k)=v(i,j,k,2)*(tem3_0(i,j-1,k)+tem3_0(i,j,k)) &
*mapfct(i,j,6)*0.5
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
tem11(i,j,k)=wcont(i,j,k) &
*(tem3_0(i,j,k-1)+tem3_0(i,j,k))*0.5
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Water vapor equation. Array uforce and vforce are used as work
! arrays. Tem9,10,11 cannot be altered until after the last
! call to a solvq routine. They contain the appropriate
! ustr,vstr, and wstr.
! Water vapor equation. Array uforce and vforce are used as work
! arrays.
!
!-----------------------------------------------------------------------
CALL solvqv(nx,ny,nz, exbcbufsz, dtbig1, &
qv,u,v,wcont, tem9,tem10,tem11, &
sdteb,sdtwb,sdtnb,sdtsb, &
rhostr,rhostri,qvbar,kmh,kmv,rprntl,qvsflx,pbldpth, &
x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,j3inv, &
qcumsrc(1,1,1,1), &
usflx,vsflx,ptsflx,ptsfc,qvsfc,ptbar,ptprt, &
exbcbuf, bcrlx, &
uforce,vforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
!-----------------------------------------------------------------------
!
! Cloud water equation. Array uforce and vforce are used as work
! arrays. Tem9,10,11 cannot be altered until after the last
! call to a solvq routine. They contain the appropriate
! ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
qflag = 2
CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag, &
qc,u,v,wcont, tem9,tem10,tem11, &
sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri, &
kmh,kmv,rprntl,x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,j3inv, &
qcumsrc(1,1,1,2),qcumsrc(1,1,1,3), &
qcumsrc(1,1,1,4),qcumsrc(1,1,1,5), &
exbcbuf, bcrlx, &
uforce,vforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
IF( mphyopt /= 0) THEN
!
!-----------------------------------------------------------------------
!
! Rain water equation. Array uforce and vforce are used as work
! arrays. Tem9,10,11 cannot be altered until after the last
! call to a solvq routine. They contain the appropriate
! ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
qflag = 3
CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag, &
qr,u,v,wcont, tem9,tem10,tem11, &
sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri, &
kmh,kmv,rprntl,x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,j3inv, &
qcumsrc(1,1,1,2),qcumsrc(1,1,1,3), &
qcumsrc(1,1,1,4),qcumsrc(1,1,1,5), &
exbcbuf, bcrlx, &
uforce,vforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
END IF
IF( lvldbg >= 1 ) THEN
CALL acct_interrupt(misc_acct)
CALL checkshx(qv(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qvx', tem1)
CALL checkshy(qv(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qvy', tem1)
CALL checkshx(qc(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qcx', tem1)
CALL checkshy(qc(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qcy', tem1)
CALL checkshx(qr(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qrx', tem1)
CALL checkshy(qr(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qry', tem1)
CALL acct_stop_inter
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Integrate the frozen water substance continuity equations
! forward one timestep.
!
! Sources and sinks due to phase changes and radiative processes
! are handled separately.
!
! If the simulation does not include ice processes then, to save
! computer resources, we skip ice variable computations.
! The flags for these options are "moist" and "ice", i.e.,
!
! Moist = 0 dry run
! Moist = 1, Ice = 0, moist run without ice processes
! Moist = 1, Ice = 1, moist run with ice processes.
!
!-----------------------------------------------------------------------
!
IF( moist == 1 .AND. ice == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! Cloud ice equation. Array uforce and vforce are used as work
! arrays. Tem9,10,11 cannot be altered until after the last
! call to a solvq routine. They contain the appropriate
! ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
qflag = 4
CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag, &
qi,u,v,wcont, tem9,tem10,tem11, &
sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri, &
kmh,kmv,rprntl,x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,j3inv, &
qcumsrc(1,1,1,2),qcumsrc(1,1,1,3), &
qcumsrc(1,1,1,4),qcumsrc(1,1,1,5), &
exbcbuf, bcrlx, &
uforce,vforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
!
!-----------------------------------------------------------------------
!
! Snow equation. Array uforce and vforce are used as work
! arrays. Tem9,10,11 cannot be altered until after the last
! call to a solvq routine. They contain the appropriate
! ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
qflag = 5
CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag, &
qs,u,v,wcont, tem9,tem10,tem11, &
sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri, &
kmh,kmv,rprntl,x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,j3inv, &
qcumsrc(1,1,1,2),qcumsrc(1,1,1,3), &
qcumsrc(1,1,1,4),qcumsrc(1,1,1,5), &
exbcbuf, bcrlx, &
uforce,vforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
!
!-----------------------------------------------------------------------
!
! Hail equation. Array uforce and vforce are used as work
! arrays. Tem9,10,11 cannot be altered until after the last
! call to a solvq routine. They contain the appropriate
! ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
qflag = 6
CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag, &
qh,u,v,wcont, tem9,tem10,tem11, &
sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri, &
kmh,kmv,rprntl,x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,j3inv, &
qcumsrc(1,1,1,2),qcumsrc(1,1,1,3), &
qcumsrc(1,1,1,4),qcumsrc(1,1,1,5), &
exbcbuf, bcrlx, &
uforce,vforce, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8, &
tem1_0,tem2_0,tem3_0)
END IF
RETURN
END SUBROUTINE tinteg
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SMLSTEP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE smlstep(mptr, nx,ny,nz, exbcbufsz, dtbig1,dtsml1, & 1,21
u,v,w,wcont,pprt,ptprt, &
udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
phydro,ptdteb,ptdtwb,ptdtnb,ptdtsb, &
ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,rhostri, &
csndsq, &
uforce,vforce,wforce,pforce,ptforce, &
x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
exbcbuf,rhofct, &
rhostru,rhostrv,rhostrw, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10, &
tem11,tem12,tem13,tem14,tem15,tem16,tem17,tem18,tem19)
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Coordinate the integration of the acoustically active parts of
! the dynamic equations.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 3/21/92.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full documentation.
!
! 4/10/93 (M. Xue & Hao Jin)
! Add the terrain.
!
! 9/10/94 (D. Weber & Y. Lu)
! Cleaned up documentation.
!
! 10/31/95 (D. Weber)
! Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
! radiation condition.
!
! 1/25/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).
!
! 10/21/97 (Donghai Wang)
! Using total density (rho) in the calculation of the pressure
! gradient force terms, and added the second order terms
! in the linerized buoyancy terms.
!
! 11/05/97 (D. Weber)
! Added phydro array for use in the bottom boundary condition for
! perturbation 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.
!
! 9/18/98 (D. Weber)
! Added arrays aj3x,y,z, rhostri,ptbari,pbari, and tem9-12 for
! use in optimizing this routine.
! Tem9-19 is used to store commonly used groups of constants.
!
!-----------------------------------------------------------------------
!
! 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
!
! dtbig1 Large time step size for this call.
! dtsml1 Small time step size for this call.
!
! 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)
! wcont Contravariant vertical velocity (m/s)
! pprt Perturbation pressure at times tpast and tpresent (Pascal)
! ptprt Perturbation potential temperature at times tpast and
! tpresent (K)
!
! 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)
!
! ptdteb Time tendency of ptprt field at east boundary (K/s)
! ptdtwb Time tendency of ptprt field at west boundary (K/s)
! ptdtnb Time tendency of ptprt field at north boundary(K/s)
! ptdtsb Time tendency of ptprt field at south boundary(K/s)
!
! phydro Big time step forcing term for use in computing the
! hydrostatic pressure at k=1.
!
! 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)
! csndsq Sound wave speed squared.
!
! uforce Acoustically inactive forcing terms in u-eq. (kg/(m*s)**2)
! vforce Acoustically inactive forcing terms in v-eq. (kg/(m*s)**2)
! wforce Acoustically inactive forcing terms in w-eq. (kg/(m*s)**2)
! pforce Acoustically inactive forcing terms in p-eq. (Pascal/s)
! ptforce Gravity wave inactive forcing terms in pt-eq. (K*kg/(m**3*s))
!
! 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)
!
! 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.
! vwork2 2-D work array for fftopt = 2.
! wsave1 Work array for fftopt = 2.
! wsave2 Work array for fftopt = 2.
!
! rhostru Averaged rhostr at u points (kg/m**3).
! rhostrv Averaged rhostr at v points (kg/m**3).
! rhostrw Averaged rhostr at w points (kg/m**3).
!
! OUTPUT:
!
! u x component of velocity at time tfuture (m/s)
! v y component of velocity at time tfuture (m/s)
! w Vertical component of Cartesian velocity at tfuture (m/s)
! pprt Perturbation pressure at time tfuture (Pascal)
! ptprt Perturbation potential temperature at time tfuture (K)
!
! 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
! 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
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'timelvls.inc'
INTEGER :: mptr
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: dtbig1 ! Large time step size for this call.
REAL :: dtsml1 ! Small time step size for this call.
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 :: pprt (nx,ny,nz,nt) ! Perturbation pressure (Pascal)
REAL :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: udteb (ny,nz) ! Time tendency of u at e-boundary (m/s**2)
REAL :: udtwb (ny,nz) ! Time tendency of u at w-boundary (m/s**2)
REAL :: udtnb (nx,nz) ! Time tendency of u at n-boundary (m/s**2)
REAL :: udtsb (nx,nz) ! Time tendency of u at s-boundary (m/s**2)
REAL :: vdteb (ny,nz) ! Time tendency of v at e-boundary (m/s**2)
REAL :: vdtwb (ny,nz) ! Time tendency of v at w-boundary (m/s**2)
REAL :: vdtnb (nx,nz) ! Time tendency of v at n-boundary (m/s**2)
REAL :: vdtsb (nx,nz) ! Time tendency of v at s-boundary (m/s**2)
REAL :: wdteb (ny,nz) ! Time tendency of w at e-boundary (m/s**2)
REAL :: wdtwb (ny,nz) ! Time tendency of w at w-boundary (m/s**2)
REAL :: wdtnb (nx,nz) ! Time tendency of w at n-boundary (m/s**2)
REAL :: wdtsb (nx,nz) ! Time tendency of w at s-boundary (m/s**2)
REAL :: pdteb (ny,nz) ! Time tendency of pprt at e-boundary (Pascal/s)
REAL :: pdtwb (ny,nz) ! Time tendency of pprt at w-boundary (Pascal/s)
REAL :: pdtnb (nx,nz) ! Time tendency of pprt at n-boundary (Pascal/s)
REAL :: pdtsb (nx,nz) ! Time tendency of pprt at s-boundary (Pascal/s)
REAL :: phydro(nx,ny) ! Big time step forcing for computing
! hydrostatic pprt at k=1.
REAL :: ptdteb(ny,nz) ! Time tendency of ptprt at e-boundary (K/s)
REAL :: ptdtwb(ny,nz) ! Time tendency of ptprt at w-boundary (K/s)
REAL :: ptdtnb(nx,nz) ! Time tendency of ptprt at n-boundary (K/s)
REAL :: ptdtsb(nx,nz) ! Time tendency of ptprt at s-boundary (K/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) ! Inverse base state density rhobar times j3.
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 of the staggered grid.
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) ! Coordinate transformation Jacobian defined as
! d( zp )/d( z ).
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 :: uforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in u-momentum equation (kg/(m*s)**2)
REAL :: vforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in v-momentum equation (kg/(m*s)**2)
REAL :: wforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in w-momentum equation (kg/(m*s)**2)
REAL :: pforce(nx,ny,nz) ! Acoustically inactive forcing terms
! in pressure equation (Pascal/s)
REAL :: ptforce(nx,ny,nz) ! Gravity wave inactive forcing terms
! in potential temperature eq. (K*kg/(m**3*s))
REAL :: rhostru(nx,ny,nz) ! Averaged rhostr at u points (kg/m**3).
REAL :: rhostrv(nx,ny,nz) ! Averaged rhostr at v points (kg/m**3).
REAL :: rhostrw(nx,ny,nz) ! Averaged rhostr at w points (kg/m**3).
INTEGER :: exbcbufsz ! EXBC buffer size
REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
REAL :: rhofct(nx,ny,nz) ! rho-factor: rhobar/rho
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
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: ismstp
REAL :: curtsml ! Current time during small time step
! integration.
REAL :: tem,tema,temb
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'bndry.inc'
INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Integrate the momentum equations using a forward-in-time (relative
! to the pressure gradient terms) scheme and the pressure equation
! using a backward-in-time (relative to the mass divergence term)
! scheme.
!
! dtsml = (2*dtbig)/nsmstp is used so that after nsmstp small time
! step integrations, the fields are brought from time t-dtbig
! to t+dtbig, i.e. from time tpast to time future.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! The first time through the small timestep loop, set the values
! at time tfuture equal to those at tpast. By doing this, we
! simplify the forward-in-time integration because the equations
! take the form:
!
! u(i,j,k,future) = u(i,j,k,future)
! : + dtsml*( acoustic + non-acoustic forcing terms )
!
! i.e., we only use a single time level in this two-time level
! scheme and automatically finish with u at time tfuture.
!
!-----------------------------------------------------------------------
!
!call test_dump(w(1,1,tfuture),"XXXsml_wfut",nx,ny,nz,3,1)
!call test_dump(w(1,1,tpast),"XXXsml_wpast",nx,ny,nz,3,1)
!call test_dump(w(1,1,tpresent),"XXXsml_wpres",nx,ny,nz,3,1)
DO k=1,nz
DO j=1,ny
DO i=1,nx
u (i,j,k,tfuture)=u (i,j,k,tpast)
v (i,j,k,tfuture)=v (i,j,k,tpast)
w (i,j,k,tfuture)=w (i,j,k,tpast)
pprt (i,j,k,tfuture)=pprt (i,j,k,tpast)
END DO
END DO
END DO
IF(ptsmlstp == 1) THEN
IF( sadvopt == 4 ) THEN
!-----------------------------------------------------------------------
!
! When FCT is used to advect ptprt, all (e.g., advection,
! mixing) except for ptbar advection terms were evaluated at tpresent,
! and the large time step integration is forward in time.
! We have to do somthing special here, so that the same ptprt(tfuture)
! results had the ptbar advection not been included when ptprt is
! is updated in small steps. Therefore, we require
!
! (ptprt(tfuture)-ptprt(tpast))/(2*dtbig1)
! = (ptprt(tfuture)-ptprt(tpresent))/dtbig
!
!-----------------------------------------------------------------------
tem = 2*dtbig1/dtbig
IF(nestgrd == 1.AND.mgrid /= 1) THEN
! Don't touch the boundary values at tpast for nested grids.
DO k=1,nz-1
DO j=2,ny-2
DO i=2,nx-2
ptprt(i,j,k,tpast)=ptprt(i,j,k,tfuture) &
-(ptprt(i,j,k,tfuture)-ptprt(i,j,k,tpresent))*tem
END DO
END DO
END DO
ELSE
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
ptprt(i,j,k,tpast)=ptprt(i,j,k,tfuture) &
-(ptprt(i,j,k,tfuture)-ptprt(i,j,k,tpresent))*tem
END DO
END DO
END DO
END IF
END IF
DO k=1,nz
DO j=1,ny
DO i=1,nx
ptprt(i,j,k,tfuture)=ptprt(i,j,k,tpast)
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Calculate wcont at future time level.
!
!-----------------------------------------------------------------------
!
CALL rhouvw(nx,ny,nz,rhostr,rhostru,rhostrv,rhostrw)
!call test_dump (w(1,1,1,tfuture),"XXX3Bwc_w1",nx,ny,nz,1,0)
CALL wcontra(nx,ny,nz, &
u(1,1,1,tfuture),v(1,1,1,tfuture), &
w(1,1,1,tfuture),mapfct,j1,j2,j3,aj3z, &
rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)
IF( lvldbg >= 3 ) THEN
CALL checkwhx(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wcxsml', tem1)
CALL checkwhy(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1, &
'wcysml', tem1)
END IF
!-----------------------------------------------------------------------
!
! Compute a number of static variables (wrt small time step) and
! pass into the appropriate subroutines.
! NOTE: ALL OF THE COMPUTATIONS ASSUME THAT THE VARIABLES USED
! ARE CONSTANT DURING THE SMALL TIME STEP!!!!!!!
! THEY CANNOT BE OVERWRITTEN UNTIL AFTER THE W-P SOLVER.
!
!-----------------------------------------------------------------------
temb = 0.5*dtsml1
DO k=2,nz-2
DO j=1,ny-1
DO i=2,nx-1
tema = 1.0/rhostru(i,j,k)
tem13(i,j,k) =temb*(rhofct(i,j,k)+rhofct(i-1,j,k)) &
*mapfct(i,j,2)*tema
uforce(i,j,k) = dtsml1*uforce(i,j,k)*tema
END DO
END DO
END DO
!call test_dump (uforce,"XXX1uforce",nx,ny,nz,1,0)
!call test_dump (rhostru,"XXXrhostru",nx,ny,nz,1,0)
!call test_dump (tem13,"XXXdtmrxrsu_tem13",nx,ny,nz,1,0)
!call test_dump2(rhofct,"XXXdtmrxrsu_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
DO k=2,nz-2
DO j=2,ny-1
DO i=1,nx-1
tema = 1.0/rhostrv(i,j,k)
tem14(i,j,k) =temb*(rhofct(i,j,k)+rhofct(i,j-1,k)) &
*mapfct(i,j,3)*tema
vforce(i,j,k) = dtsml1*vforce(i,j,k)*tema
END DO
END DO
END DO
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem17(i,j,k)=1.0/(pbar(i,j,k)+pbar(i,j,k-1)) ! used in solvwp
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem12(i,j,k)=1.0/rhostrw(i,j,k) ! used in solvwp
tem9 (i,j,k)=rhostr(i,j,k)*j3inv(i,j,k) ! used in solvwp
tem10(i,j,k)=csndsq(i,j,k)*j3inv(i,j,k)*tem9(i,j,k) ! solvwp
tem11(i,j,k)=rhostr(i,j,k)*pbari(i,j,k) ! used in solvwp
tem18(i,j,k)=rhostr(i,j,k)*ptbari(i,j,k) ! used in solvwp
END DO
END DO
END DO
IF(vimplct == 0)THEN ! pre-compute wforce, etc.. for wp-ex.
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem19(i,j,k) =temb*(rhofct(i,j,k)+rhofct(i,j,k-1)) &
*tem12(i,j,k)
wforce(i,j,k) = dtsml1*wforce(i,j,k)*tem12(i,j,k)
END DO
END DO
END DO
END IF
tema = dtsml1*dxinv
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx
tem15(i,j,k)=tema*aj3x(i,j,k)*mapfct(i,j,5) ! used in solvwp's
END DO
END DO
END DO
temb = dtsml1*dyinv
DO k=2,nz-2
DO j=1,ny
DO i=1,nx-1
tem16(i,j,k)=temb*aj3y(i,j,k)*mapfct(i,j,6) ! used in solvwp's
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! For peqopt=2, set
!
! pforce=pforce (containging exbc damping only)
! + rhostr*cs**2/ptbar*d(pt)/dt
! =pforce (containging exbc damping only)
! + cs**2/ptbar * (ptforce + buoyancy if not already included).
!
! Note that adiabet heating/colling is not included in ptforce,
! therefore is not accounted for in p equation.
!
!-----------------------------------------------------------------------
!
IF(peqopt == 2) THEN
IF( ptsmlstp == 0 ) THEN
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pforce(i,j,k)=pforce(i,j,k) &
+csndsq(i,j,k)*ptforce(i,j,k)*ptbari(i,j,k)
END DO
END DO
END DO
ELSE
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
pforce(i,j,k)=pforce(i,j,k)+csndsq(i,j,k)*ptbari(i,j,k) &
*(ptforce(i,j,k)-tem9(i,j,k)* &
0.25*(ptbar(i,j,k+1)-ptbar(i,j,k-1))*dzinv* &
(w(i,j,k+1,2)+w(i,j,k,2)) )
END DO
END DO
END DO
END IF
END IF
DO k=2,nz-2 ! preparing pforce...for all solvwp codes...
DO j=1,ny-1
DO i=1,nx-1
pforce(i,j,k) = dtsml1*pforce(i,j,k)*j3inv(i,j,k)
END DO
END DO
END DO
DO ismstp =1, nsmstp
!
!-----------------------------------------------------------------------
!
! Calculate the current time within small time step. It will be
! used to calculate the external boundary fields.
!
!-----------------------------------------------------------------------
!
curtsml = curtim + dtbig - 2*dtbig1 + ismstp*dtsml1
!
!-----------------------------------------------------------------------
!
! Advance the three momentum equations one small timestep.
! (Arguments such as u(1,1,1,tfuture) indicate the passing of
! array values at time tfuture into the subroutine. Inside the
! subroutine, these arrays are defined at a single time level.
! u, v, w at tfuture and wcont are updated in time exiting
! this routine.
!
! Advance the pressure (p at tfuture) one small timestep.
!
! tem7 is used to transfer wpgrad between subroutine SOLVUV and
! subroutine SOLVWPIM or SOLVWPEX.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Calculate the boundary time tendencies of u and v using Klemp &
! Wilhelmson type open boundary condition (rbcopt=1).
!
!-----------------------------------------------------------------------
!
IF( rbcopt == 1 .AND. mgrid == 1 ) THEN
CALL bkwsmldt(nx,ny,nz, u(1,1,1,tfuture),v(1,1,1,tfuture), &
ubar,vbar,udteb,udtwb,vdtnb,vdtsb)
END IF
CALL solvuv(mptr, nx,ny,nz, exbcbufsz, dtsml1, curtsml, &
u(1,1,1,tfuture),v(1,1,1,tfuture), &
wcont,pprt(1,1,1,tfuture), &
udteb,udtwb,udtnb,udtsb, &
vdteb,vdtwb,vdtnb,vdtsb, &
rhostr, uforce,vforce,ubar,vbar, &
x,y,z,zp,mapfct, j1,j2,j3,j3inv, &
exbcbuf,rhofct, &
rhostru,rhostrv,rhostrw,tem13,tem14,tem7, &
tem1,tem2,tem3,tem4,tem5)
!
!-----------------------------------------------------------------------
!
! Warning: tem7 is carrying wpgrad can not be overwritten before
! the subroutine for w and p-equations are called.
!
!-----------------------------------------------------------------------
!
IF( lvldbg >= 3 ) THEN
CALL checkuhx( nx,ny,nz,1,nx, 1,ny-1,1,nz-1, 'uxsml', tem1)
CALL checkuhy( nx,ny,nz,1,nx, 1,ny-1,1,nz-1, 'uysml', tem1)
CALL checkvhx( nx,ny,nz,1,nx-1,1,ny, 1,nz-1, 'vxsml', tem1)
CALL checkvhy( nx,ny,nz,1,nx-1,1,ny, 1,nz-1, 'vysml', tem1)
END IF
!
!-----------------------------------------------------------------------
!
! Advance the w and pressure (at tfuture) one small timestep.
!
!-----------------------------------------------------------------------
!
IF( vimplct == 0 ) THEN ! Not implicit
!
!-----------------------------------------------------------------------
!
! Integrate the w and pressure equations one small time step
! using explicit scheme.
!
! NOTE: tem11,tem18,tem12,tem15,tem16,tem9 are assumed to be
! constant during the small time step!!!
!
!-----------------------------------------------------------------------
!
CALL solvwpex(mptr, nx,ny,nz, exbcbufsz, dtsml1, curtsml, &
u(1,1,1,tfuture),v(1,1,1,tfuture), &
w(1,1,1,tfuture),wcont, &
ptprt(1,1,1,tfuture),pprt(1,1,1,tfuture),phydro, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
rhostr,ptbar,ptbari,pbari,csndsq, &
wforce,tem7,pforce, &
x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
rhostru,rhostrv,rhostrw, &
exbcbuf,rhofct, &
tem1,tem2,tem3,tem4,tem5, &
tem11,tem18,tem19,tem15,tem16,tem9)
!call test_dump(w(1,1,tfuture),"XXXAsolvwpex_w",nx,ny,nz,3,1)
ELSE ! Implicit
!
!-----------------------------------------------------------------------
!
! Integrate the w and pressure equations one small time step
! using vertically implicit scheme.
!
! NOTE: tem9,tem10,tem11,tem12,tem15,tem16,tem17,tem18
! are assumed to be constant during the small time step!!!
!
!-----------------------------------------------------------------------
!
IF(peqopt == 1) THEN
CALL solvwpim(mptr,nx,ny,nz,exbcbufsz, dtsml1, curtsml, &
u(1,1,1,tfuture),v(1,1,1,tfuture), &
w(1,1,1,tfuture),wcont, &
ptprt(1,1,1,tfuture),pprt(1,1,1,tfuture),phydro, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
rhostr,ptbar,ptbari,pbari,csndsq, &
wforce,tem7,pforce, &
x,y,z,zp,mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
rhostru,rhostrv,rhostrw, &
exbcbuf,rhofct, &
tem1,tem2,tem3,tem4,tem5,tem6,tem8, &
tem9,tem10,tem11,tem12,tem15,tem16,tem17)
!call test_dump(w(1,1,tfuture),"XXXAsolvwpim_w",nx,ny,nz,3,1)
ELSE
CALL solvwpim1(nx,ny,nz,exbcbufsz, dtsml1, curtsml, &
u(1,1,1,tfuture),v(1,1,1,tfuture), &
w(1,1,1,tfuture),wcont, &
ptprt(1,1,1,tfuture),pprt(1,1,1,tfuture),phydro, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
rhostr,ptbar,ptbari,pbari,csndsq, &
wforce,tem7,pforce, &
x,y,z,zp,mapfct, &
j1,j2,j3,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
rhostru,rhostrv,rhostrw, &
exbcbuf,rhofct, &
tem1,tem2,tem3,tem4,tem5,tem6,tem8, &
tem11,tem18,tem12,tem17,tem9)
!call test_dump(w(1,1,tfuture),"XXXAsolvwpim_w",nx,ny,nz,3,1)
END IF
END IF
IF( lvldbg >= 3 ) THEN
CALL checkwhx( nx,ny,nz,1,nx-1,1,ny-1,1,nz, 'wxsml', tem1)
CALL checkwhy( nx,ny,nz,1,nx-1,1,ny-1,1,nz, 'wysml', tem1)
CALL checkshx(pprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pxsml', tem1)
CALL checkshy(pprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pysml', tem1)
END IF
!-----------------------------------------------------------------------
!
! Integrate the potential temperature equation one small time step
! when ptsmlstp=1.
!
!-----------------------------------------------------------------------
!
IF( ptsmlstp == 1 ) THEN
CALL solvpt_sml(nx,ny,nz, exbcbufsz, dtbig1,dtsml1,curtsml, &
ptprt,w, ptdteb,ptdtwb,ptdtnb,ptdtsb, &
ptbar,rhostr,rhostri,rhostrw,j3,aj3z, &
ptforce, exbcbuf, &
tem1,tem2)
IF( lvldbg >= 3) THEN
CALL checkshx(ptprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptsmlx',tem1)
CALL checkshy(ptprt(1,1,1,tfuture), &
nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptsmly',tem1)
END IF
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Integrate the potential temperature equation one large time step
! when ptsmlstp=0.
!
!-----------------------------------------------------------------------
!
IF( ptsmlstp == 0 ) THEN
CALL solvpt_lrg(nx,ny,nz, exbcbufsz, dtbig1, ptprt, &
ptdteb,ptdtwb,ptdtnb,ptdtsb,rhostr,rhostri, &
ptforce,exbcbuf,tem1)
END IF
RETURN
END SUBROUTINE smlstep
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TFILT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE tfilt(nx,ny,nz, & 1,12
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Apply the Asselin time filter to all time-dependent variables.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full documentation.
!
! 9/10/94 (D. Weber & 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
!
! 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)
! 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)
!
! OUTPUT:
!
! u x component of velocity at time tpresent (m/s)
! v y component of velocity at time tpresent(m/s)
! w Vertical component of Cartesian velocity at tpresent (m/s)
! ptprt Perturbation potential temperature at time tpresent (K)
! pprt Perturbation pressure at time tpresent (Pascal)
! qv Water vapor specific humidity at time tpresent (kg/kg)
! qc Cloud water mixing ratio at time tpresent (kg/kg)
! qr Rainwater mixing ratio at time tpresent (kg/kg)
! qi Cloud ice mixing ratio at time tpresent (kg/kg)
! qs Snow mixing ratio at time tpresent (kg/kg)
! qh Hail mixing ratio at time tpresent (kg/kg)
! tke Turbulent Kinetic Energy ((m/s)**2)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx, ny, nz ! 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 :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (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
!
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
CALL aselin(u, nx,ny,nz, 1,nx,1,ny-1,1,nz-1, flteps)
CALL aselin(v, nx,ny,nz, 1,nx-1,1,ny,1,nz-1, flteps)
CALL aselin(w, nx,ny,nz, 1,nx-1,1,ny-1,1,nz, flteps)
IF(sadvopt /= 4) CALL aselin(ptprt, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
CALL aselin(pprt, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
IF( tmixopt == 4) THEN
IF(sadvopt /= 4) CALL aselin(tke, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
END IF
IF( moist /= 0 .AND. sadvopt /= 4 ) THEN
CALL aselin(qv, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
CALL aselin(qc, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
CALL aselin(qr, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
IF( ice /= 0 ) THEN
CALL aselin(qi, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
CALL aselin(qs, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
CALL aselin(qh, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
END IF
END IF
RETURN
END SUBROUTINE tfilt
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE ASELIN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE aselin(a, nx,ny,nz, nx1,nx2,ny1,ny2,nz1,nz2, flteps) 12
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Apply Asselin time filter.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full documentation.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! a A variable at three time levels.
!
! nx1,nx2 i-index array bounds where the time filter is applied
! ny1,ny2 j-index array bounds where the time filter is applied
! nz1,nz2 k-index array bounds where the time filter is applied
!
! flteps The asselin time filter coefficient.
! Defined in globcst.inc
!
! OUTPUT:
!
! a The new array values at time tpresent after the time
! filter is applied.
!
!-----------------------------------------------------------------------
!
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: a(nx,ny,nz,nt)
INTEGER :: nx1, nx2
INTEGER :: ny1, ny2
INTEGER :: nz1, nz2
REAL :: flteps ! Coefficient of the Asselin time filter
!
!-----------------------------------------------------------------------
!
! Misc. local variable:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO k=nz1,nz2
DO j=ny1,ny2
DO i=nx1,nx2
a(i,j,k,tpresent)=a(i,j,k,tpresent)+flteps* &
(a(i,j,k,tfuture)-2*a(i,j,k,tpresent)+a(i,j,k,tpast))
END DO
END DO
END DO
RETURN
END SUBROUTINE aselin
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TFLIP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE tflip(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke) 1,12
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Update the fields u, v, w, ptprt ,pprt ,qv ,qc ,qr ,qi ,qs ,qh and
! tke in time. The fields at tpast are replaced by those at tpresent.
! The fields at tpresent are replaced by those at tfuture. The fields
! at tfuture are not changed.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full 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
!
! 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)
! 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)
!
! OUTPUT:
!
! u x-velocity at tpast and tpresent updated in time (m/s)
! v y-velocity at tpast and tpresent updated in time (m/s)
! w Vertical component of Cartesian velocity at tpast
! and tpresent updated in time (m/s)
! ptprt Perturbation potential temperature at tpast and tpresent
! updated in time (K)
! pprt Perturbation pressure at tpast and tpresent
! updated in time (Pascal)
! qv Water vapor specific humidity at tpast and tpresent
! updated in time (kg/kg)
! qc Cloud water mixing ratio at tpast and tpresent
! updated in time (kg/kg)
! qr Rainwater mixing ratio at tpast and tpresent
! updated in time (kg/kg)
! qi Cloud ice mixing ratio at tpast and tpresent
! updated in time (kg/kg)
! qs Snow mixing ratio at tpast and tpresent
! updated in time (kg/kg)
! qh Hail mixing ratio at tpast and tpresent
! updated in time (kg/kg)
! tke Turbulent Kinetic Energy ((m/s)**2)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx, ny, nz ! 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 :: ptprt (nx,ny,nz,nt) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz,nt) ! Perturbation pressure (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
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
CALL tswap(nx,ny,nz, u)
CALL tswap(nx,ny,nz, v)
CALL tswap(nx,ny,nz, w)
CALL tswap(nx,ny,nz, ptprt)
CALL tswap(nx,ny,nz, pprt)
IF( tmixopt == 4 ) THEN
CALL tswap(nx,ny,nz,tke )
END IF
IF( moist /= 0 ) THEN
CALL tswap(nx,ny,nz, qv)
CALL tswap(nx,ny,nz, qc)
CALL tswap(nx,ny,nz, qr)
IF( ice /= 0 ) THEN
CALL tswap(nx,ny,nz, qi)
CALL tswap(nx,ny,nz, qs)
CALL tswap(nx,ny,nz, qh)
END IF
END IF
RETURN
END SUBROUTINE tflip
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TSWAP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE tswap(nx,ny,nz, a) 12
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Update the time levels of a time-dependent array.
!
! a(*,*,*,tpast ) = a(*,*,*,tpresent)
! a(*,*,*,tpresent) = a(*,*,*,tfuture )
! a(*,*,*,tfuture ) is not changed.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 10/10/91.
!
! MODIFICATION HISTORY:
!
! 5/20/92 (K. Droegemeier and M. Xue)
! Added full 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
!
! a A 3-D array at three time levels that will be updated
! in time.
!
! OUPTUT:
!
! a Its new values at time tpast and tfuture.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'timelvls.inc'
INTEGER :: nx, ny, nz ! Number of grid points in 3 directions
REAL :: a(nx,ny,nz,nt) ! Array whose values will be updated in time
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
a(i,j,k,tpast )=a(i,j,k,tpresent)
a(i,j,k,tpresent)=a(i,j,k,tfuture )
END DO
END DO
END DO
RETURN
END SUBROUTINE tswap