!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSINTRP_LS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
PROGRAM arpsintrp_ls,119
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This program interpolates gridded data from one ARPS grid to another.
! It can be used to prepare data for running ARPS in a one-way nested
! mode. It's exepcted to replace ARPSR2H in this capacity.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
!
! 3/27/1997. Written based on ARPSR2H and ARPSCVT.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
! nx, ny, nz: Dimensions of input grid.
! nx1, ny1, nz1: Dimensions of output grid.
!-----------------------------------------------------------------------
!
INTEGER :: nx ! Number of input grid points in the x-direction
INTEGER :: ny ! Number of input grid points in the y-direction
INTEGER :: nz ! Number of input grid points in the z-direction
INTEGER :: nx1 ! Number of output grid points in the x-direction
INTEGER :: ny1 ! Number of output grid points in the y-direction
INTEGER :: nz1 ! Number of output grid points in the z-direction
INTEGER :: nxyz,nxy,nxyz1,nxy1
!
!-----------------------------------------------------------------------
!
! Parameters for the utput grid.
!
!-----------------------------------------------------------------------
!
INTEGER :: nstyps
PARAMETER (nstyps=1)
INTEGER :: samgrd ! Are the output and the input grids the same?
! =1, the grids are the same
! =0, the grids are different
REAL :: xorig1, yorig1, zorig1
REAL :: xctr1 , yctr1
REAL :: dx1,dy1 ! Grid intervals of the refined grid.
INTEGER :: strhopt1 ! Vertical grid stretching option.
! = 0, no stretching in vertical.
! >= 1, with stretching in vertical.
REAL :: dz1 ! Average grid spacing in vertical direction in
! transformed computational space (m).
REAL :: dzmin1 ! Minimun grid spacing in vertical direction in
! physcal space (m).
REAL :: zrefsfc1 ! The reference height of the surface
! (ground level) (m)
REAL :: dlayer11 ! The depth of the lower layer with uniform
! (dz=dzmin) vertical spacing (m)
REAL :: dlayer21 ! The depth of the mid layer with stetched
! vertical spacing (m)
REAL :: strhtune1 ! Tuning parameter for stretching option 2
! A Value between 0.2 and 5.0 is appropriate.
! A larger value gives a more linear stretching.
REAL :: zflat1 ! The height at which the grid levels
! becomes flat in the terrain-following
! coordinate transformation (m).
!
!-----------------------------------------------------------------------
!
! Note:
!
! Given nx1, ny1 and nz1, the physical doma3din size of the refined
! grid will be xl1=(nx1-3)*dx1 by yl1=(ny1-3)*dy1 by zh1=(nz1-3)*dz1.
! Dx1, dy1 and dz1 are the grid intervals of the refined grid.
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: a3din(:,:,:) ! 3-D real array read in
INTEGER, ALLOCATABLE :: i2din(:,:) ! 2-D integer array read in
REAL, ALLOCATABLE :: a2din(:,:) ! 2-D real array 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 on the staggered grid.
REAL, ALLOCATABLE :: hterain(:,:) ! The height of terrain.
REAL, ALLOCATABLE :: tem3d(:,:,:) ! Work array
!
!-----------------------------------------------------------------------
!
! Arrays on the new grid:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: a3dout(:,:,:) ! 3-D real array interpolated from a3din.
INTEGER, ALLOCATABLE :: i2dout(:,:)! 2-D integer array derived from i2din
REAL, ALLOCATABLE :: a2dout(:,:) ! 2-D real array interpolated from a2din
REAL, ALLOCATABLE :: x1 (:) ! New grid x-coord. on the original grid
REAL, ALLOCATABLE :: x11 (:) ! New grid x-coord. set for new grid
REAL, ALLOCATABLE :: y1 (:) ! New grid y-coord. on the original grid
REAL, ALLOCATABLE :: y11 (:) ! New grid y-coord. set for new grid
REAL, ALLOCATABLE :: z1 (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, ALLOCATABLE :: zp1 (:,:,:) ! The physical height coordinate defined at
! w-point on the staggered grid.
REAL, ALLOCATABLE :: hterain1(:,:) ! Terrain height (m)
REAL, ALLOCATABLE :: tem3d1(:,:,:) ! Work array
REAL, ALLOCATABLE :: zp1d1 (:) ! Temporary array
REAL, ALLOCATABLE :: dzp1d1(:) ! Temporary array
REAL :: zflat11,za,zb
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: is, js
REAL :: s1,s2,s3,s4
REAL :: amin, amax
REAL :: xs1,ys1
CHARACTER (LEN=80) :: basdmpfn
INTEGER :: lbasdmpf
CHARACTER (LEN=80) :: ternfn,sfcoutfl,soiloutfl,temchar
INTEGER :: lternfn,lfn
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'bndry.inc'
INCLUDE 'indtflg.inc'
INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,nf,lenfil
PARAMETER (nhisfile_max=200)
CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
INTEGER :: ireturn
INTEGER :: houtfmt
CHARACTER (LEN=132) :: filename
INTEGER :: grdbas
REAL :: time
INTEGER :: nfilemax
PARAMETER (nfilemax=200)
CHARACTER (LEN=80) :: exbcfn
CHARACTER (LEN=80) :: timsnd
CHARACTER (LEN=80) :: new_runname
INTEGER :: tmstrln
CHARACTER (LEN=15) :: ctime
INTEGER :: nfile, length
REAL :: xeps, yeps
REAL :: ctrx,ctry,swx,swy,alatpro(2),sclf,dxscl,dyscl
REAL :: ctrlat1,ctrlon1,latitud1
CHARACTER (LEN=40) :: fmtver
PARAMETER (fmtver='004.10 Binary Data')
CHARACTER (LEN=40) :: fmtverin
CHARACTER (LEN=10) :: tmunit
CHARACTER (LEN=12) :: label
INTEGER :: nxin,nyin,nzin
INTEGER :: stgr,oldver
INTEGER :: inch,nchanl,exbchanl,sfchanl,soilchanl,trnchanl
INTEGER :: iyr, imon, idy, ihr, imin, isec
INTEGER :: ierr, itema
INTEGER :: varout1,mstout1,rainout1,prcout1,iceout1,tkeout1,trbout1
INTEGER :: sfcout1,landout1
REAL :: xgrdorg1,ygrdorg1
INTEGER :: idummy
REAL :: rdummy
INTEGER :: istat
LOGICAL :: fexist, dmpexbc, lsfcdmp, lsoildmp
!
!-----------------------------------------------------------------------
!
! namelist Declarations:
!
!-----------------------------------------------------------------------
!
NAMELIST /INPUT/ hinfmt,nhisfile, grdbasfn, hisfile
NAMELIST /output_dims/ nx1,ny1,nz1
NAMELIST /jobname/ runname
NAMELIST /newgrid/ samgrd,strhopt1,xctr1,yctr1, &
dx1,dy1,dz1,dzmin1, &
zrefsfc1,dlayer11,dlayer21,strhtune1,zflat1
NAMELIST /output/ dirname,exbcdmp,hdmpfmt,grbpkbit, &
grdout,basout,varout,mstout,rainout,prcout,iceout, &
tkeout, trbout,sfcout,landout, &
qcexout,qrexout,qiexout,qsexout,qhexout, &
totout,filcmprs
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
mgrid = 1
nestgrd = 0
WRITE(6,'(/9(/5x,a)/)') &
'###############################################################', &
'###############################################################', &
'### ###', &
'### Welcome to ARPSINTRP ###', &
'### This program converts the history dump data ###', &
'### sets generated by ARPS, between various formats. ###', &
'### ###', &
'###############################################################', &
'###############################################################'
!
!-----------------------------------------------------------------------
! Get the names of the input data files.
!-----------------------------------------------------------------------
!
CALL get_input_file_names
(hinfmt,grdbasfn,hisfile,nhisfile)
lengbf = len_trim(grdbasfn)
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nstyps, ireturn)
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz
IF( nhisfile > nhisfile_max) THEN
WRITE(6,'(1x,a,i3,/1x,a,/1x,a,/1x,a)') &
'The number of history files to be processed exceeded ', &
nhisfile_max,' please reduce the number of files to be processed',&
'in a single job, or edit the program and re-set parameter ', &
'nhisfile_max to a larger value. '
STOP 9010
END IF
length = LEN_trim( grdbasfn )
WRITE(6,'(1x,a,a)') ' grdbasfn =', grdbasfn(1:length)
DO i=1,nhisfile
length = LEN_trim( hisfile(i) )
WRITE(6,'(1x,a,i3,a,a)') ' hisfile(',i,')=',hisfile(i)(1:length)
END DO
!-----------------------------------------------------------------------
! Read output grid dimensions
!-----------------------------------------------------------------------
READ(5,output_dims, END=105)
WRITE(6,'(/a,a/)') &
'NAMELIST block intrp_dims sucessfully read.'
WRITE(6,'(3(a,i5))') 'nx1=',nx1,', ny1=',ny1,', nz1=',nz1
GO TO 10
105 WRITE(6,'(/a,a/)') &
'Error reading NAMELIST block intrp_dims. ', &
'Program ARPSINTRP stopped.'
STOP 1
10 CONTINUE
READ (5,jobname,ERR=9000)
WRITE(6,'(/5x,a,a)') 'The name of this run is: ', runname
new_runname = runname
CALL gtlfnkey(new_runname, lfnkey)
!
!-----------------------------------------------------------------------
!
! Set the output grid and the variable control parameters
!
!-----------------------------------------------------------------------
!
READ (5,newgrid,ERR=9000)
PRINT*
PRINT*,' Input parameters for the new refined grid:'
PRINT*
PRINT*,' Input samgrd:'
PRINT*,' Input was ',samgrd
PRINT*,' Input dx1:'
PRINT*,' Input was ',dx1
PRINT*,' Input dy1:'
PRINT*,' Input was ',dy1
PRINT*,' Input strhopt1:'
PRINT*,' Input was ',strhopt1
PRINT*,' Input dz1:'
PRINT*,' Input was ',dz1
PRINT*,' Input dzmin1:'
PRINT*,' Input was ',dzmin1
PRINT*,' Input xctr1:'
PRINT*,' Input was ',xctr1
PRINT*,' Input yctr1:'
PRINT*,' Input was ',yctr1
!
!-----------------------------------------------------------------------
!
! Set the control parameters for output:
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/a/)') &
' Reading in control parameters for the output data files..'
READ (5,output,ERR=9000)
houtfmt = hdmpfmt
IF( houtfmt /= 1 ) THEN
WRITE(6,'(/1x,a,a/)') 'Output format is not 1. Reset it to 1.'
houtfmt = 1
END IF
totout = 1
basout = 0
ldirnam=LEN_trim(dirname)
lengbf=LEN_trim(grdbasfn)
WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)
INQUIRE(FILE=grdbasfn(1:lengbf), EXIST = fexist )
IF( fexist ) GO TO 110
INQUIRE(FILE=grdbasfn(1:lengbf)//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( grdbasfn(1:lengbf)//'.Z' )
GO TO 110
END IF
INQUIRE(FILE=grdbasfn(1:lengbf)//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( grdbasfn(1:lengbf)//'.gz' )
GO TO 110
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//grdbasfn(1:lengbf)// &
' or its compressed version not found.', &
'Program stopped in ARPSINTRP.'
STOP 9011
110 CONTINUE
ALLOCATE(a3din(nx,ny,nz))
a3din = 0.0
ALLOCATE(i2din(nx,ny))
i2din = 0
ALLOCATE(a2din(nx,ny))
a2din =0.0
ALLOCATE(x (nx))
ALLOCATE(y (ny))
ALLOCATE(z (nz))
x = 0.0
y = 0.0
z = 0.0
ALLOCATE(zp (nx,ny,nz))
zp = 0.0
ALLOCATE(hterain(nx,ny))
hterain=0.0
ALLOCATE(tem3d(nx,ny,nz))
tem3d=0.0
!
ALLOCATE(a3dout(nx1,ny1,nz1))
a3dout=0.0
ALLOCATE(i2dout(nx1,ny1))
i2dout=0.0
ALLOCATE(a2dout(nx1,ny1))
a2dout=0.0
ALLOCATE(x1 (nx1))
x1=0.0
ALLOCATE(x11 (nx1))
x11=0.0
ALLOCATE(y1 (ny1))
y1=0.0
ALLOCATE(y11 (ny1))
y11=0.0
ALLOCATE(z1 (nz1))
z1=0.0
ALLOCATE(zp1 (nx1,ny1,nz1))
zp1=0.0
ALLOCATE(hterain1(nx1,ny1))
hterain1=0.0
ALLOCATE(tem3d1(nx1,ny1,nz1))
tem3d1=0.0
ALLOCATE(zp1d1 (nz1))
zp1d1=0.0
ALLOCATE(dzp1d1(nz1))
dzp1d1=0.0
!
!-----------------------------------------------------------------------
!
! Get the IO unit numbers for input grid and base state file
!
!-----------------------------------------------------------------------
!
CALL getunit( inch )
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(grdbasfn(1:lengbf), '-F f77 -N ieee', ierr)
!
!-----------------------------------------------------------------------
!
! Open the input grdbas file
!
!-----------------------------------------------------------------------
!
grdbas = 1
OPEN(UNIT=inch,FILE=grdbasfn(1:lengbf), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 9001
READ(inch,ERR=9110,END=9120) fmtverin
READ(inch,ERR=9110,END=9120) runname
READ(inch,ERR=9110,END=9120) nocmnt
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
READ(inch,ERR=9110,END=9120) cmnt(i)
END DO
END IF
READ(inch,ERR=9110,END=9120) time,tmunit
curtim = time
!
!-----------------------------------------------------------------------
!
! Get dimensions of data in binary file and check against
! the dimensions defined in dims.inc
!
!-----------------------------------------------------------------------
!
READ(inch,ERR=9110,END=9120) nxin, nyin, nzin
IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE(6,'(1x,a)') &
' Dimensions in ARPSINTRP inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
WRITE(6,'(1x,a)') &
' Program aborted in ARPSINTRP.'
STOP
END IF
WRITE(6,'(1x,a,f8.1,a,f8.3,a/)') &
'To read grid and base state data at time ', time, &
' secs = ',(time/60.),' mins.'
READ(inch,ERR=9110,END=9120) &
grdin,basin,varin,mstin,icein, &
trbin,idummy,idummy,landin,totin, &
tkein,idummy,idummy,mapproj,month, &
day, year, hour, minute, second
IF ( grdin /= 1 .OR. basin /= 1 ) THEN
WRITE(6,'(1x,a,a,a/a)') &
'File '//grdbasfn(1:lengbf)//' is not .bingrdbas file', &
'A .bingrdbas file is required. Program stoped in ARPSINTRP'
STOP 9012
END IF
IF ( varin == 1 ) THEN
varout1 = varout
ELSE
varout1 = 0
END IF
IF ( mstin == 1 ) THEN
mstout1 = mstout
ELSE
mstout1 = 0
END IF
IF ( icein == 1 ) THEN
iceout1 = iceout
ELSE
iceout1 = 0
END IF
IF ( tkein == 1 ) THEN
tkeout1 = tkeout
ELSE
tkeout1 = 0
END IF
IF ( trbin == 1 ) THEN
trbout1 = trbout
ELSE
trbout1 = 0
END IF
IF ( sfcin == 1 ) THEN
sfcout1 = sfcout
ELSE
sfcout1 = 0
END IF
IF ( landin == 1 ) THEN
landout1 = landout
ELSE
landout1 = 0
END IF
READ(inch,ERR=9110,END=9120) &
umove,vmove,xgrdorg,ygrdorg,trulat1, &
trulat2,trulon,sclfct,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
tstop,thisdmp,latitud,ctrlat,ctrlon
IF ( totin /= 0 ) THEN
!
!-----------------------------------------------------------------------
!
! Read in additional parameters for ARPS history dump 4.0 or later
! version.
!
!-----------------------------------------------------------------------
!
READ(inch,ERR=9110,END=9120) &
nstyp, prcout,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
IF ( nstyp < 1 ) THEN
nstyp = 1
END IF
IF ( prcin == 1 ) THEN
prcout1 = prcout
ELSE
prcout1 = 0
END IF
READ(inch,ERR=9110,END=9120) &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
END IF
!
!-----------------------------------------------------------------------
!
! Read in x,y and z at grid cell centers (scalar points).
!
!----------------------------------------------------------------------
!
IF( grdin == 1 .OR. grdbas == 1 ) THEN
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) x
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array x.'
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) y
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array y.'
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) z
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array z.'
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) zp
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array zp.'
END IF ! grdin
!
!-----------------------------------------------------------------------
!
! Set up the map projection of the input grid.
!
!-----------------------------------------------------------------------
!
alatpro(1)=trulat1
alatpro(2)=trulat2
dx = x(2)-x(1)
dy = y(2)-y(1)
IF( sclfct /= 1.0) THEN
sclf = 1.0/sclfct
dxscl = dx*sclf
dyscl = dy*sclf
ELSE
sclf = 1.0
dxscl = dx
dyscl = dy
END IF
PRINT*,mapproj,sclf,alatpro,trulon,ctrlat,ctrlon
CALL setmapr( mapproj,sclf,alatpro,trulon )
CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
swx = ctrx - (FLOAT(nx-3)/2.)*dxscl
swy = ctry - (FLOAT(ny-3)/2.)*dyscl
CALL setorig( 1, swx, swy) ! set up the model origin to the coord.
DO j=1,ny
DO i=1,nx
hterain(i,j) = zp(i,j,2)
END DO
END DO
dx = x(2) - x(1)
dy = y(2) - y(1)
dz = z(2) - z(1)
xorig = x(2)
yorig = y(2)
zorig = z(2)
WRITE(6,'(1x,a,3f15.3)') 'xorig, yorig, zorig =', &
xorig, yorig, zorig
!
!-----------------------------------------------------------------------
!
! If the new grid is same as the original one, no need to
! interpolate
!
!-----------------------------------------------------------------------
!
IF ( samgrd == 1 ) THEN
DO i=1,nx1
x1(i)=x(i)
x11(i)=x(i)
END DO
DO j=1,ny1
y1(j)=y(j)
y11(j)=y(j)
END DO
DO k=1,nz1
z1(k)=z(k)
END DO
DO j=1,ny1
DO i=1,nx1
hterain1(i,j) = hterain(i,j)
END DO
END DO
DO k=1,nz1
DO j=1,ny1
DO i=1,nx1
zp1(i,j,k) = zp(i,j,k)
END DO
END DO
END DO
ctrlat1 = ctrlat
ctrlon1 = ctrlon
latitud1 = latitud
xgrdorg1 = xgrdorg
ygrdorg1 = ygrdorg
ELSE
ebc = 3
wbc = 3
sbc = 3
nbc = 3
!
!-----------------------------------------------------------------------
!
! If output grid differs from the input grid, perform spatial
! interpolations.
!
!-----------------------------------------------------------------------
!
xorig1 = xctr1 - (nx1-3)*dx1*0.5
yorig1 = yctr1 - (ny1-3)*dy1*0.5
zorig1 = zorig
DO i=1,nx1
x1(i)=xorig1+(i-2)*dx1
END DO
DO j=1,ny1
y1(j)=xorig1+(j-2)*dy1
END DO
DO k=1,nz1
z1(k)=zorig1+(k-2)*dz1
END DO
xeps = 0.01*dx
yeps = 0.01*dy
IF(x1(1) < x(1)-xeps.OR.x1(nx1) > x(nx)+xeps.OR. &
y1(1) < y(1)-yeps.OR.y1(ny1) > y(ny)+yeps) THEN
WRITE(6,'(3(/5x,a),/5x,2(a,f12.4),2(a,i6),2(a,f12.4)/5x,a)') &
'Sorry, at least part of your new grid is outside the border of', &
'the original grid, please check input parameters', &
'dx1,dy1, nx1, ny1, xctr1 and yctr1. Currently,', &
'dx1=',dx1,', dy1=',dy1,', nx1=',nx1,', ny1=',ny1, &
', xctr1=',xctr1,', yctr1=',yctr1, &
'Job stopped in ARPSINTRP.'
WRITE(6,'(1x,4(a,f12.4)/1x,4(a,f12.4))') &
'x1(1) =',x1(1),' x1(nx1) =',x1(nx1), &
'y1(1) =',y1(1),' y1(ny1) =',y1(ny1), &
'x (1) =',x (1),' x (nx) =',x (nx), &
'y (1) =',y (1),' y (ny) =',y (ny)
STOP 9013
END IF
CALL xytoll(1,1,xctr1,yctr1,ctrlat1,ctrlon1)
PRINT*,'ctrlat1=',ctrlat1,', ctrlon1=',ctrlon1
latitud1 = ctrlat1
CALL xytoll(1,1,xorig1,yorig1,swx,swy) ! Find the lat/lon of
! the new grid origin
CALL setorig( 2, swx,swy ) ! Set the new origin
xgrdorg1 = 0.0
ygrdorg1 = 0.0
DO i=1,nx1
x11(i)=x1(i)-xorig1
END DO
DO j=1,ny1
y11(j)=y1(j)-yorig1
END DO
!
!-----------------------------------------------------------------------
!
! Intepolate terrain to the new grid
!
!-----------------------------------------------------------------------
!
DO i=1,nx1-1
DO j=1,ny1-1
xs1= (x1(i)+x1(i+1))*0.5
ys1= (y1(j)+y1(j+1))*0.5
is = MAX(1, MIN(nx-2, INT((xs1-(x(1)+x(2))*0.5)/dx)+1 ))
js = MAX(1, MIN(ny-2, INT((ys1-(y(1)+y(2))*0.5)/dy)+1 ))
s1=ABS((xs1-(x(is )+x(is+1))*0.5)*(ys1-(y(js )+y(js+1))*0.5))
s2=ABS((xs1-(x(is+1)+x(is+2))*0.5)*(ys1-(y(js )+y(js+1))*0.5))
s3=ABS((xs1-(x(is+1)+x(is+2))*0.5)*(ys1-(y(js+1)+y(js+2))*0.5))
s4=ABS((xs1-(x(is )+x(is+1))*0.5)*(ys1-(y(js+1)+y(js+2))*0.5))
hterain1(i,j) = &
(hterain(is ,js )*s3+hterain(is+1,js )*s4 &
+hterain(is+1,js+1)*s1+hterain(is ,js+1)*s2) &
/(s1+s2+s3+s4)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set up a stretched vertical grid for the output grid.
!
! For strhopt1=1, function y = a+b*x**3 is used to specify dz as a
! function of k.
! For strhopt1=2, function y = c + a*tanh(b*x) is used to specify dz
! as a function of k.
!
!-----------------------------------------------------------------------
!
IF ( strhopt1 == 0 ) THEN
DO k=1,nz1
zp1d1(k) = z1(k)
END DO
ELSE IF ( strhopt1 == 1 .OR.strhopt1 == 2 ) THEN
za = zrefsfc1 + MAX(0.0, MIN(dlayer11, z1(nz1-2)-zrefsfc1 ))
zb = za + MAX(0.0, MIN(dlayer21, z1(nz1-1)-za ))
IF( dlayer11 >= (nz1-3)*dzmin1 ) THEN
WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)') &
'Can not setup a vertical grid with uniform dz=',dzmin1, &
' over the depth of ',dlayer11,' please specify a smaller ', &
'value of dlayer11. Program stopped ARPSINTRP.'
STOP 9014
END IF
CALL strhgrd(nz1,strhopt1,zrefsfc1,za,zb,z1(nz1-1), &
dzmin1,strhtune1, zp1d1,dzp1d1)
ELSE
WRITE(6,'(1x,a,i3,a/)') &
'Invalid vertical grid stretching option, strhopt1 was ',strhopt1, &
'. Program stopped in ARPSINTRP.'
STOP 9015
END IF
!
!-----------------------------------------------------------------------
!
! Physical height of computational grid defined as
!
! Zp=(z-zrefsfc)*(Zm-hterain)/(Zm-zrefsfc)+hterain for z=<Zm.
! ZP=z for z>Zm
!
! where Zm the height at which the grid levels becomes flat.
! Hm < Zm =< Ztop, hm is the height of mounta3din and Ztop the height
! of model top.
!
!-----------------------------------------------------------------------
!
DO k=nz1-1,2,-1
IF(zp1d1(k) <= zflat1) THEN
zflat11 = zp1d1(k)
EXIT
END IF
END DO
300 CONTINUE
zflat11=MAX(MIN(z1(nz1-1),zflat11),zrefsfc1)
DO k=2,nz1-1
IF(zp1d1(k) > zflat11) THEN
DO j=1,ny1-1
DO i=1,nx1-1
zp1(i,j,k)=zp1d1(k)
END DO
END DO
ELSE
DO j=1,ny1-1
DO i=1,nx1-1
zp1(i,j,k)=(zp1d1(k)-zrefsfc1)*(zflat11-hterain1(i,j)) &
/(zflat11-zrefsfc1)+hterain1(i,j)
END DO
END DO
END IF
END DO
DO j=1,ny1-1
DO i=1,nx1-1
zp1(i,j,2)=hterain1(i,j)
zp1(i,j,1)=2.0*zp1(i,j,2)-zp1(i,j,3)
zp1(i,j,nz1)=2.0*zp1(i,j,nz1-1)-zp1(i,j,nz1-2)
END DO
END DO
! CALL jacob(nx1,ny1,nz1,x11,y11,z1,zp1,j11,j21,j31)
END IF ! samgrd
!
!-----------------------------------------------------------------------
!
! Write out terrain data
!
!-----------------------------------------------------------------------
!
CALL getunit( trnchanl )
ternfn = new_runname(1:lfnkey)//".trndata"
lternfn = lfnkey + 8
IF( dirname /= ' ' ) THEN
temchar = ternfn
ternfn = dirname(1:ldirnam)//'/'//temchar
lternfn = lternfn + ldirnam + 1
END IF
CALL fnversn(ternfn, lternfn )
PRINT *, 'Write terrain data to ',ternfn(1:lternfn)
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(ternfn, '-F f77 -N ieee', ierr)
OPEN(trnchanl,FILE=ternfn,FORM='unformatted',STATUS='new')
WRITE(trnchanl) nx1,ny1
WRITE(trnchanl) idummy,mapproj,idummy,idummy,idummy, &
idummy,idummy, idummy,idummy,idummy, &
idummy,idummy, idummy,idummy,idummy, &
idummy,idummy, idummy,idummy,idummy
WRITE(trnchanl) dx1, dy1,ctrlat1,ctrlon1,rdummy, &
rdummy,trulat1,trulat2, trulon,sclfct, &
rdummy, rdummy, rdummy, rdummy,rdummy, &
rdummy, rdummy, rdummy, rdummy,rdummy
WRITE(trnchanl) hterain1
CLOSE ( trnchanl )
CALL retunit ( trnchanl )
!
!-----------------------------------------------------------------------
!
! Open the surface property data file
!
!-----------------------------------------------------------------------
!
lsfcdmp = .false.
IF ( landin == 1 ) THEN ! take care of soil and veg data
sfcoutfl = new_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)
PRINT *, 'Write surface property data in ',sfcoutfl(1:lfn)
CALL getunit ( sfchanl )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(sfcoutfl, '-F f77 -N ieee', ierr)
OPEN (UNIT=sfchanl,FILE=sfcoutfl(1:lfn), STATUS = 'new', &
FORM='unformatted', ACCESS = 'sequential')
WRITE(sfchanl) nx1,ny1
WRITE(sfchanl) mapproj, 1, 1, 1, 1, &
1,idummy, nstyp,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
WRITE(sfchanl) dx1, dy1,ctrlat1,ctrlon1,trulat1, &
trulat2,trulon, sclfct, rdummy, rdummy, &
rdummy,rdummy, rdummy, rdummy, rdummy, &
rdummy,rdummy, rdummy, rdummy, rdummy
lsfcdmp = .true.
END IF
!
!-----------------------------------------------------------------------
!
! Open the grid base state data file for the new grid
!
!-----------------------------------------------------------------------
!
CALL gtbasfn(new_runname(1:lfnkey),dirname,ldirnam,hdmpfmt, &
1,0,basdmpfn,lbasdmpf)
PRINT*
PRINT*,'Output grid/base state file is ', basdmpfn(1:lbasdmpf)
!
!-----------------------------------------------------------------------
!
! Get the IO unit numbers for input grid and base state file
!
!-----------------------------------------------------------------------
!
CALL getunit( nchanl )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(basdmpfn(1:lbasdmpf), '-F f77 -N ieee', ierr)
!
!-----------------------------------------------------------------------
!
! Open the output grdbas file
!
!-----------------------------------------------------------------------
!
OPEN(UNIT=nchanl,FILE=basdmpfn(1:lbasdmpf), &
STATUS='new',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 9002
!
!-----------------------------------------------------------------------
!
! Read in other base state arrays and interpolate them into new grid
!
!-----------------------------------------------------------------------
!
WRITE (cmnt(nocmnt),'(a,i4,a,i4,a,i4)') &
' nx =',nx1,', ny =',ny1,', nz =',nz1
WRITE(nchanl) fmtver
WRITE(nchanl) new_runname
WRITE(nchanl) nocmnt
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
WRITE(nchanl) cmnt(i)
WRITE(6,'(1x,a)') cmnt(i)
END DO
END IF
WRITE(nchanl) time,tmunit
WRITE(nchanl) nx1,ny1,nz1
WRITE(nchanl) 1, 1, 0,mstout1, 0, &
0, 0, 0,landout1,totout, &
idummy, idummy, idummy, mapproj, month, &
day, year, hour, minute, second
WRITE(nchanl) umove,vmove,xgrdorg1,ygrdorg1,trulat1, &
trulat2,trulon,sclfct,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
tstop,thisdmp,latitud1,ctrlat1,ctrlon1
IF ( totout /= 0 ) THEN
WRITE(nchanl) nstyp, idummy, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy
WRITE(nchanl) rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy
END IF
WRITE(6,'(/1x,a/)') &
'Min. and max. of input data interpolated to the new grid:'
WRITE(nchanl) 'x coord r1d1'
WRITE(nchanl) x11
CALL a3dmax0(x11,1,nx1,1,nx1,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(/1x,2(a,e13.6))') 'xmin = ', amin,', xmax =',amax
WRITE(nchanl) 'y coord r1d2'
WRITE(nchanl) y11
CALL a3dmax0(y11,1,ny1,1,ny1,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ymin = ', amin,', ymax =',amax
WRITE(nchanl) 'z coord r1d3'
WRITE(nchanl) z1
CALL a3dmax0(z1,1,nz1,1,nz1,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'zmin = ', amin,', zmax =',amax
WRITE(nchanl) 'zp coor r3d0'
WRITE(nchanl) zp1
CALL a3dmax0(zp1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'zpmin = ', amin,', zpmax =',amax
!
!-----------------------------------------------------------------------
!
! Read in base state arrays, interpolate them into new grid, and
! dump them into the new history file
!
!-----------------------------------------------------------------------
!
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 1 ! U-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'ubar r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,', ubarmax =',amax
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 2 ! V-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'vbar r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,', vbarmax =',amax
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 3 ! W-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'wbar r3d1'
WRITE(nchanl) a3dout
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'ptbar r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,', ptbarmax=',amax
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'pbar r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,', pbarmax =',amax
IF ( mstin == 1 ) THEN
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout1 == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'qvbar r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qvbarmin= ', amin,', qvbarmax=',amax
END IF
END IF
IF ( landin == 1 ) THEN
IF (nstyp == 1) THEN
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) ((i2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array i2din'
CALL dist2d(nx,ny,nx1,ny1,x,y,x1,y1,samgrd,i2din,i2dout)
WRITE(nchanl) 'soiltyp i2d0'
WRITE(nchanl) ((i2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsfcdmp ) WRITE(sfchanl) ((i2dout(i,j),i=1,nx1),j=1,ny1)
ELSE
DO is=1,nstyp
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) &
((i2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array i2din'
CALL dist2d(nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,i2din,i2dout)
WRITE(nchanl) 'soiltyp i2d0'
WRITE(nchanl) ((i2dout(i,j),i=1,nx1),j=1,ny1)
IF( lsfcdmp ) WRITE(sfchanl) ((i2dout(i,j),i=1,nx1),j=1,ny1)
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) &
((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,a2din,a2dout )
WRITE(nchanl) 'stypfrct r2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF( lsfcdmp ) WRITE(sfchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
END DO
END IF
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) i2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array i2din'
CALL dist2d(nx,ny,nx1,ny1,x,y,x1,y1,samgrd,i2din,i2dout)
WRITE(nchanl) 'vegtyp i2d0'
WRITE(nchanl) i2dout
IF ( lsfcdmp ) WRITE(sfchanl) i2dout
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout )
WRITE(nchanl) 'lai r2d0'
WRITE(nchanl) a2dout
IF ( lsfcdmp ) WRITE(sfchanl) a2dout
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout )
WRITE(nchanl) 'roufns r2d0'
WRITE(nchanl) a2dout
IF ( lsfcdmp ) WRITE(sfchanl) a2dout
READ(inch,ERR=9110,END=9120) label
READ(inch,ERR=9110,END=9120) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout )
WRITE(nchanl) 'veg r2d0'
WRITE(nchanl) a2dout
IF ( lsfcdmp ) WRITE(sfchanl) a2dout
END IF
CLOSE (inch)
CLOSE (nchanl)
IF ( lsfcdmp ) CLOSE (sfchanl)
CALL retunit( inch )
CALL retunit( nchanl )
IF ( lsfcdmp ) CALL retunit (sfchanl)
!
!-----------------------------------------------------------------------
!
! Loop over time dependent data files
!
!-----------------------------------------------------------------------
!
grdbas = 0
DO nfile = 1,nhisfile
filename = hisfile(nfile)
lenfil=LEN_trim(filename)
WRITE(6,'(/a,a,a)') &
' Data set ', filename(1:lenfil) ,' to be processed.'
INQUIRE(FILE=filename(1:lenfil), EXIST = fexist )
IF( fexist ) THEN
GO TO 370
END IF
WRITE(6,'(a/a)') &
'File '//filename(1:lenfil)//' does not exist.', &
'Check if compressed file '//filename(1:lenfil)//'.Z' &
//' exists.'
INQUIRE(FILE=filename(1:lenfil)//'.Z',EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( filename(1:lenfil)//'.Z' )
GO TO 370
END IF
WRITE(6,'(a/a)') &
'File '//filename(1:lenfil)//'.Z'//' does not exist.', &
'Check if compressed file '//filename(1:lenfil)//'.gz' &
//' exists.'
INQUIRE(FILE=filename(1:lenfil)//'.gz',EXIST = fexist )
IF( .NOT.fexist ) THEN
CALL uncmprs( filename(1:lenfil)//'.gz' )
GO TO 370
END IF
WRITE(6,'(a/a)') &
'File '//filename(1:lenfil) &
//' or its compressed version not found.', &
'Program stopped.'
STOP
370 CONTINUE ! also continue to read another time recode
! from GrADS file
!
!-----------------------------------------------------------------------
!
! Get the IO unit numbers for input files
!
!-----------------------------------------------------------------------
!
CALL getunit( inch )
!
!-----------------------------------------------------------------------
!
! Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(filename(1:lenfil), '-F f77 -N ieee', ierr)
!
!-----------------------------------------------------------------------
!
! Open the input data file
!
!-----------------------------------------------------------------------
!
OPEN(UNIT=inch,FILE=filename(1:lenfil), &
STATUS='old',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 9003
!
!-----------------------------------------------------------------------
!
! Read all input data header
!
!-----------------------------------------------------------------------
!
READ(inch,ERR=9115,END=9125) fmtverin
IF ( fmtverin == fmtver ) THEN
oldver = 0
ELSE
oldver = 1
END IF
READ(inch,ERR=9115,END=9125) runname
READ(inch,ERR=9115,END=9125) nocmnt
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
READ(inch,ERR=9115,END=9125) cmnt(i)
END DO
END IF
READ(inch,ERR=9115,END=9125) time,tmunit
READ(inch,ERR=9115,END=9125) nxin, nyin, nzin
IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE(6,'(1x,a)') &
' Dimensions in ARPSINTRP inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
WRITE(6,'(1x,a)') &
' Program aborted in ARPSINTRP.'
STOP
END IF
WRITE(6,'(1x,a,f8.1,a,f8.3,a/)') &
'To read data at time ', time, &
' secs = ',(time/60.),' mins.'
READ(inch,ERR=9115,END=9125) &
grdin,basin,varin,mstin,icein, &
trbin, sfcin,rainin,landin,totin, &
tkein,idummy,idummy,mapproj, month, &
day, year, hour, minute, second
IF ( varin == 1 ) THEN
varout1 = varout
ELSE
varout1 = 0
END IF
IF ( mstin == 1 ) THEN
mstout1 = mstout
ELSE
mstout1 = 0
END IF
IF ( icein == 1 ) THEN
iceout1 = iceout
ELSE
iceout1 = 0
END IF
IF ( tkein == 1 ) THEN
tkeout1 = tkeout
ELSE
tkeout1 = 0
END IF
IF ( trbin == 1 ) THEN
trbout1 = trbout
ELSE
trbout1 = 0
END IF
IF ( rainin == 1 ) THEN
rainout1 = rainout
ELSE
rainout1 = 0
END IF
IF ( sfcin == 1 ) THEN
sfcout1 = sfcout
ELSE
sfcout1 = 0
END IF
IF ( landin == 1 ) THEN
landout1 = landout
ELSE
landout1 = 0
END IF
READ(inch,ERR=9115,END=9125) &
umove,vmove,xgrdorg,ygrdorg,trulat1, &
trulat2,trulon,sclfct,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
tstop,thisdmp,latitud,ctrlat,ctrlon
IF ( totin /= 0 ) THEN
!
!-----------------------------------------------------------------------
!
! Read in additional parameters for ARPS history dump 4.0 or later
! version.
!
!-----------------------------------------------------------------------
!
READ(inch,ERR=9115,END=9125) &
nstyp, prcin,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
IF ( nstyp < 1 ) THEN
nstyp = 1
END IF
IF ( prcin == 1 ) THEN
prcout1 = prcout
ELSE
prcout1 = 0
END IF
READ(inch,ERR=9115,END=9125) &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy, &
rdummy,rdummy,rdummy,rdummy,rdummy
END IF
curtim = time
!
!-----------------------------------------------------------------------
!
! Get the history file name for the new grid at time curtim
!
!-----------------------------------------------------------------------
!
runname = new_runname
CALL gtdmpfn(runname(1:lfnkey),dirname, &
ldirnam,curtim,houtfmt,1,0, hdmpfn, ldmpf)
CALL getunit( nchanl )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(hdmpfn(1:ldmpf), '-F f77 -N ieee', ierr)
OPEN(UNIT=nchanl,FILE=hdmpfn(1:ldmpf), &
STATUS='new',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 9004
!
!-----------------------------------------------------------------------
!
! Open the surface property data file
!
!-----------------------------------------------------------------------
!
lsoildmp = .false.
IF ( sfcin == 1 ) THEN ! take care of soil variables
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)
PRINT *, 'Write soil initial data to ',soiloutfl(1:lfn)
CALL getunit ( soilchanl )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(soiloutfl, '-F f77 -N ieee', ierr)
OPEN (UNIT=soilchanl,FILE=soiloutfl(1:lfn), STATUS = 'new', &
FORM='unformatted', ACCESS = 'sequential')
WRITE(soilchanl) nx1,ny1
WRITE(soilchanl) mapproj, 1, 1, 1, 1, &
1,idummy, nstyp,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
WRITE(soilchanl) dx1, dy1,ctrlat1,ctrlon1,trulat1, &
trulat2,trulon, sclfct, rdummy, rdummy, &
rdummy,rdummy, rdummy, rdummy, rdummy, &
rdummy,rdummy, rdummy, rdummy, rdummy
lsoildmp = .true.
END IF
IF ( exbcdmp == 1 ) THEN
!
!-----------------------------------------------------------------------
!
! Open the EXBC file to dump the data
!
!-----------------------------------------------------------------------
!
dmpexbc = .true.
IF ( varin == 0 ) THEN
WRITE(6,'(1x,a,a)') &
'The data file does not contain the variables needed ', &
'for EXBC data file'
dmpexbc = .false.
GO TO 385
END IF
IF ( mstout1 == 0 ) THEN
qcexout = 0
qrexout = 0
qiexout = 0
qsexout = 0
qhexout = 0
ELSE IF ( iceout1 == 0 ) THEN
qiexout = 0
qsexout = 0
qhexout = 0
END IF
CALL ctim2abss( year,month,day,hour,minute,second, itema )
itema = itema + INT(curtim)
CALL abss2ctim( itema,iyr,imon,idy,ihr,imin,isec )
WRITE (ctime,'(i4.4,2i2.2,a,3i2.2)') &
iyr,imon,idy,'.',ihr,imin,isec
exbcfn = runname(1:lfnkey)//'.'//ctime
length = lfnkey + 16
IF( dirname /= ' ' ) THEN
temchar = exbcfn
exbcfn = dirname(1:ldirnam)//'/'//temchar
length = length + ldirnam + 1
END IF
CALL fnversn(exbcfn,length)
WRITE(6,'(1x,a,a)') &
'The external boundary data file is ',exbcfn(1:length)
CALL getunit( exbchanl )
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(exbcfn(1:length), '-F f77 -N ieee', ierr)
OPEN(UNIT=exbchanl,FILE=exbcfn(1:length), &
STATUS='new',FORM='unformatted',IOSTAT=istat)
IF( istat /= 0 ) GO TO 9005
WRITE(exbchanl) nx1,ny1,nz1,dx1,dy1,dz1,ctrlat1,ctrlon1
WRITE(exbchanl) ctime
WRITE(exbchanl) varin,varin,varin,varin,varin, &
mstout1,qcexout,qrexout,qiexout,qsexout, &
qhexout,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy
END IF
385 CONTINUE
!
!-----------------------------------------------------------------------
!
! Write out the header information for new grid
!
!-----------------------------------------------------------------------
!
WRITE(nchanl) fmtver
WRITE(nchanl) new_runname
WRITE(nchanl) nocmnt
WRITE (cmnt(nocmnt),'(a,i4,a,i4,a,i4)') &
' nx =',nx1,', ny =',ny1,', nz =',nz1
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
WRITE(nchanl) cmnt(i)
WRITE(6,'(1x,a)') cmnt(i)
END DO
END IF
WRITE(nchanl) time,tmunit
WRITE(nchanl) nx1,ny1,nz1
!
!-----------------------------------------------------------------------
!
! Write the flags for different data groups.
!
!-----------------------------------------------------------------------
!
WRITE(nchanl) grdout, 0, varout1, mstout1, iceout1, &
trbout1, sfcout1,rainout1,landout1,totout, &
tkeout1, idummy, idummy, mapproj, month, &
day, year, hour, minute, second
rdummy = 0.0
WRITE(nchanl) umove, vmove,xgrdorg1,ygrdorg1, trulat1, &
trulat2, trulon, sclfct, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
tstop, thisdmp,latitud1, ctrlat1, ctrlon1
!
!-----------------------------------------------------------------------
!
! If totout=1, write additional parameters to history dump files.
! This is for ARPS version 4.1.2 or later.
!
!-----------------------------------------------------------------------
!
IF ( totout == 1 ) THEN
WRITE(nchanl) nstyp, prcout1, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy
WRITE(nchanl) rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy
END IF
!
!-----------------------------------------------------------------------
!
! Read in x,y and z at grid cell centers (scalar points).
!
!----------------------------------------------------------------------
!
IF( grdin == 1 .OR. grdbas == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) x
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array x.'
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) y
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array y.'
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) z
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array z.'
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) zp
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array zp.'
END IF ! grdin
!
!-----------------------------------------------------------------------
!
! If grdout=1 or grdbas=1, write out grid variables
!
!-----------------------------------------------------------------------
!
IF(grdout == 1 .OR. grdbas == 1 ) THEN
WRITE(nchanl) 'x coord r1d1'
WRITE(nchanl) x11
WRITE(nchanl) 'y coord r1d2'
WRITE(nchanl) y11
WRITE(nchanl) 'z coord r1d3'
WRITE(nchanl) z1
WRITE(nchanl) 'zp coor r3d0'
WRITE(nchanl) zp1
END IF ! grdout
!
!-----------------------------------------------------------------------
!
! Read in base state fields
!
!----------------------------------------------------------------------
!
IF( basin == 1 .OR. grdbas == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF( mstin == 1) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
END IF
IF (landin == 1) THEN
IF (nstyp == 1) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) &
((i2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array soiltyp for is=',is
ELSE
DO is=1,nstyp
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) &
((i2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array i2din for is=',is
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) &
((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
END DO
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) i2din
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array i2din for is=',is
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Read in arrays for interpolations
!
!----------------------------------------------------------------------
!
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 1 ! U-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'u r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'umin = ', amin,', umax =',amax
IF ( dmpexbc ) WRITE(exbchanl) a3dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 2 ! V-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'v r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vmin = ', amin,', vmax =',amax
IF ( dmpexbc ) WRITE(exbchanl) a3dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 3 ! W-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'w r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'wmin = ', amin,', wmax =',amax
IF ( dmpexbc ) WRITE(exbchanl) a3dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'pt r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ptmin = ', amin,', ptmax =',amax
IF ( dmpexbc ) WRITE(exbchanl) a3dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'p r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'pmin = ', amin,', pmax =',amax
IF ( dmpexbc ) WRITE(exbchanl) a3dout
IF( mstin == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'qv r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qvmin = ', amin,', qvmax =',amax
IF ( dmpexbc ) WRITE(exbchanl) a3dout
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'qc r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qcmin = ', amin,', qcmax =',amax
IF ( dmpexbc .AND. qcexout == 1 ) WRITE(exbchanl) a3dout
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'qr r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qrmin = ', amin,', qrmax =',amax
IF ( dmpexbc .AND. qrexout == 1 ) WRITE(exbchanl) a3dout
END IF
IF( rainin == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
IF ( mstout == 1 .AND. rainout == 1 ) THEN
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,a2din,a2dout )
WRITE(nchanl) 'raing r3d1'
WRITE(nchanl) a2dout
CALL a3dmax0(a2dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'raingmin= ', amin,', raingmax=',amax
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
IF ( mstout == 1 .AND. rainout == 1 ) THEN
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,a2din,a2dout )
WRITE(nchanl) 'rainc r3d1'
WRITE(nchanl) a2dout
CALL a3dmax0(a2dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'raincmin= ', amin,', raincmax=',amax
END IF
END IF
IF( prcin == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
IF ( mstout == 1 .AND. prcout == 1 ) THEN
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,a2din,a2dout )
WRITE(nchanl) 'prcrt1 r3d1'
WRITE(nchanl) a2dout
CALL a3dmax0(a2dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'prcr1min= ', amin,', prcr1max=',amax
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
IF ( mstout == 1 .AND. prcout == 1 ) THEN
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,a2din,a2dout )
WRITE(nchanl) 'prcrt2 r3d1'
WRITE(nchanl) a2dout
CALL a3dmax0(a2dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'prcr2min= ', amin,', prcr1max=',amax
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
IF ( mstout == 1 .AND. prcout == 1 ) THEN
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,a2din,a2dout )
WRITE(nchanl) 'prcrt3 r3d1'
WRITE(nchanl) a2dout
CALL a3dmax0(a2dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'prcr3min= ', amin,', prcr1max=',amax
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a2din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
IF ( mstout == 1 .AND. prcout == 1 ) THEN
CALL intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, &
samgrd,a2din,a2dout )
WRITE(nchanl) 'prcrt4 r3d1'
WRITE(nchanl) a2dout
CALL a3dmax0(a2dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'prcr4min= ', amin,', prcr1max=',amax
END IF
END IF
IF( icein == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 .AND. iceout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'qi r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qimin = ', amin,', qimax =',amax
IF ( dmpexbc .AND. qiexout == 1 ) WRITE(exbchanl) a3dout
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 .AND. iceout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'qs r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qsmin = ', amin,', qsmax =',amax
IF ( dmpexbc .AND. qsexout == 1 ) WRITE(exbchanl) a3dout
END IF
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 .AND. iceout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'qh r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qhmin = ', amin,', qhmax =',amax
IF ( dmpexbc .AND. qhexout == 1 ) WRITE(exbchanl) a3dout
END IF
END IF
END IF
IF( tkein == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 .AND. iceout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'tke r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'tkemin = ', amin,', tkemax =',amax
END IF
END IF
IF( trbin == 1 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 .AND. iceout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'kmh r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'kmhmin = ', amin,', kmhmax =',amax
END IF
IF ( oldver == 0 ) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) a3din
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a3din.'
IF ( mstout == 1 .AND. iceout == 1 ) THEN
stgr = 4 ! S-points
CALL intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, &
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
WRITE(nchanl) 'kmv r3d1'
WRITE(nchanl) a3dout
CALL a3dmax0(a3dout,1,nx1,1,nx1-1,1,ny1,1,ny1-1, &
1,nz1,1,nz1-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'kmvmin = ', amin,', kmvmax =',amax
END IF
END IF
END IF
IF( sfcin == 1 ) THEN
IF (nstyp == 1) THEN
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1,samgrd,a2din,a2dout)
WRITE(nchanl) 'tsfc i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1,samgrd,a2din,a2dout)
WRITE(nchanl) 'tsoil i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1,samgrd,a2din,a2dout)
WRITE(nchanl) 'wetsfc i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1,samgrd,a2din,a2dout)
WRITE(nchanl) 'wetdp i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a)') &
'Field ',label,' was read into array a2din.'
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1,samgrd,a2din,a2dout)
WRITE(nchanl) 'wetcanp i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp ) WRITE(soilchanl) a2dout
ELSE
DO is=0,nstyp
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout)
WRITE(nchanl) 'tsfc i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp .AND. is == 0 ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout)
WRITE(nchanl) 'tsoil i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp .AND. is == 0 ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout)
WRITE(nchanl) 'wetsfc i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp .AND. is == 0 ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout)
WRITE(nchanl) 'wetdp i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp .AND. is == 0 ) WRITE(soilchanl) a2dout
READ(inch,ERR=9115,END=9125) label
READ(inch,ERR=9115,END=9125) ((a2din(i,j),i=1,nx),j=1,ny)
WRITE(6,'(1x,a,a12,a,i2)') &
'Field ',label,' was read into array a2din for is=',is
CALL intrp2d(nx,ny,nx1,ny1,x,y,x1,y1, samgrd,a2din,a2dout)
WRITE(nchanl) 'wetcanp i2d0'
WRITE(nchanl) ((a2dout(i,j),i=1,nx1),j=1,ny1)
IF ( lsoildmp .AND. is == 0 ) WRITE(soilchanl) a2dout
END DO
END IF
END IF
CLOSE ( inch )
CLOSE ( nchanl )
IF ( dmpexbc ) CLOSE ( exbchanl )
IF ( lsoildmp ) CLOSE ( soilchanl )
CALL retunit ( inch )
CALL retunit ( nchanl )
IF ( dmpexbc ) CALL retunit ( exbchanl )
IF ( lsoildmp ) CALL retunit ( soilchanl )
END DO
STOP
9000 CONTINUE
WRITE(6,'(1x,a,i2,/1x,a)') &
'Namelist read error. Job stopped in ARPSINTRP.'
STOP 9000
!
!-----------------------------------------------------------------------
!
! Error during read
!
!----------------------------------------------------------------------
!
9001 CONTINUE
WRITE(6,'(/a,a/)') &
' Error in openning file ',grdbasfn(1:lengbf)
STOP 9001
9002 CONTINUE
WRITE(6,'(/a,a/)') &
' Error in openning file ',basdmpfn(1:lbasdmpf)
STOP 9002
9003 CONTINUE
WRITE(6,'(/a,a/)') &
' Error in openning file ',filename(1:lenfil)
STOP 9003
9004 CONTINUE
WRITE(6,'(/a,a/)') &
' Error in openning file ',hdmpfn(1:ldmpf)
STOP 9004
9005 CONTINUE
WRITE(6,'(/a,a/)') &
' Error in openning file ',exbcfn(1:length)
STOP 9005
9110 CONTINUE
WRITE(6,'(/a,a/)') ' Error in reading file ',basdmpfn(1:lbasdmpf)
STOP 9110
9115 CONTINUE
WRITE(6,'(/a,a/)') ' Error in reading file ',filename(1:lenfil)
STOP 9115
!
!-----------------------------------------------------------------------
!
! End-of-file during read
!
!----------------------------------------------------------------------
!
9120 CONTINUE
WRITE(6,'(/a,a/)') &
' Unexpected End of File reached in reading file ', &
basdmpfn(1:lbasdmpf)
STOP 9120
9125 CONTINUE
WRITE(6,'(/a,a/)') &
' Unexpected End of File reached in reading file ', &
filename(1:lenfil)
STOP 9125
END PROGRAM arpsintrp_ls
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INTRP3D ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE intrp3d( nx,ny,nz,nx1,ny1,nz1,x,y,z,zp,x1,y1,z1,zp1, & 20,4
samgrd,stgr,a3din,a3dout,tem3d,tem3d1 )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Interpolate a 3-d array from an ARPS grid into a new ARPS grid.
!
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 08/05/1997
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Parameters for the utput grid.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz
INTEGER :: nx1,ny1,nz1
REAL :: x(nx)
REAL :: y(ny)
REAL :: z(nz)
REAL :: zp(nx,ny,nz)
REAL :: x1(nx1)
REAL :: y1(ny1)
REAL :: z1(nz1)
REAL :: zp1(nx1,ny1,nz1)
REAL :: a3din(nx,ny,nz)
REAL :: a3dout(nx1,ny1,nz1)
REAL :: tem3d(nx,ny,nz)
REAL :: tem3d1(nx1,ny1,nz1)
INTEGER :: samgrd,stgr
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
REAL :: s1,s2,s3,s4,sgrdinv
REAL :: vgridinv,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8
REAL :: xu1,yu1,zu1,xv1,yv1,zv1,xw1,yw1,zw1,xs1,ys1,zs1
INTEGER :: iu,ju,ku,iv,jv,kv,iw,jw,kw,is,js,ks
INTEGER :: kk
REAL :: zupper,zlower
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'bndry.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( samgrd == 1 ) THEN
DO k=1,nz1
DO j=1,ny1
DO i=1,nx1
a3dout(i,j,k) = a3din(i,j,k)
END DO
END DO
END DO
RETURN
END IF
IF ( stgr == 1 ) THEN ! for u-points
DO k=1,nz-1
DO j=1,ny-1
DO i=2,nx-1
tem3d(i,j,k) = 0.25*(zp(i ,j,k)+zp(i ,j,k+1) &
+zp(i-1,j,k)+zp(i-1,j,k+1))
END DO
END DO
END DO
CALL bcsu(nx,ny,nz,1,ny,1,nz,ebc,wbc,tem3d)
DO k=1,nz1-1
DO j=1,ny1-1
DO i=2,nx1-1
tem3d1(i,j,k) = 0.25*(zp1(i ,j,k)+zp1(i ,j,k+1) &
+zp1(i-1,j,k)+zp1(i-1,j,k+1))
END DO
END DO
END DO
CALL bcsu(nx1,ny1,nz1,1,ny1,1,nz1,ebc,wbc,tem3d1)
DO i=1,nx1
DO j=1,ny1-1
DO k=1,nz1-1
xu1=x1(i)
yu1=(y1(j+1)+y1(j))*0.5
zu1= tem3d1(i,j,k)
iu = MAX(1, MIN(nx-1, INT((xu1-x(1))/dx)+1 ))
ju = MAX(1, MIN(ny-2, INT((yu1-(y(1)+y(2))*0.5)/dy)+1 ))
s1 = ABS((xu1-x(iu ))*(yu1-(y(ju )+y(ju+1))*0.5))
s2 = ABS((xu1-x(iu+1))*(yu1-(y(ju )+y(ju+1))*0.5))
s3 = ABS((xu1-x(iu+1))*(yu1-(y(ju+1)+y(ju+2))*0.5))
s4 = ABS((xu1-x(iu ))*(yu1-(y(ju+1)+y(ju+2))*0.5))
sgrdinv = 1.0/(s1+s2+s3+s4)
ku = 1
DO kk=nz-2,1,-1
IF(zu1 >= (tem3d(iu ,ju ,kk)*s3+tem3d(iu+1,ju ,kk)*s4 &
+tem3d(iu+1,ju+1,kk)*s1+tem3d(iu ,ju+1,kk)*s2) &
*sgrdinv) THEN
ku = kk
EXIT
END IF
END DO
130 CONTINUE
zlower = (tem3d(iu ,ju ,ku)*s3+tem3d(iu+1,ju ,ku)*s4 &
+tem3d(iu+1,ju+1,ku)*s1+tem3d(iu ,ju+1,ku)*s2) &
*sgrdinv
zupper = (tem3d(iu ,ju ,ku+1)*s3+tem3d(iu+1,ju ,ku+1)*s4 &
+tem3d(iu+1,ju+1,ku+1)*s1+tem3d(iu ,ju+1,ku+1)*s2) &
*sgrdinv
tmp1 = ABS(s1*(zu1-zlower))
tmp2 = ABS(s2*(zu1-zlower))
tmp3 = ABS(s3*(zu1-zlower))
tmp4 = ABS(s4*(zu1-zlower))
tmp5 = ABS(s1*(zu1-zupper))
tmp6 = ABS(s2*(zu1-zupper))
tmp7 = ABS(s3*(zu1-zupper))
tmp8 = ABS(s4*(zu1-zupper))
vgridinv=1.0/(tmp1+tmp2+tmp3+tmp4+tmp5+tmp6+tmp7+tmp8)
a3dout(i,j,k) = &
(a3din(iu ,ju ,ku )*tmp7+a3din(iu+1,ju ,ku )*tmp8 &
+a3din(iu+1,ju+1,ku )*tmp5+a3din(iu ,ju+1,ku )*tmp6 &
+a3din(iu ,ju ,ku+1)*tmp3+a3din(iu+1,ju ,ku+1)*tmp4 &
+a3din(iu+1,ju+1,ku+1)*tmp1+a3din(iu ,ju+1,ku+1)*tmp2) &
*vgridinv
END DO
END DO
END DO
ELSE IF ( stgr == 2 ) THEN
DO k=1,nz-1
DO j=2,ny-1
DO i=1,nx-1
tem3d(i,j,k) = 0.25*(zp(i,j ,k)+zp(i,j ,k+1) &
+zp(i,j-1,k)+zp(i,j-1,k+1))
END DO
END DO
END DO
CALL bcsu(nx,ny,nz,1,nx,1,nz,nbc,sbc,tem3d)
DO k=1,nz1-1
DO j=2,ny1-1
DO i=1,nx1-1
tem3d1(i,j,k) = 0.25*(zp1(i,j ,k)+zp1(i,j ,k+1) &
+zp1(i,j-1,k)+zp1(i,j-1,k+1))
END DO
END DO
END DO
CALL bcsv(nx1,ny1,nz1,1,nx1,1,nz1,nbc,sbc,tem3d1)
DO i=1,nx1-1
DO j=1,ny1
DO k=1,nz1-1
xv1=(x1(i+1)+x1(i))*0.5
yv1= y1(j)
zv1= tem3d1(i,j,k)
iv = MAX(1, MIN(nx-2, INT((xv1-(x(1)+x(2))*0.5)/dx)+1 ))
jv = MAX(1, MIN(ny-1, INT((yv1- y(1))/dy)+1 ))
s1 = ABS((xv1-(x(iv )+x(iv+1))*0.5)*(yv1-y(jv )))
s2 = ABS((xv1-(x(iv+1)+x(iv+2))*0.5)*(yv1-y(jv )))
s3 = ABS((xv1-(x(iv+1)+x(iv+2))*0.5)*(yv1-y(jv+1)))
s4 = ABS((xv1-(x(iv )+x(iv+1))*0.5)*(yv1-y(jv+1)))
sgrdinv = 1.0/(s1+s2+s3+s4)
kv = 1
DO kk=nz-2,1,-1
IF(zv1 >= (tem3d(iv ,jv ,kk)*s3+tem3d(iv+1,jv ,kk)*s4 &
+tem3d(iv+1,jv+1,kk)*s1+tem3d(iv ,jv+1,kk)*s2) &
*sgrdinv) THEN
kv = kk
EXIT
END IF
END DO
230 CONTINUE
zlower = (tem3d(iv ,jv ,kv)*s3+tem3d(iv+1,jv ,kv)*s4 &
+tem3d(iv+1,jv+1,kv)*s1+tem3d(iv ,jv+1,kv)*s2) &
*sgrdinv
zupper = (tem3d(iv ,jv ,kv+1)*s3+tem3d(iv+1,jv ,kv+1)*s4 &
+tem3d(iv+1,jv+1,kv+1)*s1+tem3d(iv ,jv+1,kv+1)*s2) &
*sgrdinv
tmp1 = ABS(s1*(zv1-zlower))
tmp2 = ABS(s2*(zv1-zlower))
tmp3 = ABS(s3*(zv1-zlower))
tmp4 = ABS(s4*(zv1-zlower))
tmp5 = ABS(s1*(zv1-zupper))
tmp6 = ABS(s2*(zv1-zupper))
tmp7 = ABS(s3*(zv1-zupper))
tmp8 = ABS(s4*(zv1-zupper))
vgridinv=1.0/(tmp1+tmp2+tmp3+tmp4+tmp5+tmp6+tmp7+tmp8)
a3dout(i,j,k) = &
(a3din(iv ,jv ,kv )*tmp7+a3din(iv+1,jv ,kv )*tmp8 &
+a3din(iv+1,jv+1,kv )*tmp5+a3din(iv ,jv+1,kv )*tmp6 &
+a3din(iv ,jv ,kv+1)*tmp3+a3din(iv+1,jv ,kv+1)*tmp4 &
+a3din(iv+1,jv+1,kv+1)*tmp1+a3din(iv ,jv+1,kv+1)*tmp2) &
*vgridinv
END DO
END DO
END DO
ELSE IF ( stgr == 3 ) THEN
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,nz1
xw1= (x1(i)+x1(i+1))*0.5
yw1= (y1(j)+y1(j+1))*0.5
zw1= zp1(i,j,k)
iw = MAX(1, MIN(nx-2, INT((xw1-(x(1)+x(2))*0.5)/dx)+1 ))
jw = MAX(1, MIN(ny-2, INT((yw1-(y(1)+y(2))*0.5)/dy)+1 ))
s1=ABS((xw1-(x(iw )+x(iw+1))*0.5) &
*(yw1-(y(jw )+y(jw+1))*0.5))
s2=ABS((xw1-(x(iw+1)+x(iw+2))*0.5) &
*(yw1-(y(jw )+y(jw+1))*0.5))
s3=ABS((xw1-(x(iw+1)+x(iw+2))*0.5) &
*(yw1-(y(jw+1)+y(jw+2))*0.5))
s4=ABS((xw1-(x(iw )+x(iw+1))*0.5) &
*(yw1-(y(jw+1)+y(jw+2))*0.5))
sgrdinv = 1.0/(s1+s2+s3+s4)
kw = 1
DO kk=nz-1,1,-1
IF(zw1 >= (zp(iw ,jw ,kk)*s3+zp(iw+1,jw ,kk)*s4 &
+zp(iw+1,jw+1,kk)*s1+zp(iw ,jw+1,kk)*s2) &
*sgrdinv) THEN
kw = kk
EXIT
END IF
END DO
320 CONTINUE
zlower = (zp(iw ,jw ,kw)*s3+zp(iw+1,jw ,kw)*s4 &
+zp(iw+1,jw+1,kw)*s1+zp(iw ,jw+1,kw)*s2) &
*sgrdinv
zupper = (zp(iw ,jw ,kw+1)*s3+zp(iw+1,jw ,kw+1)*s4 &
+zp(iw+1,jw+1,kw+1)*s1+zp(iw ,jw+1,kw+1)*s2) &
*sgrdinv
tmp1 = ABS(s1*(zw1-zlower))
tmp2 = ABS(s2*(zw1-zlower))
tmp3 = ABS(s3*(zw1-zlower))
tmp4 = ABS(s4*(zw1-zlower))
tmp5 = ABS(s1*(zw1-zupper))
tmp6 = ABS(s2*(zw1-zupper))
tmp7 = ABS(s3*(zw1-zupper))
tmp8 = ABS(s4*(zw1-zupper))
vgridinv = 1.0/(tmp1+tmp2+tmp3+tmp4+tmp5+tmp6+tmp7+tmp8)
a3dout(i,j,k) = &
(a3din(iw ,jw ,kw )*tmp7+a3din(iw+1,jw ,kw )*tmp8 &
+a3din(iw+1,jw+1,kw )*tmp5+a3din(iw ,jw+1,kw )*tmp6 &
+a3din(iw ,jw ,kw+1)*tmp3+a3din(iw+1,jw ,kw+1)*tmp4 &
+a3din(iw+1,jw+1,kw+1)*tmp1+a3din(iw ,jw+1,kw+1)*tmp2) &
*vgridinv
END DO
END DO
END DO
ELSE IF ( stgr == 4 ) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem3d(i,j,k) = 0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
END DO
END DO
DO k=1,nz1-1
DO j=1,ny1-1
DO i=1,nx1-1
tem3d1(i,j,k) = 0.5*(zp1(i,j,k)+zp1(i,j,k+1))
END DO
END DO
END DO
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,nz1-1
xs1= (x1(i)+x1(i+1))*0.5
ys1= (y1(j)+y1(j+1))*0.5
zs1= tem3d1(i,j,k)
is = MAX(1, MIN(nx-2, INT((xs1-(x(1)+x(2))*0.5)/dx)+1 ))
js = MAX(1, MIN(ny-2, INT((ys1-(y(1)+y(2))*0.5)/dy)+1 ))
s1 =ABS((xs1-(x(is )+x(is+1))*0.5) &
*(ys1-(y(js )+y(js+1))*0.5))
s2 =ABS((xs1-(x(is+1)+x(is+2))*0.5) &
*(ys1-(y(js )+y(js+1))*0.5))
s3 =ABS((xs1-(x(is+1)+x(is+2))*0.5) &
*(ys1-(y(js+1)+y(js+2))*0.5))
s4 =ABS((xs1-(x(is )+x(is+1))*0.5) &
*(ys1-(y(js+1)+y(js+2))*0.5))
sgrdinv = 1.0/(s1+s2+s3+s4)
ks = 1
DO kk=nz-2,1,-1
IF(zs1 >= (tem3d(is ,js ,kk)*s3+tem3d(is+1,js ,kk)*s4 &
+tem3d(is+1,js+1,kk)*s1+tem3d(is ,js+1,kk)*s2) &
*sgrdinv ) THEN
ks = kk
EXIT
END IF
END DO
430 CONTINUE
zlower = (tem3d(is ,js ,ks)*s3+tem3d(is+1,js ,ks)*s4 &
+tem3d(is+1,js+1,ks)*s1+tem3d(is ,js+1,ks)*s2) &
*sgrdinv
zupper = (tem3d(is ,js ,ks+1)*s3+tem3d(is+1,js ,ks+1)*s4 &
+tem3d(is+1,js+1,ks+1)*s1+tem3d(is ,js+1,ks+1)*s2) &
*sgrdinv
tmp1 = ABS(s1*(zs1-zlower))
tmp2 = ABS(s2*(zs1-zlower))
tmp3 = ABS(s3*(zs1-zlower))
tmp4 = ABS(s4*(zs1-zlower))
tmp5 = ABS(s1*(zs1-zupper))
tmp6 = ABS(s2*(zs1-zupper))
tmp7 = ABS(s3*(zs1-zupper))
tmp8 = ABS(s4*(zs1-zupper))
vgridinv=1.0/(tmp1+tmp2+tmp3+tmp4+tmp5+tmp6+tmp7+tmp8)
a3dout(i,j,k) = &
(a3din(is ,js ,ks )*tmp7+a3din(is+1,js ,ks )*tmp8 &
+a3din(is+1,js+1,ks )*tmp5+a3din(is ,js+1,ks )*tmp6 &
+a3din(is ,js ,ks+1)*tmp3+a3din(is+1,js ,ks+1)*tmp4 &
+a3din(is+1,js+1,ks+1)*tmp1+a3din(is ,js+1,ks+1)*tmp2) &
*vgridinv
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE intrp3d
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INTRP2D ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE intrp2d( nx,ny,nx1,ny1,x,y,x1,y1, & 20
samgrd,a2din,a2dout )
!
!-----------------------------------------------------------------------
!
! Interpolate a 2-d array from an ARPS grid into a new ARPS grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 08/05/1997
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Parameters for the utput grid.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny
INTEGER :: nx1,ny1
REAL :: x(nx)
REAL :: y(ny)
REAL :: x1(nx1)
REAL :: y1(ny1)
REAL :: a2din(nx,ny)
REAL :: a2dout(nx1,ny1)
INTEGER :: samgrd
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j
REAL :: s1,s2,s3,s4,sgrdinv
REAL :: xs1,ys1
INTEGER :: is,js
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( samgrd == 1 ) THEN
DO j=1,ny1
DO i=1,nx1
a2dout(i,j) = a2din(i,j)
END DO
END DO
RETURN
END IF
DO i=1,nx1-1
DO j=1,ny1-1
xs1= (x1(i)+x1(i+1))*0.5
ys1= (y1(j)+y1(j+1))*0.5
is = MAX(1, MIN(nx-2, INT((xs1-(x(1)+x(2))*0.5)/dx)+1 ))
js = MAX(1, MIN(ny-2, INT((ys1-(y(1)+y(2))*0.5)/dy)+1 ))
s1=ABS((xs1-(x(is )+x(is+1))*0.5)*(ys1-(y(js )+y(js+1))*0.5))
s2=ABS((xs1-(x(is+1)+x(is+2))*0.5)*(ys1-(y(js )+y(js+1))*0.5))
s3=ABS((xs1-(x(is+1)+x(is+2))*0.5)*(ys1-(y(js+1)+y(js+2))*0.5))
s4=ABS((xs1-(x(is )+x(is+1))*0.5)*(ys1-(y(js+1)+y(js+2))*0.5))
sgrdinv = 1.0/(s1+s2+s3+s4)
a2dout(i,j) = &
(a2din(is ,js )*s3+a2din(is+1,js )*s4 &
+a2din(is+1,js+1)*s1+a2din(is ,js+1)*s2) &
*sgrdinv
END DO
END DO
RETURN
END SUBROUTINE intrp2d
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DIST2D ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE dist2d( nx,ny,nx1,ny1,x,y,x1,y1, & 3
samgrd,i2din,i2dout )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Re-distribute a 2-d integer array from an ARPS grid into another
! grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 08/05/1997
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Parameters for the utput grid.
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny
INTEGER :: nx1,ny1
REAL :: x(nx)
REAL :: y(ny)
REAL :: x1(nx1)
REAL :: y1(ny1)
INTEGER :: i2din(nx,ny)
INTEGER :: i2dout(nx1,ny1)
INTEGER :: samgrd
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j
REAL :: xs1,ys1
INTEGER :: is,js
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF ( samgrd == 1 ) THEN
DO j=1,ny1
DO i=1,nx1
i2dout(i,j) = i2din(i,j)
END DO
END DO
RETURN
END IF
DO i=1,nx1-1
DO j=1,ny1-1
xs1= (x1(i)+x1(i+1))*0.5
ys1= (y1(j)+y1(j+1))*0.5
is = MAX(1, MIN(nx-2, INT((xs1-(x(1)+x(2))*0.5)/dx)+1 ))
js = MAX(1, MIN(ny-2, INT((ys1-(y(1)+y(2))*0.5)/dy)+1 ))
i2dout(i,j) = i2din(is,js)
END DO
END DO
RETURN
END SUBROUTINE dist2d