PROGRAM ARPSSUBDOMAIN,10
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSSUBDOMAIN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! (12/1/2002)
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'exbc.inc'
REAL, allocatable :: x (:),y (:),z (:)
REAL, allocatable :: x1(:),y1(:),z1(:)
REAL, allocatable :: xs (:), ys (:) ! x,y coord for scalar points
REAL, allocatable :: xs1(:), ys1(:) ! x,y coord for scalar points
INTEGER, allocatable :: isx(:),jsy(:),iux(:),jvy(:)
REAL, allocatable:: wgtsx(:),wgtsy(:),wgtux(:),wgtvy(:),temx1yz(:,:)
REAL, allocatable :: stypfrct(:,:,:)
INTEGER, allocatable :: soiltyp(:,:,:)
REAL, allocatable :: tem1(:,:,:)
REAL, allocatable :: tem11(:,:,:)
INTEGER(2), allocatable :: tem2(:,:,:)
INTEGER(2), allocatable :: tem21(:,:,:)
INTEGER :: nx,ny,nz,nzsoil ! Grid dimensions.
INTEGER :: nstyps ! Maximum number of soil types.
INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf
PARAMETER (nhisfile_max=1000)
CHARACTER (LEN=256) :: grdbasfn,hisfile(nhisfile_max)
INTEGER :: nxbgn, nxend, nybgn, nyend, nzbgn, nzend, intrp_opt
INTEGER :: nx1,ny1,nz1,nzsoil1,nstyp1
integer :: savespace
REAL :: xctr1,yctr1,dx1,dy1,xorig1,yorig1
NAMELIST /output_dims/ intrp_opt, &
nxbgn, nxend, nybgn, nyend, nzbgn, nzend, &
nx1,ny1,xctr1,yctr1,dx1,dy1, savespace, &
dz,strhopt,dzmin,zrefsfc,dlayer1,dlayer2,strhtune,zflat
NAMELIST /jobname/ runname
INTEGER :: istatus
NAMELIST /output/ dirname,hdmpfmt,grbpkbit,hdfcompr, &
grdout,basout,varout,mstout,rainout,prcout,iceout, &
tkeout,trbout,sfcout,landout,totout,filcmprs,readyfl, &
exbcdmp,qcexout,qrexout,qiexout,qsexout,qhexout,ngbrz,rayklow, &
terndmp
INTEGER :: i,j,k,ireturn
REAL :: time
LOGICAL :: iexist
INTEGER :: nfile
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
WRITE(6,'(/9(/5x,a)/)') &
'###############################################################', &
'###############################################################', &
'### ###', &
'### Welcome to ARPSSUBDOMAIN ###', &
'### This program read, interpolate then write out ###', &
'### variables in ARPS history for files one variable ###', &
'### at a time to minimize memory requirement. ###', &
'### ###', &
'###############################################################', &
'###############################################################'
!
!-----------------------------------------------------------------------
! Get the names of the input data files.
!-----------------------------------------------------------------------
!
hinfmt = 3
CALL get_input_file_names
(hinfmt,grdbasfn,hisfile,nhisfile)
lengbf = len_trim(grdbasfn)
WRITE(6,'(/a,a)')' The grid/base name is ', trim(grdbasfn)
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nzsoil,nstyps, ireturn)
Print*,'nx,ny,nz of input data were ', nx,ny,nz
nzsoil1 = nzsoil
nstyp1 = nstyp
allocate(x(nx ),stat=istatus)
x = 0.0
allocate(y(ny ),stat=istatus)
y = 0.0
allocate(xs(nx ),stat=istatus)
xs = 0.0
allocate(ys(ny ),stat=istatus)
ys = 0.0
allocate(z(nz ),stat=istatus)
y = 0.0
! allocate(zpsoil(nx,ny,nzsoil),stat=istatus)
! zpsoil = 0.0
CALL get_gridinfo_from_hdf(grdbasfn(1:lengbf),nx,ny,nz,x,y,z,ireturn)
! Print*,'x,y,z,zp of input data read in.'
! print*,'x(1 )=',x(1)
! print*,'x(nx)=',x(nx)
! print*,'y(1 )=',y(1)
! print*,'y(ny)=',y(ny)
dx = x(2) - x(1)
dy = y(2) - y(1)
IF (nstyps <= 0) nstyps = 1
nstyp = nstyps
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz , &
'nzsoil=',nzsoil
intrp_opt = 0
savespace = 1
READ(5,output_dims, END=100)
WRITE(6,'(/a,a/)') &
'NAMELIST block output_dims successfully read.'
IF( strhopt == 0.AND.dzmin /= dz ) THEN
WRITE(6,'(5x,a)') &
'For non-stretched case, dzmin was reset to dz.'
dzmin = dz
END IF
if( intrp_opt == 0 ) then
nx1 = nxend - nxbgn +1
ny1 = nyend - nybgn +1
else
nxbgn = 1
nxend = nx1
nybgn = 1
nyend = ny1
endif
nz1 = nzend - nzbgn +1
WRITE(6,'(4(a,i5))') 'nx1 =',nx1,', ny1=',ny1,', nz1=',nz1
allocate(x1(nx1),stat=istatus)
x1 = 0.0
allocate(y1(ny1),stat=istatus)
y1 = 0.0
allocate(xs1(nx1),stat=istatus)
xs1 = 0.0
allocate(ys1(ny1),stat=istatus)
ys1 = 0.0
allocate(z1(nz1),stat=istatus)
y1 = 0.0
! allocate(zpsoil1(nx1,ny1,nzsoil),stat=istatus)
! zpsoil1 = 0.0
allocate(stypfrct(nx,ny,nstyps),stat=istatus)
allocate(soiltyp(nx,ny,nstyps),stat=istatus)
READ (5,jobname,END=100)
WRITE(6,'(/a/)') 'Sucessfully read namelist block JOBNAME.'
WRITE(6,'(/2x,a,a)') 'The name of this run is: ', runname
!
!-----------------------------------------------------------------------
! Set the control parameters for output:
!-----------------------------------------------------------------------
!
READ (5,output,END=100)
ldirnam=LEN_trim(dirname)
print*,'ldirnam, dirname', ldirnam, trim(dirname)
IF( dirname(1:ldirnam) /= ' ') THEN
CALL inquiredir(dirname(1:ldirnam),iexist)
IF( .NOT.iexist ) THEN
WRITE(6,'(5x,a,2(/5x,a))') &
'Specified output directory '//dirname(1:ldirnam)// &
' not found.','It will be created by the program.'
CALL unixcmd( 'mkdir -p '//dirname(1:ldirnam) )
END IF
END IF
WRITE(6,'(5x,a)') &
'Output files will be in directory '//dirname(1:ldirnam)//'.'
WRITE(6,'(/a/)') 'Output control parameters read in are:'
WRITE(6,output)
ALLOCATE(tem1(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nx*ny*nz,'tem1')
tem1 = 0.0
ALLOCATE(tem2(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nx*ny*nz,'tem2')
tem2 = 0.0
ALLOCATE(tem11(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nx1*ny1*nz1,'tem11')
tem11 = 0.0
ALLOCATE(tem21(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nx1*ny1*nz1,'tem21')
tem21 = 0.0
if( intrp_opt == 0 ) then
DO i=1,nx1
x1(i) = x(nxbgn+i-1)
END DO
DO j=1,ny1
y1(j) =y(nybgn+j-1)
END DO
else
xorig1 = xctr1 - (nx1-3)*dx1*0.5
yorig1 = yctr1 - (ny1-3)*dy1*0.5
DO i=1,nx1
x1(i) =xorig1+(i-2.0)*dx1
xs1(i)=xorig1+(i-1.5)*dx1
END DO
DO j=1,ny1
y1 (j) =yorig1+(j-2.0)*dy1
ys1(j) =yorig1+(j-1.5)*dy1
END DO
DO i=1,nx-1
xs(i)=0.5*(x(i)+x(i+1))
END DO
DO j=1,ny-1
ys(j)=0.5*(y(j)+y(j+1))
END DO
print*,'xorig1,xs1(1),xs1(nx1-1)=',xorig1,xs1(1),xs1(nx1-1)
print*,'yorig1,ys1(1),ys1(ny1-1)=',yorig1,ys1(1),ys1(ny1-1)
allocate(isx(nx1),stat=istatus)
isx = 0.0
allocate(jsy(ny1),stat=istatus)
jsy = 0.0
allocate(iux(nx1),stat=istatus)
iux = 0.0
allocate(jvy(ny1),stat=istatus)
jvy = 0.0
allocate(wgtsx(nx1),stat=istatus)
wgtsx = 0.0
allocate(wgtsy(ny1),stat=istatus)
wgtsy = 0.0
allocate(wgtux(nx1),stat=istatus)
wgtux = 0.0
allocate(wgtvy(ny1),stat=istatus)
wgtvy = 0.0
allocate(temx1yz(nx1,ny),stat=istatus)
temx1yz =0.0
DO i=1,nx1-1
isx(i) = MAX(1, MIN(nx-2, INT((xs1(i)-xs(1))/dx)+1 ))
wgtsx(i)= (xs(isx(i)+1)-xs1(i))/(xs(isx(i)+1)-xs(isx(i)))
END DO
DO j=1,ny1-1
jsy(j) = MAX(1, MIN(ny-2, INT((ys1(j)-ys(1))/dy)+1 ))
wgtsy(j)= (ys(jsy(j)+1)-ys1(j))/(ys(jsy(j)+1)-ys(jsy(j)))
END DO
print*,'isx(1),jsy(1)=', isx(1),jsy(1)
DO i=1,nx1
iux(i) = MAX(1, MIN(nx-1, INT((x1 (i)-x (1))/dx)+1 ))
wgtux(i)= (x (iux(i)+1)-x1 (i))/(x (iux(i)+1)-x (iux(i)))
END DO
DO j=1,ny1
jvy(j) = MAX(1, MIN(ny-1, INT((y1 (j)-y (1))/dy)+1 ))
wgtvy(j)= (y (jvy(j)+1)-y1 (j))/(y (jvy(j)+1)-y (jvy(j)))
END DO
endif
print*,'nx1, ny1 =', nx1, ny1
print*,'x1(1), x1(nx1), y1(1), y1(ny1) for output grid are :'
print*,x1(1), x1(nx1), y1(1), y1(ny1)
DO k=1,nz1
z1(k)=z(nzbgn+k-1)
END DO
DO nfile = 1,nhisfile
WRITE(6,'(/a,a,a)') ' Data set ', trim(hisfile(nfile)),' to be subsampled.'
CALL dtareaddump(nx,ny,nz,nzsoil,x,y,z,nstyps, &
trim(grdbasfn),trim(hisfile(nfile)), &
nx1,ny1,nz1,x1,y1,z1,soiltyp,stypfrct,xctr1,yctr1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend, &
tem1,tem2,tem11,tem21,intrp_opt,savespace, &
isx,jsy,iux,jvy,wgtsx,wgtsy,wgtux,wgtvy,temx1yz)
ENDDO
STOP
100 WRITE(6,'(a)') 'Error reading NAMELIST file. Program ARPSINTRP_LS stopped.'
STOP
END PROGRAM ARPSSUBDOMAIN
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE DTAREADDUMP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE dtareaddump(nx,ny,nz,nzsoil,x,y,z,nstyps, & 1,8
grdbasfn,datafn, &
nx1,ny1,nz1,x1,y1,z1,soiltyp,stypfrct,xctr1,yctr1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend, &
tem1,tem2,tem11,tem21,intrp_opt,savespace, &
isx,jsy,iux,jvy,wgtsx,wgtsy,wgtux,wgtvy,temx1yz)
!
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
INTEGER :: nzsoil ! Number of grid points in the soil
INTEGER :: nstyps ! Number of soil type
INTEGER :: hinfmt ! The format of the history data dump
CHARACTER(LEN=*) :: grdbasfn ! Name of the grid/base state array file
CHARACTER(LEN=*) :: datafn ! Name of the other time dependent data file
REAL :: x(nx),y(ny),z(nz)
INTEGER :: nx1,ny1,nz1
REAL :: x1(nx1),y1(ny1),z1(nz1)
REAL :: stypfrct(nx,ny,nstyps)
INTEGER :: soiltyp(nx,ny,nstyps)
INTEGER :: nxbgn,nxend,nybgn,nyend,nzbgn,nzend
REAL :: tem1(nx,ny,nz) ! Temporary array
INTEGER(2) :: tem2(nx,ny,nz) ! Temporary array
REAL :: tem11(nx1,ny1,nz1) ! Temporary array
INTEGER(2) :: tem21(nx1,ny1,nz1) ! Temporary array
integer :: savespace,INTRP_OPT
integer :: isx(nx1),jsy(ny1),iux(nx1),jvy(ny1)
real :: wgtsx(nx1),wgtsy(ny1),wgtux(nx1),wgtvy(ny1)
real :: temx1yz(nx1,ny)
INTEGER :: grdread,iread
SAVE grdread
INTEGER :: istat
INTEGER :: ireturn ! Return status indicator
INTEGER :: grdbas ! Wether this is a grid/base state
! array dump
INTEGER :: i,j,k
LOGICAL :: fexist
INTEGER :: is
REAL :: xctr1,yctr1
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'indtflg.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid & map parameters.
INCLUDE 'mp.inc' ! mpi parameters.
INCLUDE 'exbc.inc'
INCLUDE 'phycst.inc'
DATA grdread /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ireturn = 0
hinfmt = 3
!
!-----------------------------------------------------------------------
!
! Open and read grid and base state data file depending on the
! values of parameters grdin and basin, which are read in from the
! time dependent data set. If grdin or basin is zero, the grid and
! base state arrays have to be read in from a separate file.
!
!-----------------------------------------------------------------------
!
IF( grdread == 0 ) THEN
! print*,'grdread inside if block=', grdread
grdbas = 1
INQUIRE(FILE=grdbasfn, EXIST = fexist )
IF( fexist ) GO TO 200
INQUIRE(FILE=grdbasfn//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( grdbasfn//'.Z' )
GO TO 200
END IF
INQUIRE(FILE=grdbasfn//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( grdbasfn//'.gz' )
GO TO 200
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//grdbasfn// &
' or its compressed version not found.', &
'Program stopped in DTAREAD.'
CALL arpsstop('arpsstop called from dtareaddump during base state read',1)
200 CONTINUE
!
!-----------------------------------------------------------------------
! Read grid and base state fields.
!-----------------------------------------------------------------------
!
CALL hdfreaddump(nx,ny,nz,nzsoil,x,y,z,nstyps,grdbas,trim(grdbasfn), &
nx1,ny1,nz1,x1,y1,z1,soiltyp,stypfrct,xctr1,yctr1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend, &
tem1,tem2,tem11,tem21,intrp_opt,savespace, &
isx,jsy,iux,jvy,wgtsx,wgtsy,wgtux,wgtvy,temx1yz,&
0, 1,1,1,1,1,1, qcexout,qrexout,qiexout,qsexout,qhexout,ngbrz,rayklow)
grdread = 1
END IF
!
!
!-----------------------------------------------------------------------
! Read time dependent data fields.
!-----------------------------------------------------------------------
!
!
grdbas = 0
INQUIRE(FILE=trim(datafn), EXIST = fexist )
IF( fexist ) GO TO 100
INQUIRE(FILE=trim(datafn)//'.Z', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(datafn)//'.Z' )
GO TO 100
END IF
INQUIRE(FILE=trim(datafn)//'.gz', EXIST = fexist )
IF( fexist ) THEN
CALL uncmprs( trim(datafn)//'.gz' )
GO TO 100
END IF
WRITE(6,'(/1x,a,/1x,a/)') &
'File '//trim(datafn) &
//' or its compressed version not found.', &
'Program stopped in DTAREADDUMP.'
CALL arpsstop('arpsstop called from dtareaddump during base read-2',1)
100 CONTINUE
CALL hdfreaddump(nx,ny,nz,nzsoil,x,y,z,nstyps,grdbas, trim(datafn), &
nx1,ny1,nz1,x1,y1,z1,soiltyp,stypfrct,xctr1,yctr1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend, &
tem1,tem2,tem11,tem21,intrp_opt,savespace, &
isx,jsy,iux,jvy,wgtsx,wgtsy,wgtux,wgtvy,temx1yz,&
exbcdmp, 1,1,1,1,1,1, qcexout,qrexout,qiexout,qsexout,qhexout,ngbrz,rayklow)
RETURN
END SUBROUTINE dtareaddump
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HDFREAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE hdfreaddump(nx,ny,nz,nzsoil,x,y,z,nstyps,grdbas,filename, & 2,333
nx1,ny1,nz1,x1,y1,z1,soiltyp,stypfrct,xctr1,yctr1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend, &
tem1,itmp,tem1o,itmpo,intrp_opt,savespace, &
isx,jsy,iux,jvy,wgtsx,wgtsy,wgtux,wgtvy,temx1yz,&
exbcdmpopt, ubcdmp,vbcdmp,wbcdmp,ptbcdmp,prbcdmp, &
qvbcdmp,qcbcdmp,qrbcdmp,qibcdmp,qsbcdmp,qhbcdmp,ngbrz,rayklowest)
!-----------------------------------------------------------------------
! PURPOSE:
! Read in history data in the NCSA HDF4 format.
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 2000/04/15
!
! MODIFICATION HISTORY:
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nzsoil Number of grid points in the soil
!
! grdbas Data read flag.
! =1, only grid and base state arrays will be read
! =0, all arrays will be read based on data
! parameter setting.
! filename Character variable nhming the input HDF file
!-----------------------------------------------------------------------
! Variable Declarations.
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: nzsoil
INTEGER :: grdbas
CHARACTER (LEN=*) :: filename
CHARACTER (LEN=256) :: filename_out
INTEGER :: lfilename_out
INTEGER :: itema, lenstr
INTEGER :: iyr, imon, idy, ihr, imin, isec
CHARACTER (LEN=256) :: exbcfn
CHARACTER (LEN=256) :: temchar
CHARACTER (LEN=15) :: ctime
INTEGER :: exbcdmpopt
INTEGER :: ubcdmp,vbcdmp,wbcdmp,ptbcdmp,prbcdmp, &
qvbcdmp,qcbcdmp,qrbcdmp,qibcdmp,qsbcdmp,qhbcdmp
INTEGER :: ngbrz, rayklowest
INTEGER :: exbccompr
REAL :: x (nx) ! x coord.
REAL :: y (ny) ! y coord.
REAL :: z (nz) ! z coord.
! REAL :: zpsoil(nx,ny,nzsoil) ! physical x coord. for soil (m)
REAL, ALLOCATABLE :: zpsoil(:,:,:) ! physical x coord. for soil (m)
INTEGER :: nstyps
INTEGER :: nx1,ny1,nz1
REAL :: x1(nx1),y1(ny1),z1(nz1)
REAL :: x1_out(nx1),y1_out(ny1)
! REAL :: zpsoil1(nx1,ny1,nzsoil)
REAL, ALLOCATABLE :: zpsoil1(:,:,:) ! physical x coord. for soil (m)
REAL :: stypfrct(nx,ny,nstyps)
INTEGER :: soiltyp(nx,ny,nstyps)
INTEGER :: nxbgn,nxend,nybgn,nyend,nzbgn,nzend
integer :: intrp_opt
integer :: isx(nx1),jsy(ny1),iux(nx1),jvy(ny1)
real :: wgtsx(nx1),wgtsy(ny1),wgtux(nx1),wgtvy(ny1)
real :: temx1yz(nx1,ny)
REAL :: tem1(nx,ny,nz) ! Temporary array
INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz) ! Temporary array
REAL, ALLOCATABLE :: hmax(:), hmin(:) ! Temporary array
REAL, ALLOCATABLE :: hmaxsoil(:), hminsoil(:) ! Temporary array
INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmpsoil(:,:,:)
INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmpsoil4(:,:,:,:)
REAL :: tem1o(nx1,ny1,nz1) ! Temporary work array
INTEGER(2) :: itmpo(nx1,ny1,nz1) ! Temporary array
REAL, allocatable :: hmaxo(:), hmino(:) ! Temporary array
INTEGER(2), allocatable :: itmpsoilo(:,:,:) ! Temporary array
INTEGER(2), allocatable :: itmpsoilo4(:,:,:,:) ! Temporary array
REAL, allocatable :: hmaxsoilo(:), hminsoilo(:)! Temporary array
INTEGER :: ireturn
REAL :: time, dx1, dy1, temvar
REAL :: ptmin, ptmax
! surface and soil variables
! by mhu
! INTEGER, allocatable :: soiltyp(:,:,:)
! REAL, allocatable :: stypfrct(:,:,:)
INTEGER, allocatable :: vegtyp(:,:)
INTEGER, allocatable :: soiltyp1(:,:,:)
REAL, allocatable :: stypfrct1(:,:,:)
INTEGER, allocatable :: vegtyp1(:,:)
REAL, allocatable :: rsoil(:,:,:,:)
REAL, allocatable :: rsoil1(:,:,:,:)
REAL, allocatable :: wetcsoil(:,:,:)
REAL, allocatable :: wetcsoil1(:,:,:)
!-----------------------------------------------------------------------
! Parameters describing routine that wrote the gridded data
!-----------------------------------------------------------------------
!
! 06/28/2002 Zuwen He
!
! fmtver??: to label each data a version.
! intver??: an integer to allow faster comparison than fmtver??,
! which are strings.
!
! Verion 5.00: significant change in soil variables since version 4.10.
!
!-----------------------------------------------------------------------
CHARACTER (LEN=40) :: fmtver410,fmtver500
INTEGER :: intver,intver410,intver500
PARAMETER (fmtver410='004.10 HDF4 Coded Data',intver410=410)
PARAMETER (fmtver500='005.00 HDF4 Coded Data',intver500=500)
CHARACTER (LEN=40) :: fmtverin
CHARACTER (LEN=10) :: tmunit
!-----------------------------------------------------------------------
! Map
!-----------------------------------------------------------------------
REAL :: alatpro(2)
REAL :: ctrlat1,ctrlon1
REAL :: ctrx, ctry, dxscl, dyscl
REAL :: xctr1,yctr1
REAL :: swx,swy
!-----------------------------------------------------------------------
! Terrain
!-----------------------------------------------------------------------
REAL, allocatable :: hterain(:,:)
CHARACTER(len=256) :: ternfn
INTEGER :: nunit,lternfn
!-----------------------------------------------------------------------
! Misc. local variables
!-----------------------------------------------------------------------
INTEGER :: lchanl
PARAMETER (lchanl=6) ! Channel number for formatted printing.
INTEGER :: i,j,k,is
INTEGER :: nxin,nyin,nzin,nzsoilin
INTEGER :: bgrdin,bbasin,bvarin,bicein,btrbin,btkein
INTEGER :: istat, sd_id,sd_id1, sd_id2
INTEGER :: nstyp1,nstypin
CHARACTER (LEN=30) :: substring
INTEGER :: savespace
INTEGER :: varflg
CHARACTER (LEN=40) :: runname_old
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc' ! Grid parameters
INCLUDE 'indtflg.inc'
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF (exbcdmpopt /= 0) THEN
IF (exbcdmpopt == 3) THEN
exbccompr = hdfcompr
ELSE IF (exbcdmpopt == 4) THEN
exbccompr = 5
ELSE
exbccompr = 0
END IF
ENDIF
! if( nxbgn>999.or.nxend>999.or.nybgn>999.or.nyend>999 &
! .or.nzbgn>999.or.nzend>999) then
! write(substring,'(6(a1,i4.4))') &
! '.',nxbgn,'-',nxend,'x',nybgn,'-',nyend,'x',nzbgn,'-',nzend
! else
! write(substring,'(6(a1,i3.3))') &
! '.',nxbgn,'-',nxend,'x',nybgn,'-',nyend,'x',nzbgn,'-',nzend
! endif
if( nxbgn>999.or.nxend>999.or.nybgn>999.or.nyend>999 &
.or.nzbgn>999.or.nzend>999) then
write(substring,'(6(a1,i4.4))') &
'.',nxend,'x',nyend,'x',nzend
else
write(substring,'(6(a1,i3.3))') &
'.',nxend,'x',nyend,'x',nzend
endif
WRITE(*,*) 'HDFREAD: Reading HDF file: ', trim(filename)
ALLOCATE (hmax(nz),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFREAD: ERROR allocating hmax, returning"
RETURN
END IF
ALLOCATE (hmin(nz),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFREAD: ERROR allocating hmin, returning"
RETURN
END IF
!-----------------------------------------------------------------------
! Open file for reading
!-----------------------------------------------------------------------
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "HDFREAD: ERROR opening ", &
trim(filename)," for reading."
GO TO 110
END IF
!-----------------------------------------------------------------------
! Create file for writing
!-----------------------------------------------------------------------
lfilename_out = len_trim( filename )
DO i=lfilename_out,1,-1
if( filename(i:i) .eq.'/') exit
enddo
filename_out = trim(dirname)//filename(i+1:lfilename_out)//trim(substring)
lfilename_out = len_trim( filename_out )
CALL fnversn( filename_out, lfilename_out )
if( hdmpfmt /= 0 ) then
CALL hdfopen(filename_out(1:lfilename_out),2,sd_id1)
IF (sd_id1 < 0) THEN
WRITE (6,*) "HDFDUMP: ERROR creating HDF4 file: ",filename_out(1:lfilename_out)
WRITE (6,*) "Program stopped."
STOP
END IF
endif
fmtverin = fmtver500
intver = intver500
WRITE(6,'(/1x,a,a/)') &
'Incoming data format, fmtverin=',fmtverin
CALL hdfrdc(sd_id,40,"runname",runname_old,istat)
runname=runname_old
CALL hdfrdi(sd_id,"nocmnt",nocmnt,istat)
IF( nocmnt > 0 ) THEN
CALL hdfrdc(sd_id,80*nocmnt,"cmnt",cmnt,istat)
END IF
WRITE(6,'(//'' THE NAME OF THIS RUN IS: '',A//)') trim(runname)
WRITE (6,*) "Comments:"
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
WRITE(6,'(1x,a)') cmnt(i)
END DO
END IF
if( hdmpfmt /= 0 ) then
CALL hdfwrtc(sd_id1, 40, 'fmtver', fmtverin, istat)
CALL hdfwrtc(sd_id1, 40, 'runname', runname, istat)
CALL hdfwrti(sd_id1, 'nocmnt', nocmnt, istat)
IF( nocmnt > 0 ) THEN
CALL hdfwrtc(sd_id1, 80*nocmnt, 'cmnt', cmnt, istat)
END IF
endif
WRITE (6,*) " "
CALL hdfrdc(sd_id,10,"tmunit",tmunit,istat)
CALL hdfrdr(sd_id,"time",time,istat)
if( hdmpfmt /= 0 ) then
WRITE(6,'(1x,a,f13.3,a,a,a/)') &
'Writing HDF4 data at time=', curtim,' into file ',filename_out(1:lfilename_out)
else
WRITE(6,'(1x,a)') 'hdmpfmt =0, no history dump will be created.'
endif
if( hdmpfmt /= 0 ) then
CALL hdfwrtc(sd_id1, 7, 'tmunit', 'seconds', istat)
CALL hdfwrtr(sd_id1, 'time', time, istat)
endif
dx1 = x1(2)-x1(1)
dy1 = y1(2)-y1(1)
DO i=1,nx1
x1_out(i) = x1(i)-x1(2)
END DO
DO j=1,ny1
y1_out(j) = y1(j)-y1(2)
END DO
!-----------------------------------------------------------------------
! Get dimensions of data in binary file and check against
! the dimensions passed to HDFREAD
!-----------------------------------------------------------------------
CALL hdfrdi(sd_id,"nx",nxin,istat)
CALL hdfrdi(sd_id,"ny",nyin,istat)
CALL hdfrdi(sd_id,"nz",nzin,istat)
IF ( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
WRITE(6,'(1x,a)') ' Dimensions in HDFREAD inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
WRITE(6,'(1x,a,3I15)') ' Expected: ', nx, ny, nz
WRITE(6,'(1x,a)') ' Program aborted in HDFREAD.'
CALL arpsstop('arpsstop called from HDFREAD due to nxin...',1)
END IF
!-----------------------------------------------------------------------
! Read in x,y and z at grid cell centers (scalar points).
!-----------------------------------------------------------------------
IF( grdin == 1 .OR. grdbas == 1 ) THEN
CALL hdfrd1d(sd_id,"x",nx,x,istat)
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"y",ny,y,istat)
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"z",nz,z,istat)
IF (istat /= 0) GO TO 110
! CALL hdfrd3d(sd_id,"zpsoil",nx,ny,nzsoil,zpsoil,istat, &
! itmpsoil,hmaxsoil,hminsoil)
! IF (istat /= 0) GO TO 110
END IF ! grdin
if( hdmpfmt /= 0 ) then
CALL hdfwrti(sd_id1, 'nx', nx1, istat)
CALL hdfwrti(sd_id1, 'ny', ny1, istat)
CALL hdfwrti(sd_id1, 'nz', nz1, istat)
endif
CALL hdfrdi(sd_id,"nzsoil",nzsoilin,istat)
if( hdmpfmt /= 0 ) &
CALL hdfwrti(sd_id1, 'nzsoil', nzsoilin, istat)
IF (nzsoilin /= nzsoil) THEN
WRITE(6,'(1x,a)') &
' Dimensions in HDFREADDUMP inconsistent with data.'
WRITE(6,'(1x,a,I15)') ' Read were: ', nzsoilin
WRITE(6,'(1x,a,I15)') ' Expected: ', nzsoil
WRITE(6,'(1x,a)') ' Program aborted in HDFREADDUMP.'
CALL arpsstop('arpsstop called from HDFREADDUMP due to nzsoilin...',1)
END IF
IF (hdfcompr > 3 .or. exbccompr > 3 ) THEN
ALLOCATE (hmaxo(nz1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmaxo, returning"
RETURN
END IF
ALLOCATE (hmino(nz1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmino, returning"
RETURN
END IF
END IF
!-----------------------------------------------------------------------
!
! Read in flags for different data groups
!
!-----------------------------------------------------------------------
IF ( grdbas == 1 ) THEN ! Read grid and base state arrays
WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)') &
'To read grid and base state data at time ', time, &
' secs = ',(time/60.),' mins.'
CALL hdfrdi(sd_id,"grdflg",bgrdin,istat)
CALL hdfrdi(sd_id,"basflg",bbasin,istat)
CALL hdfrdi(sd_id,"varflg",bvarin,istat)
CALL hdfrdi(sd_id,"mstflg",mstin,istat)
CALL hdfrdi(sd_id,"iceflg",bicein,istat)
CALL hdfrdi(sd_id,"trbflg",btrbin,istat)
CALL hdfrdi(sd_id,"landflg",landin,istat)
CALL hdfrdi(sd_id,"totflg",totin,istat)
CALL hdfrdi(sd_id,"tkeflg",btkein,istat)
if( hdmpfmt /= 0 ) then
! mhu landout = 0 ! writing code not implemented in this version
CALL hdfwrti(sd_id1, 'grdflg', 1, istat)
CALL hdfwrti(sd_id1, 'basflg', 1, istat)
CALL hdfwrti(sd_id1, 'varflg', 0, istat)
CALL hdfwrti(sd_id1, 'mstflg', 1, istat)
CALL hdfwrti(sd_id1, 'iceflg', 0, istat)
CALL hdfwrti(sd_id1, 'trbflg', 0, istat)
CALL hdfwrti(sd_id1, 'sfcflg', 0, istat)
CALL hdfwrti(sd_id1, 'rainflg', 0, istat)
CALL hdfwrti(sd_id1, 'landflg', landout, istat)
CALL hdfwrti(sd_id1, 'totflg', 1, istat)
CALL hdfwrti(sd_id1, 'tkeflg', 0, istat)
endif
ELSE ! Normal data reading
WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:', &
time,' secs = ',(time/60.),' mins.'
CALL hdfrdi(sd_id,"grdflg",grdin,istat)
CALL hdfrdi(sd_id,"basflg",basin,istat)
CALL hdfrdi(sd_id,"varflg",varin,istat)
CALL hdfrdi(sd_id,"mstflg",mstin,istat)
CALL hdfrdi(sd_id,"iceflg",icein,istat)
CALL hdfrdi(sd_id,"trbflg",trbin,istat)
CALL hdfrdi(sd_id,"sfcflg",sfcin,istat)
CALL hdfrdi(sd_id,"rainflg",rainin,istat)
CALL hdfrdi(sd_id,"landflg",landin,istat)
CALL hdfrdi(sd_id,"totflg",totin,istat)
CALL hdfrdi(sd_id,"tkeflg",tkein,istat)
if( hdmpfmt /= 0 ) then
landout=0 ! no land data for normal dump file
CALL hdfwrti(sd_id1, 'grdflg', grdin*grdout, istat)
CALL hdfwrti(sd_id1, 'basflg', basin*basout, istat)
if( varin == 0 ) then
varflg = 0
elseif( varin == 1 .and. varout == 1 ) then
varflg = 1
elseif( varin == 1 .and. (varout == 2 .or. varout == 3) ) then
varflg = varout
elseif( varin == 2 .and. varout == 2 ) then
varflg = 2
elseif( varin == 3 .and. varout == 3 ) then
varflg = 3
else
varflg = 0
endif
CALL hdfwrti(sd_id1, 'varflg', varflg , istat)
CALL hdfwrti(sd_id1, 'mstflg', mstin*mstout, istat)
CALL hdfwrti(sd_id1, 'iceflg', icein*iceout, istat)
CALL hdfwrti(sd_id1, 'trbflg', trbin*trbout, istat)
CALL hdfwrti(sd_id1, 'sfcflg', sfcin*sfcout, istat)
CALL hdfwrti(sd_id1, 'rainflg',rainin*rainout, istat)
CALL hdfwrti(sd_id1, 'landflg',landin*landout, istat)
CALL hdfwrti(sd_id1, 'totflg', 1, istat)
CALL hdfwrti(sd_id1, 'tkeflg', tkein*tkeout, istat)
endif
END IF
CALL hdfrdi(sd_id,"nstyp",nstyp1,istat)
IF ( nstyp1 < 1 ) THEN
nstyp1 = 1
END IF
if( hdmpfmt /= 0 ) &
CALL hdfwrti(sd_id1, 'nstyp', nstyp1, istat)
IF (nstyp1 > nstyp) THEN
WRITE (6,*) "HDFREAD: WARNING, nstyp in file (",nstyp1, &
") greater than that specified in input file (",nstyp, &
"), using only",nstyp
nstypin = nstyp
ELSE
nstypin = nstyp1
ENDIF
CALL hdfrdi(sd_id,"prcflg",prcin,istat)
CALL hdfrdi(sd_id,"radflg",radin,istat)
CALL hdfrdi(sd_id,"flxflg",flxin,istat)
CALL hdfrdi(sd_id,"snowflg",snowin,istat)
CALL hdfrdi(sd_id,"month",month,istat)
CALL hdfrdi(sd_id,"day",day,istat)
CALL hdfrdi(sd_id,"year",year,istat)
CALL hdfrdi(sd_id,"hour",hour,istat)
CALL hdfrdi(sd_id,"minute",minute,istat)
CALL hdfrdi(sd_id,"second",second,istat)
CALL hdfrdr(sd_id,"umove",umove,istat)
CALL hdfrdr(sd_id,"vmove",vmove,istat)
CALL hdfrdr(sd_id,"xgrdorg",xgrdorg,istat)
CALL hdfrdr(sd_id,"ygrdorg",ygrdorg,istat)
CALL hdfrdi(sd_id,"mapproj",mapproj,istat)
CALL hdfrdr(sd_id,"trulat1",trulat1,istat)
CALL hdfrdr(sd_id,"trulat2",trulat2,istat)
CALL hdfrdr(sd_id,"trulon",trulon,istat)
CALL hdfrdr(sd_id,"sclfct",sclfct,istat)
CALL hdfrdr(sd_id,"tstop",tstop,istat)
CALL hdfrdr(sd_id,"thisdmp",thisdmp,istat)
CALL hdfrdr(sd_id,"latitud",latitud,istat)
CALL hdfrdr(sd_id,"ctrlat",ctrlat,istat)
CALL hdfrdr(sd_id,"ctrlon",ctrlon,istat)
alatpro(1)=trulat1
alatpro(2)=trulat2
CALL setmapr( mapproj,sclfct,alatpro,trulon )
CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
dxscl=x(2)-x(1)
dyscl=y(2)-y(1)
swx = ctrx - (FLOAT(nx-3)/2.)*dxscl
swy = ctry - (FLOAT(ny-3)/2.)*dyscl
CALL setorig( 1, swx, swy)
CALL xytoll(1,1,xctr1,yctr1, ctrlat1,ctrlon1)
write(*,*) 'The new domain has a center point at:'
write(*,*) 'lat=',ctrlat1,'lon=',ctrlon1
if( hdmpfmt /= 0 ) then
CALL hdfwrti(sd_id1, 'prcflg', prcin*prcout, istat)
CALL hdfwrti(sd_id1, 'radflg', radin*radout, istat)
CALL hdfwrti(sd_id1, 'flxflg', flxout*flxout, istat)
CALL hdfwrti(sd_id1, 'snowflg', snowin*snowout, istat)
CALL hdfwrti(sd_id1, 'day', day, istat)
CALL hdfwrti(sd_id1, 'year', year, istat)
CALL hdfwrti(sd_id1, 'month', month, istat)
CALL hdfwrti(sd_id1, 'hour', hour, istat)
CALL hdfwrti(sd_id1, 'minute', minute, istat)
CALL hdfwrti(sd_id1, 'second', second, istat)
CALL hdfwrtr(sd_id1, 'umove', umove, istat)
CALL hdfwrtr(sd_id1, 'vmove', vmove, istat)
CALL hdfwrtr(sd_id1, 'xgrdorg', xgrdorg, istat)
CALL hdfwrtr(sd_id1, 'ygrdorg', ygrdorg, istat)
CALL hdfwrti(sd_id1, 'mapproj', mapproj, istat)
CALL hdfwrtr(sd_id1, 'trulat1', trulat1, istat)
CALL hdfwrtr(sd_id1, 'trulat2', trulat2, istat)
CALL hdfwrtr(sd_id1, 'trulon', trulon, istat)
CALL hdfwrtr(sd_id1, 'sclfct', sclfct, istat)
CALL hdfwrtr(sd_id1, 'tstop', tstop, istat)
CALL hdfwrtr(sd_id1, 'thisdmp', thisdmp, istat)
CALL hdfwrtr(sd_id1, 'latitud', latitud, istat)
CALL hdfwrtr(sd_id1, 'ctrlat', ctrlat1, istat)
CALL hdfwrtr(sd_id1, 'ctrlon', ctrlon1, istat)
CALL hdfwrtr(sd_id1, 'dx', dx1, istat)
CALL hdfwrtr(sd_id1, 'dy', dy1, istat)
endif
!-----------------------------------------------------------------------
! Write exbc file
!-----------------------------------------------------------------------
IF (exbcdmpopt /= 0) THEN
CALL ctim2abss( year,month,day,hour,minute,second, itema )
itema = itema + INT(time)
CALL abss2ctim( itema,iyr,imon,idy,ihr,imin,isec )
WRITE (ctime,'(i4.4,2i2.2,a,3i2.2)')iyr,imon,idy,'.',ihr,imin,isec
lfnkey = len_trim(runname)
CALL gtlfnkey( runname, lfnkey )
exbcfn = runname(1:lfnkey)//'.'//ctime
lenstr = lfnkey + 16
IF( dirname /= ' ' ) THEN
temchar = exbcfn
exbcfn = trim(dirname)//'/'//temchar
lenstr = len_trim(exbcfn)
END IF
CALL fnversn(exbcfn,lenstr)
WRITE(6,'(1x,a,a)') &
'Dumping to the external boundary format file: ',exbcfn(1:lenstr)
CALL hdfopen(exbcfn(1:lenstr), 2, sd_id2)
IF (sd_id2 < 0) THEN
WRITE (6,*) "WRITEXBC: ERROR creating HDF4 file: ",exbcfn(1:lenstr)
STOP
END IF
CALL hdfwrtc(sd_id2, 40, 'fmtver', fmtver500, istat)
CALL hdfwrtc(sd_id2, 15, 'ctime', ctime, istat)
CALL hdfwrti(sd_id2, 'nx', nx1, istat)
CALL hdfwrti(sd_id2, 'ny', ny1, istat)
CALL hdfwrti(sd_id2, 'nz', nz1, istat)
CALL hdfwrtr(sd_id2, 'dx', dx1, istat)
CALL hdfwrtr(sd_id2, 'dy', dy1, istat)
CALL hdfwrtr(sd_id2, 'dz', dz, istat)
CALL hdfwrtr(sd_id2, 'dzmin', dzmin, istat) ! no known yet
CALL hdfwrti(sd_id2, 'strhopt', strhopt, istat)
CALL hdfwrtr(sd_id2, 'zrefsfc', zrefsfc, istat)
CALL hdfwrtr(sd_id2, 'dlayer1', dlayer1, istat)
CALL hdfwrtr(sd_id2, 'dlayer2', dlayer2, istat)
CALL hdfwrtr(sd_id2, 'zflat', zflat, istat)
CALL hdfwrtr(sd_id2, 'strhtune', strhtune, istat)
CALL hdfwrti(sd_id2, 'mapproj', mapproj, istat)
CALL hdfwrtr(sd_id2, 'trulat1', trulat1, istat)
CALL hdfwrtr(sd_id2, 'trulat2', trulat2, istat)
CALL hdfwrtr(sd_id2, 'trulon', trulon, istat)
CALL hdfwrtr(sd_id2, 'sclfct', sclfct, istat)
CALL hdfwrtr(sd_id2, 'ctrlat', ctrlat1, istat)
CALL hdfwrtr(sd_id2, 'ctrlon', ctrlon1, istat)
CALL hdfwrti(sd_id2,"ubcflg",ubcdmp,istat)
CALL hdfwrti(sd_id2,"vbcflg",vbcdmp,istat)
CALL hdfwrti(sd_id2,"wbcflg",wbcdmp,istat)
CALL hdfwrti(sd_id2,"ptbcflg",ptbcdmp,istat)
CALL hdfwrti(sd_id2,"prbcflg",prbcdmp,istat)
CALL hdfwrti(sd_id2,"qvbcflg",qvbcdmp,istat)
CALL hdfwrti(sd_id2,"qcbcflg",qcbcdmp,istat)
CALL hdfwrti(sd_id2,"qrbcflg",qrbcdmp,istat)
CALL hdfwrti(sd_id2,"qibcflg",qibcdmp,istat)
CALL hdfwrti(sd_id2,"qsbcflg",qsbcdmp,istat)
CALL hdfwrti(sd_id2,"qhbcflg",qhbcdmp,istat)
IF (exbcdmpopt == 4) THEN
CALL hdfwrti(sd_id2, 'clipxy', ngbrz, istat)
CALL hdfwrti(sd_id2, 'clipz', rayklowest, istat)
END IF
ENDIF
!-----------------------------------------------------------------------
! If grdout=1 or grdbas=1, write out grid variables
!-----------------------------------------------------------------------
IF((grdout == 1 .OR. grdbas == 1).and.(hdmpfmt /= 0) ) THEN
CALL hdfwrt1d(x1_out,nx1,sd_id1,'x','x coordinate','m')
CALL hdfwrt1d(y1_out,ny1,sd_id1,'y','y coordinate','m')
CALL hdfwrt1d(z1,nz1,sd_id1,'z','z coordinate','m')
IF( intver >= intver500 .and. grdbas == 1 ) then
ALLOCATE (zpsoil(nx,ny,nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating zpsoil, returning"
RETURN
END IF
ALLOCATE (zpsoil1(nx1,ny1,nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating zpsoil1, returning"
RETURN
END IF
ALLOCATE (itmpsoil(nx,ny,nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating itmpsoil, returning"
RETURN
END IF
ALLOCATE (hmaxsoil(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmaxsoil, returning"
RETURN
END IF
ALLOCATE (hminsoil(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminisoil, returning"
RETURN
END IF
CALL hdfrd3d(sd_id,"zpsoil",nx,ny,nzsoil,zpsoil,istat, &
itmpsoil,hmaxsoil,hminsoil)
DEALLOCATE(itmpsoil)
DEALLOCATE(hmaxsoil)
DEALLOCATE(hminsoil)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(zpsoil,nx,ny,nzsoil, zpsoil1, nx1,ny1,nzsoil, &
nxbgn,nxend,nybgn,nyend,1,nzsoil)
else
call intrpxy3d(zpsoil,nx,1,nx-1,ny,1,ny-1,nzsoil,1,nzsoil, &
wgtsx,isx,wgtsy,jsy, intrp_opt, &
zpsoil1,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
DEALLOCATE(zpsoil)
ALLOCATE (itmpsoilo(nx1,ny1,nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating itmpsoilo, returning"
RETURN
END IF
ALLOCATE (hmaxsoilo(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmaxsoilo, returning"
RETURN
END IF
ALLOCATE (hminsoilo(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminisoilo, returning"
RETURN
END IF
CALL hdfwrt3d(zpsoil1,nx1,ny1,nzsoil,sd_id1,0,hdfcompr, &
'zpsoil','Physical height coordinate (soil)','m', &
itmpsoilo,hmaxsoilo,hminsoilo)
DEALLOCATE(zpsoil1)
DEALLOCATE(itmpsoilo)
DEALLOCATE(hmaxsoilo)
DEALLOCATE(hminsoilo)
ENDIF
if(grdbas == 1) then
CALL hdfrd3d(sd_id,"zp",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt, &
tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,0,hdfcompr, &
'zp','Physical height coordinate','m', &
itmpo,hmaxo,hmino)
endif
if(grdbas == 1 .and. terndmp >= 1 ) then
!
!-----------------------------------------------------------------------
!
! Write out terrain data
!
!-----------------------------------------------------------------------
!
ALLOCATE (hterain(nx1,ny1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hterain, returning"
RETURN
END IF
hterain(:,:)=tem1o(:,:,2)
CALL getunit( nunit )
lfilename_out = len_trim( runname )
DO i=1,lfilename_out
if( runname(i:i) .eq.',' .or. runname(i:i) .eq.' ' ) exit
enddo
filename_out = trim(dirname)//runname(1:i-1)
ternfn = trim(filename_out)//".trndata"
lternfn = len_trim( ternfn )
CALL fnversn(ternfn, lternfn )
PRINT *, 'Write terrain data to ',trim(ternfn)
CALL writtrn(nx1,ny1,ternfn(1:lternfn), dx1,dy1, &
mapproj,trulat1,trulat2,trulon,sclfct, &
ctrlat1,ctrlon1,hterain)
endif ! terrain
ENDIF
!-----------------------------------------------------------------------
! Read in base state fields
!-----------------------------------------------------------------------
! Print*,'start doing 3d arrays'
IF( (basout== 1 .OR. grdbas == 1) .and. hdmpfmt /= 0) THEN
CALL hdfrd3d(sd_id,"ubar",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx,ny,1,ny-1,nz,nzbgn,nzend, &
wgtux,iux,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'ubar','Base state u-velocity','m/s', &
itmpo,hmaxo,hmino)
CALL hdfrd3d(sd_id,"vbar",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny,nz,nzbgn,nzend, &
wgtsx,isx,wgtvy,jvy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'vbar','Base state v-velocity','m/s', &
itmpo,hmaxo,hmino)
CALL hdfrd3d(sd_id,"wbar",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
tem1o = 0.0
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'wbar','Base state w-velocity','m/s', &
itmpo,hmaxo,hmino)
CALL hdfrd3d(sd_id,"ptbar",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'ptbar','Base state potential temperature','K', &
itmpo,hmaxo,hmino)
print*,'ptbar(1,1,1)=', tem1o(1,1,1)
CALL hdfrd3d(sd_id,"pbar",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'pbar','Base state pressure','Pascal', &
itmpo,hmaxo,hmino)
IF( mstin == 1 ) THEN
CALL hdfrd3d(sd_id,"qvbar",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
ELSE
tem1 = 0.0
END IF
! IF(mstout == 1 ) THEN
! always output qvbar
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'qvbar','Base state water vapor specific humidity','kg/kg', &
itmpo,hmaxo,hmino)
! END IF
IF(landin == 1 .AND. landout == 1 .and. hdmpfmt /= 0 ) THEN
! allocate(soiltyp(nx,ny,nstyp),stat=istat)
! IF (istat /= 0) THEN
! WRITE (6,*) "HDFDUMP: ERROR allocating soiltyp, returning"
! RETURN
! END IF
allocate(soiltyp1(nx1,ny1,nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating soiltyp1, returning"
RETURN
END IF
!! allocate(stypfrct(nx,ny,nstyp),stat=istat)
! IF (istat /= 0) THEN
! WRITE (6,*) "HDFDUMP: ERROR allocating stypfrct, returning"
! RETURN
! END IF
allocate(stypfrct1(nx1,ny1,nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating stypfrct1, returning"
RETURN
END IF
allocate(vegtyp(nx,ny),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating vegtyp, returning"
RETURN
END IF
allocate(vegtyp1(nx1,ny1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating vegtyp1, returning"
RETURN
END IF
allocate(itmpsoil(nx,ny,nstyp),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating itmpsoil, returning"
RETURN
END IF
allocate(hmaxsoil(nstyp),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmaxsoil, returning"
RETURN
END IF
allocate(hminsoil(nstyp),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoil, returning"
RETURN
END IF
CALL hdfrd3di(sd_id,"soiltyp",nx,ny,nstyp,soiltyp,istat)
CALL hdfrd3d(sd_id,"stypfrct",nx,ny,nstyp,stypfrct,istat, &
itmpsoil,hmaxsoil,hminsoil)
CALL hdfrd2di(sd_id,"vegtyp",nx,ny,vegtyp,istat)
deallocate(itmpsoil)
deallocate(hmaxsoil)
deallocate(hminsoil)
if( intrp_opt == 0 ) then
call copyarrayi(soiltyp,nx,ny,nstyp,soiltyp1,nx1,ny1,nstyp1, &
nxbgn,nxend,nybgn,nyend,1,nstyp1)
call copyarray(stypfrct,nx,ny,nstyp,stypfrct1,nx1,ny1,nstyp1, &
nxbgn,nxend,nybgn,nyend,1,nstyp1)
call copyarrayi(vegtyp,nx,ny,1,vegtyp1,nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
CALL intrp_soil_int(nx,ny,nx1,ny1,nstyp,nstyp1,wgtsx,wgtsy,isx,jsy, &
soiltyp,stypfrct,vegtyp, &
soiltyp1,stypfrct1,vegtyp1)
endif
allocate(itmpsoilo(nx1,ny1,nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
allocate(hmaxsoilo(nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
allocate(hminsoilo(nstyp1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
CALL hdfwrt3di(soiltyp1,nx1,ny1,nstyp1,sd_id1,0,0, &
'soiltyp','Soil type','index')
CALL hdfwrt3d(stypfrct1,nx1,ny1,nstyp1,sd_id1,0,hdfcompr, &
'stypfrct','Soil type fractional coverage','fraction', &
itmpsoilo,hmaxsoilo,hminsoilo)
CALL hdfwrt2di(vegtyp1,nx1,ny1,sd_id1,0,0,'vegtyp', &
'Vegetation type','index')
deallocate(itmpsoilo,stat=istat)
deallocate(hmaxsoilo,stat=istat)
deallocate(hminsoilo,stat=istat)
deallocate(soiltyp1,stat=istat)
deallocate(stypfrct1,stat=istat)
deallocate(vegtyp,stat=istat)
deallocate(vegtyp1,stat=istat)
! lai
CALL hdfrd2d(sd_id,"lai",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'lai','Leaf Area Index','index',itmpo)
! roufns
CALL hdfrd2d(sd_id,"roufns",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'roufns','Surface roughness','0-1',itmpo)
! veg
CALL hdfrd2d(sd_id,"veg",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'veg','Vegetation fraction','fraction',itmpo)
END IF ! landin ==1
END IF ! grdbas == 1
IF( grdbas == 1 ) GO TO 930
IF( varflg == 1 .or. varflg == 2 ) then
!-----------------------------------------------------------------------
! Read in total values of variables from history dump
!-----------------------------------------------------------------------
CALL hdfrd3d(sd_id,"u",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx,ny,1,ny-1,nz,nzbgn,nzend, &
wgtux,iux,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'u','u-velocity','m/s', itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. ubcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1, rayklowest )
tem1o((2+ngbrz):(nx1-2-ngbrz),(1+ngbrz):(ny1-1-ngbrz),k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'u','u-velocity','m/s',itmpo,hmaxo,hmino)
END IF
CALL hdfrd3d(sd_id,"v",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny,nz,nzbgn,nzend, &
wgtsx,isx,wgtvy,jvy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'v','v-velocity','m/s', itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. vbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,2+ngbrz:ny1-2-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'v','v-velocity','m/s',itmpo,hmaxo,hmino)
END IF
CALL hdfrd3d(sd_id,"w",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:, 1)=0.0
tem1o(:,:,nz1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'w','w-velocity','m/s', itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. wbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'w','w-velocity','m/s',itmpo,hmaxo,hmino)
END IF
ENDIF
IF( varflg == 1 .or. varflg == 3 ) then
CALL hdfrd3d(sd_id,"pt",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
DO k=1,nz-1
ptmin = tem1(1,1,k)
ptmax = tem1(1,1,k)
DO i=1,nx-1
DO j=1,ny-1
ptmin=min(ptmin,tem1(i,j,k))
ptmax=max(ptmax,tem1(i,j,k))
ENDDO
ENDDO
print*,'k, ptmin, ptmax=', k, ptmin, ptmax
enddo
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'pt','Potential temperature','K', itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. ptbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'pt','Potential temperature','K', itmpo,hmaxo,hmino)
END IF
CALL hdfrd3d(sd_id,"p",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'p','Pressure','Pascal', itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. prbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'p','Pressure','Pascal', itmpo,hmaxo,hmino)
END IF
END IF
!-----------------------------------------------------------------------
! Read in moisture variables
!-----------------------------------------------------------------------
IF( mstin == 1 .AND. mstout == 1) THEN
CALL hdfrd3d(sd_id,"qv",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'qv','Water vapor specific humidity','kg/kg', itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. qvbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'qv','Water vapor specific humidity','kg/kg', itmpo,hmaxo,hmino)
END IF
CALL hdfrd3d(sd_id,"qc",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'qc','Cloud water mixing ratio','kg/kg',itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. qcbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'qc','Cloud water mixing ratio','kg/kg',itmpo,hmaxo,hmino)
END IF
CALL hdfrd3d(sd_id,"qr",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'qr','Rain water mixing ratio','kg/kg',itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. qrbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'qr','Rain water mixing ratio','kg/kg',itmpo,hmaxo,hmino)
END IF
IF( icein == 1 .AND. iceout == 1) THEN
CALL hdfrd3d(sd_id,"qi",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'qi','Cloud ice mixing ratio','kg/kg',itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. qibcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'qi','Cloud ice mixing ratio','kg/kg',itmpo,hmaxo,hmino)
END IF
CALL hdfrd3d(sd_id,"qs",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'qs','Snow mixing ratio','kg/kg',itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. qsbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'qs','Snow mixing ratio','kg/kg',itmpo,hmaxo,hmino)
END IF
CALL hdfrd3d(sd_id,"qh",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,nzbgn,nzend, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'qh','Hail mixing ratio','kg/kg',itmpo,hmaxo,hmino)
IF( exbcdmpopt /= 0 .and. qhbcdmp == 1) THEN
IF (exbcdmpopt == 4) THEN
do k=1,min(nz1,rayklowest)
tem1o(1+ngbrz:nx1-1-ngbrz,1+ngbrz:ny1-1-ngbrz,k:k)=tem1o(2,2,k)
enddo
ENDIF
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id2,1,exbccompr, &
'qh','Hail mixing ratio','kg/kg',itmpo,hmaxo,hmino)
END IF
END IF
IF( rainin == 1 .AND. rainout == 1 .and. hdmpfmt /= 0) THEN
CALL hdfrd2d(sd_id,"raing",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'raing','Grid supersaturation rain','mm',itmpo)
CALL hdfrd2d(sd_id,"rainc",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'rainc','Cumulus convective rain','mm',itmpo)
END IF
IF( prcin == 1 .AND. prcout ==1 .and. hdmpfmt /= 0) THEN
CALL hdfrd2d(sd_id,"prcrate1",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'prcrate1','Total precip. rate','kg/(m**2*s)',itmpo)
CALL hdfrd2d(sd_id,"prcrate2",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'prcrate2','Grid scale precip. rate','kg/(m**2*s)',itmpo)
CALL hdfrd2d(sd_id,"prcrate3",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'prcrate3','Cumulative precip. rate','kg/(m**2*s)',itmpo)
CALL hdfrd2d(sd_id,"prcrate4",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'prcrate4','Microphysics precip. rate','kg/(m**2*s)',itmpo)
END IF
END IF
IF( tkein == 1 .AND. tkeout == 1 .and. hdmpfmt /= 0) THEN
CALL hdfrd3d(sd_id,"tke",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'tke','Turbulent Kinetic Energy','(m/s)**2', &
itmpo,hmaxo,hmino)
END IF
IF( trbin == 1 .AND. trbout == 1 .and. hdmpfmt /= 0 ) THEN
CALL hdfrd3d(sd_id,"kmh",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'kmh','Hori. turb. mixing coef. for momentum','m**2/s', &
itmpo,hmaxo,hmino)
CALL hdfrd3d(sd_id,"kmv",nx,ny,nz,tem1,istat,itmp,hmax,hmin)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,nz, tem1o, nx1,ny1,nz1, &
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt, tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( savespace == 1 ) then
tem1o(:,:,nz1)=0.0
tem1o(:,:,1)=0.0
tem1o(:,:,nz1-1)=0.0
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt3d(tem1o,nx1,ny1,nz1,sd_id1,1,hdfcompr, &
'kmv','Vert. turb. mixing coef. for momentum','m**2/s', &
itmpo,hmaxo,hmino)
END IF
!
!-----------------------------------------------------------------------
!
! If sfcout = 1, write out the surface variables,
! tsoil, qsoil, and wetcanp.
! working here
!-----------------------------------------------------------------------
IF( sfcin == 1 .AND. sfcout == 1 .and. hdmpfmt /= 0 ) THEN
ALLOCATE (itmpsoil4(nx,ny,nzsoil,nstyp+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating itmpsoil, returning"
RETURN
END IF
ALLOCATE (hmaxsoil(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmaxsoil, returning"
RETURN
END IF
ALLOCATE (hminsoil(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminisoil, returning"
RETURN
END IF
allocate(itmpsoilo4(nx1,ny1,nzsoil,nstyp1+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
allocate(hmaxsoilo(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
allocate(hminsoilo(nzsoil),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
ALLOCATE (rsoil(nx,ny,nzsoil,nstyp+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating rsoil, returning"
RETURN
END IF
ALLOCATE (rsoil1(nx1,ny1,nzsoil,nstyp1+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating rsoil1, returning"
RETURN
END IF
CALL hdfrd4d(sd_id,"tsoil",nx,ny,nzsoil,nstyp1+1,rsoil,istat,&
itmpsoil4,hmaxsoil,hminsoil)
if( intrp_opt == 0 ) then
call copyarray4(rsoil,nx,ny,nstyp,rsoil1,nx1,ny1,nstyp1,nzsoil, &
nxbgn,nxend,nybgn,nyend,1,nstyp1)
else
call intrp_soil_real(nx,ny,nx1,ny1,nstyp,nstyp1,nzsoil,wgtsx,wgtsy,isx,jsy, &
rsoil,soiltyp,stypfrct,rsoil1)
endif
CALL hdfwrt4d(rsoil1,nx1,ny1,nzsoil,nstyp1+1,sd_id1,0, &
hdfcompr,'tsoil','Soil temperature','K', &
itmpsoilo4,hmaxsoilo,hminsoilo)
!
CALL hdfrd4d(sd_id,"qsoil",nx,ny,nzsoil,nstyp+1,rsoil,istat,&
itmpsoil4,hmaxsoil,hminsoil)
if( intrp_opt == 0 ) then
call copyarray4(rsoil,nx,ny,nstyp,rsoil1,nx1,ny1,nstyp1,nzsoil, &
nxbgn,nxend,nybgn,nyend,1,nstyp1)
else
call intrp_soil_real(nx,ny,nx1,ny1,nstyp,nstyp1,nzsoil,wgtsx,wgtsy,isx,jsy, &
rsoil,soiltyp,stypfrct,rsoil1)
endif
CALL hdfwrt4d(rsoil1,nx1,ny1,nzsoil,nstyp1+1,sd_id1,0, &
hdfcompr,'qsoil','Soil moisture','fraction', &
itmpsoilo4,hmaxsoilo,hminsoilo)
DEALLOCATE(itmpsoil4)
DEALLOCATE(hmaxsoil)
DEALLOCATE(hminsoil)
DEALLOCATE(itmpsoilo4)
DEALLOCATE(hmaxsoilo)
DEALLOCATE(hminsoilo)
DEALLOCATE(rsoil)
DEALLOCATE(rsoil1)
ALLOCATE (itmpsoil(nx,ny,nstyp+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating itmpsoil, returning"
RETURN
END IF
ALLOCATE (hmaxsoil(nstyp+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hmaxsoil, returning"
RETURN
END IF
ALLOCATE (hminsoil(nstyp+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminisoil, returning"
RETURN
END IF
allocate(itmpsoilo(nx1,ny1,nstyp1+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
allocate(hmaxsoilo(nstyp1+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
allocate(hminsoilo(nstyp1+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating hminsoilo, returning"
RETURN
END IF
ALLOCATE (wetcsoil(nx,ny,nstyp+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating wetcsoil, returning"
RETURN
END IF
ALLOCATE (wetcsoil1(nx1,ny1,nstyp1+1),stat=istat)
IF (istat /= 0) THEN
WRITE (6,*) "HDFDUMP: ERROR allocating wetcsoil1, returning"
RETURN
END IF
CALL hdfrd3d(sd_id,"wetcanp",nx,ny,nstyp+1,wetcsoil,istat, &
itmpsoil,hmaxsoil,hminsoil)
if( intrp_opt == 0 ) then
call copyarray4(wetcsoil,nx,ny,nstyp,wetcsoil1,nx1,ny1,nstyp1,1, &
nxbgn,nxend,nybgn,nyend,1,nstyp1)
else
call intrp_soil_real(nx,ny,nx1,ny1,nstyp,nstyp1,1,wgtsx,wgtsy,isx,jsy, &
wetcsoil,soiltyp,stypfrct,wetcsoil1)
endif
CALL hdfwrt3d(wetcsoil1,nx1,ny1,nstyp1+1,sd_id1,0,hdfcompr, &
'wetcanp','Canopy water amount','fraction', &
itmpsoilo,hmaxsoilo,hminsoilo)
DEALLOCATE(itmpsoil)
DEALLOCATE(hmaxsoil)
DEALLOCATE(hminsoil)
DEALLOCATE(itmpsoilo)
DEALLOCATE(hmaxsoilo)
DEALLOCATE(hminsoilo)
DEALLOCATE(wetcsoil)
DEALLOCATE(wetcsoil1)
IF (snowout == 1) THEN
CALL hdfrd2d(sd_id,"snowdpth",nx,ny,tem1,istat,itmp)
IF (istat /= 0) GO TO 110
if( intrp_opt == 0 ) then
call copyarray(tem1,nx,ny,1, tem1o, nx1,ny1,1, &
nxbgn,nxend,nybgn,nyend,1,1)
else
call intrpxy3d(tem1,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrp_opt,tem1o,nx1,1,nx1-1,ny1,1,ny1-1, temx1yz)
endif
if( hdmpfmt /= 0 ) &
CALL hdfwrt2d(tem1o,nx1,ny1,sd_id1,0,hdfcompr, &
'snowdpth','Snow depth','m',itmpo)
END IF
END IF
!-----------------------------------------------------------------------
!
! Friendly exit message
!
!-----------------------------------------------------------------------
930 CONTINUE
WRITE(6,'(/a,F8.1,a/)') &
' Data at time=', time/60,' (min) were successfully read.'
ireturn = 0
GO TO 130
!-----------------------------------------------------------------------
!
! Error during read
!
!-----------------------------------------------------------------------
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in HDFREAD'
ireturn=1
GO TO 130
!-----------------------------------------------------------------------
!
! End-of-file during read
!
!-----------------------------------------------------------------------
! 120 CONTINUE
! WRITE(6,'(/a/)') ' End of file reached in HDFREAD'
! ireturn=2
130 CONTINUE
CALL hdfclose(sd_id,istat)
IF (ireturn == 0) THEN
IF (istat == 0) THEN
WRITE(6,'(/a/a)') &
"HDFREADDUMP: Successfully read ", trim(filename)
ELSE
WRITE(6,'(/a,i3,a/,a)') &
"HDFREADDUMP: ERROR (status=", istat, ") closing ", trim(filename)
END IF
END IF
if( hdmpfmt /= 0 ) then
CALL hdfclose(sd_id1,istat)
IF (istat == 0) THEN
WRITE(6,'(/a/,a,a)') &
"HDFREADDUMP: Successfully wrote ", filename_out(1:lfilename_out)
endif
endif
IF (exbcdmpopt /= 0) THEN
CALL hdfclose(sd_id2,istat)
ENDIF
! DEALLOCATE (itmp,stat=istat)
DEALLOCATE (hmax,stat=istat)
DEALLOCATE (hmin,stat=istat)
IF( hdfcompr > 3 .or. exbccompr > 3) then
! DEALLOCATE (itmpo,stat=istat)
DEALLOCATE (hmaxo,stat=istat)
DEALLOCATE (hmino,stat=istat)
ENDIF
RETURN
END SUBROUTINE hdfreaddump
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GET_GRIDINFO_FROM_HDF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
SUBROUTINE get_gridinfo_from_hdf(filename,nx,ny,nz,x,y,z,ireturn) 2,15
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 7/17/2000.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! filename Channel number for binary reading.
!
! OUTPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: stat, sd_id
CHARACTER (LEN=*) :: filename
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: x(nx),y(ny),z(nz)
INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmp(:,:,:) ! Temporary array
REAL, ALLOCATABLE :: hmax(:), hmin(:) ! Temporary array
INTEGER :: ireturn ! Return status indicator
INTEGER istat
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "get_gridinfo_from_hdf: ERROR opening ", &
trim(filename)," for reading."
GO TO 110
ELSE
WRITE(6,*) 'File ',filename,' openned.'
END IF
ALLOCATE (itmp(nx,ny,nz),stat=istat)
ALLOCATE (hmax(nz),stat=istat)
ALLOCATE (hmin(nz),stat=istat)
! print*,'sd_id, nx =', sd_id, nx
CALL hdfrd1d(sd_id,"x",nx,x,istat)
! print*,'istat after reading x =', istat
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"y",ny,y,istat)
IF (istat /= 0) GO TO 110
CALL hdfrd1d(sd_id,"z",nz,z,istat)
IF (istat /= 0) GO TO 110
ireturn = 0
GO TO 130
!-----------------------------------------------------------------------
!
! Error during read
!
!-----------------------------------------------------------------------
110 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in GET_GRIDINFO_FROM_HDF.'
ireturn=1
130 CONTINUE
!tmp stat = sfendacc(sd_id) ! is this necessary?
CALL hdfclose(sd_id,stat)
DEALLOCATE (itmp)
DEALLOCATE (hmax)
DEALLOCATE (hmin)
RETURN
END SUBROUTINE get_gridinfo_from_hdf
SUBROUTINE copyarray(arrayin,nx,ny,nz, arrayout, nx1,ny1,nz1, & 32
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
INTEGER :: nx,ny,nz, nx1,ny1,nz1,nxbgn,nxend,nybgn,nyend,nzbgn,nzend
REAL :: arrayin(nx,ny,nz)
REAL :: arrayout(nx1,ny1,nz1)
DO k=1,nz1
DO j=1,ny1
DO i=1,nx1
arrayout(i,j,k)=arrayin(nxbgn+i-1,nybgn+j-1,nzbgn+k-1)
ENDDO
ENDDO
ENDDO
return
END SUBROUTINE copyarray
SUBROUTINE copyarray4(arrayin,nx,ny,nz, arrayout, nx1,ny1,nz1,nzsoil, & 3
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
INTEGER :: nx,ny,nz, nx1,ny1,nz1,nxbgn,nxend,nybgn,nyend,nzbgn,nzend
INTEGER :: nzsoil
REAL :: arrayin(nx,ny,nzsoil,nz)
REAL :: arrayout(nx1,ny1,nzsoil,nz1)
DO k=1,nz1
DO n=1,nzsoil
DO j=1,ny1
DO i=1,nx1
arrayout(i,j,n,k)=arrayin(nxbgn+i-1,nybgn+j-1,n,nzbgn+k-1)
ENDDO
ENDDO
ENDDO
ENDDO
return
END SUBROUTINE copyarray4
SUBROUTINE copyarrayi(arrayin,nx,ny,nz, arrayout, nx1,ny1,nz1, & 2
nxbgn,nxend,nybgn,nyend,nzbgn,nzend)
INTEGER :: nx,ny,nz, nx1,ny1,nz1,nxbgn,nxend,nybgn,nyend,nzbgn,nzend
INTEGER :: arrayin(nx,ny,nz)
INTEGER :: arrayout(nx1,ny1,nz1)
DO k=1,nz1
DO j=1,ny1
DO i=1,nx1
arrayout(i,j,k)=arrayin(nxbgn+i-1,nybgn+j-1,nzbgn+k-1)
ENDDO
ENDDO
ENDDO
return
END SUBROUTINE copyarrayi