!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSTINTRP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
PROGRAM arpstintrp,99
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This program interpolates two ARPS history data on grid of the same
! size to a time inbetween them. The output will be written into a new
! history dump file.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
!
! 2/25/1999. Written based on ARPSTINTRP.
!
! MODIFICATION HISTORY:
!
! 1999/10/13 (Gene Bassett)
! Corrected a roundoff error problem for the history dump reference
! times. Made the history dump output characterists similar
! to arpsintrp.
!
! 2001/06/18 (Gene Bassett)
! Corrected error with absolute time (iabstinit variables).
!
! 2002/03/19 (Keith Brewster)
! Corrected time calculations for the case when the user chooses
! to have curtim be relative to initime of the input file.
!
! 1 June 2002 Eric Kemp
! Soil variable updates.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Dimension of the base grid (input data).
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz
INTEGER :: nzsoil, nstyps
!
!-----------------------------------------------------------------------
!
! ARPS arrays. The last dimension is for the two sets of variables.
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: u (:,:,:,:) ! Total u-velocity (m/s).
REAL, ALLOCATABLE :: v (:,:,:,:) ! Total v-velocity (m/s).
REAL, ALLOCATABLE :: w (:,:,:,:) ! Total w-velocity (m/s).
REAL, ALLOCATABLE :: ptprt (:,:,:,:) ! Perturbation potential temperature
! from that of base state atmosphere (Kelvin).
REAL, ALLOCATABLE :: pprt (:,:,:,:) ! Perturbation pressure from that
! of base state atmosphere (Pascal).
REAL, ALLOCATABLE :: qv (:,:,:,:) ! Water vapor mixing ratio (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 )
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Fraction of soil type
INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsoil (:,:,:,:,:) ! Deep soil temperature (K)
REAL, ALLOCATABLE :: qsoil (:,:,:,:,:) ! Deep soil temperature (K)
REAL, ALLOCATABLE :: wetcanp(:,:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:,:) ! Snow cover
REAL, ALLOCATABLE :: raing(:,:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: prcrate(:,:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(:,:,:,:) = total precip. rate
! prcrate(:,:,:,:) = grid scale precip. rate
! prcrate(:,:,:,:) = cumulus precip. rate
! prcrate(:,:,:,:) = 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 :: 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,2).
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 on the staggered grid.
REAL, ALLOCATABLE :: zpsoil (:,:,:) ! The physical height coordinate defined at
REAL, ALLOCATABLE :: uprt (:,:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vprt (:,:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: qvprt (:,:,:,:) ! Perturbation water vapor specific humidity (kg/kg)
REAL, ALLOCATABLE :: tem1 (:,:,:) ! Temporary array
REAL, ALLOCATABLE :: tem2 (:,:,:) ! Temporary array
REAL, ALLOCATABLE :: tem3 (:,:,:) ! Temporary array
REAL, ALLOCATABLE :: tem4 (:,:,:) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
REAL :: amin, amax
CHARACTER (LEN=256) :: basdmpfn
INTEGER :: lbasdmpf
CHARACTER (LEN=256) :: ternfn,sfcoutfl,soiloutfl,temchar
INTEGER :: lternfn,lfn
INTEGER :: iss,is
REAL :: zpmax
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'indtflg.inc'
INCLUDE 'grid.inc'
INCLUDE 'exbc.inc'
INTEGER :: hinfmt,houtfmt, nchin, nchout
CHARACTER (LEN=256) :: filename
CHARACTER (LEN=256) :: grdbasfn
INTEGER :: lenfil,lengbf
INTEGER :: grdbas
INTEGER :: ireturn
REAL :: time
INTEGER :: gboutcnt, vroutcnt
DATA gboutcnt, vroutcnt /0,0/
INTEGER :: nfilemax
PARAMETER (nfilemax=2)
CHARACTER (LEN=256) :: hisfile(nfilemax)
INTEGER :: nhisfile,nd, length, lenstr
CHARACTER (LEN=80) :: timsnd
CHARACTER (LEN=80) :: new_runname
INTEGER :: tmstrln
REAL :: times(nfilemax), outtime, alpha, beta
INTEGER :: iabstinit,iabstinit1
INTEGER :: ioffset
INTEGER :: year1,month1,day1,hour1,minute1,second1,ioutabst, &
ioutabstinit
!
!-----------------------------------------------------------------------
!
! namelist Declarations:
!
!-----------------------------------------------------------------------
!
INTEGER :: use_data_t
CHARACTER (LEN=19) :: initime ! Real time in form of 'year-mo-dy:hr:mn:ss'
NAMELIST /INPUT/hinfmt,nhisfile,grdbasfn,hisfile
NAMELIST /output/ runname,use_data_t,initime,outtime, &
dirname,exbcdmp,exbchdfcompr,hdmpfmt,grbpkbit,hdfcompr, &
grdout,basout,varout,mstout,rainout,prcout,iceout, &
tkeout, trbout,sfcout,landout,radout,flxout, &
qcexout,qrexout,qiexout,qsexout,qhexout, &
totout,filcmprs,sfcdmp,soildmp,ngbrz,zbrdmp
INTEGER :: nsize, nxy,nxyz
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
mgrid = 1
nestgrd = 0
!
!-----------------------------------------------------------------------
!
! Set the default parameters
!
!-----------------------------------------------------------------------
!
hinfmt = 10
grdbasfn = 'X'
nhisfile = 1
hisfile(1) = 'X'
use_data_t = 1
initime ='0000-00-00:00:00:00'
runname = 'runname_not_set'
dirname = './'
exbcdmp = 1
exbchdfcompr = 0
hdmpfmt = 1
grbpkbit = 16
hdfcompr = 0
filcmprs = 0
basout = 0
grdout = 0
varout = 1
mstout = 1
iceout = 1
tkeout = 1
trbout = 0
rainout = 0
sfcout = 0
snowout = 0
landout = 0
qcexout = 0
qrexout = 0
qiexout = 0
qsexout = 0
qhexout = 0
sfcdmp = 1
soildmp = 1
ngbrz = 5
zbrdmp = 10000.0
WRITE(6,'(/9(/2x,a)/)') &
'###############################################################', &
'###############################################################', &
'### ###', &
'### Welcome to ARPSTINTRP ###', &
'### ###', &
'###############################################################', &
'###############################################################'
READ(5,INPUT, END=100)
nhisfile = 2
WRITE(6,'(/a/)') ' Input control parameters read in are:'
IF( hinfmt == 5) THEN
WRITE(6,'(2(2x,a))') &
'The Savi3D not supported as an INPUT file format', &
'Job stopped in ARPSTINTRP.'
STOP 9102
END IF
WRITE(6,'(1x,a,i3)') ' hinfmt =', hinfmt
WRITE(6,'(1x,a,i3)') ' nhisfile =', nhisfile
length = LEN( grdbasfn )
CALL strlnth
( grdbasfn, length )
WRITE(6,'(1x,a,a)') ' grdbasfn =', grdbasfn(1:length)
DO i=1,nhisfile
length = LEN( hisfile(i) )
CALL strlnth( hisfile(i), length )
WRITE(6,'(1x,a,i3,a,a)') ' hisfile(',i,')=',hisfile(i)(1:length)
END DO
!
!-----------------------------------------------------------------------
!
! Set the control parameters for output:
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/a/)') &
' Reading in control parameters for the output data files..'
READ (5,output,END=100)
IF ( exbchdfcompr > 4 ) rayklow = -1
WRITE(6,'(/2x,a,a)') &
'The run name to be used for constructing output file name is ', &
runname
new_runname = runname
totout = 1
CALL get_dims_from_data(hinfmt,trim(grdbasfn), &
nx,ny,nz,nzsoil,nstyps, ireturn)
IF (nstyps <= 0) nstyps = 1
nstyp = nstyps
!
!-----------------------------------------------------------------------
!
! Allocate the arrays.
!
!-----------------------------------------------------------------------
!
ALLOCATE(u(nx,ny,nz,2))
ALLOCATE(v(nx,ny,nz,2))
ALLOCATE(w(nx,ny,nz,2))
ALLOCATE(ptprt(nx,ny,nz,2))
ALLOCATE(pprt(nx,ny,nz,2))
ALLOCATE(qv(nx,ny,nz,2))
ALLOCATE(qc(nx,ny,nz,2))
ALLOCATE(qr(nx,ny,nz,2))
ALLOCATE(qi(nx,ny,nz,2))
ALLOCATE(qs(nx,ny,nz,2))
ALLOCATE(qh(nx,ny,nz,2))
ALLOCATE(tke(nx,ny,nz,2))
ALLOCATE(kmh(nx,ny,nz,2))
ALLOCATE(kmv(nx,ny,nz,2))
ALLOCATE(soiltyp(nx,ny,nstyps))
ALLOCATE(stypfrct(nx,ny,nstyps))
ALLOCATE(vegtyp(nx,ny))
ALLOCATE(roufns(nx,ny))
ALLOCATE(lai(nx,ny))
ALLOCATE(veg(nx,ny))
ALLOCATE(tsoil(nx,ny,nzsoil,0:nstyps,2))
ALLOCATE(qsoil(nx,ny,nzsoil,0:nstyps,2))
ALLOCATE(wetcanp(nx,ny,0:nstyps,2))
ALLOCATE(snowdpth(nx,ny,2))
ALLOCATE(raing(nx,ny,2))
ALLOCATE(rainc(nx,ny,2))
ALLOCATE(prcrate(nx,ny,4,2))
ALLOCATE(radfrc(nx,ny,nz,2))
ALLOCATE(radsw(nx,ny,2))
ALLOCATE(rnflx(nx,ny,2))
ALLOCATE(radswnet(nx,ny,2))
ALLOCATE(radlwin(nx,ny,2))
ALLOCATE(usflx(nx,ny,2))
ALLOCATE(vsflx(nx,ny,2))
ALLOCATE(ptsflx(nx,ny,2))
ALLOCATE(qvsflx(nx,ny,2))
ALLOCATE(ubar(nx,ny,nz))
ALLOCATE(vbar(nx,ny,nz))
ALLOCATE(wbar(nx,ny,nz))
ALLOCATE(ptbar(nx,ny,nz))
ALLOCATE(pbar(nx,ny,nz))
ALLOCATE(rhobar(nx,ny,nz))
ALLOCATE(qvbar(nx,ny,nz))
ALLOCATE(x(nx))
ALLOCATE(y(ny))
ALLOCATE(z(nz))
ALLOCATE(zp(nx,ny,nz))
ALLOCATE(zpsoil(nx,ny,nzsoil))
ALLOCATE(uprt(nx,ny,nz,2))
ALLOCATE(vprt(nx,ny,nz,2))
ALLOCATE(qvprt(nx,ny,nz,2))
ALLOCATE(tem1(nx,ny,nz))
ALLOCATE(tem2(nx,ny,nz))
ALLOCATE(tem3(nx,ny,nz))
ALLOCATE(tem4(nx,ny,nz))
nxyz = nx*ny*nz
nxy = nx*ny
CALL flzero(u , nxyz*2)
CALL flzero(v , nxyz*2)
CALL flzero(uprt , nxyz*2)
CALL flzero(vprt , nxyz*2)
CALL flzero(w , nxyz*2)
CALL flzero(ptprt , nxyz*2)
CALL flzero(pprt , nxyz*2)
CALL flzero(qvprt , nxyz*2)
CALL flzero(qc , nxyz*2)
CALL flzero(qr , nxyz*2)
CALL flzero(qi , nxyz*2)
CALL flzero(qs , nxyz*2)
CALL flzero(qh , nxyz*2)
CALL flzero(tke , nxyz*2)
CALL flzero(kmh , nxyz*2)
CALL flzero(kmv , nxyz*2)
CALL flzero(radfrc, nxyz*2)
CALL flzero(ubar , nxyz)
CALL flzero(vbar , nxyz)
CALL flzero(wbar , nxyz)
CALL flzero(ptbar , nxyz)
CALL flzero(pbar , nxyz)
CALL flzero(rhobar, nxyz)
CALL flzero(qvbar , nxyz)
CALL flzero(x, nx)
CALL flzero(y, ny)
CALL flzero(z, nz)
CALL flzero(zp , nxyz)
nsize = nx*ny*nzsoil
CALL flzero(zpsoil, nsize)
nsize = nx*ny*nzsoil*(1+nstyps)*2
CALL flzero(tsoil,nsize)
CALL flzero(qsoil,nsize)
nsize = nx*ny*(1+nstyps)*2
CALL flzero(wetcanp,nsize)
DO iss=1,nstyps
DO j=1,ny
DO i=1,nx
soiltyp (i,j,iss) = 0
stypfrct(i,j,iss) = 0.0
END DO
END DO
END DO
DO is=1,2
DO j=1,ny
DO i=1,nx
snowdpth(i,j,is) = 0
vegtyp (i,j) = 0
lai (i,j) = 0.0
roufns (i,j) = 0.0
veg (i,j) = 0.0
raing (i,j,is) = 0.0
rainc (i,j,is) = 0.0
prcrate(i,j,1,is) = 0.0
prcrate(i,j,2,is) = 0.0
prcrate(i,j,3,is) = 0.0
prcrate(i,j,4,is) = 0.0
radsw (i,j,is) = 0.0
rnflx (i,j,is) = 0.0
radswnet(i,j,is) = 0.0
radlwin(i,j,is) = 0.0
usflx (i,j,is) = 0.0
vsflx (i,j,is) = 0.0
ptsflx (i,j,is) = 0.0
qvsflx (i,j,is) = 0.0
END DO
END DO
END DO
ldirnam=LEN(dirname)
CALL strlnth( dirname , ldirnam)
lengbf=LEN(grdbasfn)
CALL strlnth( grdbasfn, lengbf)
WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)
!
!-----------------------------------------------------------------------
!
! Loop over data files
!
!-----------------------------------------------------------------------
!
ireturn = 0
DO nd=1,2
filename = hisfile(nd)
lenfil=LEN(filename)
CALL strlnth( filename, lenfil)
WRITE(6,'(/a,a,a)') &
' Data set ', filename(1:lenfil) ,' to be processed.'
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
CALL dtaread(nx,ny,nz,nzsoil,nstyps, &
hinfmt,nchin,grdbasfn(1:lengbf),lengbf, &
filename(1:lenfil),lenfil,time, x,y,z,zp,zpsoil, &
uprt(1,1,1,nd),vprt (1,1,1,nd),w (1,1,1,nd),ptprt(1,1,1,nd), &
pprt(1,1,1,nd),qvprt(1,1,1,nd),qc (1,1,1,nd),qr (1,1,1,nd), &
qi (1,1,1,nd),qs (1,1,1,nd),qh (1,1,1,nd),tke (1,1,1,nd), &
kmh (1,1,1,nd),kmv (1,1,1,nd), &
ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil(1,1,1,0,nd),qsoil(1,1,1,0,nd), &
wetcanp(1,1,0,nd),snowdpth(1,1,nd), &
raing(1,1,nd),rainc(1,1,nd),prcrate(1,1,1,nd), &
radfrc(1,1,1,nd),radsw(1,1,nd),rnflx(1,1,nd), &
radswnet(1,1,nd),radlwin(1,1,nd), &
usflx(1,1,nd),vsflx(1,1,nd),ptsflx(1,1,nd),qvsflx(1,1,nd), &
ireturn, tem1, tem2, tem3)
WRITE(6,'(/a,/2(a,i2),a,i4,a,/3(a,i2),a,f13.3,a/)') &
'History data read in for time: ', &
'month=',month,', day=', day,', year=',year,',', &
'hour =',hour ,',minute=',minute,', second=',second, &
', time=',time,'(s)'
CALL ctim2abss(year,month,day,hour,minute,second,iabstinit)
IF(nd == 1) THEN ! Save the values for data set 1
year1 = year
month1 = month
day1 = day
hour1 = hour
minute1= minute
second1= second
iabstinit1 = iabstinit
END IF
times(nd) = time + int(iabstinit - iabstinit1)
END DO
IF( hinfmt == 9 .AND. ireturn == 2 ) THEN
WRITE(6,'(/1x,a/)') 'The end of GrADS file was reached.'
CLOSE ( nchin )
CALL retunit( nchin )
GO TO 9001
END IF
IF( ireturn /= 0 ) GO TO 9002 ! Read was unsuccessful
!
IF( use_data_t == 1) THEN ! Use init time in input file
ioutabstinit=iabstinit1
year = year1
month = month1
day = day1
hour = hour1
minute= minute1
second= second1
ELSE
READ(initime,'(i4.4,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2)') &
year,month,day,hour,minute,second
CALL ctim2abss(year,month,day,hour,minute,second,ioutabstinit)
END IF
PRINT*,'ioutabstinit=',ioutabstinit
PRINT*,'iabstinit1=',iabstinit1
ioffset = ioutabstinit - iabstinit1
times(1) = times(1) - ioffset
times(2) = times(2) - ioffset
curtim = outtime
WRITE(6,'(/a,/2(a,i2),a,i4,a,/3(a,i2),a,f13.3,a/)') &
'In output file, the reference time is', &
'month=',month,', day=', day,', year=',year,',', &
'hour =',hour ,',minute=',minute,', second=',second, &
', & the time relative to this reference =',curtim
IF ( curtim > MAX(times(1),times(2)) .OR. &
curtim < MIN(times(1),times(2))) THEN
WRITE (*,*) "WARNING: Performing extrapolation. Desired time ", &
"is outside the range of the reference files."
END IF
IF( times(2) == times(1)) THEN
WRITE (*,*) "ERROR: times in reference files are the same, ", &
"can't perform interpolation."
STOP 1
END IF
alpha = (times(2)-curtim)/(times(2)-times(1))
beta = 1.0-alpha
WRITE (*,*)
WRITE (*,*) "Relative weights: file 1",alpha," file 2",beta
!-----------------------------------------------------------------------
!
! Calculate total fields from that for base state and perturbations
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
uprt (i,j,k,1)=alpha*uprt (i,j,k,1)+beta*uprt (i,j,k,2)
vprt (i,j,k,1)=alpha*vprt (i,j,k,1)+beta*vprt (i,j,k,2)
w (i,j,k,1)=alpha*w (i,j,k,1)+beta*w (i,j,k,2)
ptprt (i,j,k,1)=alpha*ptprt (i,j,k,1)+beta*ptprt (i,j,k,2)
pprt (i,j,k,1)=alpha*pprt (i,j,k,1)+beta*pprt (i,j,k,2)
qvprt (i,j,k,1)=alpha*qvprt (i,j,k,1)+beta*qvprt (i,j,k,2)
qc (i,j,k,1)=alpha*qc (i,j,k,1)+beta*qc (i,j,k,2)
qr (i,j,k,1)=alpha*qr (i,j,k,1)+beta*qr (i,j,k,2)
qi (i,j,k,1)=alpha*qi (i,j,k,1)+beta*qi (i,j,k,2)
qs (i,j,k,1)=alpha*qs (i,j,k,1)+beta*qs (i,j,k,2)
qh (i,j,k,1)=alpha*qh (i,j,k,1)+beta*qh (i,j,k,2)
tke (i,j,k,1)=alpha*tke (i,j,k,1)+beta*tke (i,j,k,2)
kmh (i,j,k,1)=alpha*kmh (i,j,k,1)+beta*kmh (i,j,k,2)
kmv (i,j,k,1)=alpha*kmv (i,j,k,1)+beta*kmv (i,j,k,2)
u (i,j,k,1)=uprt (i,j,k,1)+ubar (i,j,k)
v (i,j,k,1)=vprt (i,j,k,1)+vbar (i,j,k)
qv (i,j,k,1)=qvprt(i,j,k,1)+qvbar(i,j,k)
radfrc(i,j,k,1)=alpha*radfrc(i,j,k,1)+beta*radfrc(i,j,k,2)
END DO
END DO
END DO
DO is=0,nstyp
DO k=1,nzsoil
DO j=1,ny-1
DO i=1,nx-1
tsoil(i,j,k,is,1)=alpha*tsoil(i,j,k,is,1) &
+beta*tsoil(i,j,k,is,2)
qsoil(i,j,k,is,1)=alpha*qsoil(i,j,k,is,1) &
+beta*qsoil(i,j,k,is,2)
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
wetcanp(i,j,is,1)=alpha*wetcanp(i,j,is,1) &
+beta*wetcanp(i,j,is,2)
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
snowdpth(i,j,1)=alpha*snowdpth(i,j,1)+beta*snowdpth(i,j,2)
raing (i,j,1)=alpha*raing (i,j,1)+beta*raing (i,j,2)
rainc (i,j,1)=alpha*rainc (i,j,1)+beta*rainc (i,j,2)
prcrate(i,j,1,1)=alpha*prcrate(i,j,1,1)+beta*prcrate(i,j,1,2)
prcrate(i,j,2,1)=alpha*prcrate(i,j,2,1)+beta*prcrate(i,j,2,2)
prcrate(i,j,3,1)=alpha*prcrate(i,j,3,1)+beta*prcrate(i,j,3,2)
prcrate(i,j,4,1)=alpha*prcrate(i,j,4,1)+beta*prcrate(i,j,4,2)
radsw (i,j,1)=alpha*radsw (i,j,1)+beta*radsw (i,j,2)
rnflx (i,j,1)=alpha*rnflx (i,j,1)+beta*rnflx (i,j,2)
radswnet(i,j,1)=alpha*radswnet(i,j,1)+beta*radswnet(i,j,2)
radlwin(i,j,1)=alpha*radlwin(i,j,1)+beta*radlwin(i,j,2)
usflx (i,j,1)=alpha*usflx (i,j,1)+beta*usflx (i,j,2)
vsflx (i,j,1)=alpha*vsflx (i,j,1)+beta*vsflx (i,j,2)
ptsflx (i,j,1)=alpha*ptsflx (i,j,1)+beta*ptsflx (i,j,2)
qvsflx (i,j,1)=alpha*qvsflx (i,j,1)+beta*qvsflx (i,j,2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Print out the max/min of output varaibles.
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/1x,a/)') &
'Min. and max. of data interpolated to the new time:'
CALL a3dmax0(x,1,nx,1,nx,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(/1x,2(a,e13.6))') 'xmin = ', amin,', xmax =',amax
CALL a3dmax0(y,1,ny,1,ny,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ymin = ', amin,', ymax =',amax
CALL a3dmax0(z,1,nz,1,nz,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'zmin = ', amin,', zmax =',amax
CALL a3dmax0(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'zpmin = ', amin,', zpmax =',amax
CALL a3dmax0(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,', ubarmax =',amax
CALL a3dmax0(vbar,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,', vbarmax =',amax
CALL a3dmax0(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,', ptbarmax=',amax
CALL a3dmax0(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,', pbarmax =',amax
CALL a3dmax0(rhobar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'rhobarmin=', amin,', rhobarmax=',amax
CALL a3dmax0(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin,', qvbarmax=',amax
CALL a3dmax0(uprt,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'uprtmin = ', amin,', uprtmax =',amax
CALL a3dmax0(vprt,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vprtmin = ', amin,', vprtmax =',amax
CALL a3dmax0(w,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'wmin = ', amin,', wmax =',amax
CALL a3dmax0(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,', ptprtmax=',amax
CALL a3dmax0(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,', pprtmax =',amax
CALL a3dmax0(qvprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qvprtmin= ', amin,', qvprtmax=',amax
CALL a3dmax0(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qcmin = ', amin,', qcmax =',amax
CALL a3dmax0(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qrmin = ', amin,', qrmax =',amax
CALL a3dmax0(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qimin = ', amin,', qimax =',amax
CALL a3dmax0(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qsmin = ', amin,', qsmax =',amax
CALL a3dmax0(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qhmin = ', amin,', qhmax =',amax
CALL a3dmax0(tke,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'tkemin = ', amin,', tkemax =',amax
CALL a3dmax0(kmh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'kmhmin = ', amin,', kmhmax =',amax
CALL a3dmax0(kmv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'kmvmin = ', amin,', kmvmax =',amax
CALL a3dmax0(raing,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'raingmin= ', amin,', raingmax=',amax
CALL a3dmax0(rainc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'raincmin= ', amin,', raincmax=',amax
CALL a3dmax0(prcrate(1,1,1,1),1,nx,1,nx-1,1,ny,1,ny-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'prcr1min= ', amin,', prcr1max=',amax
CALL a3dmax0(prcrate(1,1,2,1),1,nx,1,nx-1,1,ny,1,ny-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'prcr2min= ', amin,', prcr2max=',amax
CALL a3dmax0(prcrate(1,1,3,1),1,nx,1,nx-1,1,ny,1,ny-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'prcr3min= ', amin,', prcr3max=',amax
CALL a3dmax0(prcrate(1,1,4,1),1,nx,1,nx-1,1,ny,1,ny-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'prcr4min= ', amin,', prcr4max=',amax
DO iss = 0, nstyp
CALL a3dmax0(tsoil(1,1,1,iss,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'tsoil_sfcmin = ', amin,', tsoil_sfcmax =',amax,' for soil type=',iss
CALL a3dmax0(tsoil(1,1,2,iss,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'tsoil_dpmin= ', amin,', tsoil_dpmax=',amax,' for soil type=',iss
CALL a3dmax0(qsoil(1,1,1,iss,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'qsoil_sfcmin = ', amin,', qsoil_sfcmax =',amax,' for soil type=',iss
CALL a3dmax0(qsoil(1,1,2,iss,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'qsoil_dpmin = ', amin,', qsoil_dpmax =',amax,' for soil type=',iss
CALL a3dmax0(wetcanp(1,1,iss,1), &
1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'wetcmin = ', amin,', wetcmax =',amax,' for soil type=',iss
END DO
CALL a3dmax0(roufns,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'roufnmin =', amin,', roufnmax =',amax
CALL a3dmax0(veg,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vegmin = ', amin,', vegmax =',amax
CALL a3dmax0(radfrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'radfnmin =', amin,', radfnmax =',amax
CALL a3dmax0(radsw,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'radswmin =', amin,', radswmax =',amax
CALL a3dmax0(rnflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'rnflxmin =', amin,', rnflxmax =',amax
CALL a3dmax0(radswnet,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'radswnetmin =', amin,', radswnetmax =',amax
CALL a3dmax0(radlwin,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'radlwinmin =', amin,', radlwinmax =',amax
CALL a3dmax0(usflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'usflxnmin =', amin,', usflxmax =',amax
CALL a3dmax0(vsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vsflxmin =', amin,', vsflxmax =',amax
CALL a3dmax0(ptsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ptflxmin =', amin,', ptflxmax =',amax
CALL a3dmax0(qvsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qvflxmin =', amin,', qvflxmax =',amax
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1 (i,j,k)= 0.0 ! To be put in place of wbar
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Data dump of the model grid and base state arrays:
!
! First find a unique name basdmpfn(1:lbasdmpf) for the grid and
! base state array dump file
!
! If grid/base state data has been written out once, skip
! the following writing block. Also no need to write out
! separate data for Savi3D dump. The same for GrADS dump.
!
!-----------------------------------------------------------------------
!
WRITE (cmnt(nocmnt),'(a,i4,a,i4,a,i4)') &
' nx =',nx,', ny =',ny,', nz =',nz
runname = new_runname
houtfmt = hdmpfmt
grbpkbit = 16
CALL gtlfnkey(runname, lfnkey)
IF(houtfmt /= 9 ) THEN
IF( gboutcnt == 1 ) GO TO 500 ! If done already, skip this part.
CALL gtbasfn(runname(1:lfnkey),dirname,ldirnam,hdmpfmt, &
1,0,basdmpfn,lbasdmpf)
PRINT*
PRINT*,'Output grid/base state file is ', basdmpfn(1:lbasdmpf)
grdbas = 1 ! Dump out grd and base state arrays only
CALL dtadump(nx,ny,nz,nzsoil,nstyps,hdmpfmt,nchout, &
basdmpfn(1:lbasdmpf),grdbas,filcmprs, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
tke,kmh,kmv, &
ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,zpsoil, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
tem2,tem3,tem4)
gboutcnt = 1
500 CONTINUE
END IF
!
!-----------------------------------------------------------------------
!
! Then the time dependent fields:
!
!-----------------------------------------------------------------------
!
IF( .NOT. (houtfmt == 9 .AND. vroutcnt == 1) ) THEN
!
!-----------------------------------------------------------------------
!
! Reconstruct the file name using the specified directory name
!
!-----------------------------------------------------------------------
!
CALL gtdmpfn(runname(1:lfnkey),dirname, &
ldirnam,curtim,hdmpfmt,1,0, hdmpfn, ldmpf)
END IF
WRITE(6,'(a,a)') 'Writing t-dependent variable history dump ', &
hdmpfn(1:ldmpf)
grdbas = 0
CALL dtadump(nx,ny,nz,nzsoil,nstyps,hdmpfmt,nchout, &
hdmpfn(1:ldmpf),grdbas,filcmprs, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
tke,kmh,kmv, &
ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,zpsoil, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
tem2,tem3,tem4)
!
!-----------------------------------------------------------------------
!
! Write out soil model variable file
!
!-----------------------------------------------------------------------
!
IF ( sfcin == 1 ) THEN
CALL cvttsnd( curtim, timsnd, tmstrln )
soiloutfl = runname(1:lfnkey)//".soilvar."//timsnd(1:tmstrln)
lfn = lfnkey + 9 + tmstrln
IF( dirname /= ' ' ) THEN
temchar = soiloutfl
soiloutfl = dirname(1:ldirnam)//'/'//temchar
lfn = lfn + ldirnam + 1
END IF
CALL fnversn(soiloutfl, lfn)
IF (soildmp > 0) THEN
PRINT *, 'Writing soil data to ',soiloutfl(1:lfn)
CALL wrtsoil(nx,ny,nzsoil,nstyps,soiloutfl(1:lfn),dx,dy,zpsoil, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
1,1,1,1,1, &
tsoil,qsoil,wetcanp,snowdpth,soiltyp)
IF (soildmp == 1) CALL soilcntl(nx,ny,nzsoil,zpsoil,soiloutfl(1:lfn), &
1,1,1,1,1,x,y)
END IF
END IF ! sfcin.eq.1
!-----------------------------------------------------------------------
!
! Write out surface property data file: sfcoutfl .
!
!-----------------------------------------------------------------------
!
IF ( landin == 1 ) THEN
sfcoutfl = runname(1:lfnkey)//".sfcdata"
lfn = lfnkey + 8
IF( dirname /= ' ' ) THEN
temchar = sfcoutfl
sfcoutfl = dirname(1:ldirnam)//'/'//temchar
lfn = lfn + ldirnam + 1
END IF
CALL fnversn(sfcoutfl, lfn)
IF (sfcdmp > 0) THEN
PRINT *, 'Write surface property data in ',sfcoutfl(1:lfn)
CALL wrtsfcdt(nx,ny,nstyps,sfcoutfl(1:lfn), dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
1,1,1,1,1,0, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg,veg)
IF (sfcdmp == 1) CALL sfccntl(nx,ny, sfcoutfl(1:lfn), &
1,1,1,1,1,0, x,y, tem1,tem2)
END IF
END IF ! landin.eq.1
STOP 0
100 WRITE(6,'(a)') &
'Error reading NAMELIST file. Program ARPSTINTRP stopped.'
STOP 9104
9001 CONTINUE
WRITE(6,'(/2x,a)')'For the output grid:'
WRITE(6,'(2x,a,f12.4)') &
'The latitude of the output grid center, ctrlat=',ctrlat
WRITE(6,'(2x,a,f12.4/)') &
'The longitude of the output grid center, ctrlon=',ctrlon
WRITE(6,'(2x,a/2x,a,2f15.4,a)') &
'The SW corner (i,j)=(2,2) of the grid is located at ', &
'(',xgrdorg,ygrdorg,') of the input grid.'
STOP 9001
9002 CONTINUE
WRITE(6,'(1x,a,i2,/1x,a)') &
'Data read was unsuccessful. ireturn =', ireturn, &
'Job stopped in ARPSTINTRP.'
STOP 9002
END PROGRAM arpstintrp