!
!##################################################################
!##################################################################
!###### ######
!###### Advanced Regional Prediction System (ARPS) ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
PROGRAM arps,173
!
!-----------------------------------------------------------------------
!
! CONTACT:
!
! For further information, contact:
!
! Center for Analysis and Prediction of Storms
! University of Oklahoma
! Sarkeys Energy Center,
! 100 East Boyd
! Norman, OK 73019 USA
! Phone: (405) 325-0453
! FAX : (405) 325-7614
! E-mail: kkd@ou.edu or: mxue@ou.edu
!
! User support: arpssupport@ou.edu
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This program is the driver for the Advanced Regional Prediction
! System (ARPS) model.
!
! The ARPS is a fully three-dimensional, compressible,
! nonhydrostatic model based on the Arakawa C grid. The model
! solves three momentum equations, a thermodynamic energy
! equation, a pressure equation, six moisture equations (water
! vapor, cloud water, rainwater, cloud ice, hail, and snow),
! and also provides for the statistical representation of
! unresolvable processes via a turbulence parameterization scheme
! (user option). The model is written in a terrain-following
! coordinate system, and has provisions for several boundary
! condition options.
!
! For a full description of ARPS, please consult the ARPS 4.0
! User's Guide.
!
!-----------------------------------------------------------------------
!
! MODIFICATION HISTORY:
!
! The modification history can be found in file HISTORY.
!
!-----------------------------------------------------------------------
!
! GRID STRUCTURE AND STENCIL NUMBERING CONVENTION
!
! Each of the three velocity components is located at a physical
! boundary in the model. The numbering convention for the
! staggered grid is shown below for each of the coordinate
! directions, where S denotes a scalar variable.
!
! Arakawa C-grid is used by this model.
!
!
! X-Direction (East.......West):
!
! West boundary East boundary
! | |
! <-------- Physical Domain ------>
! +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
! U S U S U S U S U S U
!
! i= 1 1 2 2 . . . NX-2 NX-2 NX-1 NX-1 NX
!
!
! Y-Direction (North......South):
!
! South boundary North boundary
! | |
! <------- Physical Domain ------->
! +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
! V S V S V S V S V S V
!
! j= 1 1 2 2 . . . NY-2 NY-2 NY-1 NY-1 NY
!
!
! Z-Direction:
!
! NZ + W
! |
! |
! NZ-1 + S
! |
! |
! NZ-1 + W -- Physical Top boundary.
! |
! .
! .
!
! |
! 2 + S
! |
! |
! 2 + W -- Physical Ground Level
! |
! |
! 1 + S
! |
! |
! k= 1 + W
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE ! Force explicit declarations
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc' ! Global constants that control model
! execution
INCLUDE 'bndry.inc'
INCLUDE 'timelvls.inc'
INCLUDE 'radcst.inc' ! Includes radiation grid size information.
INCLUDE 'mp.inc' ! Message passing parameters.
INCLUDE 'alloc.inc' ! Memory allocation declaration and interfaces
INCLUDE 'nudging.inc' ! Memory allocation declaration and interfaces
!wdt tmp buffer
!INCLUDE "mpif.h"
!
!-----------------------------------------------------------------------
!
! Time dependent variables:
!
!-----------------------------------------------------------------------
REAL, ALLOCATABLE :: u (:,:,:,:) ! Total u-velocity (m/s)
REAL, ALLOCATABLE :: v (:,:,:,:) ! Total v-velocity (m/s)
REAL, ALLOCATABLE :: w (:,:,:,:) ! Total w-velocity (m/s)
REAL, ALLOCATABLE :: wcont (:,:,:) ! Contravariant vertical velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:,:) ! Perturbation potential temperature
! from that of base state atmosphere (K)
REAL, ALLOCATABLE :: pprt (:,:,:,:) ! Perturbation pressure from that
! of base state atmosphere (Pascal)
REAL, ALLOCATABLE :: qv (:,:,:,:) ! Water vapor specific humidity (kg/kg)
REAL, ALLOCATABLE :: qc (:,:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qr (:,:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qi (:,:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qs (:,:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qh (:,:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: tke (:,:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, ALLOCATABLE :: pbldpth(:,:,:) ! Planetary boundary layer depth (m)
REAL, ALLOCATABLE :: kmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: rprntl(:,:,:) ! Reciprocal of Prandtl number
!
!-----------------------------------------------------------------------
!
! Time tendencies of certain time dependent variables at the lateral
! boundaries, which are used to set the lateral boundary conditions
! for options 4 and 5.
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: udteb (:,:) ! Time tendency of u at e-boundary (m/s**2)
REAL, ALLOCATABLE :: udtwb (:,:) ! Time tendency of u at w-boundary (m/s**2)
REAL, ALLOCATABLE :: udtnb (:,:) ! Time tendency of u at n-boundary (m/s**2)
REAL, ALLOCATABLE :: udtsb (:,:) ! Time tendency of u at s-boundary (m/s**2)
REAL, ALLOCATABLE :: vdteb (:,:) ! Time tendency of v at e-boundary (m/s**2)
REAL, ALLOCATABLE :: vdtwb (:,:) ! Time tendency of v at w-boundary (m/s**2)
REAL, ALLOCATABLE :: vdtnb (:,:) ! Time tendency of v at n-boundary (m/s**2)
REAL, ALLOCATABLE :: vdtsb (:,:) ! Time tendency of v at s-boundary (m/s**2)
REAL, ALLOCATABLE :: wdteb (:,:) ! Time tendency of w at e-boundary (m/s**2)
REAL, ALLOCATABLE :: wdtwb (:,:) ! Time tendency of w at w-boundary (m/s**2)
REAL, ALLOCATABLE :: wdtnb (:,:) ! Time tendency of w at n-boundary (m/s**2)
REAL, ALLOCATABLE :: wdtsb (:,:) ! Time tendency of w at s-boundary (m/s**2)
REAL, ALLOCATABLE :: pdteb (:,:) ! Time tendency of pprt at e-boundary (Pascal/s)
REAL, ALLOCATABLE :: pdtwb (:,:) ! Time tendency of pprt at w-boundary (Pascal/s)
REAL, ALLOCATABLE :: pdtnb (:,:) ! Time tendency of pprt at n-boundary (Pascal/s)
REAL, ALLOCATABLE :: pdtsb (:,:) ! Time tendency of pprt at s-boundary (Pascal/s)
REAL, ALLOCATABLE :: sdteb (:,:) ! Time tendency of a scalar at e-boundary (m/s**2)
REAL, ALLOCATABLE :: sdtwb (:,:) ! Time tendency of a scalar at w-boundary (m/s**2)
REAL, ALLOCATABLE :: sdtnb (:,:) ! Time tendency of a scalar at n-boundary (m/s**2)
REAL, ALLOCATABLE :: sdtsb (:,:) ! Time tendency of a scalar at s-boundary (m/s**2)
!
!-----------------------------------------------------------------------
!
! Base state variables:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: ptbari(:,:,:) ! Inverse Base state pot. temperature (K)
REAL, ALLOCATABLE :: pbari (:,:,:) ! Inverse Base state pressure (Pascal)
REAL, ALLOCATABLE :: rhostr(:,:,:) ! Base state density rhobar times j3.
REAL, ALLOCATABLE :: rhostri(:,:,:)! Inv. base state density rhobar times j3.
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific humidity (kg/kg)
REAL, ALLOCATABLE :: ppi (:,:,:) ! Exner function.
REAL, ALLOCATABLE :: csndsq(:,:,:) ! Sound wave speed squared.
!
!-----------------------------------------------------------------------
!
! Arrays related to model grid definition:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x (:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, ALLOCATABLE :: y (:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, ALLOCATABLE :: z (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, ALLOCATABLE :: zp (:,:,:) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, ALLOCATABLE :: hterain(:,:) ! The height of the terrain.
REAL, ALLOCATABLE :: mapfct(:,:,:) ! Map factors at scalar, u and v points
REAL, ALLOCATABLE :: j1 (:,:,:) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( x ).
REAL, ALLOCATABLE :: j2 (:,:,:) ! Coordinate transformation Jacobian defined
! as - d( zp )/d( y ).
REAL, ALLOCATABLE :: j3 (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ).
REAL, ALLOCATABLE :: aj3x (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
REAL, ALLOCATABLE :: aj3y (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
REAL, ALLOCATABLE :: aj3z (:,:,:) ! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL, ALLOCATABLE :: j3inv (:,:,:) ! Inverse of J3.
REAL, ALLOCATABLE :: trigs1(:) ! Array containing pre-computed trig
! function for fftopt=1.
REAL, ALLOCATABLE :: trigs2(:) ! Array containing pre-computed trig
! function for fftopt=1.
INTEGER, ALLOCATABLE :: ifax1(:) ! Array containing the factors of nx for
! fftopt=1.
INTEGER, ALLOCATABLE :: ifax2(:) ! Array containing the factors of ny for
! fftopt=1.
REAL, ALLOCATABLE :: vwork1 (:,:) ! 2-D work array for fftopt=2.
REAL, ALLOCATABLE :: vwork2 (:,:) ! 2-D work array for fftopt=2.
REAL, ALLOCATABLE :: wsave1 (:) ! Work array for fftopt =2.
REAL, ALLOCATABLE :: wsave2 (:) ! Work array for fftopt =2.
REAL, ALLOCATABLE :: sinlat(:,:) ! Sin of latitude at each grid point
REAL, ALLOCATABLE :: ptcumsrc(:,:,:) ! Source term in pt-equation due
! to cumulus parameterization
REAL, ALLOCATABLE :: qcumsrc(:,:,:,:)! 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, ALLOCATABLE :: w0avg(:,:,:) ! a closing running average vertical
! velocity in 10min for K-F scheme
!
!-----------------------------------------------------------------------
!
! Arrays for soil model.
!
!-----------------------------------------------------------------------
!
INTEGER, ALLOCATABLE :: soiltyp(:,:,:) ! Soil type at each point
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Fraction of soil type
INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type at each point
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsfc (:,:,:) ! Ground temperature at surface (K)
REAL, ALLOCATABLE :: tsoil (:,:,:) ! Deep soil temperature (K)
REAL, ALLOCATABLE :: wetsfc (:,:,:) ! Surface soil moisture
REAL, ALLOCATABLE :: wetdp (:,:,:) ! Deep soil moisture
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: qvsfc (:,:,:) ! Effective sfc. qv
REAL, ALLOCATABLE :: ptsfc (:,:) ! Ground surface potential temperature (K)
REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rates (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulative precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, ALLOCATABLE :: kfraincv(:,:) ! K-F convective rainfall (cm)
INTEGER, ALLOCATABLE :: nca(:,:) ! K-F counter for CAPE release
!EMK BMJ
REAL, ALLOCATABLE :: cldefi(:,:) ! BMJ cloud efficiency
REAL, ALLOCATABLE :: xland(:,:) ! BMJ land mask
! (1.0 = land, 2.0 = sea)
REAL, ALLOCATABLE :: bmjraincv(:,:) ! BMJ convective rainfall (cm)
!EMK END
!
!-----------------------------------------------------------------------
!
! Arrays for radiation
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface
REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface
REAL, ALLOCATABLE :: rad2d (:,:,:) ! 2-D arrays for radiation (see radcst.inc)
REAL, ALLOCATABLE :: radbuf(:) ! Buffer to carry working arrays for
! radiation computing (see radcst.inc)
!
!-----------------------------------------------------------------------
!
! Arrays that carry surface fluxes
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: ptsflx(:,:) ! Surface heat flux (K*kg/(m**2*s))
REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
! Buffer arrays that carry external boundary 3-D arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: exbcbuf( : ) ! EXBC buffer array
REAL, ALLOCATABLE :: bcrlx(:,:) ! EXBC relaxation function coefficients
!
!-----------------------------------------------------------------------
!
! Arrays for analysis increment updating (a type of nudging) output
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: uincr(:,:,:)
REAL, ALLOCATABLE :: vincr(:,:,:)
REAL, ALLOCATABLE :: wincr(:,:,:)
REAL, ALLOCATABLE :: pincr(:,:,:)
REAL, ALLOCATABLE :: ptincr(:,:,:)
REAL, ALLOCATABLE :: qvincr(:,:,:)
REAL, ALLOCATABLE :: qcincr(:,:,:)
REAL, ALLOCATABLE :: qrincr(:,:,:)
REAL, ALLOCATABLE :: qiincr(:,:,:)
REAL, ALLOCATABLE :: qsincr(:,:,:)
REAL, ALLOCATABLE :: qhincr(:,:,:)
!
!-----------------------------------------------------------------------
!
! Work arrays that do not carry physical meaning in the code
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: temxy1(:,:) ! Temporary work array, used as phydro
REAL, ALLOCATABLE :: tem1 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem2 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem3 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem4 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem5 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem6 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem7 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem8 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem9 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem10 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem11 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem12 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem13 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem14 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem15 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem16 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem17 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem18 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem19 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem20 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem21 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem22 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem23 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem24 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem25 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem26 (:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem1_0(:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem2_0(:,:,:) ! Temporary work array.
REAL, ALLOCATABLE :: tem3_0(:,:,:) ! Temporary work array.
!
!-----------------------------------------------------------------------
!
! ARPS dimensions:
!
! nx, ny, nz: Dimensions of computatioal grid. When run on MPP
! with PVM or MPI, they represent of the size of the
! subdomains. See below.
!
! Given nx, ny and nz, the physical domain size will be
! xl=(nx-3)*dx by yl=(ny-3)*dy by zh=(nz-3)*dz. The
! variables nx, ny, nz, dx, dy and dz are read in from
! the input file by subroutine INITPARA.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx ! Number of grid points in the x-direction
INTEGER :: ny ! Number of grid points in the y-direction
INTEGER :: nz ! Number of grid points in the z-direction
INTEGER :: nxndg ! Number of x grid points for nudging (1 or nx)
INTEGER :: nyndg ! Number of y grid points for nudging (1 or ny)
INTEGER :: nzndg ! Number of z grid points for nudging (1 or nz)
!
!-----------------------------------------------------------------------
!
! Soil types.
!
!-----------------------------------------------------------------------
!
INTEGER :: nstyps ! Number of soil types
!
!-----------------------------------------------------------------------
!
! Radiation buffer size. Set equal to
! n2d_radiat*nx*ny + n3d_radiat*nx*ny*nz
! in subroutine inipara if radopt=2, otherwise set to 1.
!
!-----------------------------------------------------------------------
!
INTEGER :: rbufsz
!
!-----------------------------------------------------------------------
!
! External boundary buffer size. Set equal to
! 22*nx*ny*nz
! in subroutine inipara if lbcopt=2, otherwise set to 1.
!
!-----------------------------------------------------------------------
!
INTEGER :: exbcbufsz
INTEGER :: frstep ! Flag for the initial time step
INTEGER :: mptr ! Grid identifier.
INTEGER :: ierr, i,j,k
REAL :: eps ! Small value to compensate the roundoff
! error in checking of the end of time
! integration.
DATA eps /0.01/
INTEGER istatus, nxy, nxyz
REAL :: twall0, t_walltime
!wdt tmp buffer
!INTEGER :: n_buffer
!REAL, ALLOCATABLE :: mpi_buffer(:)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! First, initialize the model parameters. Most of them are
! global parameters passed among subroutines through common blocks.
!
! These parameters are declared in the include files:
! globcst.inc, bndry.inc, exbc.inc, radcst.inc etc.
!
!-----------------------------------------------------------------------
!
CALL acct_init
CALL set_acct(init_acct)
CALL initpara(nx,ny,nz,nstyps)
! Set radiation and external boundary work array buffer sizes:
IF (radopt == 2) THEN
rbufsz = n2d_radiat*nx*ny + n3d_radiat*nx*ny*nz
ELSE
rbufsz = 1
END IF
IF (lbcopt == 2) THEN
exbcbufsz = 22*nx*ny*nz
ELSE
exbcbufsz = 1
END IF
IF ( nudgopt > 0 ) THEN
nxndg=nx
nyndg=ny
nzndg=nz
ELSE
nxndg=1
nyndg=1
nzndg=1
END IF
!
!-----------------------------------------------------------------------
!
! Print a log file in namelist format.
!
!-----------------------------------------------------------------------
!
CALL acct_interrupt(output_acct)
IF (myproc == 0) THEN
IF (mp_opt == 0) THEN
CALL prtlog(nx,ny,nz,6)
CALL prtlog(nx,ny,nz,0)
ELSE
CALL prtlog(nproc_x*(nx-3)+3,nproc_y*(ny-3)+3,nz,6)
CALL prtlog(nproc_x*(nx-3)+3,nproc_y*(ny-3)+3,nz,0)
ENDIF
END IF
CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
! Allocate arrays etc.
!
!-----------------------------------------------------------------------
!
mptr = 1 ! Set the grid number to 1 for a single grid run.
nestgrd = 0 ! Set the grid nesting flag to zero, i.e. no nesting.
nxy = nx*ny
nxyz = nxy*nz
IF (myproc == 0) WRITE(6,*) 'Allocating arrays'
!wdt tmp buffer
!n_buffer = 100*(max(nx,ny)*nz+MPI_BSEND_OVERHEAD)
!ALLOCATE(mpi_buffer(n_buffer))
!CALL MPI_BUFFER_ATTACH(mpi_buffer,n_buffer*4,istatus)
ALLOCATE(u (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:u")
u = 0
ALLOCATE(v (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:v")
v = 0
ALLOCATE(w (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:w")
w = 0
ALLOCATE(wcont(nx,ny,nz ),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wcont")
wcont = 0
ALLOCATE(ptprt(nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptprt")
ptprt = 0
ALLOCATE(pprt (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pprt")
pprt = 0
ALLOCATE(qv (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qv")
qv = 0
ALLOCATE(qc (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qc")
qc = 0
ALLOCATE(qr (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qr")
qr = 0
ALLOCATE(qi (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qi")
qi = 0
ALLOCATE(qs (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qs")
qs = 0
ALLOCATE(qh (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qh")
qh = 0
ALLOCATE(tke (nx,ny,nz,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tke")
tke = 0
ALLOCATE(pbldpth(nx,ny,nt),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pbldpth")
pbldpth = 0
ALLOCATE(kmh(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:kmh")
kmh = 0
ALLOCATE(kmv(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:kmv")
kmv = 0
ALLOCATE(rprntl(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rprntl")
rprntl = 0
ALLOCATE(udteb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:udteb")
udteb = 0
ALLOCATE(udtwb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:udtwb")
udtwb = 0
ALLOCATE(udtnb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:udtnb")
udtnb = 0
ALLOCATE(udtsb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:udtsb")
udtsb = 0
ALLOCATE(vdteb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vdteb")
vdteb = 0
ALLOCATE(vdtwb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vdtwb")
vdtwb = 0
ALLOCATE(vdtnb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vdtnb")
vdtnb = 0
ALLOCATE(vdtsb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vdtsb")
vdtsb = 0
ALLOCATE(wdteb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wdteb")
wdteb = 0
ALLOCATE(wdtwb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wdtwb")
wdtwb = 0
ALLOCATE(wdtnb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wdtnb")
wdtnb = 0
ALLOCATE(wdtsb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wdtsb")
wdtsb = 0
ALLOCATE(pdteb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pdteb")
pdteb = 0
ALLOCATE(pdtwb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pdteb")
pdtwb = 0
ALLOCATE(pdtnb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pdtnb")
pdtnb = 0
ALLOCATE(pdtsb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pdtsb")
pdtsb = 0
ALLOCATE(sdteb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:sdteb")
sdteb = 0
ALLOCATE(sdtwb(ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:sdtwb")
sdtwb = 0
ALLOCATE(sdtnb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:sdtnb")
sdtnb = 0
ALLOCATE(sdtsb(nx,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:sdtsb")
sdtsb = 0
ALLOCATE(ubar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ubar")
ubar = 0
ALLOCATE(vbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vbar")
vbar = 0
ALLOCATE(ptbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptbar")
ptbar = 0
ALLOCATE(pbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pbar")
pbar = 0
ALLOCATE(ptbari(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptbari")
ptbari = 0
ALLOCATE(pbari (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pbari")
pbari = 0
ALLOCATE(rhostr(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rhostr")
rhostr = 0
ALLOCATE(rhostri(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rhostri")
rhostri = 0
ALLOCATE(qvbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qvbar")
qvbar = 0
ALLOCATE(ppi (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ppi")
ppi = 0
ALLOCATE(csndsq(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:csndsq")
csndsq = 0
ALLOCATE( CALL check_alloc_status(istatus, "arps:x")
x = 0
ALLOCATE( CALL check_alloc_status(istatus, "arps:y")
y = 0
ALLOCATE( CALL check_alloc_status(istatus, "arps:z")
z = 0
ALLOCATE(zp(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:zp")
zp = 0
ALLOCATE(hterain(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:hterain")
hterain = 0
ALLOCATE(mapfct(nx,ny,8),STAT=istatus)
CALL check_alloc_status(istatus, "arps:mapfct")
mapfct = 0
ALLOCATE(j1 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:j1")
j1 = 0
ALLOCATE(j2 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:j2")
j2 = 0
ALLOCATE(j3 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:j3")
j3 = 0
ALLOCATE(aj3x (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:aj3x")
aj3x = 0
ALLOCATE(aj3y (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:aj3y")
aj3y = 0
ALLOCATE(aj3z (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:aj3z")
aj3z = 0
ALLOCATE(j3inv(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:j3inv")
j3inv = 0
ALLOCATE(trigs1(3*(nx-1)/2+1),STAT=istatus)
CALL check_alloc_status(istatus, "arps:trigs1")
trigs1 = 0
ALLOCATE(trigs2(3*(ny-1)/2+1),STAT=istatus)
CALL check_alloc_status(istatus, "arps:trigs2")
trigs2 = 0
ALLOCATE(ifax1(13),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ifax1")
ifax1 = 0
ALLOCATE(ifax2(13),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ifax2")
ifax2 = 0
ALLOCATE(vwork1(nx+1,ny+1),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vwork1")
vwork1 = 0
ALLOCATE(vwork2(ny,nx+1),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vwork2")
vwork2 = 0
ALLOCATE(wsave1(3*(ny-1)+15),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wsave1")
wsave1 = 0
ALLOCATE(wsave2(3*(nx-1)+15),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wsave2")
wsave2 = 0
ALLOCATE(sinlat(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:sinlat")
sinlat = 0
ALLOCATE(ptcumsrc(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptcumsrc")
ptcumsrc = 0
ALLOCATE(qcumsrc(nx,ny,nz,5),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qcumsrc")
qcumsrc = 0
ALLOCATE(w0avg(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:w0avg")
w0avg = 0
ALLOCATE(soiltyp (nx,ny,nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:soiltyp")
soiltyp = 0
ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:stypfrct")
stypfrct = 0
ALLOCATE(vegtyp(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vegtyp")
vegtyp = 0
ALLOCATE(lai (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:lai")
lai = 0
ALLOCATE(roufns(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:roufns")
roufns = 0
ALLOCATE(veg (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:veg")
veg = 0
ALLOCATE(tsfc (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tsfc")
tsfc = 0
ALLOCATE(tsoil (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tsoil")
tsoil = 0
ALLOCATE(wetsfc (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wetsfc")
wetsfc = 0
ALLOCATE(wetdp (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wetdp")
wetdp = 0
ALLOCATE(wetcanp (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wetcanp")
wetcanp = 0
ALLOCATE(qvsfc (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qvsfc")
qvsfc = 0
ALLOCATE(ptsfc (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptsfc")
ptsfc = 0
ALLOCATE(snowdpth(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:snowdpth")
snowdpth = 0
ALLOCATE(raing (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:raing")
raing = 0
ALLOCATE(rainc (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rainc")
rainc = 0
ALLOCATE(kfraincv (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:kfraincv")
kfraincv = 0
ALLOCATE(nca (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:nca")
nca = 0
ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
CALL check_alloc_status(istatus, "arps:prcrate")
prcrate = 0
!EMK BMJ
ALLOCATE(cldefi(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:cldefi")
cldefi(:,:) = 1
ALLOCATE(xland(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:xland")
xland(:,:) = 0
ALLOCATE(bmjraincv (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:bmjraincv")
bmjraincv = 0
! EMK END
ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radfrc")
radfrc = 0
ALLOCATE(radsw (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radsw")
radsw = 0
ALLOCATE(rnflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rnflx")
rnflx = 0
ALLOCATE(rad2d (nx,ny,nrad2d),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rad2d")
rad2d = 0
ALLOCATE(radbuf(rbufsz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radbuf")
radbuf = 0
ALLOCATE(usflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:usflx")
usflx = 0
ALLOCATE(vsflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vsflx")
vsflx = 0
ALLOCATE(ptsflx(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptsflx")
ptsflx = 0
ALLOCATE(qvsflx(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qvsflx")
qvsflx = 0
ALLOCATE(exbcbuf(exbcbufsz ),STAT=istatus)
CALL check_alloc_status(istatus, "arps:exbcbuf")
exbcbuf = 0
ALLOCATE(bcrlx(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:bcrlx")
bcrlx = 0
IF (myproc == 0) print *, ' nxndg,nyndg,nzndg = ',nxndg,nyndg,nzndg
ALLOCATE(uincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:uincr")
uincr = 0.
ALLOCATE(vincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vincr")
vincr = 0.
ALLOCATE(wincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wincr")
wincr = 0.
ALLOCATE(pincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pincr")
pincr = 0.
ALLOCATE(ptincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptincr")
ptincr = 0.
ALLOCATE(qvincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qvincr")
qvincr = 0.
ALLOCATE(qcincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qcincr")
qcincr = 0.
ALLOCATE(qrincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qrincr")
qrincr = 0.
ALLOCATE(qiincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qiincr")
qiincr = 0.
ALLOCATE(qsincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qsincr")
qsincr = 0.
ALLOCATE(qhincr(nxndg,nyndg,nzndg),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qhincr")
qhincr = 0.
ALLOCATE(temxy1(nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:temxy1")
temxy1 = 0
ALLOCATE(tem1 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem1")
tem1 = 0
ALLOCATE(tem2 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem2")
tem2 = 0
ALLOCATE(tem3 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem3")
tem3 = 0
ALLOCATE(tem4 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem4")
tem4 = 0
ALLOCATE(tem5 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem5")
tem5 = 0
ALLOCATE(tem6 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem6")
tem6 = 0
ALLOCATE(tem7 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem7")
tem7 = 0
ALLOCATE(tem8 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem8")
tem8 = 0
ALLOCATE(tem9 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem9")
tem9 = 0
ALLOCATE(tem10(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem10")
tem10 = 0
ALLOCATE(tem11(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem11")
tem11 = 0
ALLOCATE(tem12(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem12")
tem12 = 0
ALLOCATE(tem13(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem13")
tem13 = 0
ALLOCATE(tem14(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem14")
tem14 = 0
ALLOCATE(tem15(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem15")
tem15 = 0
ALLOCATE(tem16(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem16")
tem16 = 0
ALLOCATE(tem17(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem17")
tem17 = 0
ALLOCATE(tem18(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem18")
tem18 = 0
ALLOCATE(tem19(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem19")
tem19 = 0
ALLOCATE(tem20(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem20")
tem20 = 0
ALLOCATE(tem21(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem21")
tem21 = 0
ALLOCATE(tem22(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem22")
tem22 = 0
ALLOCATE(tem23(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem23")
tem23 = 0
ALLOCATE(tem24(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem24")
tem24 = 0
ALLOCATE(tem25(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem25")
tem25 = 0
ALLOCATE(tem26(nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem26")
tem26 = 0
ALLOCATE(tem1_0(0:nx,0:ny,0:nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem1_0")
tem1_0 = 0
ALLOCATE(tem2_0(0:nx,0:ny,0:nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem2_0")
tem2_0 = 0
ALLOCATE(tem3_0(0:nx,0:ny,0:nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem3_0")
tem3_0 = 0
IF (myproc == 0) WRITE(6,*) 'Done allocating arrays'
!
!-----------------------------------------------------------------------
!
! Set up exbcbuf pointers.
!
!-----------------------------------------------------------------------
!
IF( lbcopt == 2 .AND. mptr == 1 ) CALL setexbcptr(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! Initialize all arrays and variables, and set initial
! condition for thermal disturbance if desired.
!
!-----------------------------------------------------------------------
!
!EMK BMJ
CALL initial(mptr,nx,ny,nz,nxndg,nyndg,nzndg,nstyps,exbcbufsz, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb, &
wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
sdteb,sdtwb,sdtnb,sdtsb, &
ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,rhostri, &
qvbar,ppi,csndsq, &
x,y,z,zp,hterain,mapfct, &
j1,j2,j3,aj3x,aj3y,aj3z,j3inv, &
trigs1,trigs2,ifax1,ifax2, &
wsave1,wsave2,vwork1,vwork2, &
sinlat,kmh,kmv,rprntl, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,ptsfc,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
raing,rainc,prcrate,exbcbuf,bcrlx,radfrc,radsw, &
rnflx,usflx,vsflx,ptsflx,qvsflx, &
uincr,vincr,wincr,pincr,ptincr,qvincr, &
qcincr,qrincr,qiincr,qsincr,qhincr, &
temxy1,tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11,tem12,tem13, &
tem14,tem15,tem16,tem17,tem18,tem19, &
tem20,tem21,tem22,tem23,tem24,tem25,tem26, &
tem1_0,tem2_0,tem3_0)
!EMK END
!call test_dump(u,"XXXAinitial_u",nx,ny,3*nz,1,1)
!call test_dump(v,"XXXAinitial_v",nx,ny,3*nz,2,1)
!call test_dump(w,"XXXAinitial_w",nx,ny,3*nz,3,1)
!call test_dump(wcont,"XXXAinitial_wcont",nx,ny,nz,0,1)
!call test_dump(ptprt,"XXXAinitial_ptprt",nx,ny,3*nz,0,1)
!call test_dump(pprt,"XXXAinitial_pprt",nx,ny,3*nz,0,1)
!call test_dump(qv,"XXXAinitial_qv",nx,ny,3*nz,0,1)
!call test_dump(qc,"XXXAinitial_qc",nx,ny,3*nz,0,1)
!call test_dump(qr,"XXXAinitial_qr",nx,ny,3*nz,0,1)
!call test_dump(qi,"XXXAinitial_qi",nx,ny,3*nz,0,1)
!call test_dump(qs,"XXXAinitial_qs",nx,ny,3*nz,0,1)
!call test_dump(qh,"XXXAinitial_qh",nx,ny,3*nz,0,1)
!call test_dump(tke,"XXXAinitial_tke",nx,ny,3*nz,0,1)
!call test_dump(udteb,"XXXAinitial_udteb",1,ny,nz,0,1)
!call test_dump(udtwb,"XXXAinitial_udtwb",1,ny,nz,0,1)
!call test_dump(udtnb,"XXXAinitial_udtnb",nx,1,nz,0,1)
!call test_dump(udtsb,"XXXAinitial_udtsb",nx,1,nz,0,1)
!call test_dump(vdteb,"XXXAinitial_vdteb",1,ny,nz,0,1)
!call test_dump(vdtwb,"XXXAinitial_vdtwb",1,ny,nz,0,1)
!call test_dump(vdtnb,"XXXAinitial_vdtnb",nx,1,nz,0,1)
!call test_dump(vdtsb,"XXXAinitial_vdtsb",nx,1,nz,0,1)
!call test_dump(wdteb,"XXXAinitial_wdteb",1,ny,nz,0,1)
!call test_dump(wdtwb,"XXXAinitial_wdtwb",1,ny,nz,0,1)
!call test_dump(wdtnb,"XXXAinitial_wdtnb",nx,1,nz,0,1)
!call test_dump(wdtsb,"XXXAinitial_wdtsb",nx,1,nz,0,1)
!call test_dump(pdteb,"XXXAinitial_pdteb",1,ny,nz,0,1)
!call test_dump(pdtwb,"XXXAinitial_pdtwb",1,ny,nz,0,1)
!call test_dump(pdtnb,"XXXAinitial_pdtnb",nx,1,nz,0,1)
!call test_dump(pdtsb,"XXXAinitial_pdtsb",nx,1,nz,0,1)
!call test_dump(sdteb,"XXXAinitial_sdteb",1,ny,nz,0,1)
!call test_dump(sdtwb,"XXXAinitial_sdtwb",1,ny,nz,0,1)
!call test_dump(sdtnb,"XXXAinitial_sdtnb",nx,1,nz,0,1)
!call test_dump(sdtsb,"XXXAinitial_sdtsb",nx,1,nz,0,1)
!call test_dump(ubar,"XXXAinitial_ubar",nx,ny,nz,1,1)
!call test_dump(vbar,"XXXAinitial_vbar",nx,ny,nz,2,1)
!call test_dump(ptbar,"XXXAinitial_ptbar",nx,ny,nz,0,1)
!call test_dump(pbar,"XXXAinitial_pbar",nx,ny,nz,0,1)
!call test_dump(ptbari,"XXXAinitial_ptbari",nx,ny,nz,0,1)
!call test_dump(pbari,"XXXAinitial_pbari",nx,ny,nz,0,1)
!call test_dump(rhostr,"XXXAinitial_rhostr",nx,ny,nz,0,1)
!call test_dump(rhostri,"XXXAinitial_rhostri",nx,ny,nz,0,1)
!call test_dump(qvbar,"XXXAinitial_qvbar",nx,ny,nz,0,1)
!call test_dump(ppi,"XXXAinitial_ppi",nx,ny,nz,0,1)
!call test_dump(csndsq,"XXXAinitial_csndsq",nx,ny,nz,0,1)
!call test_dump(x,"XXXAinitial_x",nx,1,1,1,1)
!call test_dump(y,"XXXAinitial_y",1,ny,1,2,1)
!call test_dump(z,"XXXAinitial_z",1,1,nz,3,1)
!call test_dump(zp,"XXXAinitial_zp",nx,ny,nz,3,1)
!call test_dump(hterain,"XXXAinitial_hterain",nx,ny,1,0,1)
!call test_dump(mapfct,"XXXAinitial_mapfct",nx,ny,8,0,1)
!call test_dump(j1,"XXXAinitial_j1",nx,ny,nz,1,1)
!call test_dump(j2,"XXXAinitial_j2",nx,ny,nz,2,1)
!call test_dump(j3,"XXXAinitial_j3",nx,ny,nz,3,1)
!call test_dump(aj3x,"XXXAinitial_aj3x",nx,ny,nz,1,1)
!call test_dump(aj3y,"XXXAinitial_aj3y",nx,ny,nz,2,1)
!call test_dump(aj3z,"XXXAinitial_aj3z",nx,ny,nz,3,1)
!call test_dump(j3inv,"XXXAinitial_j3inv",nx,ny,nz,3,1)
!call test_dump(sinlat,"XXXAinitial_sinlat",nx,ny,1,0,1)
!call test_dump(kmh,"XXXAinitial_kmh",nx,ny,nz,0,1)
!call test_dump(kmv,"XXXAinitial_kmv",nx,ny,nz,0,1)
!call test_dump(rprntl,"XXXAinitial_rprntl",nx,ny,nz,0,1)
!call test_dump(soiltyp,"XXXAinitial_soiltyp",nx,ny,nstyps,0,1)
!call test_dump(stypfrct,"XXXAinitial_stypfrct",nx,ny,nstyps,0,1)
!call test_dump(vegtyp,"XXXAinitial_vegtyp",nx,ny,1,0,1)
!call test_dump(lai,"XXXAinitial_lai",nx,ny,1,0,1)
!call test_dump(roufns,"XXXAinitial_roufns",nx,ny,1,0,1)
!call test_dump(veg,"XXXAinitial_veg",nx,ny,1,0,1)
!call test_dump(tsfc,"XXXAinitial_tsfc",nx,ny,nstyps+1,0,1)
!call test_dump(tsoil,"XXXAinitial_tsoil",nx,ny,nstyps+1,0,1)
!call test_dump(wetsfc,"XXXAinitial_wetsfc",nx,ny,nstyps+1,0,1)
!call test_dump(wetdp,"XXXAinitial_wetdp",nx,ny,nstyps+1,0,1)
!call test_dump(wetcanp,"XXXAinitial_wetcanp",nx,ny,nstyps+1,0,1)
!call test_dump(snowdpth,"XXXAinitial_snowdpth",nx,ny,1,0,1)
!call test_dump(ptsfc,"XXXAinitial_ptsfc",nx,ny,1,0,1)
!call test_dump(qvsfc,"XXXAinitial_qvsfc",nx,ny,1,0,1)
!call test_dump(bcrlx,"XXXAinitial_bcrlx",nx,ny,1,0,1)
!
!-----------------------------------------------------------------------
!
! The current model time (initial time):
!
!-----------------------------------------------------------------------
!
curtim = tstart
IF( myproc == 0 ) WRITE(6,'(1x,a,f13.3,a)') &
'The initial model time is at ', curtim,' seconds.'
!
!-----------------------------------------------------------------------
!
! Output the initial fields
!
!-----------------------------------------------------------------------
!
CALL set_acct(output_acct)
CALL initout(mptr,nx,ny,nz,nstyps,exbcbufsz, &
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,hterain,mapfct,j1,j2,j3,j3inv, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
raing,rainc,prcrate,exbcbuf, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11)
IF (mp_opt > 0) THEN
CALL mpbarrier ()
END IF
CALL set_acct(misc_acct)
nstep = nint( curtim/dtbig )
!
!-----------------------------------------------------------------------
!
! Time integration loop begins ----------------------------------->
!
!-----------------------------------------------------------------------
!
900 CONTINUE ! Entry point of time step integration.
!
!-----------------------------------------------------------------------
!
! Define the time step counter nstep reference to time zero
! (curtim=0.0). This means for tstart.ne.0 case, nstep for the
! first step of integration is not 1.
!
!-----------------------------------------------------------------------
!
nstep = nstep + 1
IF( (ABS(curtim-tstart) <= 1.0E-10) .AND. (restrt /= 1) ) THEN
frstep=1 ! Indicate that this is the initial step of
! model integration. For this step, a forward
! time integration scheme will be used.
ELSE ! For non-first step or restart run
frstep=0 ! When frstep=0, leapfrog scheme is used.
END IF
!
!-----------------------------------------------------------------------
!
! Perform one time step integration for all equations:
! On exit of this routine, all time dependent fields are
! advanced by one time step.
!
!-----------------------------------------------------------------------
!
!EMK
CALL cordintg(mptr,frstep, nx,ny,nz,nxndg,nyndg,nzndg, &
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, &
temxy1,tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11,tem12,tem13, &
tem14,tem15,tem16,tem17,tem18,tem19, &
tem20,tem21,tem22,tem23,tem24,tem25,tem26, &
tem1_0,tem2_0,tem3_0)
!EMK BMJ
CALL set_acct(misc_acct)
!
!-----------------------------------------------------------------------
!
! Update physical time and generate output files
!
!-----------------------------------------------------------------------
!
curtim = curtim + dtbig
!
!-----------------------------------------------------------------------
!
! Print timestep marker
!
!-----------------------------------------------------------------------
!
!wdt
IF (myproc == 0) &
WRITE(6,'(a,i8,a,f10.2,a,f8.3,a)') &
' End of integration time step ', nstep, &
', Model time=',curtim,' s (', curtim/3600.0, ' h)'
CALL flush (6) ! GMB wdt
CALL set_acct(output_acct)
CALL output(mptr,nx,ny,nz,nstyps,exbcbufsz, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
udteb,udtwb,vdtnb,vdtsb,pdteb,pdtwb,pdtnb,pdtsb, &
ubar,vbar,ptbar,pbar,rhostr,qvbar,kmh,kmv, &
x,y,z,zp,hterain, mapfct, j1,j2,j3,j3inv, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc, &
ptcumsrc,qcumsrc,w0avg,nca,kfraincv, &
cldefi,xland,bmjraincv, &
raing,rainc,prcrate,exbcbuf, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
tem1, tem2,tem3,tem4,tem5, tem6,tem7, &
tem8,tem9,tem10,tem11)
CALL set_acct(misc_acct)
!
!-----------------------------------------------------------------------
!
! 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.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem11(i,j,k) = rhostr(i,j,k)*j3inv(i,j,k)
END DO
END DO
END DO
CALL chkstab(mptr,nx,ny,nz,nstyps, &
u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke, &
ubar,vbar,ptbar,pbar,tem11,qvbar,kmh,kmv, &
x,y,z,zp,hterain,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)
!
!-----------------------------------------------------------------------
!
! Calculate and apply the adjustment of domain translation
! speed using either cell-tracking or optimal mean speed algorithm.
!
!-----------------------------------------------------------------------
!
CALL grdtran(nx,ny,nz,ubar,vbar,u,v,w,ptprt,pprt, &
qv,qc,qr,qi,qs,qh,qvbar,rhostr,x,y,zp,j3,j3inv, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7, &
tem8,tem9,tem10,tem11)
!
!-----------------------------------------------------------------------
!
! Time integration loop ends, go back to the beginning of the
! time step loop if stop time is not reached.
!
!-----------------------------------------------------------------------
!
IF( curtim+eps*dtbig < tstop ) GO TO 900
!
!-----------------------------------------------------------------------
!
! End of entire model time integration. The program stops.
!
! Close opened files.
!
!-----------------------------------------------------------------------
!
IF (hdmpfmt == 5) THEN
CALL mclosescheme (gridid, ierr)
CALL mclosedataset (dsindex, ierr)
ELSE IF (hdmpfmt == 9) THEN
CLOSE (nchdmp)
END IF
CALL acct_finish
CALL acct_report_arps
IF (myproc == 0) THEN
WRITE(6,'(//1x,a,/1x,a,f13.3,a,/1x,a)') &
'ARPS stopped normally in the main program. ', &
'The ending time was ', curtim,' seconds.', &
'Thanks for using ARPS.'
WRITE (6,*) "Maxumum memory allocation (in words):", &
max_memory_use
END IF
IF (mp_opt > 0) THEN
CALL mpexit(0)
END IF
STOP
END PROGRAM ARPS