PROGRAM arpsensic,110
!
!##################################################################
!##################################################################
!###### ######
!###### 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 real variable iorder, a real factor
!
! 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
!
! 02/01/2005 (F. KONG)
! Change iorder to real factor (originally interger)
!
! 04/10/2005 (F. KONG)
! Change random perturbation to normal distribution
! with zero mean and specified variance (add subroutine
! normalrand for this purpose). Specify variance via
! arpsensic.input; Add computer assigning seed value
!
! 02/15/2007 (F. KONG)
! Major rewrite in this revision to:
! - Changed to MPI version.
! - Remove domain matching check to streamline the code
! (Always assuming all read in data have the same ARPS
! domain seeting)
!
!-----------------------------------------------------------------------
!
! DATA ARRAYS READ IN:
!
! x x coordinate of grid points in computational space (m)
! y y coordinate of grid points in computational space (m)
! z z coordinate of grid points in computational space (m)
! zp z coordinate of grid points in physical space (m)
! zpsoil z coordinate of grid points in the soil (m)
!
! uprt perturbation x component of velocity (m/s)
! vprt perturbation y component of velocity (m/s)
! wprt perturbation z component of velocity (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
! kmh
! kmv
!
! soiltyp Soil type
! stypfrct
! 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
! snowdpth
!
! raing Grid supersaturation rain
! rainc Cumulus convective rain
!
!
! 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'
INCLUDE 'mp.inc'
INTEGER :: nx,ny,nz,nzsoil
INTEGER :: nstyps ! Maximum number of soil types in each
! grid box
!
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
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 (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 :: 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 (kg/kg)
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil fraction
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(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
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)
!
!-----------------------------------------------------------------------
!
! 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 (kg/kg)
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 (kg/kg)
INTEGER, ALLOCATABLE :: vsoiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: vstypfrct(:,:,:) ! Soil fraction
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 :: 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 :: lat(:,:),lon(:,:)
REAL, ALLOCATABLE :: vxs(:)
REAL, ALLOCATABLE :: vys(:)
REAL, ALLOCATABLE :: vzps(:,:,:)
REAL, ALLOCATABLE :: dxfld(:)
REAL, ALLOCATABLE :: dyfld(:)
REAL, ALLOCATABLE :: rdxfld(:)
REAL, ALLOCATABLE :: rdyfld(:)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: fcrnam,runnmin
CHARACTER (LEN=256) :: filename(3),grdbasfn(3)
CHARACTER (LEN=256) :: fbasfn,filnam
CHARACTER (LEN=40) :: q3dname,q3dunit
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 :: hinfmt,iread,isread,iensopt,inibred,iseed
REAL :: iorder,uvar,vvar,wvar,ptvar,qvvar
INTEGER :: ireturn
LOGICAL :: comcoord
INTEGER :: nchin
INTEGER :: zpsoilin,tsoilin,qsoilin,wcanpin,snowdin
INTEGER :: dmp_out_joined, pertout
!
INTEGER :: istatus
INTEGER :: idate(8)
REAL :: utot,utot2,usd,uscl
REAL :: vtot,vtot2,vsd,vscl
REAL :: wtot,wtot2,wsd,wscl
REAL :: pttot,pttot2,ptsd,ptscl
REAL :: qvtot,qvtot2,qvsd,qvscl
REAL :: ptot,ptot2,psd,pscl
REAL :: totalpoint
REAL :: ampu,ampv,ampw,amppt,ampqv,ampp,ampke,enorm,escl
REAL :: rateu,ratev,ratew,ratept,rateqv,ratep,rateke
INTEGER :: nprocx_in, nprocy_in
NAMELIST /message_passing/nproc_x, nproc_y, readsplit, nprocx_in, nprocy_in
NAMELIST /prtbpara/ iensopt,inibred,iseed,iorder,isread, &
soilinfl,soilfmt,iread, &
uvar,vvar,wvar,ptvar,qvvar
NAMELIST /input_data/ hinfmt,grdbasfn,filename
NAMELIST /outpt_data/ dmp_out_joined,hdmpfmt,runnmin,grdout,basout, &
varout,mstout,iceout,trbout,sfcout, &
rainout,snowout,filcmprs,pertout
!-----------------------------------------------------------------------
!
! Variables for mpi jobs
!
!-----------------------------------------------------------------------
INTEGER, PARAMETER :: fzone = 3
INTEGER :: nxlg, nylg ! global domain
INTEGER :: ii,jj,ia,ja
CHARACTER(LEN=256) :: tmpstr
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Initializations
!
!-----------------------------------------------------------------------
!
CALL mpinit_proc
IF(myproc == 0) WRITE(6,'(/6(/5x,a)/)') &
'###############################################################', &
'# #', &
'# Welcome to ARPSENSIC, a program that reads in history files #', &
'# generated by ARPS and produces perturbation grids variables #', &
'# #', &
'###############################################################'
soilfmt = 1
!------------------------------------------------------------------------
! Read namelist
!------------------------------------------------------------------------
IF(myproc == 0) THEN
READ(5,message_passing)
WRITE(6,'(a)')'Namelist message_passing was successfully read.'
write(6,message_passing)
END IF
CALL mpupdatei(nproc_x,1)
CALL mpupdatei(nproc_y,1)
CALL mpupdatei(readsplit,1)
IF (mp_opt == 0) THEN
nproc_x = 1
nproc_y = 1
nprocx_in = 1
nprocy_in = 1
readsplit = 0
nprocs = 1
max_fopen = 1
ELSE
max_fopen = nproc_x * nproc_y
END IF
CALL mpinit_var
IF(myproc == 0) THEN
READ(5,prtbpara,END=999)
write(6,prtbpara)
END IF
CALL mpupdatei(iensopt,1)
CALL mpupdatei(inibred,1)
CALL mpupdatei(iseed,1)
CALL mpupdater(uvar,1)
CALL mpupdater(vvar,1)
CALL mpupdater(wvar,1)
CALL mpupdater(ptvar,1)
CALL mpupdater(qvvar,1)
CALL mpupdater(iorder,1)
CALL mpupdatec(soilinfl,256)
CALL mpupdatei(isread,1)
CALL mpupdatei(iread,1)
IF(myproc == 0) THEN
READ(5,input_data,END=999)
write(6,input_data)
END IF
CALL mpupdatei(hinfmt,1)
CALL mpupdatec(grdbasfn,3*256)
CALL mpupdatec(filename,3*256)
IF(myproc == 0) THEN
READ(5,outpt_data,END=999)
write(6,outpt_data)
END IF
CALL mpupdatei(dmp_out_joined,1)
CALL mpupdatei(hdmpfmt,1)
CALL mpupdater(runnmin,80)
CALL mpupdatei(grdout,1)
CALL mpupdatei(basout,1)
CALL mpupdatei(varout,1)
CALL mpupdatei(mstout,1)
CALL mpupdatei(iceout,1)
CALL mpupdatei(trbout,1)
CALL mpupdatei(rainout,1)
CALL mpupdatei(sfcout,1)
CALL mpupdatei(snowout,1)
CALL mpupdatei(filcmprs,1)
CALL mpupdatei(pertout,1)
joindmp = dmp_out_joined
IF (mp_opt > 0) THEN ! should moved into mpinit_var later
dumpstride = max_fopen
IF (joindmp > 0) dumpstride = nprocs ! join and dump
END IF
GO TO 1001
999 IF(myproc == 0) PRINT *, 'ERROR in reading from input file, Program Terminated'
STOP
1001 CONTINUE
!------------------------------------------------------------------------
!
! Get the dimensions of the input data files
!
!------------------------------------------------------------------------
lengbf(1) = LEN_TRIM(grdbasfn(1))
IF(mp_opt > 0 .AND. readsplit <= 0) THEN
tmpstr = grdbasfn(1)
WRITE(grdbasfn(1),'(a,a,2i2.2)') tmpstr(1:lengbf(1)),'_',loc_x,loc_y
lengbf(1) = lengbf(1) + 5
END IF
IF(myproc == 0) THEN
PRINT *, 'Get dimension from:',trim(grdbasfn(1))
CALL get_dims_from_data(hinfmt,grdbasfn(1)(1:lengbf(1)), &
nx,ny,nz,nzsoil,nstyps, ireturn)
IF (mp_opt > 0 .AND. readsplit > 0) THEN
!
! Fiddle with nx/ny, which apparently are wrong.
!
nx = (nx - 3) / nproc_x + 3
ny = (ny - 3) / nproc_y + 3
END IF
IF (nstyps <= 0) nstyps = 1
nstyp = nstyps ! Copy to global variabl
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
CALL arpsstop('Problem with data.',1)
END IF
END IF
CALL mpupdatei(nx,1)
CALL mpupdatei(ny,1)
CALL mpupdatei(nz,1)
CALL mpupdatei(nzsoil,1)
CALL mpupdatei(nstyps,1)
CALL mpupdatei(nstyp,1)
IF(myproc == 0) &
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,' nzsoil =',nzsoil
IF(myproc == 0) print*,'nstyps =', nstyps
!
!-----------------------------------------------------------------------
!
! Allocate the variables and initialize the them to zero
!
!-----------------------------------------------------------------------
ALLOCATE(x(nx),STAT=istatus)
x=0
ALLOCATE(y(ny),STAT=istatus)
y=0
ALLOCATE(z(nz),STAT=istatus)
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(nx),STAT=istatus)
vx=0
ALLOCATE(vy(ny),STAT=istatus)
vy=0
ALLOCATE(vz(nz),STAT=istatus)
vz=0
ALLOCATE(vzp(nx,ny,nz),STAT=istatus)
vzp=0
ALLOCATE(vzpsoil(nx,ny,nzsoil),STAT=istatus)
vzpsoil=0
ALLOCATE(vuprt(nx,ny,nz),STAT=istatus)
vuprt=0
ALLOCATE(vvprt(nx,ny,nz),STAT=istatus)
vvprt=0
ALLOCATE(vwprt(nx,ny,nz),STAT=istatus)
vwprt=0
ALLOCATE(vptprt(nx,ny,nz),STAT=istatus)
vptprt=0
ALLOCATE(vpprt(nx,ny,nz),STAT=istatus)
vpprt=0
ALLOCATE(vqvprt(nx,ny,nz),STAT=istatus)
vqvprt=0
ALLOCATE(vqc(nx,ny,nz),STAT=istatus)
vqc=0
ALLOCATE(vqr(nx,ny,nz),STAT=istatus)
vqr=0
ALLOCATE(vqi(nx,ny,nz),STAT=istatus)
vqi=0
ALLOCATE(vqs(nx,ny,nz),STAT=istatus)
vqs=0
ALLOCATE(vqh(nx,ny,nz),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(nx,ny,nz),STAT=istatus)
vubar=0
ALLOCATE(vvbar(nx,ny,nz),STAT=istatus)
vvbar=0
ALLOCATE(vwbar(nx,ny,nz),STAT=istatus)
vwbar=0
ALLOCATE(vptbar(nx,ny,nz),STAT=istatus)
vptbar=0
ALLOCATE(vpbar(nx,ny,nz),STAT=istatus)
vpbar=0
ALLOCATE(vrhobar(nx,ny,nz),STAT=istatus)
vrhobar=0
ALLOCATE(vqvbar(nx,ny,nz),STAT=istatus)
vqvbar=0
ALLOCATE(vsoiltyp(nx,ny,nstyps),STAT=istatus)
vsoiltyp=0
ALLOCATE(vstypfrct(nx,ny,nstyps),STAT=istatus)
vstypfrct=0
ALLOCATE(vvegtyp(nx,ny),STAT=istatus)
vvegtyp=0
ALLOCATE(vlai(nx,ny),STAT=istatus)
vlai=0
ALLOCATE(vroufns(nx,ny),STAT=istatus)
vroufns=0
ALLOCATE(vveg(nx,ny),STAT=istatus)
vveg=0
ALLOCATE(vtsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus)
vtsoil=0
ALLOCATE(vqsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus)
vqsoil=0
ALLOCATE(vwetcanp(nx,ny,0:nstyps),STAT=istatus)
vwetcanp=0
ALLOCATE(vsnowdpth(nx,ny),STAT=istatus)
vsnowdpth=0
ALLOCATE(vraing(nx,ny),STAT=istatus)
vraing=0
ALLOCATE(vrainc(nx,ny),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(nx,ny,nz),STAT=istatus)
vtem1=0
ALLOCATE(vtem2(nx,ny,nz),STAT=istatus)
vtem2=0
ALLOCATE(vtem3(nx,ny,nz),STAT=istatus)
vtem3=0
ALLOCATE(vtem3dsoil(nx,ny,nzsoil),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(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(nx,ny,nz),STAT=istatus)
vzps=0
ALLOCATE(dxfld(nx),STAT=istatus)
dxfld=0
ALLOCATE(dyfld(ny),STAT=istatus)
dyfld=0
ALLOCATE(rdxfld(nx),STAT=istatus)
rdxfld=0
ALLOCATE(rdyfld(ny),STAT=istatus)
rdyfld=0
!
! mpi variables
!
nxlg = (nx-fzone)*nproc_x + fzone
nylg = (ny-fzone)*nproc_y + fzone
totalpoint = (nxlg-3)*(nylg-3)*(nz-1)
IF(myproc == 0) print *, 'totalpoint =',totalpoint
101 CONTINUE
!
!-----------------------------------------------------------------------
!
! Get the name of the grid/base data set.
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) THEN
DO i=1,3
WRITE(6,'(/a,i5,a)') ' For data set ',i,':'
WRITE(6,'(/a,/1x,a)')' The base set name is ',trim(grdbasfn(i))
WRITE(6,'(/a,/1x,a)')' The data set name is ',trim(filename(i))
END DO
WRITE(6,'(/5x,a,a)') 'The output run name is: ', trim(runnmin)
WRITE(6,'(a,i5)') ' The output data format flag is: ',hdmpfmt
END IF
!-----------------------------------------------------------------------
!
! Read all input data arrays (a - forecast data)
!
!-----------------------------------------------------------------------
!
lenfil(1) = LEN_TRIM(filename(1))
IF (mp_opt > 0 .AND. readsplit <= 0) THEN
tmpstr = filename(1)
WRITE(filename(1),'(2a,2i2.2)') tmpstr(1:lenfil(1)),'_',loc_x,loc_y
lenfil(1) = lenfil(1) +5
END IF
IF (myproc == 0) THEN
WRITE(6,'(/2a/)') ' Reading file: ',filename(1)(1:lenfil(1))
WRITE(6,'(/2a/)') ' Reading file: ',grdbasfn(1)(1:lengbf(1))
END IF
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt,nchin,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 ! currently, no MPI function
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
!
IF( ireturn /= 0 ) CALL arpsstop('dtaread errors.',1)
!
IF (inibred == 1) THEN ! initial bred, preparing sd(f)
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-1
do i=2,nx-2
do j=2,ny-2
utot=utot+(ubar(i,j,k)+uprt(i,j,k))
vtot=vtot+(vbar(i,j,k)+vprt(i,j,k))
wtot=wtot+(wbar(i,j,k)+wprt(i,j,k))
pttot=pttot+ptprt(i,j,k)
! pttot=pttot+ptprt(i,j,k)+ptbar(i,j,k)
qvtot=qvtot+(qvbar(i,j,k)+qvprt(i,j,k))
ptot=ptot+pprt(i,j,k)
! ptot=ptot+(pprt(i,j,k)+pbar(i,j,k))
utot2=utot2+(ubar(i,j,k)+uprt(i,j,k))*(ubar(i,j,k)+uprt(i,j,k))
vtot2=vtot2+(vbar(i,j,k)+vprt(i,j,k))*(vbar(i,j,k)+vprt(i,j,k))
wtot2=wtot2+(wbar(i,j,k)+wprt(i,j,k))*(wbar(i,j,k)+wprt(i,j,k))
pttot2=pttot2+(ptprt(i,j,k)*ptprt(i,j,k))
! pttot2=pttot2+(ptbar(i,j,k)+ptprt(i,j,k))*(ptbar(i,j,k)+ptprt(i,j,k))
qvtot2=qvtot2+(qvbar(i,j,k)+qvprt(i,j,k))*(qvbar(i,j,k)+qvprt(i,j,k))
ptot2=ptot2+(pprt(i,j,k)*pprt(i,j,k))
! ptot2=ptot2+(pbar(i,j,k)+pprt(i,j,k))*(pbar(i,j,k)+pprt(i,j,k))
enddo
enddo
enddo
IF (mp_opt > 0) THEN
CALL mpsumr(utot, 1)
CALL mpsumr(vtot, 1)
CALL mpsumr(wtot, 1)
CALL mpsumr(pttot, 1)
CALL mpsumr(qvtot, 1)
CALL mpsumr(ptot, 1)
CALL mpsumr(utot2, 1)
CALL mpsumr(vtot2, 1)
CALL mpsumr(wtot2, 1)
CALL mpsumr(pttot2, 1)
CALL mpsumr(qvtot2, 1)
CALL mpsumr(ptot2, 1)
END IF
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)
IF (myproc == 0) THEN
print *,'totalpoint=',totalpoint
print *,'variance info:'
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
write(12,*) 'utot2,vtot2,wtot2,pttot2,qvtot2,ptot2,enorm - total mean'
write(12,'(7e11.4)') sqrt(utot2/totalpoint),sqrt(vtot2/totalpoint), &
sqrt(wtot2/totalpoint),sqrt(pttot2/totalpoint), &
sqrt(qvtot2/totalpoint),sqrt(ptot2/totalpoint), &
sqrt(0.5*(utot2+vtot2+wtot2)/totalpoint)
endif
! set initial (random) amplitude for each variable
ampu = uvar ! provided from input file
ampv = vvar
ampw = wvar
amppt= ptvar
ampqv= qvvar
ampp = 0.0
IF (myproc == 0) THEN
print *,'Initial amplitudes (for random perturbation):'
print *,'ampu =',ampu
print *,'ampv =',ampv
print *,'ampw =',ampw
print *,'amppt=',amppt
print *,'ampqv=',ampqv
print *,'ampp=',ampp
endif
ELSE
ampu = uvar ! provided from input file
ampv = vvar
ampw = wvar
amppt= ptvar
ampqv= qvvar
ampp = 0.0
IF (myproc == 0) THEN
print *,'Perturbation switches ( >0 on; =0 off):'
print *,'pert_u =',ampu
print *,'pert_v =',ampv
print *,'pert_w =',ampw
print *,'pert_pt=',amppt
print *,'pert_qv=',ampqv
print *,'pert_p =',ampp
endif
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*((nxlg-3)/2)
fy0=yctr-fdy*((nylg-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.
!
!-----------------------------------------------------------------------
!
lenfil(2) = len_trim(filename(2))
lengbf(2) = len_trim(grdbasfn(2))
IF (mp_opt > 0 .AND. readsplit <= 0) THEN
tmpstr = filename(2)
WRITE(filename(2),'(2a,2i2.2)') tmpstr(1:lenfil(2)),'_',loc_x,loc_y
lenfil(2) = lenfil(2) + 5
tmpstr = grdbasfn(2)
WRITE(grdbasfn(2),'(2a,2i2.2)') tmpstr(1:lengbf(2)),'_',loc_x,loc_y
lengbf(2) = lengbf(2) + 5
END IF
IF (myproc == 0) THEN
WRITE(6,'(/2a/)') ' Reading file: ',filename(2)(1:lenfil(2))
WRITE(6,'(/2a/)') ' Reading file: ',grdbasfn(2)(1:lengbf(2))
END IF
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt,nchin,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(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 ) CALL arpsstop('dtaread errors.',1)
! 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
! IF(myproc == 0) PRINT *, 'Grids 1 and 2 shares a common coordinate system'
! ELSE
! IF(myproc == 0) THEN
! PRINT *, 'Grids 1/2 are different, CHECK the PROGRAM or data'
! PRINT *, 'Forced to STOP'
! ENDIF
! 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 )
! write out 3D perturbation fields for plotting purpose
IF(pertout > 0) THEN
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
filnam = filename(1)(1:lenfil(1)-9)//'ptpert'// &
filename(1)(lenfil(1)-5:lenfil(1))
! Here, further action is needed for split file name (KONG)
q3dname='pt_pert'
q3dunit='K'
CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,ptprt, &
1,1,1)
filnam = filename(1)(1:lenfil(1)-9)//'qvpert'// &
filename(1)(lenfil(1)-5:lenfil(1))
q3dname='qv_pert'
q3dunit='g/kg'
CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,1e3*qvprt, &
1,1,1)
filnam = filename(1)(1:lenfil(1)-9)//'u_pert'// &
filename(1)(lenfil(1)-5:lenfil(1))
q3dname='u_pert'
q3dunit='m/s'
CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,uprt, &
0,1,1)
filnam = filename(1)(1:lenfil(1)-9)//'v_pert'// &
filename(1)(lenfil(1)-5:lenfil(1))
q3dname='v_pert'
q3dunit='m/s'
CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,vprt, &
1,0,1)
filnam = filename(1)(1:lenfil(1)-9)//'w_pert'// &
filename(1)(lenfil(1)-5:lenfil(1))
q3dname='w_pe'
CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,wprt, &
1,1,0)
filnam = filename(1)(1:lenfil(1)-9)//'p_pert'// &
filename(1)(lenfil(1)-5:lenfil(1))
q3dname='p_pert'
q3dunit='Pa'
CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,pprt, &
1,1,1)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
END IF ! End writing 3D perturbation fields
!
! Select perturbation fields based on switches
IF(ampu < 1e-5) uprt = 0.0
IF(ampv < 1e-5) vprt = 0.0
IF(ampw < 1e-5) wprt = 0.0
IF(amppt< 1e-5) ptprt= 0.0
IF(ampqv< 1e-5) qvprt= 0.0
IF(ampp < 1e-5) pprt = 0.0
utot2=0.0
vtot2=0.0
wtot2=0.0
pttot2=0.0
qvtot2=0.0
ptot2=0.0
do k=1,nz-1
do i=2,nx-2
do j=2,ny-2
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
IF (mp_opt > 0) THEN
CALL mpsumr(utot2, 1)
CALL mpsumr(vtot2, 1)
CALL mpsumr(wtot2, 1)
CALL mpsumr(pttot2, 1)
CALL mpsumr(qvtot2, 1)
CALL mpsumr(ptot2, 1)
END IF
usd=sqrt(utot2/totalpoint)
vsd=sqrt(vtot2/totalpoint)
wsd=sqrt(wtot2/totalpoint)
ptsd=sqrt(pttot2/totalpoint)
qvsd=sqrt(qvtot2/totalpoint)
psd=sqrt(ptot2/totalpoint)
enorm=sqrt(0.5*(utot2+vtot2+wtot2)/totalpoint)
IF (myproc == 0) THEN
print *,'Perturbation norms before scaling:'
print *,'unorm =',usd
print *,'vnorm =',vsd
print *,'wnorm =',wsd
print *,'ptnorm=',ptsd
print *,'qvnorm=',qvsd
print *,'Pnorm =',psd
print *,'KEnorm =',enorm
write(13,*) 'usd,vsd,wsd,ptsd,qvsd,psd,enorm'
write(13,'(7e11.4)') usd,vsd,wsd,ptsd,qvsd,psd,enorm
endif
IF (iensopt == 1) THEN ! scaling the perturbation in breeding IC
IF (myproc == 0) THEN
read(10,*) ampu,ampv,ampw,amppt,ampqv,ampp,ampke ! now hold initial norms
endif
CALL mpupdater(ampu,1)
CALL mpupdater(ampv,1)
CALL mpupdater(ampw,1)
CALL mpupdater(amppt,1)
CALL mpupdater(ampqv,1)
CALL mpupdater(ampp,1)
CALL mpupdater(ampke,1)
! for rescaling
! or assign valuess here manually for rescaling
IF (myproc == 0) THEN
print *,'Initial norms read back:'
print *,ampu,ampv,ampw,amppt,ampqv,ampp,ampke
endif
if(usd > 0.0) uscl = ampu/usd
if(vsd > 0.0) vscl = ampv/vsd
if(wsd > 0.0) wscl = ampw/wsd
if(ptsd > 0.0) ptscl = amppt/ptsd
if(qvsd > 0.0) qvscl = ampqv/qvsd
if(psd > 0.0) pscl = ampp/psd
if(enorm > 0.0) escl=ampke/enorm
rateu=-1.0 ! special value (no meaning)
ratev=-1.0
ratew=-1.0
ratept=-1.0
rateqv=-1.0
ratep=-1.0
rateke=-1.0
if(uscl > 0.0) rateu=1.0/uscl
if(vscl > 0.0) ratev=1.0/vscl
if(wscl > 0.0) ratew=1.0/wscl
if(ptscl > 0.0) ratept=1.0/ptscl
if(qvscl > 0.0) rateqv=1.0/qvscl
if(pscl > 0.0) ratep=1.0/pscl
if(escl > 0.0) rateke=1.0/escl
IF(myproc == 0) THEN
write(11,'(7f11.4)') usd,vsd,wsd,ptsd,qvsd,psd,enorm
write(11,'(7f11.4)') uscl,vscl,wscl,ptscl,qvscl,pscl,escl
write(11,'(7f11.4)') rateu,ratev,ratew,ratept,rateqv,ratep,rateke
print *,'totalpoint=',totalpoint
print *,'usd=',usd,' uscl=',uscl
print *,'vsd=',vsd,' vscl=',vscl
print *,'wsd=',wsd,' wscl=',wscl
print *,'ptsd=',ptsd,' ptscl=',ptscl
print *,'qvsd=',qvsd,' qvscl=',qvscl
print *,'psd=',psd,' pscl=',pscl
endif
uprt = uscl*uprt
vprt = vscl*vprt
wprt = wscl*wprt
ptprt= ptscl*ptprt
qvprt= qvscl*qvprt
pprt = pscl*pprt
ELSE
IF(myproc == 0) THEN
print *,'Perturbation norms after scaling:'
print *,'unorm =',abs(iorder)*usd
print *,'vnorm =',abs(iorder)*vsd
print *,'wnorm =',abs(iorder)*wsd
print *,'ptnorm=',abs(iorder)*ptsd
print *,'qvnorm=',abs(iorder)*qvsd
print *,'Pnorm =',abs(iorder)*psd
print *,'KEnorm =',abs(iorder)*enorm
write(14,*) 'usd,vsd,wsd,ptsd,qvsd,psd,enorm - after scaling'
write(14,'(7e11.4)') abs(iorder)*usd,abs(iorder)*vsd,abs(iorder)*wsd, &
abs(iorder)*ptsd,abs(iorder)*qvsd,abs(iorder)*psd,abs(iorder)*enorm
endif
END IF
ELSE IF (inibred ==1) THEN ! generate random perturbations
uprt=0.0
vprt=0.0
wprt=0.0
ptprt=0.0
qvprt=0.0
pprt=0.0
if(iseed < 0) then ! computer select iseed
CALL DATE_AND_TIME(VALUES=idate)
iseed=idate(8)
IF(myproc == 0) print *,'Machine decided ISEED=',iseed
IF(myproc == 0) write(98,*) iseed
else if (iseed == 9999) then
IF(myproc == 0) read(98,*) iseed
IF(myproc == 0) print *,'Read in ISEED: ',iseed
CALL mpupdatei(iseed,1)
endif
IF(ampu > 1e-10) THEN
IF(myproc == 0) print *,'U-PERT initial iseed =',iseed
call normalrand(iseed,nx*ny*nz,tem1)
do k=1,nz-1
do i=1,nx
do j=1,ny-1
uprt(i,j,k) = ampu*tem1(i,j,k)
enddo
enddo
enddo
END IF
IF(ampv > 1e-10) THEN
IF(myproc == 0) print *,'V-PERT initial iseed =',iseed
call normalrand(iseed,nx*ny*nz,tem1)
do k=1,nz-1
do i=1,nx-1
do j=1,ny
vprt(i,j,k) = ampv*tem1(i,j,k)
enddo
enddo
enddo
END IF
IF(ampw > 1e-10) THEN
IF(myproc == 0) print *,'W-PERT initial iseed =',iseed
call normalrand(iseed,nx*ny*nz,tem1)
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
wprt(i,j,k) = ampw*tem1(i,j,k)
enddo
enddo
enddo
END IF
IF(amppt > 1e-10) THEN
IF(myproc == 0) print *,'PT-PERT initial iseed =',iseed
call normalrand(iseed,nx*ny*nz,tem1)
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
ptprt(i,j,k) = amppt*tem1(i,j,k)
enddo
enddo
enddo
END IF
IF(ampqv > 1e-10) THEN
IF(myproc == 0) print *,'QV-PERT initial iseed =',iseed
call normalrand(iseed,nx*ny*nz,tem1)
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
qvprt(i,j,k) = ampqv*tem1(i,j,k)
enddo
enddo
enddo
END IF
IF(ampp > 1e-10) THEN
IF(myproc == 0) print *,'P-PERT initial iseed =',iseed
call normalrand(iseed,nx*ny*nz,tem1)
do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
pprt(i,j,k) = ampp*tem1(i,j,k)
enddo
enddo
enddo
END IF
tem1=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
! Calculate and save initial perturbations (norm)
utot2=0.0
vtot2=0.0
wtot2=0.0
pttot2=0.0
qvtot2=0.0
ptot2=0.0
do k=1,nz-1
do i=2,nx-2
do j=2,ny-2
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
IF (mp_opt > 0) THEN
CALL mpsumr(utot2, 1)
CALL mpsumr(vtot2, 1)
CALL mpsumr(wtot2, 1)
CALL mpsumr(pttot2, 1)
CALL mpsumr(qvtot2, 1)
CALL mpsumr(ptot2, 1)
END IF
usd=sqrt(utot2/totalpoint)
vsd=sqrt(vtot2/totalpoint)
wsd=sqrt(wtot2/totalpoint)
ptsd=sqrt(pttot2/totalpoint)
qvsd=sqrt(qvtot2/totalpoint)
psd=sqrt(ptot2/totalpoint)
enorm=sqrt(0.5*(utot2+vtot2+wtot2)/totalpoint)
IF(myproc == 0) THEN
print *,'Initial (random) perturbation norms:'
print *,'unorm =',usd
print *,'vnorm =',vsd
print *,'wnorm =',wsd
print *,'ptnorm=',ptsd
print *,'qvnorm=',qvsd
print *,'Pnorm =',psd
print *,'KEnorm =',enorm
write(10,'(7f11.4)') usd,vsd,wsd,ptsd,qvsd,psd,enorm
endif
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)
lenfil(3) = len_trim(filename(3))
lengbf(3) = len_trim(grdbasfn(3))
IF (mp_opt > 0 .AND. readsplit <= 0) THEN
tmpstr = filename(3)
WRITE(filename(3),'(2a,2i2.2)') tmpstr(1:lenfil(3)),'_',loc_x,loc_y
lenfil(3) = lenfil(3) + 5
tmpstr = grdbasfn(3)
WRITE(grdbasfn(3),'(2a,2i2.2)') tmpstr(1:lengbf(3)),'_',loc_x,loc_y
lengbf(3) = lengbf(3) + 5
END IF
IF (myproc == 0) THEN
WRITE(6,'(/2a/)') ' Reading file: ',filename(3)(1:lenfil(3))
WRITE(6,'(/2a/)') ' Reading file: ',grdbasfn(3)(1:lengbf(3))
END IF
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt,nchin,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 ) CALL arpsstop('dtaread errors.',1)
! 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 (abs(fx0-vx0) <=1.0 .AND.abs(fy0-vy0) <=1.0 .AND. &
! 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 (abs(iorder) > 1e-5) 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))
IF(myproc == 0) &
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 (abs(iorder) > 1e-5) 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
IF (myproc == 0) &
PRINT*,'Writing base state history dump ', grdbasfn(1)(1:lengbf(1))
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL 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)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
!
!-----------------------------------------------------------------------
!
! 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)
IF(myproc == 0) &
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
! blocking inserted for ordering i/o for message passing
DO i=0,nprocs-1,max_fopen
IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
CALL 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)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
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,6
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
INCLUDE 'mp.inc'
INTEGER :: nx,ny,nz,nzsoil ! 4 dimensions of array
REAL :: 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
!
!-----------------------------------------------------------------------
tem1=0.0
tem3dsoil = 0.0
IF(myproc == 0) PRINT *, ' ptprt: '
CALL subtrprt(nx,ny,nz, ptprt,ptbar,aptprt,aptbar,dptprt,nscale, &
is,js,ks,ie,je,ke)
IF(myproc == 0) PRINT *, ' pprt: '
CALL subtrprt(nx,ny,nz, pprt, pbar, apprt, apbar, dpprt,nscale, &
is,js,ks,ie,je,ke)
IF(myproc == 0) 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)
dqc=qc
dqr=qr
dqi=qi
dqs=qs
dqh=qh
dtke=tke
dkmh=kmh
dkmv=kmv
!-----------------------------------------------------------------------
!
! u wind components
!
!-----------------------------------------------------------------------
ie=nx
IF(myproc == 0) 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
IF(myproc == 0) 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
IF(myproc == 0) 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)
dtsoil=tsoil
dqsoil=qsoil
dwetcanp=wetcanp
draing=0.0
drainc=0.0
!
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, & 6
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
INCLUDE 'mp.inc'
INTEGER :: nx,ny,nz ! 3 dimensions of array
REAL :: 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
IF(myproc == 0) PRINT *, 'imid,jmid,kmid: ',imid,jmid,kmid
!
!-----------------------------------------------------------------------
!
! Tell us about a sample input point
!
!-----------------------------------------------------------------------
!
IF (abs(nscale) < 1e-5) THEN
IF(myproc == 0) &
PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), &
' b= ',(b(imid,jmid,kmid)+bbar(imid,jmid,kmid))
ELSE
IF(myproc == 0) &
PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), &
' b= ',b(imid,jmid,kmid)*abs(nscale)
END IF
!
!-----------------------------------------------------------------------
!
! Subtraction
!
!-----------------------------------------------------------------------
!
DO k=kstr,kend
DO j=jstr,jend
DO i=istr,iend
IF (abs(nscale) < 1e-5) 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)*nscale
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Tell us about a sample output point
!
!-----------------------------------------------------------------------
!
IF (abs(nscale) < 1e-5) THEN
IF(myproc == 0) &
PRINT *, ' c= ',c(imid,jmid,kmid)
ELSE
IF(myproc == 0) &
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
SUBROUTINE normalrand(iseed,n,rndn) 6,3
IMPLICIT NONE
INCLUDE 'mp.inc'
INTEGER :: iseed ! random seed
INTEGER :: n ! dimension size
INTEGER :: i, n0
REAL :: rndn (n) ! return standard normal distri.
REAL, ALLOCATABLE :: work1(:),work2(:)
REAL, PARAMETER :: pi=3.141592653589
REAL :: ave,var
ALLOCATE (work1(n),work2(n))
n0 = n
do i=1,n
iseed = mod(iseed*7141+54773,259200)
work1(i)=float(iseed)/259199.
enddo
do i=1,n
iseed = mod(iseed*7141+54773,259200)
work2(i)=float(iseed)/259199.
enddo
rndn = sqrt(-2.0*log(work1+tiny(1.0)))*cos(2.0*pi*work2)
! N(0,1) standard normal distribution
! statistics
! ave = sum(rndn)/float(n)
ave = sum(rndn)
IF (mp_opt > 0) THEN
CALL mpsumr(ave, 1)
CALL mptotali(n0)
END IF
ave = ave/float(n0)
work1 = rndn-ave
! var = dot_product(work1,work1)/float(n-1)
var = dot_product(work1,work1)
IF (mp_opt > 0) CALL mpsumr(var, 1)
var = var/float(n0-1)
IF(myproc == 0) print *,' AVE,VAR,SQRT(AVR):',ave,var,sqrt(var)
DEALLOCATE (work1,work2)
RETURN
END SUBROUTINE normalrand