PROGRAM arpsensic,25
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSENSIC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Generate one perturbation initial condition files from two
! sets of ARPS output files and write it out for the use in
! ENSEMBLE forecast.
! The idea is based on the SLAF (Scaled Lagged Average Forecast).
! The procedure is as follows:
! 1, read file a;
! 2, read in file b;
! 3, find the difference bwtween a and b, store it in a. a=a-b
! 4, read in the base file c or use b as the control file (iread=0).
! 5, generate the perturbation c=c+a/n (b is used as c in code)
! 6, n is input as the variable iorder; it can be ... -2,-1,+1,+2 ...
!
! It shares with the model the include file dims.inc for
! definition of dimensions and domain size. a,b,c have the same
! dimensions and the same grid structure.
!
! Parameters grdout,varout,mstout,iceout and trbout should be input
! with the same values as in the data dump subroutines in the model.
!
! AUTHOR: Dingchen Hou
! History: Apr. 30, 1998: developed from the framework of ARPSDIFF.
! Sep. 15, 1999: modified to include soil variables in
! perturbation change input to namelist format.
! Feb-Apr, 2002: (F.KONG) Major modifications to:
! - accommodate BGM (Breeding Fast Growing Mode) IC
! generation, including the initial random perturbation
! procedure (with inibred = 1 and specified iseed).
! During the regular breeding cycles, certain scale
! factors are calculated to control the amplitude of
! the growing perturbations (with iensopt = 1). When
! generating BGM IC, the two 12h forecasts from the
! paired breeding members are specified as a and b,
! and the analysis data (c) must be read in (iread=1).
!
! When generating initial BGM IC, the only one data
! needed is the analysis (both a and b)
!
! - read in domain config directly from data
!
! MODIFIED:
!
! 05/28/2002 (J. Brotzge)
! Added tsoil/qsoil to accomodate new soil scheme.
!
! 1 June 2002 Eric Kemp
! Soil variable updates
!
!-----------------------------------------------------------------------
!
! DATA ARRAYS READ IN:
!
! 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 z coordinate of grid points in physical space (m)
! zpsoil z coordinate of grid points in the soil (m)
!
! u x component of velocity (m/s)
! v y component of velocity (m/s)
! w vertical component of velocity in Cartesian
! coordinates (m/s).
!
! ptprt perturbation potential temperature (K)
! pprt perturbation pressure (Pascal)
!
! qv water vapor mixing ratio (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! wbar Base state z velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhobar Base state density (kg/m**3)
! qvbar Base state water vapor mixing ratio (kg/kg)
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m**3/m**3)
! wetcanp Canopy water amount
!
! raing Grid supersaturation rain
! rainc Cumulus convective rain
!
! CALCULATED DATA ARRAYS:
!
! uprt perturbation x component of velocity (m/s)
! vprt perturbation y component of velocity (m/s)
! wprt perturbation z component of velocity (m/s)
!
! WORK ARRAYS:
!
! tem1 Temporary work array.
! tem2 Temporary work array.
! tem3dsoil Work array
! vtem3dsoil Work array
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'indtflg.inc'
INTEGER :: nx,ny,nz,nzsoil
INTEGER :: vnx,vny,vnz,vnzsoil
! PARAMETER (vnx=nx,vny=ny,vnz=nz,vnzsoil=nzsoil)
!
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
INTEGER :: nstyps ! Maximum number of soil types in each
! grid box
! PARAMETER (nstyps=4)
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 :: zpsoil(:,:,:) ! The physical height coordinate defined at
! w-point of the soil grid.
REAL, ALLOCATABLE :: uprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: wprt (:,:,:) ! Perturbation w-velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal)
REAL, ALLOCATABLE :: qvprt (:,:,:) ! Perturbation water vapor specific humidity
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 :: kmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: wbar (:,:,:) ! Base state w-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: rhobar (:,:,:) ! Base state air density (kg/m**3)
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific humidity
INTEGER, ALLOCATABLE :: soiltyp (:,:,:)! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type
INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! Soil temperature (K)
REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! Soil moisture (m**3/m**3)
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: raing (:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: rainc (:,:) ! Cumulus convective rain
!
!-----------------------------------------------------------------------
!
! Verification Arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: vx (:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, ALLOCATABLE :: vy (:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, ALLOCATABLE :: vz (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, ALLOCATABLE :: vzp (:,:,:) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, ALLOCATABLE :: vzpsoil(:,:,:) ! The physical height coordinate defined at
! w-point of the soil grid.
REAL, ALLOCATABLE :: vuprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vvprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: vwprt (:,:,:) ! Perturbation w-velocity (m/s)
REAL, ALLOCATABLE :: vptprt (:,:,:) ! Perturbation potential temperature (K)
REAL, ALLOCATABLE :: vpprt (:,:,:) ! Perturbation pressure (Pascal)
REAL, ALLOCATABLE :: vqvprt (:,:,:) ! Perturbation water vapor specific humidity
REAL, ALLOCATABLE :: vqc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: vqr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: vqi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: vqs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: vqh (:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: vtke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, ALLOCATABLE :: vkmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: vkmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: vubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vvbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: vwbar (:,:,:) ! Base state w-velocity (m/s)
REAL, ALLOCATABLE :: vptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: vpbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: vrhobar(:,:,:) ! Base state air density (kg/m**3)
REAL, ALLOCATABLE :: vqvbar (:,:,:) ! Base state water vapor specific humidity
INTEGER, ALLOCATABLE :: vsoiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: vstypfrct(:,:,:) ! Soil type
INTEGER, ALLOCATABLE :: vvegtyp (:,:) ! Vegetation type
REAL, ALLOCATABLE :: vlai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: vroufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: vveg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: vtsoil (:,:,:,:) ! Soil temperature (K)
REAL, ALLOCATABLE :: vqsoil (:,:,:,:) ! Soil moisture (m**3/m**3)
REAL, ALLOCATABLE :: vwetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: vsnowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: vraing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: vrainc(:,:) ! Cumulus convective rain
!
!-----------------------------------------------------------------------
!
! Work Arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, 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 :: radswnet(:,:) ! Net shortwave radiation
REAL, ALLOCATABLE :: radlwin(:,:) ! Incoming longwave radiation
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*s**2))
REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s))
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
REAL, ALLOCATABLE :: tem3dsoil(:,:,:)
REAL, ALLOCATABLE :: vtem1(:,:,:)
REAL, ALLOCATABLE :: vtem2(:,:,:)
REAL, ALLOCATABLE :: vtem3(:,:,:)
REAL, ALLOCATABLE :: vtem3dsoil(:,:,:)
REAL, ALLOCATABLE :: xs(:)
REAL, ALLOCATABLE :: ys(:)
REAL, ALLOCATABLE :: zps(:,:,:)
REAL, ALLOCATABLE :: x2d(:,:)
REAL, ALLOCATABLE :: y2d(:,:)
REAL, ALLOCATABLE :: lat(:,:),lon(:,:)
REAL, ALLOCATABLE :: vxs(:)
REAL, ALLOCATABLE :: vys(:)
REAL, ALLOCATABLE :: vzps(:,:,:)
REAL, ALLOCATABLE :: dxfld(:)
REAL, ALLOCATABLE :: dyfld(:)
REAL, ALLOCATABLE :: rdxfld(:)
REAL, ALLOCATABLE :: rdyfld(:)
INTEGER, ALLOCATABLE :: iloc(:,:),jloc(:,:)
REAL, ALLOCATABLE :: zpver(:,:,:)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: fcrnam,runnmin
CHARACTER (LEN=132) :: filename(3),grdbasfn(3)
INTEGER :: lengbf(3),lenfil(3)
INTEGER :: ifproj,ivproj
REAL :: flatnot(2),vlatnot(2)
REAL :: fscale,ftrulon,fdx,fdy,fx0,fy0
REAL :: fctrlat,fctrlon
REAL :: vscale,vtrulon,vdx,vdy,vx0,vy0
REAL :: vctrlat,vctrlon
REAL :: time,xctr,yctr
INTEGER :: i,j,k
INTEGER :: grdbas
INTEGER :: iorder,hinfmt,iread,isread,iensopt,inibred,iseed
INTEGER :: ireturn
LOGICAL :: comcoord
INTEGER :: nch
INTEGER :: zpsoilin,tsoilin,qsoilin,wcanpin,snowdin
!
INTEGER :: istatus
REAL :: utot,utot2,usd,uscl
REAL :: vtot,vtot2,vsd,vscl
REAL :: wtot,wtot2,wsd
REAL :: pttot,pttot2,ptsd,ptscl
REAL :: qvtot,qvtot2,qvsd,qvscl
REAL :: ptot,ptot2,psd,pscl
REAL :: totalpoint
NAMELIST /prtbpara/ iensopt,inibred,iseed,iorder,isread,soilinfl,iread
NAMELIST /input_data/ hinfmt,grdbasfn,filename
NAMELIST /outpt_data/ hdmpfmt,runnmin,grdout,basout,varout,mstout, &
iceout,trbout,sfcout,rainout,snowout,filcmprs
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
WRITE(6,'(/6(/5x,a)/)') &
'###############################################################', &
'# #', &
'# Welcome to ARPSENSIC, a program that reads in history files #', &
'# generated by ARPS and produces perturbation grids variables #', &
'# #', &
'###############################################################'
!------------------------------------------------------------------------
! Read namelist
!------------------------------------------------------------------------
READ(5,prtbpara,END=999)
! PRINT *,iensopt,inibred,iseed,iorder,isread,soilinfl,iread
write(6,prtbpara)
READ(5,input_data,END=999)
! PRINT *,hinfmt, grdbasfn, filename
write(6,input_data)
READ(5,outpt_data,END=999)
! PRINT *,hdmpfmt, runnmin, &
! grdout,basout,varout,mstout, &
! iceout,trbout,sfcout,rainout,snowout,filcmprs
write(6,outpt_data)
GO TO 1001
999 PRINT *, 'ERROR in reading from input file, Program Terminated'
STOP
1001 CONTINUE
! WRITE(6,'(/a,i5/a,i5)') &
! ' The scaling factor is :',iorder, &
! ' and the soil field reading flag is :',isread, &
! ' The 3rd set data reading flag is: ',iread
!
! WRITE(6,'(a,i5)')' The data format flag value is: ',hinfmt
!------------------------------------------------------------------------
!
! Get the dimensions of the input data files
!
!------------------------------------------------------------------------
lengbf(1) = len_trim(grdbasfn(1))
CALL get_dims_from_data
(hinfmt,grdbasfn(1)(1:lengbf(1)), &
nx,ny,nz,nzsoil,nstyps, ireturn)
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,' nzsoil =',nzsoil
print*,'nstyps =', nstyps
vnx=nx
vny=ny
vnz=nz
vnzsoil = nzsoil
totalpoint = nx*ny*nz
!-----------------------------------------------------------------------
!
! Allocate the variables and initialize the them to zero
!
!-----------------------------------------------------------------------
ALLOCATE( x=0
ALLOCATE( y=0
ALLOCATE( z=0
ALLOCATE(zp(nx,ny,nz),STAT=istatus)
zp=0
ALLOCATE(zpsoil(nx,ny,nzsoil),STAT=istatus)
zpsoil=0
ALLOCATE(uprt(nx,ny,nz),STAT=istatus)
uprt=0
ALLOCATE(vprt(nx,ny,nz),STAT=istatus)
vprt=0
ALLOCATE(wprt(nx,ny,nz),STAT=istatus)
wprt=0
ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
ptprt=0
ALLOCATE(pprt(nx,ny,nz),STAT=istatus)
pprt=0
ALLOCATE(qvprt(nx,ny,nz),STAT=istatus)
qvprt=0
ALLOCATE(qc(nx,ny,nz),STAT=istatus)
qc=0
ALLOCATE(qr(nx,ny,nz),STAT=istatus)
qr=0
ALLOCATE(qi(nx,ny,nz),STAT=istatus)
qi=0
ALLOCATE(qs(nx,ny,nz),STAT=istatus)
qs=0
ALLOCATE(qh(nx,ny,nz),STAT=istatus)
qh=0
ALLOCATE(tke(nx,ny,nz),STAT=istatus)
tke=0
ALLOCATE(kmh(nx,ny,nz),STAT=istatus)
kmh=0
ALLOCATE(kmv(nx,ny,nz),STAT=istatus)
kmv=0
ALLOCATE(ubar(nx,ny,nz),STAT=istatus)
ubar=0
ALLOCATE(vbar(nx,ny,nz),STAT=istatus)
vbar=0
ALLOCATE(wbar(nx,ny,nz),STAT=istatus)
wbar=0
ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
ptbar=0
ALLOCATE(pbar(nx,ny,nz),STAT=istatus)
pbar=0
ALLOCATE(rhobar(nx,ny,nz),STAT=istatus)
rhobar=0
ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
qvbar=0
ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus)
soiltyp=0
ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
stypfrct=0
ALLOCATE(vegtyp(nx,ny),STAT=istatus)
vegtyp=0
ALLOCATE(lai(nx,ny),STAT=istatus)
lai=0
ALLOCATE(roufns(nx,ny),STAT=istatus)
roufns=0
ALLOCATE(veg(nx,ny),STAT=istatus)
veg=0
ALLOCATE(tsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus)
tsoil=0
ALLOCATE(qsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus)
qsoil=0
ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus)
wetcanp=0
ALLOCATE(snowdpth(nx,ny),STAT=istatus)
snowdpth=0
ALLOCATE(raing(nx,ny),STAT=istatus)
raing=0
ALLOCATE(rainc(nx,ny),STAT=istatus)
rainc=0
ALLOCATE(vx(vnx),STAT=istatus)
vx=0
ALLOCATE(vy(vny),STAT=istatus)
vy=0
ALLOCATE(vz(vnz),STAT=istatus)
vz=0
ALLOCATE(vzp(vnx,vny,vnz),STAT=istatus)
vzp=0
ALLOCATE(vzpsoil(vnx,vny,vnzsoil),STAT=istatus)
vzpsoil=0
ALLOCATE(vuprt(vnx,vny,vnz),STAT=istatus)
vuprt=0
ALLOCATE(vvprt(vnx,vny,vnz),STAT=istatus)
vvprt=0
ALLOCATE(vwprt(vnx,vny,vnz),STAT=istatus)
vwprt=0
ALLOCATE(vptprt(vnx,vny,vnz),STAT=istatus)
vptprt=0
ALLOCATE(vpprt(vnx,vny,vnz),STAT=istatus)
vpprt=0
ALLOCATE(vqvprt(vnx,vny,vnz),STAT=istatus)
vqvprt=0
ALLOCATE(vqc(vnx,vny,vnz),STAT=istatus)
vqc=0
ALLOCATE(vqr(vnx,vny,vnz),STAT=istatus)
vqr=0
ALLOCATE(vqi(vnx,vny,vnz),STAT=istatus)
vqi=0
ALLOCATE(vqs(vnx,vny,vnz),STAT=istatus)
vqs=0
ALLOCATE(vqh(vnx,vny,vnz),STAT=istatus)
vqh=0
ALLOCATE(vtke(nx,ny,nz),STAT=istatus)
vtke=0
ALLOCATE(vkmh(nx,ny,nz),STAT=istatus)
vkmh=0
ALLOCATE(vkmv(nx,ny,nz),STAT=istatus)
vkmv=0
ALLOCATE(vubar(vnx,vny,vnz),STAT=istatus)
vubar=0
ALLOCATE(vvbar(vnx,vny,vnz),STAT=istatus)
vvbar=0
ALLOCATE(vwbar(vnx,vny,vnz),STAT=istatus)
vwbar=0
ALLOCATE(vptbar(vnx,vny,vnz),STAT=istatus)
vptbar=0
ALLOCATE(vpbar(vnx,vny,vnz),STAT=istatus)
vpbar=0
ALLOCATE(vrhobar(vnx,vny,vnz),STAT=istatus)
vrhobar=0
ALLOCATE(vqvbar(vnx,vny,vnz),STAT=istatus)
vqvbar=0
ALLOCATE(vsoiltyp(vnx,vny,nstyps),STAT=istatus)
vsoiltyp=0
ALLOCATE(vstypfrct(vnx,vny,nstyps),STAT=istatus)
vstypfrct=0
ALLOCATE(vvegtyp(vnx,vny),STAT=istatus)
vvegtyp=0
ALLOCATE(vlai(vnx,vny),STAT=istatus)
vlai=0
ALLOCATE(vroufns(vnx,vny),STAT=istatus)
vroufns=0
ALLOCATE(vveg(vnx,vny),STAT=istatus)
vveg=0
ALLOCATE(vtsoil(vnx,vny,vnzsoil,0:nstyps),STAT=istatus)
vtsoil=0
ALLOCATE(vqsoil(vnx,vny,vnzsoil,0:nstyps),STAT=istatus)
vqsoil=0
ALLOCATE(vwetcanp(vnx,vny,0:nstyps),STAT=istatus)
vwetcanp=0
ALLOCATE(vsnowdpth(nx,ny),STAT=istatus)
vsnowdpth=0
ALLOCATE(vraing(vnx,vny),STAT=istatus)
vraing=0
ALLOCATE(vrainc(vnx,vny),STAT=istatus)
vrainc=0
ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
prcrate=0
ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
radfrc=0
ALLOCATE(radsw(nx,ny),STAT=istatus)
radsw=0
ALLOCATE(rnflx(nx,ny),STAT=istatus)
rnflx=0
ALLOCATE(radswnet(nx,ny),STAT=istatus)
radswnet=0
ALLOCATE(radlwin(nx,ny),STAT=istatus)
radlwin=0
ALLOCATE(usflx(nx,ny),STAT=istatus)
usflx=0
ALLOCATE(vsflx(nx,ny),STAT=istatus)
vsflx=0
ALLOCATE(ptsflx(nx,ny),STAT=istatus)
ptsflx=0
ALLOCATE(qvsflx(nx,ny),STAT=istatus)
qvsflx=0
ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
tem1=0
ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
tem2=0
ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
tem3=0
ALLOCATE(tem3dsoil(nx,ny,nzsoil),STAT=istatus)
tem3dsoil=0
ALLOCATE(vtem1(vnx,vny,vnz),STAT=istatus)
vtem1=0
ALLOCATE(vtem2(vnx,vny,vnz),STAT=istatus)
vtem2=0
ALLOCATE(vtem3(vnx,vny,vnz),STAT=istatus)
vtem3=0
ALLOCATE(vtem3dsoil(vnx,vny,vnzsoil),STAT=istatus)
vtem3dsoil=0
ALLOCATE(xs(nx),STAT=istatus)
xs=0
ALLOCATE(ys(ny),STAT=istatus)
ys=0
ALLOCATE(zps(nx,ny,nz),STAT=istatus)
zps=0
ALLOCATE(x2d(nx,ny),STAT=istatus)
x2d=0
ALLOCATE(y2d(nx,ny),STAT=istatus)
y2d=0
ALLOCATE(lat(nx,ny),STAT=istatus)
ALLOCATE(lon(nx,ny),STAT=istatus)
lat=0
lon=0
ALLOCATE(vxs(nx),STAT=istatus)
vxs=0
ALLOCATE(vys(ny),STAT=istatus)
vys=0
ALLOCATE(vzps(vnx,vny,vnz),STAT=istatus)
vzps=0
ALLOCATE(dxfld(vnx),STAT=istatus)
dxfld=0
ALLOCATE(dyfld(vny),STAT=istatus)
dyfld=0
ALLOCATE(rdxfld(vnx),STAT=istatus)
rdxfld=0
ALLOCATE(rdyfld(vny),STAT=istatus)
rdyfld=0
ALLOCATE(iloc(nx,ny),STAT=istatus)
ALLOCATE(jloc(nx,ny),STAT=istatus)
iloc=0
jloc=0
ALLOCATE(zpver(nx,ny,vnz),STAT=istatus)
zpver=0
101 CONTINUE
!
!-----------------------------------------------------------------------
!
! Get the name of the grid/base data set.
!
!-----------------------------------------------------------------------
!
DO i=1,3
WRITE(6,'(/a,i5,a)')' For data set ',i,':'
lengbf(i) = len_trim(grdbasfn(i))
! lengbf(i)=LEN(grdbasfn(i))
! CALL strlnth( grdbasfn(i), lengbf(i))
WRITE(6,'(/a,a)')' The grid/base name is ', &
grdbasfn(i)(1:lengbf(i))
lenfil(i) = len_trim(filename(i))
! lenfil(i) = LEN(filename(i))
! CALL strlnth( filename(i), lenfil(i))
WRITE(6,'(/a,/1x,a)')' The data set name is ', &
filename(i)(1:lenfil(i))
END DO
!
WRITE(6,'(/5x,a,a)') 'The output run name is: ', runnmin
WRITE(6,'(a,i5)')' The output data format flag is: ',hdmpfmt
!-----------------------------------------------------------------------
!
! Read all input data arrays (a - forecast data)
!
!-----------------------------------------------------------------------
!
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt,nch,grdbasfn(1)(1:lengbf(1)),lengbf(1), &
filename(1)(1:lenfil(1)),lenfil(1),time, &
x,y,z,zp,zpsoil, uprt ,vprt ,wprt ,ptprt, pprt , &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1,tem2,tem3)
IF (isread == 1) THEN
CALL readsoil(nx,ny,nzsoil,nstyps,soilinfl,dx,dy,zpsoil, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
zpsoilin,tsoilin,qsoilin,wcanpin,snowdin, &
tsoil,qsoil,wetcanp,snowdpth,soiltyp)
END IF
!
!-----------------------------------------------------------------------
!
! ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
IF( ireturn /= 0 ) THEN ! failed read
PRINT *,'Problem reading forecast data, forced to STOP'
STOP
END IF ! successful read
IF (inibred == 1) THEN
utot=0.0
vtot=0.0
wtot=0.0
pttot=0.0
qvtot=0.0
ptot=0.0
utot2=0.0
vtot2=0.0
wtot2=0.0
pttot2=0.0
qvtot2=0.0
ptot2=0.0
do k=1,nz
do i=1,nx
do j=1,ny
utot=utot+uprt(i,j,k)
vtot=vtot+vprt(i,j,k)
wtot=wtot+wprt(i,j,k)
pttot=pttot+ptprt(i,j,k)
qvtot=qvtot+qvprt(i,j,k)
ptot=ptot+pprt(i,j,k)
utot2=utot2+uprt(i,j,k)*uprt(i,j,k)
vtot2=vtot2+vprt(i,j,k)*vprt(i,j,k)
wtot2=wtot2+wprt(i,j,k)*wprt(i,j,k)
pttot2=pttot2+ptprt(i,j,k)*ptprt(i,j,k)
qvtot2=qvtot2+qvprt(i,j,k)*qvprt(i,j,k)
ptot2=ptot2+pprt(i,j,k)*pprt(i,j,k)
enddo
enddo
enddo
utot=utot/totalpoint
vtot=utot/totalpoint
wtot=wtot/totalpoint
pttot=pttot/totalpoint
qvtot=qvtot/totalpoint
ptot=ptot/totalpoint
usd=sqrt(utot2/totalpoint-utot*utot)
vsd=sqrt(vtot2/totalpoint-vtot*vtot)
wsd=sqrt(wtot2/totalpoint-wtot*wtot)
ptsd=sqrt(pttot2/totalpoint-pttot*pttot)
qvsd=sqrt(qvtot2/totalpoint-qvtot*qvtot)
psd=sqrt(ptot2/totalpoint-ptot*ptot)
print *,'totalpoint=',totalpoint
print *,'usd=',usd,' utot=',utot,' utot2=',utot2/totalpoint
print *,'vsd=',vsd,' vtot=',vtot,' vtot2=',vtot2/totalpoint
print *,'wsd=',wsd,' wtot=',wtot,' wtot2=',wtot2/totalpoint
print *,'ptsd=',ptsd,' pttot=',pttot,' pttot2=',pttot2/totalpoint
print *,'qvsd=',qvsd,' qvtot=',qvtot,' qvtot2=',qvtot2/totalpoint
print *,'psd=',psd,' ptot=',ptot,' ptot2=',ptot2/totalpoint
END IF
curtim=time
fcrnam=runname
ifproj=mapproj
fscale=sclfct
flatnot(1)=trulat1
flatnot(2)=trulat2
ftrulon=trulon
fdx=x(3)-x(2)
fdy=y(3)-y(2)
fctrlat=ctrlat
fctrlon=ctrlon
CALL setmapr(ifproj,fscale,flatnot,ftrulon)
CALL lltoxy(1,1,fctrlat,fctrlon,xctr,yctr)
fx0=xctr-fdx*((nx-3)/2)
fy0=yctr-fdy*((ny-3)/2)
CALL setorig(1,fx0,fy0)
!
!-----------------------------------------------------------------------
!
! Get verification data (b - analysis (or gridded obs.) data)
!
!-----------------------------------------------------------------------
!
! Set the gridread parameter to 0 so that the verification
! grid/base file will be read.
!
!-----------------------------------------------------------------------
!
CALL setgbrd (0)
!
!-----------------------------------------------------------------------
!
! Read in the verification data.
!
!-----------------------------------------------------------------------
!
CALL dtaread(vnx,vny,vnz,vnzsoil,nstyps, &
hinfmt,nch,grdbasfn(2)(1:lengbf(2)),lengbf(2), &
filename(2)(1:lenfil(2)),lenfil(2),time, &
vx,vy,vz,vzp,vzpsoil,vuprt,vvprt,vwprt,vptprt,vpprt, &
vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, &
vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, &
vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, &
vtsoil,vqsoil,vwetcanp,vsnowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, vtem1,vtem2,vtem3)
IF (isread == 1) THEN
CALL readsoil(vnx,vny,vnzsoil,nstyps,soilinfl,dx,dy,vzpsoil, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
zpsoilin,tsoilin,qsoilin,wcanpin,snowdin, &
vtsoil,vqsoil,vwetcanp,vsnowdpth,vsoiltyp)
END IF
IF( ireturn /= 0 ) THEN ! failed read
PRINT *,'Problem reading verification data, forced to STOP'
STOP
END IF ! successful read
! IF (vnx /= nx.OR.vny /= ny.OR.vnz /= nz) THEN
! PRINT *,'nx,ny,nz,','=/=','vnx,vny,vnz'
! PRINT *,nx,ny,nz,'=/=',vnx,vny,vnz
! PRINT *, 'forced to stop'
! STOP
! END IF ! don't need
ivproj=mapproj
vscale=sclfct
vlatnot(1)=trulat1
vlatnot(2)=trulat2
vtrulon=trulon
vdx=vx(3)-vx(2)
vdy=vy(3)-vy(2)
vctrlat=ctrlat
vctrlon=ctrlon
CALL setmapr(ivproj,vscale,vlatnot,vtrulon)
CALL lltoxy(1,1,vctrlat,vctrlon,xctr,yctr)
vx0=xctr-vdx*((vnx-3)/2)
vy0=yctr-vdy*((vny-3)/2)
CALL setorig(1,vx0,vy0)
IF (fx0 == vx0.AND.fy0 == vy0.AND. &
flatnot(1) == vlatnot(1).AND.flatnot(2) == vlatnot(2).AND. &
ftrulon == vtrulon.AND.ifproj == ivproj .AND. &
fscale == vscale ) THEN
PRINT *, 'Grids 1 and 2 shares a common coordinate system'
ELSE
PRINT *, 'Grids 1/2 are different, CHECK the PROGRAM or data'
PRINT *, 'Forced to STOP'
STOP
END IF
IF (inibred == 0) THEN
!
!-----------------------------------------------------------------------
!
! Find difference = forecast - verification
! a=a-b (note the forth parameter is 0)
! To reduce memory requirements, the difference fields are
! written to the same arrays as the interpolated fields.
!
!-----------------------------------------------------------------------
!
CALL prtfield(nx,ny,nz,nzsoil,0, &
uprt, vprt, wprt, ptprt, pprt, &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
tsoil,qsoil,wetcanp, &
raing,rainc, &
vuprt, vvprt, vwprt, vptprt, vpprt, &
vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, &
vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, &
vtsoil,vqsoil,vwetcanp, &
vraing,vrainc, &
uprt, vprt, wprt, ptprt, pprt, &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
tsoil,qsoil,wetcanp, &
raing,rainc, &
tem1,tem3dsoil,ireturn )
IF (iensopt == 1) THEN ! scaling the perturbation in breeding IC
utot2=0.0
vtot2=0.0
pttot2=0.0
qvtot2=0.0
ptot2=0.0
do k=1,nz
do i=1,nx
do j=1,ny
utot2=utot2+uprt(i,j,k)*uprt(i,j,k)
vtot2=vtot2+vprt(i,j,k)*vprt(i,j,k)
pttot2=pttot2+ptprt(i,j,k)*ptprt(i,j,k)
qvtot2=qvtot2+qvprt(i,j,k)*qvprt(i,j,k)
ptot2=ptot2+pprt(i,j,k)*pprt(i,j,k)
enddo
enddo
enddo
usd=sqrt(utot2/totalpoint)
vsd=sqrt(vtot2/totalpoint)
ptsd=sqrt(pttot2/totalpoint)
qvsd=sqrt(qvtot2/totalpoint)
psd=sqrt(ptot2/totalpoint)
if(usd > 0.0) uscl = min(0.20*10.0/usd, 1.0)
vscl = uscl
if(ptsd > 0.0) ptscl = min(0.20*11.0/ptsd, 1.0)
if(qvsd > 0.0) qvscl = min(0.20*2e-3/qvsd, 1.0)
if(psd > 0.0) pscl = min(0.20*1000.0/psd, 1.0)
print *,'totalpoint=',totalpoint
print *,'usd=',usd,' uscl=',uscl
print *,'vsd=',vsd,' vscl=',vscl
print *,'ptsd=',ptsd,' ptscl=',ptscl
print *,'qvsd=',qvsd,' qvscl=',qvscl
print *,'psd=',psd,' pscl=',pscl
uprt = uscl*uprt
vprt = vscl*vprt
wprt = 0.0
ptprt = ptscl*ptprt
qvprt = qvscl*qvprt
pprt = pscl*pprt
END IF
ELSE IF (inibred ==1) THEN ! generate random perturbations
do k=1,nz-1
do i=1,nx
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
uprt(i,j,k) = 2.*0.20*usd*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny
iseed=mod(iseed*7141+54773,259200)
vprt(i,j,k) = 2.*0.20*vsd*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
ptprt(i,j,k) = 2.*0.20*ptsd*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
qvprt(i,j,k) = 2.*0.20*qvsd*(iseed/259199.-.5)
enddo
enddo
enddo
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
iseed=mod(iseed*7141+54773,259200)
pprt(i,j,k) = 2.*0.20*psd*(iseed/259199.-.5)
enddo
enddo
enddo
wprt=0.0
qc=0.0
qr=0.0
qi=0.0
qs=0.0
qh=0.0
tke=0.0
kmh=0.0
kmv=0.0
tsoil=0.0
qsoil=0.0
wetcanp=0.0
raing=0.0
rainc=0.0
END IF
!
!-----------------------------------------------------------------------
!
! Get the name of the CONTROL data set, field c or c=b (if iread=0)
! (note: It is needed only if control .ne. verification)
!
!-----------------------------------------------------------------------
!
IF (iread /= 0) THEN
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
CALL setgbrd (0)
CALL dtaread(vnx,vny,vnz,vnzsoil,nstyps, &
hinfmt,nch,grdbasfn(3)(1:lengbf(3)),lengbf(3), &
filename(3)(1:lenfil(3)),lenfil(3),time, &
vx,vy,vz,vzp,vzpsoil,vuprt,vvprt,vwprt,vptprt,vpprt, &
vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, &
vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, &
vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, &
vtsoil,vqsoil,vwetcanp, vsnowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, vtem1,vtem2,vtem3)
IF (isread == 1) THEN
CALL readsoil(nx,ny,nzsoil,nstyps,soilinfl,dx,dy,vzpsoil, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
zpsoilin,tsoilin,qsoilin,wcanpin,snowdin, &
vtsoil, vqsoil, vwetcanp,vsnowdpth,vsoiltyp)
END IF
IF( ireturn /= 0 ) THEN ! failed read
PRINT *,'Problem reading control data, forced to STOP'
STOP
END IF ! successful read
! IF (vnx /= nx.OR.vny /= ny.OR.vnz /= nz) THEN
! PRINT *,'nx,ny,nz','=/=','vnx,vny,vnz'
! PRINT *,nx,ny,nz,'=/=',vnx,vny,vnz
! PRINT *, 'forced to stop'
! STOP
! END IF
ivproj=mapproj
vscale=sclfct
vlatnot(1)=trulat1
vlatnot(2)=trulat2
vtrulon=trulon
vdx=vx(3)-vx(2)
vdy=vy(3)-vy(2)
vctrlat=ctrlat
vctrlon=ctrlon
CALL setmapr(ivproj,vscale,vlatnot,vtrulon)
CALL lltoxy(1,1,vctrlat,vctrlon,xctr,yctr)
vx0=xctr-vdx*((vnx-3)/2)
vy0=yctr-vdy*((vny-3)/2)
CALL setorig(1,vx0,vy0)
!
IF (fx0 == vx0.AND.fy0 == vy0.AND. &
flatnot(1) == vlatnot(1).AND.flatnot(2) == vlatnot(2).AND. &
ftrulon == vtrulon.AND.ifproj == ivproj .AND. &
fscale == vscale ) THEN
PRINT *, 'Grids 2 and 3 shares a common coordinate system'
ELSE
PRINT *, 'Grids 2/3 are different, CHECK the PROFRAM or data'
PRINT *, 'Forced to STOP'
STOP
END IF
END IF
!-----------------------------------------------------------------------
!
! Set output variables to forecast coordinates
!
!-----------------------------------------------------------------------
!
curtim=time
mapproj=ivproj
sclfct=vscale
trulat1=vlatnot(1)
trulat2=vlatnot(2)
trulon=vtrulon
ctrlat=vctrlat
ctrlon=vctrlon
!
!-----------------------------------------------------------------------
!
! Find ptb field = filed c + a/n (note iorder=/=0)
! (b=b+a/n in the code)
! To reduce memory requirements, the resulted fields are
! written to the same arrays as the control fields.
!
!-----------------------------------------------------------------------
!
IF (iorder /= 0) THEN
CALL prtfield(nx,ny,nz,nzsoil,iorder, &
vuprt, vvprt, vwprt, vptprt, vpprt, &
vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, &
vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, &
vtsoil,vqsoil,vwetcanp, &
vraing,vrainc, &
uprt, vprt, wprt, ptprt, pprt, &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
tsoil,qsoil,wetcanp, &
raing,rainc, &
vuprt, vvprt, vwprt, vptprt, vpprt, &
vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, &
vtsoil,vqsoil,vwetcanp, &
vraing,vrainc, &
vtem1,vtem3dsoil,ireturn )
END IF
!-----------------------------------------------------------------------
! Assign the average of soil variables to every soil type
!-----------------------------------------------------------------------
CALL dhslini(nx,ny,nzsoil,nstyps, &
tsoil,qsoil,wetcanp)
!
!-----------------------------------------------------------------------
!
! Get output info
!
!-----------------------------------------------------------------------
!
runname=runnmin
!
!-----------------------------------------------------------------------
!
! Find out the number of characters to be used to construct file
! names.
!
!-----------------------------------------------------------------------
!
CALL gtlfnkey( runname, lfnkey )
!
!-----------------------------------------------------------------------
!
! Find out the number of characters to be used to construct file
! names.
!
!-----------------------------------------------------------------------
!
CALL gtlfnkey( runname, lfnkey )
!
IF (hdmpfmt == 10.AND.grbpkbit == 0) THEN
grbpkbit=16
END IF
!
!-----------------------------------------------------------------------
!
! Set control parameters for
! grid, base state, moisture, and ice variable dumping.
!
!-----------------------------------------------------------------------
!
varout=1
totout = totin
grdout = grdin
basout = basin
varout = varin
mstout = mstin
rainout = rainin
prcout = prcin
iceout = icein
tkeout = tkein
trbout = trbin
sfcout = sfcin
snowout = snowin
landout = landin
radout = radin
flxout = flxin
CALL gtbasfn(runname(1:lfnkey),'./',2,hdmpfmt,mgrid,nestgrd, &
grdbasfn(1), lengbf(1))
WRITE(6,'(/1x,a,a)') &
'Output grid/base state file is ', grdbasfn(1)(1:lengbf(1))
nchdmp = 80
grdbas = 1 ! Dump out grd and base state arrays only
DO k=1,nz
DO j=1,ny
DO i=1,nx
vuprt(i,j,k)=vubar(i,j,k)+vuprt(i,j,k)
vvprt(i,j,k)=vvbar(i,j,k)+vvprt(i,j,k)
vwprt(i,j,k)=vwbar(i,j,k)+vwprt(i,j,k)
vqvprt(i,j,k)=vqvbar(i,j,k)+vqvprt(i,j,k)
END DO
END DO
END DO
IF (iorder /= 0) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
IF (k /= 1) THEN
IF ((vpprt(i,j,k)+vpbar(i,j,k)) >= &
(vpprt(i,j,k-1)+vpbar(i,j,k-1)) ) THEN
vpprt(i,j,k)=vpprt(i,j,k-1)+vpbar(i,j,k-1)-vpbar(i,j,k) &
-(vpbar(i,j,k-1)-vpbar(i,j,k))*0.001
END IF
END IF
IF (vqvprt(i,j,k) <= 1.0E-8) THEN
vqvprt(i,j,k)=1.0E-8
END IF
END DO
END DO
END DO
END IF
CALL dtadump(nx,ny,nz,nzsoil,nstyps, &
hdmpfmt,nchdmp,grdbasfn (1)(1:lengbf(1)),grdbas,filcmprs, &
vuprt,vvprt,vwprt,vptprt,vpprt, &
vqvprt,vqc,vqr,vqi,vqs,vqh,vtke,vkmh,vkmv, &
vubar,vvbar,vwbar,vptbar,vpbar,vrhobar,vqvbar, &
vx,vy,vz,vzp,vzpsoil, &
vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, &
vtsoil,vqsoil,vwetcanp,vsnowdpth, &
vraing,vrainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
vtem1,vtem2,vtem3)
!
!-----------------------------------------------------------------------
!
! Find a unique name hdmpfn(1:ldmpf) for history dump data set
! at time 'curtim'.
!
!-----------------------------------------------------------------------
!
CALL gtdmpfn(runname(1:lfnkey),'./',2, &
curtim,hdmpfmt, &
mgrid,nestgrd, hdmpfn, ldmpf)
WRITE(6,'(/1x,a,f10.0,a,a)') &
'Output file at time ',curtim,' (s) is ', hdmpfn(1:ldmpf)
grdbas = 0 ! Not just dump out grd and base state arrays only
CALL dtadump(nx,ny,nz,nzsoil,nstyps, &
hdmpfmt,nchdmp,hdmpfn(1:ldmpf),grdbas,filcmprs, &
vuprt,vvprt,vwprt,vptprt,vpprt, &
vqvprt,vqc,vqr,vqi,vqs,vqh,vtke,vkmh,vkmv, &
vubar,vvbar,vwbar,vptbar,vpbar,vrhobar,vqvbar, &
vx,vy,vz,vzp,vzpsoil, &
vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, &
vtsoil,qsoil,vwetcanp,vsnowdpth, &
vraing,vrainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
vtem1,vtem2,vtem3)
STOP
END PROGRAM arpsensic
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PRTFIELD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE prtfield(nx,ny,nz,nzsoil,nscale, & 2,19
uprt, vprt, wprt, ptprt, pprt, &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
tsoil,qsoil,wetcanp, &
raing,rainc, &
auprt, avprt, awprt, aptprt, apprt, &
aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv, &
aubar, avbar, awbar, aptbar, apbar, arhobar, aqvbar, &
atsoil,aqsoil,awetcanp, &
araing,arainc, &
duprt, dvprt, dwprt, dptprt, dpprt, &
dqvprt, dqc, dqr, dqi, dqs, dqh, dtke,dkmh,dkmv, &
dtsoil,dqsoil,dwetcanp, &
draing,drainc, &
tem1,tem3dsoil,ireturn)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! (1) nscale=0, Subtract the forecast fields from the verification
! fields (names beginning with "a") and output to the difference
! fields (names beginning with "d"). d=( )bar+( )-abar-a
! (2) nscale=/=0. Generate perturbation fields d=( )+a/nscale with
! the bar arraies being neglected.
! The input ( ( ) ) and output ( d )
! fields may share the same storage location. For this subroutine
! it is assumed the forecast and corresponding verification
! data are at the same physical location, however, the physical
! location may differ between variables. That is uprt and auprt
! are at the same location, but that may differ from pprt and apprt.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster Ou School of Meteorology. April 1992
!
! MODIFICATION HISTORY:
! 14 May 1992 (KB) changed from arps2.5 to arps3.0
! 03 Aug 1992 (KB) updated to account for changes in arps3.0
!
! 09/07/1995 (KB)
! Added differencing of surface (soil) fields.
!
! 15/05/1998 (DH)
! converted from the difference scheme to the current multi-purpose
! version used in ARPSENS group.
!
! 15/12/1998 (DH)
! Added the 3-Dimensionity of surface (soil) fields.
!
! 05/28/2002 (J. Brotzge)
! Added tsoil/qsoil to accomodate new soil schemes.
!
!-----------------------------------------------------------------------
!
! INPUT :
! nx,ny,nz,nzsoil Array dimensions for forecast field.
!
! FORECAST FIELDS:
!
! uprt perturbation x component of velocity (m/s)
! vprt perturbation y component of velocity (m/s)
! wprt perturbation vertical component of velocity in Cartesian
! coordinates (m/s).
!
! ptprt perturbation potential temperature (K)
! pprt perturbation pressure (Pascal)
!
! qvprt perturbation water vapor mixing ratio (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! wbar Base state z velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhobar Base state density (kg/m**3)
! qvbar Base state water vapor mixing ratio (kg/kg)
! tke Turbulent Kinetic Energy ((m/s)**2)
! kmh Horizontal turbulent mixing coefficient (m**2/s)
! kmv Vertical turbulent mixing coefficient (m**2/s)
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m**3/m**3)
! wetcanp Canopy water amount
!
! raing Grid supersaturation rain
! rainc Cumulus convective rain
!
! INTERPOLATED VERIFICATION FIELDS:
!
! auprt perturbation x component of velocity (m/s)
! avprt perturbation y component of velocity (m/s)
! awprt perturbation vertical component of velocity in Cartesian
! coordinates (m/s).
!
! aptprt perturbation potential temperature (K)
! apprt perturbation pressure (Pascal)
!
! aqvprt perturbation water vapor mixing ratio (kg/kg)
! aqc Cloud water mixing ratio (kg/kg)
! aqr Rainwater mixing ratio (kg/kg)
! aqi Cloud ice mixing ratio (kg/kg)
! aqs Snow mixing ratio (kg/kg)
! aqh Hail mixing ratio (kg/kg)
!
! aubar Base state x velocity component (m/s)
! avbar Base state y velocity component (m/s)
! awbar Base state z velocity component (m/s)
! aptbar Base state potential temperature (K)
! apbar Base state pressure (Pascal)
! arhobar Base state density (kg/m**3)
! aqvbar Base state water vapor mixing ratio (kg/kg)
!
! atsoil Soil temperature (K)
! aqsoil Soil moisture (m**3/m**3)
! awetcanp Canopy water amount
!
! araing Grid supersaturation rain
! arainc Cumulus convective rain
!
! OUTPUT :
!
! DIFFERENCE FIELDS (may share storage with forecast fields
! or interpolated fields in calling program):
!
! duprt perturbation x component of velocity (m/s)
! dvprt perturbation y component of velocity (m/s)
! dwprt perturbation vertical component of velocity in Cartesian
! coordinates (m/s).
!
! dptprt perturbation potential temperature (K)
! dpprt perturbation pressure (Pascal)
!
! dqvprt perturbation water vapor mixing ratio (kg/kg)
! dqc Cloud water mixing ratio (kg/kg)
! dqr Rainwater mixing ratio (kg/kg)
! dqi Cloud ice mixing ratio (kg/kg)
! dqs Snow mixing ratio (kg/kg)
! dqh Hail mixing ratio (kg/kg)
!
! dtsoil Soil temperature (K)
! dqsoil Soil moisture (m**3/m**3)
! dwetcanp Canopy water amount
!
! draing Grid supersaturation rain
! drainc Cumulus convective rain
!
! tem1 Work array
! tem3dsoil Work array
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nzsoil ! 4 dimensions of array
INTEGER :: nscale
!
!-----------------------------------------------------------------------
!
! Model Arrays
!
!-----------------------------------------------------------------------
!
REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s)
REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s)
REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor specific humidity
REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2)
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 :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: rhobar (nx,ny,nz) ! Base state air density (kg/m**3)
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity
REAL :: tsoil (nx,ny,nzsoil) ! Deep soil temperature (K)
REAL :: qsoil (nx,ny,nzsoil) ! Deep soil moisture
REAL :: wetcanp(nx,ny) ! Canopy water amount
REAL :: raing (nx,ny) ! Grid supersaturation rain
REAL :: rainc (nx,ny) ! Cumulus convective rain
!
!-----------------------------------------------------------------------
!
! Verification data interpolated to model grid
!
!-----------------------------------------------------------------------
!
REAL :: auprt (nx,ny,nz) ! Perturbation u-velocity (m/s)
REAL :: avprt (nx,ny,nz) ! Perturbation v-velocity (m/s)
REAL :: awprt (nx,ny,nz) ! Perturbation w-velocity (m/s)
REAL :: aptprt (nx,ny,nz) ! Perturbation potential temperature (K)
REAL :: apprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: aqvprt (nx,ny,nz) ! Perturbation water vapor specific humidity
REAL :: aqc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: aqr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: aqi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: aqs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: aqh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: atke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2)
REAL :: akmh (nx,ny,nz) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: akmv (nx,ny,nz) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: aubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: avbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: awbar (nx,ny,nz) ! Base state w-velocity (m/s)
REAL :: aptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: arhobar(nx,ny,nz) ! Base state density (kg/m**3)
REAL :: apbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: aqvbar (nx,ny,nz) ! Base state water vapor specific humidity
REAL :: atsoil (nx,ny,nzsoil) ! Soil temperature (K)
REAL :: aqsoil (nx,ny,nzsoil) ! Soil moisture (m**3/m**3)
REAL :: awetcanp(nx,ny) ! Canopy water amount
REAL :: araing (nx,ny) ! Grid supersaturation rain
REAL :: arainc (nx,ny) ! Cumulus convective rain
!
!------------------------------------------------------------------------
!
! Difference arrays
!
!-----------------------------------------------------------------------
!
REAL :: duprt (nx,ny,nz) ! perturbation x component of velocity (m/s)
REAL :: dvprt (nx,ny,nz) ! perturbation y component of velocity (m/s)
REAL :: dwprt (nx,ny,nz) ! perturbation vertical component of
! velocity in Cartesian coordinates (m/s)
REAL :: dptprt (nx,ny,nz) ! perturbation potential temperature (K)
REAL :: dpprt (nx,ny,nz) ! perturbation pressure (Pascal)
REAL :: dqvprt (nx,ny,nz) ! perturbation water vapor mixing ratio (kg/kg)
REAL :: dqc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: dqr (nx,ny,nz) ! Rainwater mixing ratio (kg/kg)
REAL :: dqi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: dqs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: dqh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: dtke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2)
REAL :: dkmh (nx,ny,nz) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: dkmv (nx,ny,nz) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL :: dtsoil (nx,ny,nzsoil) ! Soil temperature (K)
REAL :: dqsoil (nx,ny,nzsoil) ! Soil moisture (m**3/m**3)
REAL :: dwetcanp(nx,ny) ! Canopy water amount
REAL :: draing (nx,ny) ! Grid supersaturation rain
REAL :: drainc (nx,ny) ! Cumulus convective rain
REAL :: tem1 (nx,ny,nz) ! A work array
REAL :: tem3dsoil(nx,ny,nzsoil) ! A work array
INTEGER :: ireturn, i,j,k
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: is,js,ks,ls,ie,je,ke,le
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
is=1
js=1
ks=1
ls=1
ie=nx-1
je=ny-1
ke=nz-1
le=nzsoil-1
!-----------------------------------------------------------------------
!
! Scalars
!
!-----------------------------------------------------------------------
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k)=0.0
END DO
END DO
END DO
tem3dsoil = 0.0
PRINT *, ' ptprt: '
CALL subtrprt(nx,ny,nz, ptprt,ptbar,aptprt,aptbar,dptprt,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' pprt: '
CALL subtrprt(nx,ny,nz, pprt, pbar, apprt, apbar, dpprt,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qvprt: '
CALL subtrprt(nx,ny,nz, qvprt,qvbar,aqvprt,aqvbar,dqvprt,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qc: '
CALL subtrprt(nx,ny,nz, qc, tem1, aqc, tem1, dqc,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qr: '
CALL subtrprt(nx,ny,nz, qr, tem1, aqr, tem1, dqr,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qi: '
CALL subtrprt(nx,ny,nz, qi, tem1, aqi, tem1, dqi,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qs: '
CALL subtrprt(nx,ny,nz, qs, tem1, aqs, tem1, dqs,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' qh: '
CALL subtrprt(nx,ny,nz, qh, tem1, aqh, tem1, dqh,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' tke: '
CALL subtrprt(nx,ny,nz, tke, tem1, atke, tem1, dtke,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' kmh: '
CALL subtrprt(nx,ny,nz, kmh, tem1, akmh, tem1, dkmh,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' kmv: '
CALL subtrprt(nx,ny,nz, kmv, tem1, akmv, tem1, dkmv,nscale, &
is,js,ks,ie,je,ke)
!-----------------------------------------------------------------------
!
! u wind components
!
!-----------------------------------------------------------------------
ie=nx
PRINT *, ' uprt: '
CALL subtrprt(nx,ny,nz,uprt,ubar,auprt,aubar,duprt,nscale, &
is,js,ks,ie,je,ke)
!-----------------------------------------------------------------------
!
! v wind components
!
!-----------------------------------------------------------------------
ie=nx-1
je=ny
PRINT *, ' vprt: '
CALL subtrprt(nx,ny,nz,vprt,vbar,avprt,avbar,dvprt,nscale, &
is,js,ks,ie,je,ke)
!-----------------------------------------------------------------------
!
! w wind components
!
!-----------------------------------------------------------------------
je=ny-1
ke=nz
PRINT *, ' wprt: '
CALL subtrprt(nx,ny,nz,wprt,tem1,awprt,tem1,dwprt,nscale, &
is,js,ks,ie,je,ke)
!-----------------------------------------------------------------------
!
! 2-d/3-d surface (soil) variables
!
!-----------------------------------------------------------------------
ie=nx-1
je=ny-1
le=nzsoil
ks=1
ke=1
le=1
!
PRINT *, ' tsoil:'
CALL subtrprt(nx,ny,nzsoil,tsoil,tem3dsoil,atsoil,tem3dsoil,dtsoil, &
nscale,is,js,ls,ie,je,le)
PRINT *, ' qsoil:'
CALL subtrprt(nx,ny,nzsoil,qsoil,tem3dsoil,aqsoil,tem3dsoil,dqsoil, &
nscale,is,js,ls,ie,je,le)
PRINT *, ' wetcanp:'
CALL subtrprt(nx,ny,1,wetcanp,tem1,awetcanp,tem1,dwetcanp,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' raing:'
CALL subtrprt(nx,ny,1, raing,tem1, araing,tem1, draing,nscale, &
is,js,ks,ie,je,ke)
PRINT *, ' rainc:'
CALL subtrprt(nx,ny,1, rainc,tem1, arainc,tem1, drainc,nscale, &
is,js,ks,ie,je,ke)
!
RETURN
END SUBROUTINE prtfield
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SUBTRPRT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE subtrprt(nx,ny,nz, a,abar,b,bbar,c,nscale, & 19
istr,jstr,kstr,iend,jend,kend)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Subtracts 2 three-dimensional arrays, represented by
! means plus perturbations.
!
! AUTHOR: Keith Brewster OU School of Meteorology. Feb 1992
!
! MODIFICATION HISTORY:
! 11 Aug 1992 (KB) changed from arps2.5 to arps3.0
!
!
!-----------------------------------------------------------------------
!
! INPUT:
! a perturbation data array
! abar mean data array
! b perturbation data array to subtract from a
! bbar mean data array to subtract from a
!
! OUTPUT:
! c difference array a-b
! (may share storage in calling program with array a or b)
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! 3 dimensions of array
INTEGER :: nscale
REAL :: a (nx,ny,nz) ! data array
REAL :: abar(nx,ny,nz) ! base state of data array a
REAL :: b (nx,ny,nz) ! data array to subtract from a
REAL :: bbar(nx,ny,nz) ! base state of data arrya b
REAL :: c (nx,ny,nz) ! difference array a-b
INTEGER :: istr,jstr,kstr
INTEGER :: iend,jend,kend
INTEGER :: i,j,k,imid,jmid,kmid
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
imid=nint(0.5*FLOAT(istr+iend))
jmid=nint(0.5*FLOAT(jstr+jend))
kmid=nint(0.5*FLOAT(kstr+kend))
IF (nz < 10) kmid=kstr
PRINT *, 'imid,jmid,kmid: ',imid,jmid,kmid
!
!-----------------------------------------------------------------------
!
! Tell us about a sample input point
!
!-----------------------------------------------------------------------
!
IF (nscale == 0) THEN
PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), &
' b= ',(b(imid,jmid,kmid)+bbar(imid,jmid,kmid))
ELSE
PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), &
' b= ',b(imid,jmid,kmid)
END IF
!
!-----------------------------------------------------------------------
!
! Subtraction
!
!-----------------------------------------------------------------------
!
DO k=kstr,kend
DO j=jstr,jend
DO i=istr,iend
IF (nscale == 0) THEN
c(i,j,k)=a(i,j,k)+abar(i,j,k)-(b(i,j,k)+bbar(i,j,k))
ELSE
c(i,j,k)=a(i,j,k)+b(i,j,k)/FLOAT(nscale)
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Tell us about a sample output point
!
!-----------------------------------------------------------------------
!
IF (nscale == 0) THEN
PRINT *, ' c= ',c(imid,jmid,kmid)
ELSE
PRINT *, ' c= ',c(imid,jmid,kmid)+abar(imid,jmid,kmid)
END IF
RETURN
END SUBROUTINE subtrprt
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DHSLINI ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE dhslini(nx,ny,nzsoil,nstyps, & 1
tsoil,qsoil,wetcanp)
IMPLICIT NONE
INTEGER :: nx,ny,nzsoil,nstyps
INTEGER :: i,j,k,l
REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3)
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
DO k=1,nstyps
DO j=1,ny
DO i=1,nx
DO l=1,nzsoil
tsoil(i,j,l,k)=tsoil(i,j,l,0)
qsoil(i,j,l,k)=qsoil(i,j,l,0)
END DO
wetcanp(i,j,k)=wetcanp(i,j,0)
END DO
END DO
END DO
RETURN
END SUBROUTINE dhslini