!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSINTRP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
PROGRAM arpsintrp,335
!
!-----------------------------------------------------------------------
!
! 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:
!
! 7/31/1997 (M. Xue)
! Added options for specifying the terrain for the new grid, in
! addition to interpolated terrain.
!
! 4/24/1998 (D. Weber)
! Added hydrostatic extrapolation calculations for pprt and pbar
! when the fine grid terrain is lower in elevation than the
! coarse grid terrain. (DO LOOP 1423) Note: this update is
! functional for bglopt=2 ONLY.
!
! 4/24/1998 (D. Weber)
! Added option for quadratic interpolation in the vertical.
! intrpvopt = 1 for linear
! intrpvopt = 2 for quadratic
!
! 4/1/1999 (M. Xue)
!
! Rewrote major portions of the program. Special treatment given
! to grid points below input grid ground. When output grid
! terrain falls below ground, base-state variables are reconstructed
! from averages on horizontal planes to ensure x or y independency
! below ground. This is done only for bglopt=4, however.
! The new code also runs much faster.
!
! 2000/04/05 (Gene Bassett)
! Added bglopt=5, which is similar to bglopt=4 but uses the
! method in ext2arps for constructing the base sounding.
!
! 2000/04/17 (Ming Xue)
! Added an option that allows one to specify input history data
! at a constant time interval.
!
! 2000/05/20 (Gene Bassett)
! Converted to F90, creating allocation and main subroutines.
!
! 2000/07/28 (Ming Xue)
! Converted to F90 free format. Use ALLOCATABLE instead of
! POINTER allocation to avoid double memory usage.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2000/10/31 (Gene Bassett)
! Refined handling of interpolating soil variables.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
! 2001/06/28 (Gene Bassett)
! Added ntagopt, option to set surface winds in new grid to the surface
! winds in the original grid for areas where new terrain is above the
! original terrain.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL, allocatable :: u (:,:,:) ! Total u-velocity (m/s).
REAL, allocatable :: v (:,:,:) ! Total v-velocity (m/s).
REAL, allocatable :: w (:,:,:) ! Total w-velocity (m/s).
REAL, allocatable :: ptprt (:,:,:) ! Perturbation potential temperature
! from that of base state atmosphere (Kelvin).
REAL, allocatable :: pprt (:,:,:) ! Perturbation pressure from that
! of base state atmosphere (Pascal).
REAL, allocatable :: qv (:,:,:) ! Water vapor specific humidity (kg/kg).
REAL, allocatable :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg).
REAL, allocatable :: qr (:,:,:) ! Rain water mixing ratio (kg/kg).
REAL, allocatable :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg).
REAL, allocatable :: qs (:,:,:) ! Snow mixing ratio (kg/kg).
REAL, allocatable :: qh (:,:,:) ! Hail mixing ratio (kg/kg).
REAL, allocatable :: tke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, allocatable :: kmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, allocatable :: kmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
INTEGER, allocatable :: soiltyp(:,:,:) ! Soil type
REAL, allocatable :: stypfrct(:,:,:) ! Fraction of soil type
INTEGER, allocatable :: vegtyp(:,:) ! Vegetation type
REAL, allocatable :: lai (:,:) ! Leaf Area Index
REAL, allocatable :: roufns (:,:) ! Surface roughness
REAL, allocatable :: veg (:,:) ! Vegetation fraction
REAL, allocatable :: ndvi (:,:) ! NDVI
REAL, allocatable :: tsfc (:,:,:) ! Ground sfc. temperature (K)
REAL, allocatable :: tsoil (:,:,:) ! Deep soil temperature (K)
REAL, allocatable :: wetsfc (:,:,:) ! Surface soil moisture
REAL, allocatable :: wetdp (:,:,:) ! Deep soil moisture
REAL, allocatable :: wetcanp (:,:,:) ! Canopy water amount
REAL, allocatable :: snowdpth(:,:) ! Snow depth (m)
REAL, allocatable :: raing(:,:) ! Grid supersaturation rain
REAL, allocatable :: rainc(:,:) ! Cumulus convective rain
REAL, allocatable :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, allocatable :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, allocatable :: radsw (:,:) ! Solar radiation reaching the surface
REAL, allocatable :: rnflx (:,:) ! Net radiation flux absorbed by surface
REAL, allocatable :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, allocatable :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, allocatable :: ptsflx(:,:) ! Surface heat flux (K*kg/(m*s**2))
REAL, allocatable :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s))
REAL, allocatable :: ubar (:,:,:) ! Base state u-velocity (m/s).
REAL, allocatable :: vbar (:,:,:) ! Base state v-velocity (m/s).
REAL, allocatable :: wbar (:,:,:) ! Base state w-velocity (m/s).
REAL, allocatable :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, allocatable :: pbar (:,:,:) ! Base state pressure (Pascal).
REAL, allocatable :: rhobar(:,:,:) ! Base state air density (kg/m**3)
REAL, allocatable :: qvbar (:,:,:) ! Base state water vapor specific humidity
! (kg/kg).
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 :: uprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, allocatable :: vprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, allocatable :: qvprt (:,:,:) ! Perturbation water vapor specific
! humidity (kg/kg)
REAL, allocatable :: tem1 (:,:,:) ! Temporary array
REAL, allocatable :: tem2 (:,:,:) ! Temporary array
REAL, allocatable :: tem3 (:,:,:) ! Temporary array
!wdt update
REAL, allocatable :: tem4 (:,:,:) ! Temporary array
!
!-----------------------------------------------------------------------
!
! Arrays on the new grid:
!
!-----------------------------------------------------------------------
!
REAL, allocatable :: u1 (:,:,:) ! Total u-velocity (m/s).
REAL, allocatable :: v1 (:,:,:) ! Total v-velocity (m/s).
REAL, allocatable :: w1 (:,:,:) ! Total w-velocity (m/s).
REAL, allocatable :: ptprt1(:,:,:) ! Perturbation potential temperature
! from that of base state atmosphere (Kelvin).
REAL, allocatable :: pprt1 (:,:,:) ! Perturbation pressure from that
! of base state atmosphere (Pascal).
REAL, allocatable :: qv1 (:,:,:) ! Water vapor specific humidity (kg/kg).
REAL, allocatable :: qc1 (:,:,:) ! Cloud water mixing ratio (kg/kg).
REAL, allocatable :: qr1 (:,:,:) ! Rain water mixing ratio (kg/kg).
REAL, allocatable :: qi1 (:,:,:) ! Cloud ice mixing ratio (kg/kg).
REAL, allocatable :: qs1 (:,:,:) ! Snow mixing ratio (kg/kg).
REAL, allocatable :: qh1 (:,:,:) ! Hail mixing ratio (kg/kg).
REAL, allocatable :: tke1 (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, allocatable :: kmh1 (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, allocatable :: kmv1 (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
INTEGER, allocatable :: soiltyp1(:,:,:) ! Soil type
REAL, allocatable :: stypfrct1(:,:,:) ! Frction of soil type
INTEGER, allocatable :: vegtyp1 (:,:) ! Vegetation type
REAL, allocatable :: lai1 (:,:) ! Leaf Area Index
REAL, allocatable :: roufns1 (:,:) ! Surface roughness
REAL, allocatable :: veg1 (:,:) ! Vegetation fraction
REAL, allocatable :: tsfc1 (:,:,:) ! Ground sfc. temperature (K)
REAL, allocatable :: tsoil1 (:,:,:) ! Deep soil temperature (K)
REAL, allocatable :: wetsfc1 (:,:,:) ! Surface soil moisture
REAL, allocatable :: wetdp1 (:,:,:) ! Deep soil moisture
REAL, allocatable :: wetcanp1(:,:,:) ! Canopy water amount
REAL, allocatable :: snowdpth1(:,:) ! Snow depth (m)
REAL, allocatable :: raing1(:,:) ! Grid supersaturation rain
REAL, allocatable :: rainc1(:,:) ! Cumulus convective rain
REAL, allocatable :: prcrate1(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, allocatable :: radfrc1(:,:,:) ! Radiation forcing (K/s)
REAL, allocatable :: radsw1 (:,:) ! Solar radiation reaching the surface
REAL, allocatable :: rnflx1 (:,:) ! Net radiation flux absorbed by surface
REAL, allocatable :: usflx1 (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, allocatable :: vsflx1 (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, allocatable :: ptsflx1(:,:) ! Surface heat flux (K*kg/(m*s**2))
REAL, allocatable :: qvsflx1(:,:) ! Surface moisture flux (kg/(m**2*s))
REAL, allocatable :: ubar1 (:,:,:) ! Base state u-velocity (m/s).
REAL, allocatable :: vbar1 (:,:,:) ! Base state v-velocity (m/s).
REAL, allocatable :: ptbar1(:,:,:) ! Base state potential temperature (K)
REAL, allocatable :: pbar1 (:,:,:) ! Base state pressure (Pascal).
REAL, allocatable :: rhobar1(:,:,:) ! Base state air density (kg/m**3)
REAL, allocatable :: qvbar1(:,:,:) ! Base state water vapor specific humidity
! (kg/kg).
REAL, allocatable :: x1 (:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, allocatable :: y1 (:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, allocatable :: z1 (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
!wdt update
REAL, allocatable :: x1_out(:)
REAL, allocatable :: y1_out(:)
REAL, allocatable :: zp1 (:,:,:) ! The physical height coordinate defined at
! w-point on the staggered grid.
REAL, allocatable :: hterain1(:,:) ! Terrain height (m)
!wdt update
REAL, allocatable :: htrn1orig(:,:) ! Terrain height (m)
REAL, allocatable :: j11 (:,:,:) ! Coordinate transformation Jacobian -d(zp)/dx.
REAL, allocatable :: j21 (:,:,:) ! Coordinate transformation Jacobian -d(zp)/dy.
REAL, allocatable :: j31 (:,:,:) ! Coordinate transformation Jacobian d(zp)/dz.
REAL, allocatable :: tem11 (:,:,:) ! Work array
REAL, allocatable :: tem21 (:,:,:) ! Work array
REAL, allocatable :: wgtsx(:,:) ! Weight for interpolation in x-dir for scalar points
REAL, allocatable :: wgtsy(:,:) ! Weight for interpolation in y-dir for scalar points
REAL, allocatable :: wgtux(:,:) ! Weight for interpolation in x-dir for u points
REAL, allocatable :: wgtvy(:,:) ! Weight for interpolation in y-dir for v points
REAL, allocatable :: wgtz(:,:,:,:) ! Weight for interpolation in z-dir for v points
INTEGER, allocatable :: isx(:),jsy(:),iux(:),jvy(:),kz(:,:,:)
REAL, allocatable :: zp1d1 (:) ! Temporary array
REAL, allocatable :: dzp1d1(:) ! Temporary array
!
!-----------------------------------------------------------------------
!
! More arrays.
!
!-----------------------------------------------------------------------
!
REAL, allocatable :: xs (:), ys (:) ! x,y coord for scalar points
REAL, allocatable :: xs1(:), ys1(:) ! x,y coord for scalar points
REAL, allocatable :: xs_2d(:,:), ys_2d(:,:) ! used by subroutine extmnsnd
REAL, allocatable :: temx1yz (:,:,:)
REAL, allocatable :: temx1y1z (:,:,:)
REAL, allocatable :: temx1y1zb(:,:,:)
REAL, allocatable :: temz1d1(:),temz1d2(:),temz1d3(:),temz1d4(:), &
temz1d5(:),temz1d6(:),temz1d7(:)
REAL, allocatable :: zsnd(:),ptsnd(:),qvsnd(:),psnd(:), &
ubsnd(:),vbsnd(:)
REAL, allocatable :: ptpsfc(:,:),ppsfc (:,:),qvsfc(:,:),qcsfc(:,:)
REAL, allocatable :: qrsfc (:,:),qisfc (:,:),qssfc(:,:),qhsfc(:,:)
REAL, allocatable :: tkesfc(:,:),kmhsfc(:,:),kmvsfc(:,:)
REAL, allocatable :: ptbsfc(:,:),pbsfc(:,:),qvbsfc(:,:)
REAL, allocatable :: htrnx1y1(:,:),radsfc(:,:)
INTEGER, allocatable :: ktrnx1y1(:,:)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
INTEGER, allocatable :: ix(:,:),jy(:,:)
REAL, allocatable :: xw(:,:),yw(:,:)
!
!-----------------------------------------------------------------------
! 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
!
!-----------------------------------------------------------------------
!
! Soil types.
!
!-----------------------------------------------------------------------
!
INTEGER :: nstypin ! Number of soil types in input data
INTEGER :: nstyp1 ! Number of soil types for ouput data
!
!-----------------------------------------------------------------------
!
! Misc
!
!-----------------------------------------------------------------------
!
INTEGER :: lvlprof ! Number of levels in 1-d average sounding
!
!-----------------------------------------------------------------------
!
INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,nf,lenfil
PARAMETER (nhisfile_max=200)
CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
INTEGER :: ireturn
NAMELIST /output_dims/ nx1,ny1,nz1
INTEGER :: max_plevels
PARAMETER (max_plevels=1000)
REAL :: plevels(max_plevels)
! Array contain the values of pressure levels to
! which fields will be interpolated.
!
!-----------------------------------------------------------------------
!
! Parameters for the output grid.
!
! Note:
!
! Given nx1, ny1 and nz1, the physical domain 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 :: xorig1, yorig1, zorig1
REAL :: xctr1 , yctr1
REAL :: ctrlat1, ctrlon1
REAL :: origlat, origlon
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).
!
!-----------------------------------------------------------------------
!
! Terrain option parameters for the new grid:
!
!-----------------------------------------------------------------------
!
INTEGER :: ternopt1 ! Model terrain option.
! = 0, no terrain, the ground is flat;
! = 1, bell-shaped mountain;
! = 2, terrain data read from terrain data
! base (not implemented yet)
!wdt-williams 2002-01-17 GMB ntmerge
INTEGER :: ntmerge ! Number of zones over which to merge original
! and new terrain at the bondaries
! (used when ternopt1=4).
INTEGER :: ternfmt1 ! Terrain data file format.
INTEGER :: mntopt1 ! Option for choosing idealized mountain
! type.
REAL :: hmount1 ! The mountain height (m)
REAL :: mntwidx1 ! The half-width of bell-shaped mountain
! in x-dir.
REAL :: mntwidy1 ! The half-width of bell-shaped mountain
! in y-dir.
REAL :: mntctrx1 ! The x coordinate of the bell-shaped
! mountain center.
REAL :: mntctry1 ! The y coordinate of the bell-shaped
! mountain center.
CHARACTER (LEN=132) :: terndta1 ! Name of the terrain data file
REAL :: zflat11,za,zb
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: nunit
INTEGER :: i, j, k
REAL :: amin, amax
REAL :: tmp1,tmp2
REAL :: radnd,pi2,z0
INTEGER :: k0,kk
CHARACTER (LEN=132) :: basdmpfn
INTEGER :: lbasdmpf
CHARACTER (LEN=132) :: ternfn,sfcoutfl,soiloutfl,temchar
INTEGER :: lternfn,lfn
INTEGER :: iss
!wdt-williams 2002-01-17 GMB ntmerge
DOUBLE PRECISION :: ntmergeinv, mfac
INTEGER :: idist
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'bndry.inc'
INCLUDE 'indtflg.inc'
INCLUDE 'alloc.inc'
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
INCLUDE 'exbc.inc'
INTEGER :: houtfmt, nchin, nchout
CHARACTER (LEN=132) :: filename
INTEGER :: grdbas
REAL :: time
!wdt update
INTEGER :: vroutcnt
DATA vroutcnt /0/
INTEGER :: nfile, length, lenstr
CHARACTER (LEN=132) :: timsnd
CHARACTER (LEN=80) :: new_runname
INTEGER :: tmstrln
REAL :: xeps, yeps
REAL :: ctrx,ctry,swx,swy,alatpro(2),sclf,dxscl,dyscl
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
INTEGER :: bglopt,ntagopt
REAL :: misvalue,aghght
REAL :: snddelz,delz,dtdz,tbark1,tbark,ttotk1,ttotk,lnpbar,lnptot
INTEGER :: i1mn,i1mx,j1mn,j1mx
INTEGER :: npoint_below_ground ! flag indicating if any point in the
! output grid terrain is below input grid terrain
INTEGER :: redo_base_state
REAL :: zpmin,zpmax,zpmin1,pbartop, tem,tema,a,b,c,d
!
!-----------------------------------------------------------------------
!
! namelist Declarations:
!
!-----------------------------------------------------------------------
!
NAMELIST /jobname/ runname
INTEGER :: z_or_p ! Control parameter denoting if the output grid is
! in height or pressure coordinates.
INTEGER :: xy_or_ll ! Control parameter denoting if (xctr1,yctr1) or
! (ctrlat1,ctrlon1) are to be used to specify the new
! grid center.
INTEGER :: intrphopt ! Option for horizontal interpolation
! = 1 linear, =2 quadratic
INTEGER :: intrpvopt ! Option for vertical interpolation
! = 1 linear, =2 quadratic
INTEGER :: snap_to_grid, same_res
INTEGER :: istatus, ib,jb
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
INTEGER, PARAMETER :: MAX_GRD=200
INTEGER :: noutgrds
REAL :: xctr_grd(MAX_GRD), yctr_grd(MAX_GRD)
REAL :: clat_grd(MAX_GRD), clon_grd(MAX_GRD)
CHARACTER (LEN=80 ) :: name_grd(MAX_GRD)
INTEGER :: ng
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
REAL :: ctrlat_sv, ctrlon_sv, latitud_sv, dx_sv, dy_sv, dz_sv
REAL :: dzmin_sv, zrefsfc_sv, dlayer1_sv, dlayer2_sv, zflat_sv, strhtune_sv
INTEGER :: strhopt_sv, nstyp_sv
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
NAMELIST /newgrid/ z_or_p,xy_or_ll,strhopt1,xctr1,yctr1, &
ctrlat1,ctrlon1,snap_to_grid,same_res,dx1,dy1,dz1,dzmin1, &
zrefsfc1,dlayer11,dlayer21,strhtune1,zflat1,plevels,nstyp1, &
noutgrds,xctr_grd,yctr_grd,clat_grd,clon_grd,name_grd
!wdt-williams 2002-01-17 GMB ntmerge
NAMELIST /newterrain/ ternopt1,mntopt1,hmount1,mntwidx1,mntwidy1, &
mntctrx1,mntctry1,terndta1,ternfmt1,ntmerge
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
NAMELIST /bgloption/ bglopt,misvalue,intrphopt,intrpvopt,ntagopt,aghght
NAMELIST /sfc_data/ sfcdat,styp,vtyp,lai0,roufns0,veg0, &
sfcdtfl,sfcfmt
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
NAMELIST /output/ dirname,exbcdmp,hdmpfmt,grbpkbit,hdfcompr, &
grdout,basout,varout,mstout,rainout,prcout,iceout, &
tkeout, trbout,sfcout,landout,radout,flxout, &
qcexout,qrexout,qiexout,qsexout,qhexout, &
totout,filcmprs,readyfl,sfcdmp,soildmp,terndmp,ngbrz,zbrdmp
!wdt update
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
REAL :: cpu0, cpu_time, f_cputime
REAL :: xctr1_old,yctr1_old,xsw,ysw,xsw1,ysw1,xsw1_old,ysw1_old
REAL :: ctrlat1_old, ctrlon1_old
!wdt update
REAL :: rhmax = 1.0
REAL :: qvsat
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
REAL :: hght0, dhght, fac
!wdt update
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
cpu0 = f_cputime
()
WRITE(6,'(/9(/2x,a)/)') &
'###############################################################', &
'###############################################################', &
'### ###', &
'### Welcome to ARPSINTRP ###', &
'### This program converts the history dump data ###', &
'### sets generated by ARPS, between various formats. ###', &
'### ###', &
'###############################################################', &
'###############################################################'
!
!-----------------------------------------------------------------------
!
! Read in the input parameters.
!
!-----------------------------------------------------------------------
!
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,nstypin, ireturn)
nstypin = max(1,nstypin)
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
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
WRITE (6,'(/a,g13.3/)') "Current memory allocation (in words):", &
current_memory_use
!
IF (max_plevels < nz1) THEN
WRITE (6,*) "ARPSINTRP: ERROR, nz1 < max_plevels. Stopping."
WRITE (6,*) "nz1 =",nz1," max_plevels = ",max_plevels
STOP 1
END IF
mgrid = 1
nestgrd = 0
!
!-----------------------------------------------------------------------
!
! Set the default parameters
!
!-----------------------------------------------------------------------
!
runname = 'intrp'
z_or_p = 1
xy_or_ll = 1
xctr1 = 512000.0
yctr1 = 512000.0
ctrlat1 = 34.0
ctrlon1 = -98.0
dx1 = 4000.0
dy1 = 4000.0
dz1 = 350.0
strhopt1 = 2
dzmin1 = 20.0
zrefsfc1 = 0.0
dlayer11 = 0.0
dlayer21 = 1.0E5
strhtune1 = 1.0
zflat1 = 1.0E5
nstyp1 = -1
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
noutgrds = 0
xctr_grd = 0
yctr_grd = 0
clat_grd = ctrlat1
clon_grd = ctrlon1
name_grd = 'NULL'
ternopt1 = 3
!wdt-williams 2002-01-17 GMB ntmerge
ntmerge = 1
mntopt1 = 1
hmount1 = 0.000
mntwidx1 = 10000.000
mntwidy1 = 10000.000
mntctrx1 = 10000.000
mntctry1 = 10000.000
terndta1 = 'arpstern.dat'
ternfmt1 = 0
bglopt = 1
misvalue = -9999.0
intrphopt = 1
intrpvopt = 1
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
ntagopt = 0
aghght = 500.0
sfcdat = 1
styp = 3
vtyp = 10
lai0 = 0.31
roufns0 = 0.1
veg0 = 0.0
sfcdtfl = 'arpssfc.data'
sfcfmt = 0
dirname = './'
exbcdmp = 1
hdmpfmt = 1
grbpkbit = 16
hdfcompr = 0
filcmprs = 0
readyfl = 1
basout = 0
grdout = 0
varout = 1
mstout = 1
iceout = 1
tkeout = 1
trbout = 0
rainout = 0
sfcout = 0
snowout = 0
landout = 0
qcexout = 0
qrexout = 0
qiexout = 0
qsexout = 0
qhexout = 0
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
ngbrz = 5
zbrdmp = 10000.0
sfcdmp = 1
soildmp = 1
terndmp = 1
same_res = 0
snap_to_grid = 0
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
new_runname = runname
!
!-----------------------------------------------------------------------
!
! Set the output grid and the variable control parameters
!
!-----------------------------------------------------------------------
!
READ (5,newgrid)
WRITE(6,'(/a/)') 'Sucessfully read namelist block NEWGRID.'
PRINT*
PRINT*,' Input parameters for the new refined grid:'
PRINT*
PRINT*,' z_or_p= ', z_or_p
PRINT*,' xy_or_ll= ', xy_or_ll
PRINT*,' Input dx1:'
PRINT*,' dx1 = ',dx1
PRINT*,' dy1 = ',dy1
PRINT*,' strhopt1= ',strhopt1
PRINT*,' dz1 = ',dz1
PRINT*,' dzmin1 = ',dzmin1
PRINT*,' xctr1 = ',xctr1
PRINT*,' yctr1 = ',yctr1
PRINT*,' ctrlat1 = ',ctrlat1
PRINT*,' ctrlon1 = ',ctrlon1
PRINT*,'plevels =', (plevels(k),k=1,nz1)
PRINT*,' nstyp1 = ',nstyp1
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
PRINT*,'noutgrds = ',noutgrds
IF (noutgrds > MAX_GRD) THEN
WRITE (*,*) "ERROR, noutgrds greater than MAX_GRD. ", &
"Set noutgrds to less than or equal to ",MAX_GRD," ABORTING!"
STOP
ENDIF
IF (noutgrds > 0) THEN
DO ng = 1,noutgrds
PRINT*,'xctr_grd(',ng,') = ',xctr_grd(ng)
PRINT*,'yctr_grd(',ng,') = ',yctr_grd(ng)
PRINT*,'clat_grd(',ng,') = ',clat_grd(ng)
PRINT*,'clon_grd(',ng,') = ',clon_grd(ng)
PRINT*,'name_grd(',ng,') = ',name_grd(ng)
END DO
ENDIF
!
!-----------------------------------------------------------------------
!
! Specify the terrain initialization option parameters:
!
!-----------------------------------------------------------------------
!
READ (5,newterrain,END=100)
WRITE(6,'(/a)') 'Sucessfully read namelist block NEWTERRAIN.'
WRITE(6,'(1x,a,i4)') 'ternopt1 =', ternopt1
!wdt-williams 2002-01-17 GMB ntmerge
WRITE(6,'(1x,a,i4)') 'ntmerge =', ntmerge
WRITE(6,'(1x,a,i4)') 'mntopt1 =', mntopt1
WRITE(6,'(1x,a,f10.3)') 'hmount1 =',hmount1
WRITE(6,'(1x,a,f10.3)') 'mntwidx1=',mntwidx
WRITE(6,'(1x,a,f10.3)') 'mntwidy1=',mntwidy
WRITE(6,'(1x,a,f10.3)') 'mntctrx1=',mntctrx
WRITE(6,'(1x,a,f10.3)') 'mntctry1=',mntctry
lenstr = 132
CALL strlnth( terndta1, lenstr)
WRITE(6,'(2x,a,a/)') 'terndta1 = ',terndta1(1:lenstr)
WRITE(6,'(1x,a,i4)') 'ternfmt1 =', ternfmt1
ternfmt = ternfmt1
READ(5,bgloption,END=100)
WRITE(6,'(/a)') 'Sucessfully read namelist block BLGOPTION.'
WRITE(6,'(1x,a,i4)') 'bglopt =', bglopt
WRITE(6,'(1x,a,f10.3)') 'misvalue=', misvalue
WRITE(6,'(1x,a,i4)') 'intrphopt =', intrphopt
WRITE(6,'(1x,a,i4)') 'intrpvopt =', intrpvopt
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
WRITE(6,'(1x,a,i4)') 'ntagopt =', ntagopt
WRITE(6,'(1x,a,f10.3)') 'aghght =', aghght
!
!-----------------------------------------------------------------------
!
! Set the control parameters for soil data:
!
!-----------------------------------------------------------------------
!
READ (5,sfc_data,END=100)
WRITE(6,'(/a/)') 'Sucessfully read namelist block SOIL_DATA.'
IF (sfcdat .eq. 1) THEN
nstyp = 1
nstypin = nstyp
ELSE
nstyp = nstypin
ENDIF
IF (nstyp1 .eq. -1) nstyp1 = nstyp
WRITE (6,'(a,i3)') "Number of soil types to be used for input data:",nstypin
!
!-----------------------------------------------------------------------
!
! Set the control parameters for output:
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a)') &
' Reading in control parameters for the output data files..'
READ (5,output,END=100)
WRITE(6,'(a)') 'Sucessfully read namelist block OUTPUT.'
houtfmt = hdmpfmt
! IF( houtfmt.eq.0 ) THEN
! write(6,'(/1x,a/)') 'Output format is 0, no data is dumped.'
! STOP 9103
! ENDIF
IF ( houtfmt == 10 .AND. grbpkbit <= 0 ) THEN
WRITE(6,'(a,a,i2/a)') &
'The bit width for packing GRIB data was invalid, ', &
'The old value was ', grbpkbit, 'Reset it to 16 bits'
grbpkbit = 16
END IF
totout = 1
IF (sfcout == 1) THEN
snowout = 1
END IF
GO TO 15
100 WRITE(6,'(a)') &
'Error reading NAMELIST file. Program ARPSINTRP stopped.'
STOP 9104
15 CONTINUE
ldirnam=LEN(dirname)
CALL strlnth( dirname , ldirnam)
lengbf=len_trim(grdbasfn)
WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)
!
!-----------------------------------------------------------------------
!
! Allocate arrays.
!
!-----------------------------------------------------------------------
!
WRITE (6,'(/a,g13.3/)') "Current memory allocation (in words):", &
current_memory_use
!
lvlprof=MAX(601,nz)
nxyz = nx*ny*nz
nxy = nx*ny
allocate( CALL alloc_status_accounting(istatus,nxyz,'u')
allocate( CALL alloc_status_accounting(istatus,nxyz,'v')
allocate( CALL alloc_status_accounting(istatus,nxyz,'w')
allocate(ptprt(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'ptprt')
allocate(pprt(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'pprt')
allocate(qv(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qv')
allocate(qc(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qc')
allocate(qr(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qr')
allocate(qi(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qi')
allocate(qs(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qs')
allocate(qh(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qh')
allocate(tke(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'tke')
allocate(kmh(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'kmh')
allocate(kmv(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'kmh')
allocate(soiltyp(nx,ny,nstypin),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*nstypin,'soiltyp')
allocate(stypfrct(nx,ny,nstypin),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*nstypin,'stypfrct')
allocate(vegtyp(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'vegtyp ')
allocate(lai(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'lai ')
allocate(roufns(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'roufns ')
allocate(veg(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'veg ')
allocate(ndvi(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'ndvi ')
allocate(tsfc(nx,ny,0:nstypin),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'tsfc ')
allocate(tsoil(nx,ny,0:nstypin),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'tsoil ')
allocate(wetsfc(nx,ny,0:nstypin),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'wetsfc ')
allocate(wetdp(nx,ny,0:nstypin),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'wetdp ')
allocate(wetcanp(nx,ny,0:nstypin),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'wetcanp')
allocate(snowdpth(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'snowdpth')
allocate(raing(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'raing')
allocate(rainc(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'rainc')
allocate(prcrate(nx,ny,4),stat=istatus)
CALL alloc_status_accounting(istatus,nxy*4,'prcrate')
allocate(radfrc(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'radfrc')
allocate(radsw(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'radsw ')
allocate(rnflx(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'rnflx ')
allocate(usflx(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'usflx ')
allocate(vsflx(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'vsflx ')
allocate(ptsflx(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'ptsflx')
allocate(qvsflx(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'qvsflx')
allocate(ubar(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'ubar')
allocate(vbar(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'vbar')
allocate(wbar(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'wbar')
allocate(ptbar(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'ptbar')
allocate(pbar(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'pbar')
allocate(rhobar(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'rhobar')
allocate(qvbar(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qvbar')
allocate( CALL alloc_status_accounting(istatus,nx,'x')
allocate( CALL alloc_status_accounting(istatus,ny,'y')
allocate( CALL alloc_status_accounting(istatus,nz,'z')
allocate(zp(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'zp')
allocate(hterain(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'hterain')
allocate(uprt(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'uprt')
allocate(vprt(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'vprt')
allocate(qvprt(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'qvprt')
allocate(tem1(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'tem1')
allocate(tem2(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'tem2')
allocate(tem3(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'tem3')
!wdt update
allocate(tem4(nx,ny,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'tem4')
allocate(u1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'u1')
allocate(v1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'v1')
allocate(w1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'w1')
allocate(ptprt1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'ptprt1')
allocate(pprt1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz,'pprt1')
nxyz1 = nx1*ny1*nz1
nxy1 = nx1*ny1
print*,'Start allocating arrays for output grid.'
allocate(qv1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'qv1')
qv1 = 0.0
allocate(qc1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'qc1')
qc1 = 0.0
allocate(qr1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'qr1')
qr1 = 0.0
allocate(qi1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'qi1')
qi1 = 0.0
allocate(qs1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'qs1')
qs1 = 0.0
allocate(qh1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'qh1')
qh1 = 0.0
allocate(tke1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'tke1')
tke1= 0.0
allocate(kmh1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'kmh1')
kmh1= 0.0
allocate(kmv1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'kmv1')
kmv1= 0.0
allocate(soiltyp1(nx1,ny1,nstyp1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*nstyp1,'soiltyp1')
allocate(stypfrct1(nx1,ny1,nstyp1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*nstyp1,'stypfrct1')
soiltyp1 = 0
stypfrct1 = 0.0
allocate(vegtyp1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'vegtyp1')
vegtyp1 = 0
allocate(lai1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'lai1')
lai1 = 0.0
allocate(roufns1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'roufns1')
roufns1 = 0.0
allocate(veg1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'veg1')
veg1 = 0.0
allocate(tsfc1(nx1,ny1,0:nstyp1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'tsfc1')
allocate(tsoil1(nx1,ny1,0:nstyp1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'tsoil1')
allocate(wetsfc1(nx1,ny1,0:nstyp1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'wetsfc1')
allocate(wetdp1(nx1,ny1,0:nstyp1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'wetdp1')
allocate(wetcanp1(nx1,ny1,0:nstyp1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'wetcanp1')
tsfc1 = 0.0
tsoil1 = 0.0
wetsfc1 = 0.0
wetdp1 = 0.0
wetcanp1 = 0.0
allocate(snowdpth1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'snowdpth1')
snowdpth1 = 0.0
allocate(raing1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'raing1')
raing1 = 0.0
allocate(rainc1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'rainc1')
rainc1 = 0.0
allocate(prcrate1(nx1,ny1,4),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1*4,'prcrate1')
prcrate1 = 0.0
allocate(radfrc1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'radfrc1')
radfrc1 = 0.0
allocate(radsw1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'radsw1')
radsw1 = 0.0
allocate(rnflx1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'rnflx1')
rnflx1 = 0.0
allocate(usflx1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'usflx1')
usflx1 = 0.0
allocate(vsflx1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'vsflx1')
vsflx1 = 0.0
allocate(ptsflx1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'ptsflx1')
ptsflx1 = 0.0
allocate(qvsflx1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qvsflx1')
qvsflx1 = 0.0
allocate(ubar1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'ubar1')
ubar1 = 0.0
allocate(vbar1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'vbar1')
vbar1 = 0.0
allocate(ptbar1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'pbar1')
ptbar1 = 0.0
allocate(pbar1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'ptbar1')
pbar1 = 0.0
allocate(rhobar1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'rhobar1')
rhobar1 = 0.0
allocate(qvbar1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'qvbar1')
qvbar1 = 0.0
allocate(x1(nx1),stat=istatus)
CALL alloc_status_accounting(istatus,nx1,'x1')
x1 = 0.0
allocate(y1(ny1),stat=istatus)
CALL alloc_status_accounting(istatus,ny1,'y1')
y1 = 0.0
!wdt update
allocate(x1_out(nx1),stat=istatus)
CALL alloc_status_accounting(istatus,nx1,'x1_out')
x1_out = 0.0
allocate(y1_out(ny1),stat=istatus)
CALL alloc_status_accounting(istatus,ny1,'y1_out')
y1_out = 0.0
allocate(z1(nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nz1,'z1')
z1 = 0.0
allocate(zp1(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'zp1')
zp1 = 0.0
allocate(hterain1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'hterain1')
hterain1 = 0.0
!wdt update
allocate(htrn1orig(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'htrn1orig')
htrn1orig = 0.0
allocate(j11(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'j11')
j11 = 0.0
allocate(j21(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'j21')
j21 = 0.0
allocate(j31(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'j31')
j31 = 1.0
allocate(tem11(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'tem11')
tem11 = 0.0
allocate(tem21(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'tem21')
tem21 = 0.0
allocate(wgtsx(nx1,3),stat=istatus)
CALL alloc_status_accounting(istatus,nx1*3,'wgtsx')
wgtsx = 0.0
allocate(wgtsy(ny1,3),stat=istatus)
CALL alloc_status_accounting(istatus,ny1*3,'wgtsy')
wgtsy = 0.0
allocate(wgtux(nx1,3),stat=istatus)
CALL alloc_status_accounting(istatus,nx1*3,'wgtux')
wgtux = 0.0
allocate(wgtvy(ny1,3),stat=istatus)
CALL alloc_status_accounting(istatus,ny1*3,'wgtvy')
wgtvy = 0.0
allocate(wgtz(nx1,ny1,nz1,3),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1*3,'wgtz')
wgtz = 0.0
allocate(isx(nx1),stat=istatus)
CALL alloc_status_accounting(istatus,nx1,'isx')
isx = 0.0
allocate(jsy(ny1),stat=istatus)
CALL alloc_status_accounting(istatus,ny1,'jsy')
jsy = 0.0
allocate(iux(nx1),stat=istatus)
CALL alloc_status_accounting(istatus,nx1,'iux')
iux = 0.0
allocate(jvy(ny1),stat=istatus)
CALL alloc_status_accounting(istatus,ny1,'jvy')
jvy = 0.0
allocate(kz(nx1,ny1,nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nxyz1,'kz')
kz = 0.0
allocate(zp1d1(nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nz1,'zp1d1')
zp1d1 = 0.0
allocate(dzp1d1(nz1),stat=istatus)
CALL alloc_status_accounting(istatus,nz1,'dzp1d1')
dzp1d1 = 0.0
allocate(xs(nx ),stat=istatus)
CALL alloc_status_accounting(istatus,nx,'xs')
xs = 0.0
allocate(ys(ny ),stat=istatus)
CALL alloc_status_accounting(istatus,ny,'ys')
ys = 0.0
allocate(xs1(nx1),stat=istatus)
CALL alloc_status_accounting(istatus,nx1,'xs1')
xs1 = 0.0
allocate(ys1(ny1),stat=istatus)
CALL alloc_status_accounting(istatus,ny1,'ys1')
ys1 = 0.0
allocate(xs_2d(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'xs_2d')
xs_2d = 0.0
allocate(ys_2d(nx,ny),stat=istatus)
CALL alloc_status_accounting(istatus,nxy,'ys_2d')
ys_2d = 0.0
allocate(temx1yz(nx1,ny ,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nx1*ny*nz,'temx1yz')
temx1yz = 0.0
allocate(temx1y1z(nx1,ny1,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nx1*ny1*nz,'temx1y1z')
temx1y1z = 0.0
allocate(temx1y1zb(nx1,ny1,nz),stat=istatus)
CALL alloc_status_accounting(istatus,nx1*ny1*nz,'temx1y1zb')
temx1y1zb = 0.0
allocate(temz1d1(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'temz1d1')
temz1d1 = 0.0
allocate(temz1d2(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'temz1d2')
temz1d2 = 0.0
allocate(temz1d3(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'temz1d3')
temz1d3 = 0.0
allocate(temz1d4(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'temz1d4')
temz1d4 = 0.0
allocate(temz1d5(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'temz1d5')
temz1d5 = 0.0
allocate(temz1d6(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'temz1d6')
temz1d6 = 0.0
allocate(temz1d7(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'temz1d7')
temz1d7 = 0.0
allocate(zsnd(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'zsnd')
zsnd = 0.0
allocate(ptsnd(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'ptsnd')
ptsnd = 0.0
allocate(qvsnd(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'qvsnd')
qvsnd = 0.0
allocate(psnd(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'psnd')
psnd = 0.0
allocate(ubsnd(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'ubsnd')
ubsnd = 0.0
allocate(vbsnd(lvlprof),stat=istatus)
CALL alloc_status_accounting(istatus,lvlprof,'vbsnd')
vbsnd = 0.0
allocate(ptpsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'ptpsfc')
ptpsfc = 0.0
allocate(ppsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'ppsfc')
ppsfc = 0.0
allocate(qvsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qvsfc')
qvsfc = 0.0
allocate(qcsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qcsfc')
qcsfc = 0.0
allocate(qrsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qrsfc')
qrsfc = 0.0
allocate(qisfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qisfc')
qisfc = 0.0
allocate(qssfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qssfc')
qssfc = 0.0
allocate(qhsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qhsfc')
qhsfc = 0.0
allocate(tkesfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'tkesfc')
tkesfc = 0.0
allocate(kmhsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'kmhsfc')
kmhsfc = 0.0
allocate(kmvsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'kmvsfc')
kmvsfc = 0.0
allocate(ptbsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'ptbsfc')
ptbsfc = 0.0
allocate(pbsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'pbsfc')
pbsfc = 0.0
allocate(qvbsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'qvbsfc')
qvbsfc = 0.0
allocate(htrnx1y1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'htrnx1y1')
htrnx1y1 = 0.0
allocate(radsfc(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'radsfc')
radsfc = 0.0
allocate(ktrnx1y1(nx1,ny1),stat=istatus)
CALL alloc_status_accounting(istatus,nxy1,'ktrnx1y1')
ktrnx1y1 = 0.0
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
allocate (ix(nx1,ny1),stat=istatus)
allocate (jy(nx1,ny1),stat=istatus)
allocate (xw(nx1,ny1),stat=istatus)
allocate (yw(nx1,ny1),stat=istatus)
nxyz = nx*ny*nz
CALL flzero(u , nxyz)
CALL flzero(v , nxyz)
CALL flzero(uprt , nxyz)
CALL flzero(vprt , nxyz)
CALL flzero(w , nxyz)
CALL flzero(ptprt , nxyz)
CALL flzero(pprt , nxyz)
CALL flzero(qv , nxyz)
CALL flzero(qvprt , nxyz)
CALL flzero(qc , nxyz)
CALL flzero(qr , nxyz)
CALL flzero(qi , nxyz)
CALL flzero(qs , nxyz)
CALL flzero(qh , nxyz)
CALL flzero(tke , nxyz)
CALL flzero(kmh , nxyz)
CALL flzero(kmv , nxyz)
CALL flzero(ubar , nxyz)
CALL flzero(vbar , nxyz)
CALL flzero(wbar , nxyz)
CALL flzero(ptbar , nxyz)
CALL flzero(pbar , nxyz)
CALL flzero(rhobar, nxyz)
CALL flzero(qvbar , nxyz)
CALL flzero(x, nx)
CALL flzero(y, ny)
CALL flzero(z, nz)
CALL flzero(zp , nxyz)
DO iss=0,nstyp
DO j=1,ny
DO i=1,nx
tsfc (i,j,iss) = 0.0
tsoil (i,j,iss) = 0.0
wetsfc (i,j,iss) = 0.0
wetdp (i,j,iss) = 0.0
wetcanp(i,j,iss) = 0.0
END DO
END DO
END DO
DO iss=1,nstyp
DO j=1,ny
DO i=1,nx
soiltyp (i,j,iss) = 0
stypfrct(i,j,iss) = 0.0
END DO
END DO
END DO
vegtyp = 0
snowdpth= 0.0
lai = 0.0
roufns = 0.0
veg = 0.0
ndvi = 0.0
raing = 0.0
rainc = 0.0
prcrate = 0.0
prcrate = 0.0
prcrate = 0.0
prcrate = 0.0
radsw = 0.0
rnflx = 0.0
usflx = 0.0
vsflx = 0.0
ptsflx = 0.0
qvsflx = 0.0
radfrc = 0.0
!
!-----------------------------------------------------------------------
!
! Run interpolation.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Loop over data files
!
!-----------------------------------------------------------------------
!
ireturn = 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.'
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
102 CONTINUE
CALL dtaread(nx,ny,nz,nstyp, &
hinfmt, nchin,grdbasfn(1:lengbf),lengbf, &
filename(1:lenfil),lenfil,time, &
x,y,z,zp, uprt ,vprt ,w,ptprt, pprt , &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc, tsoil, wetsfc, wetdp, wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1, tem2, tem3)
!wdt update moved code block to here
!wdt begin block
curtim = time
IF( hinfmt == 9 .AND. ireturn == 2 ) THEN
WRITE(6,'(/1x,a/)') 'The end of GrADS file was reached.'
CLOSE ( nchin )
CALL retunit( nchin )
GO TO 9001
END IF
IF( ireturn /= 0 ) GO TO 9002 ! Read was unsuccessful
!
!-----------------------------------------------------------------------
!
! Grid and base related calculations that need to be done only once.
!
!-----------------------------------------------------------------------
!
IF( nfile == 1) THEN ! Need to do this only once.
!
!-----------------------------------------------------------------------
!
! Set-up land surface property data, if requested.
!
!-----------------------------------------------------------------------
!
IF ( sfcdat == 0 ) GO TO 110
IF ( sfcdat == 1 ) THEN
nstyp = 1
DO j=1, ny-1
DO i=1, nx-1
soiltyp(i,j,1) = styp
vegtyp (i,j) = vtyp
lai (i,j) = lai0
roufns (i,j) = roufns0
veg (i,j) = veg0
END DO
END DO
ELSE IF (sfcdat == 2 .OR. (sfcdat == 3.AND.landin /= 1) ) THEN
!
!-----------------------------------------------------------------------
!
! Read the surface property data from file sfcdtfl (passed
! in globcst.inc).
!
!-----------------------------------------------------------------------
!
CALL readsfcdt( nx,ny,nstyp,sfcdtfl,dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi )
ELSE IF(sfcdat == 3 .AND. landin == 1) THEN
WRITE(6,'(1x,a/,a/)') &
'Surface property data in the history file was used.', &
'Data in ',sfcdtfl,' ignored.'
ELSE
WRITE(6,'(1x,a,i3,a/)') &
'Invalid surface data input option. sfcdat =',sfcdat, &
'. Program stopped in ARPSINTRP'
STOP
END IF ! sfcdat option
IF ( nstyp == 1 ) THEN
DO j=1,ny
DO i=1,nx
stypfrct(i,j,1) = 1.0
END DO
END DO
END IF
110 CONTINUE
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)
xorig = x(2)
yorig = y(2)
zorig = z(2)
WRITE(6,'(1x,a,3f15.3)') 'xorig, yorig, zorig =', &
xorig, yorig, zorig
ebc = 3
wbc = 3
sbc = 3
nbc = 3
dx = x(2)-x(1)
dy = y(2)-y(1)
dz = z(2)-z(1)
END IF
!wdt end block
!
!-----------------------------------------------------------------------
!
! 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
!wdt update: moved code block
IF (same_res == 1) THEN
IF(dx /= dx1 .OR. dy /= dy1) THEN
WRITE (*,*) "WARNING: resetting dx1 and/or dy1 to match input grid."
ENDIF
dx1 = dx
dy1 = dy
ENDIF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
ng = 0
500 CONTINUE ! Loop here for each output grid (noutgrds)
IF (noutgrds > 0) THEN
ng = ng + 1
ctrlat1 = clat_grd(ng)
ctrlon1 = clon_grd(ng)
xctr1 = xctr_grd(ng)
yctr1 = yctr_grd(ng)
ENDIF
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.
CALL xytoll(1,1,(FLOAT(nx-3)/2.)*dx,(FLOAT(ny-3)/2.)*dy, &
tmp1,tmp2)
PRINT*,'ctrlat=',tmp1,', ctrlon=',tmp2
IF (xy_or_ll == 1) THEN
CALL xytoll(1,1,xctr1,yctr1, ctrlat1,ctrlon1)
ELSE IF (xy_or_ll == 2) THEN
CALL lltoxy(1,1, ctrlat1,ctrlon1, xctr1,yctr1)
ELSE
WRITE (*,*) "ERROR: xy_or_ll =",xy_or_ll," not a valid value."
STOP
END IF
IF (same_res == 1) THEN
IF(dx /= dx1 .OR. dy /= dy1) THEN
WRITE (*,*) "WARNING: resetting dx1 and/or dy1 to match input grid."
ENDIF
dx1 = dx
dy1 = dy
ENDIF
ctrlat1_old = ctrlat1
ctrlon1_old = ctrlon1
IF( snap_to_grid /= 0 .OR. same_res == 1) THEN
xsw = x(1)
ysw = y(1)
xsw1 = xctr1-(nx1-1)*dx1*0.5
ysw1 = yctr1-(ny1-1)*dy1*0.5
xsw1_old = xsw1
ysw1_old = ysw1
PRINT*,'mod =', MOD(dx,dx1),MOD(dy,dy1)
IF(((dx >= dx1 .AND. &
ABS((MOD(dx,dx1)-dx1)*MOD(dx,dx1)) <= 0.001*dx1).OR. &
(dx <= dx1 .AND. &
ABS((MOD(dx1,dx)-dx )*MOD(dx1,dx)) <= 0.001*dx )).AND. &
((dy >= dy1 .AND. &
ABS((MOD(dy,dy1)-dy1)*MOD(dy,dy1)) <= 0.001*dy1).OR. &
(dy <= dy1 .AND. &
ABS((MOD(dy1,dy)-dy )*MOD(dy1,dy)) <= 0.001*dy ) )) THEN
PRINT*,'is multiple'
IF(snap_to_grid == 1) THEN
xsw1 = xsw + nint( (xsw1_old-xsw)/dx1 )*dx1
ysw1 = ysw + nint( (ysw1_old-ysw)/dy1 )*dy1
ELSE IF(snap_to_grid == 2) THEN
xsw1 = xsw + nint( (xsw1_old-xsw)/dx )*dx
ysw1 = ysw + nint( (ysw1_old-ysw)/dy )*dy
END IF
PRINT*,'xsw,ysw,xsw1,ysw1=',xsw,ysw,xsw1,ysw1
PRINT*,'xsw1_old,ysw1_old=',xsw1_old,ysw1_old
IF((snap_to_grid /= 0) .AND. ( &
(ABS((xsw1_old - xsw1)) > 0.00001*dx1).OR. &
(ABS((ysw1_old - ysw1)) > 0.00001*dy1) )) THEN
xctr1_old = xctr1
yctr1_old = yctr1
xctr1 = xsw1+(nx1-1)*dx1*0.5
yctr1 = ysw1+(ny1-1)*dy1*0.5
!
! Note there is a possibility that the shifted grid is out of bound
! The checking is done later.
!
WRITE(6,'(/1x,a/)') &
'################################################################'
WRITE(6,'(/1x,a/,2(/1x,a)/)') &
' ATTENTION:', &
'xctr1 and/or yctr1 reset so that some of the output ', &
'grid lines match those of input grid.'
WRITE(6,'(1x,2(a,f17.5))')'Original xctr1_old=',xctr1_old, &
', yctr1_old=',yctr1_old
WRITE(6,'(1x,2(a,f17.5))')'New xctr1 = ',xctr1, &
', yctr1 = ',yctr1
WRITE(6,'(/1x,a/)') &
'################################################################'
CALL xytoll(1,1,xctr1,yctr1, ctrlat1,ctrlon1)
END IF
ELSE
same_res = 0 ! ???
END IF
END IF
IF( nz1 /= nz .AND. same_res == 1) THEN
WRITE (*,*) "WARNING: nz1 (",nz1,") does not match nz (",nz,")."
WRITE (*,*) " Resetting same_res to 0."
same_res = 0
ENDIF
WRITE(6,'(1x,2(a,f17.5))') &
'Original ctrlat1_old = ',ctrlat1_old,', ctrlon1_old = ', &
ctrlon1_old
WRITE(6,'(1x,2(a,f17.5))') &
'Final ctrlat1 = ',ctrlat1 ,', ctrlon1 = ',ctrlon1
WRITE(6,'(1x,2(a,f17.5))') &
'Final xctr1 = ',xctr1, ', yctr1 = ',yctr1
WRITE(6,'(a,i2/)')'same_res = ', same_res
!wdt update: moved code block to above
!
!-----------------------------------------------------------------------
!
! Perform spatial interpolations.
!
!-----------------------------------------------------------------------
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
! IF( nfile == 1) THEN
IF( nfile == 1 .or. noutgrds > 0) THEN
! Need to do this for the first file only
!wdt update moved code block here
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.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
ib= nint(x1(1)/dx)+1
jb= nint(y1(1)/dy)+1
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
DO k=1,nz1
z1(k)=zorig1+(k-2.0)*dz1
END DO
IF( nfile == 1) THEN ! Need to do this for the first file only
xeps = 0.01*dx
yeps = 0.01*dy
PRINT*,'x (2),x (nx -1),y (2),y (ny -1) =', &
x(2),x(nx-1),y(2),y(ny-1)
PRINT*,'x1(2),x1(nx1-1),y1(2),y1(ny1-1) =', &
x1(2),x1(nx1-1),y1(2),y1(ny1-1)
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(/2x,a),/2x,2(a,f12.4),2(a,i5),2(a,f12.4),/2x,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.'
GOTO 600 !wdt
STOP 1001
END IF
END IF
IF( same_res == 0 ) THEN
IF( intrphopt == 1) THEN ! Weights for 1st order interpolation
DO i=1,nx1-1
isx(i) = MAX(1, MIN(nx-2, INT((xs1(i)-xs(1))/dx)+1 ))
wgtsx(i,1)= (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,1)= (ys(jsy(j)+1)-ys1(j))/(ys(jsy(j)+1)-ys(jsy(j)))
END DO
DO i=1,nx1
iux(i) = MAX(1, MIN(nx-1, INT((x1 (i)-x (1))/dx)+1 ))
wgtux(i,1)= (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,1)= (y (jvy(j)+1)-y1 (j))/(y (jvy(j)+1)-y (jvy(j)))
END DO
ELSE ! Weights for second order interpolation
DO i=1,nx1-1
isx(i) = MAX(2, MIN(nx-2, INT((xs1(i)-xs(1))/dx)+1 ))
d=xs1(i)
a=xs(isx(i)-1)
b=xs(isx(i))
c=xs(isx(i)+1)
wgtsx(i,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtsx(i,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtsx(i,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
DO j=1,ny1-1
jsy(j) = MAX(2, MIN(ny-2, INT((ys1(j)-ys(1))/dy)+1 ))
d=ys1(j)
a=ys(jsy(j)-1)
b=ys(jsy(j))
c=ys(jsy(j)+1)
wgtsy(j,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtsy(j,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtsy(j,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
DO i=1,nx1
iux(i) = MAX(2, MIN(nx-1, INT((x1 (i)-x (1))/dx)+1 ))
d=x1(i)
a=x(iux(i)-1)
b=x(iux(i))
c=x(iux(i)+1)
wgtux(i,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtux(i,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtux(i,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
DO j=1,ny1
jvy(j) = MAX(2, MIN(ny-1, INT((y1 (j)-y (1))/dy)+1 ))
d=y1(j)
a=y(jvy(j)-1)
b=y(jvy(j))
c=y(jvy(j)+1)
wgtvy(j,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtvy(j,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtvy(j,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
END IF ! intrphopt = ?
IF( ternopt1 == 0 ) THEN
DO i=1,nx1-1
DO j=1,ny1-1
hterain1(i,j) = zrefsfc1
END DO
END DO
ELSE IF( ternopt1 == 1 ) THEN ! Bell-shaped mountain
!
!-----------------------------------------------------------------------
!
! Define the bell-shaped mountain
!
!-----------------------------------------------------------------------
!
pi2 = 1.5707963267949
DO j=1,ny1-1
DO i=1,nx1-1
IF( mntwidy1 < 0.0 ) THEN ! 2-d terrain in x-z plane.
radnd = 1.0+((xs1(i)-mntctrx1)/mntwidx1)**2
ELSE IF( mntwidx1 < 0.0 ) THEN ! 2-d terrain in y-z plane.
radnd = 1.0+((ys1(j)-mntctry1)/mntwidy1)**2
ELSE ! 3-d terrain
radnd = 1.0+((xs1(i)-mntctrx1)/mntwidx1)**2 &
+((ys1(j)-mntctry1)/mntwidy1)**2
END IF
hterain1(i,j) = zrefsfc1 + hmount1/radnd
END DO
END DO
ELSE IF( ternopt1 == 2 .or. ternopt1 == 4 ) THEN
! Read from terrain data base
!
!-----------------------------------------------------------------------
!
! Read in the terrain data.
!
!-----------------------------------------------------------------------
!
lenstr = 132
CALL strlnth( terndta1, lenstr)
CALL readtrn( nx1,ny1,dx1,dy1, terndta1(1:lenstr), &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat1,ctrlon1, &
hterain1 )
!
!-----------------------------------------------------------------------
!
! Define the new grid terrain by interpolation
!
!-----------------------------------------------------------------------
!
!wdt update
END IF
!wdt update
IF( ternopt1 == 3 .or. ternopt1 == 4 ) THEN
!wdt update htrn1orig
CALL intrpxy3d(hterain,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
htrn1orig,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!wdt-williams 2002-01-17 GMB ntmerge
!hterain1 = htrn1orig ! uncomment if ternopt1=4 (ntmerge) dissabled
IF (ternopt1 == 3) THEN
hterain1 = htrn1orig
ELSE
ntmergeinv = 1.d0/ntmerge
DO j = 1,ny1
DO i = 1,nx1
idist = max(0,min(ntmerge,i-2,nx1-2-i,j-2,ny1-2-j))
mfac = idist*ntmergeinv
hterain1(i,j) = (1.d0-mfac)*htrn1orig(i,j) &
+ mfac*hterain1(i,j)
END DO
END DO
END IF
END IF
npoint_below_ground = 0
IF( ternopt1 /= 3) THEN
CALL intrpxy3d(hterain,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
temx1y1z,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL a3dmax0(hterain,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'htmin = ', amin,', htmax =',amax
CALL a3dmax0(temx1y1z,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ht1min = ', amin,', ht1max =',amax
CALL a3dmax0(hterain1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ht1min = ', amin,', ht1max =',amax
npoint_below_ground = 0
DO i=1,nx1-1
DO j=1,ny1-1
IF( hterain1(i,j) < temx1y1z(i,j,1)-1.0E-2) THEN
npoint_below_ground = 1
GO TO 2100
END IF
END DO
END DO
2100 CONTINUE
END IF
!
!-----------------------------------------------------------------------
!
! 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 in ARPSINTRP.'
STOP 9105
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 INIGRD.'
STOP 9106
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 mountain 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
525 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,x1,y1,z1,zp1,j11,j21,j31)
END IF ! same_res = 0?
END IF ! IF (nfile.eq.1)
!
!-----------------------------------------------------------------------
!
! Calculate total fields from that for base state and perturbations
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
u (i,j,k)=uprt (i,j,k)+ubar (i,j,k)
v (i,j,k)=vprt (i,j,k)+vbar (i,j,k)
qv(i,j,k)=qvprt(i,j,k)+qvbar(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Start interpolation or sampling in the case of same_res=1.
!
!-----------------------------------------------------------------------
!
IF( same_res == 1 ) THEN
DO k=1,nz1
DO j=1,ny1
DO i=1,nx1
u1 (i,j,k)=u (i+ib,j+jb,k)
v1 (i,j,k)=v (i+ib,j+jb,k)
w1 (i,j,k)=w (i+ib,j+jb,k)
ptprt1(i,j,k)=ptprt(i+ib,j+jb,k)
pprt1 (i,j,k)=pprt (i+ib,j+jb,k)
qv1 (i,j,k)=qv (i+ib,j+jb,k)
qc1 (i,j,k)=qc (i+ib,j+jb,k)
qr1 (i,j,k)=qr (i+ib,j+jb,k)
qi1 (i,j,k)=qi (i+ib,j+jb,k)
qs1 (i,j,k)=qs (i+ib,j+jb,k)
qh1 (i,j,k)=qh (i+ib,j+jb,k)
tke1 (i,j,k)=tke (i+ib,j+jb,k)
kmh1 (i,j,k)=kmh (i+ib,j+jb,k)
kmv1 (i,j,k)=kmv (i+ib,j+jb,k)
radfrc1(i,j,k)=radfrc(i+ib,j+jb,k)
rhobar1(i,j,k)=rhobar(i+ib,j+jb,k)
END DO
END DO
END DO
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
! IF( nfile == 1) THEN
IF( nfile == 1 .or. noutgrds > 0) THEN
DO k=1,nz1
DO i=1,nx1
DO j=1,ny1
ubar1 (i,j,k)=ubar (i+ib,j+jb,k)
vbar1 (i,j,k)=vbar (i+ib,j+jb,k)
ptbar1(i,j,k)=ptbar(i+ib,j+jb,k)
pbar1 (i,j,k)=pbar (i+ib,j+jb,k)
qvbar1(i,j,k)=qvbar(i+ib,j+jb,k)
zp1 (i,j,k)=zp (i+ib,j+jb,k)
END DO
END DO
END DO
CALL jacob(nx1,ny1,nz1,x1,y1,z1,zp1,j11,j21,j31)
END IF
DO i=1,nx1
DO j=1,ny1
hterain1(i,j)=hterain(i+ib,j+jb)
soiltyp1 (i,j,:)=soiltyp (i+ib,j+jb,:)
stypfrct1(i,j,:)=stypfrct(i+ib,j+jb,:)
vegtyp1(i,j)=vegtyp(i+ib,j+jb)
lai1 (i,j)=lai (i+ib,j+jb)
roufns1(i,j)=roufns(i+ib,j+jb)
veg1 (i,j)=veg (i+ib,j+jb)
tsfc1 (i,j,:)=tsfc (i+ib,j+jb,:)
tsoil1 (i,j,:)=tsoil (i+ib,j+jb,:)
wetsfc1 (i,j,:)=wetsfc (i+ib,j+jb,:)
wetdp1 (i,j,:)=wetdp (i+ib,j+jb,:)
wetcanp1 (i,j,:)=wetcanp (i+ib,j+jb,:)
snowdpth1(i,j)=snowdpth(i+ib,j+jb)
raing1 (i,j)=raing(i+ib,j+jb)
rainc1 (i,j)=rainc(i+ib,j+jb)
prcrate1(i,j,1:4)=prcrate(i+ib,j+jb,1:4)
radsw1 (i,j)=radsw (i+ib,j+jb)
rnflx1 (i,j)=rnflx (i+ib,j+jb)
usflx1 (i,j)=usflx (i+ib,j+jb)
vsflx1 (i,j)=vsflx (i+ib,j+jb)
ptsflx1(i,j)=ptsflx(i+ib,j+jb)
qvsflx1(i,j)=qvsflx(i+ib,j+jb)
END DO
END DO
GO TO 800
END IF
!
!-----------------------------------------------------------------------
!
! Intepolate the 2-D arrays to the new grid
!
!-----------------------------------------------------------------------
!
CALL intrpxy3d(rainc,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
rainc1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(raing,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
raing1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(prcrate,nx,1,nx-1,ny,1,ny-1,4,1,4, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
prcrate1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(radsw,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
radsw1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(rnflx,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
rnflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(usflx,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
usflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(vsflx,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
vsflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(ptsflx,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
ptsflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(qvsflx,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
qvsflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(lai,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
lai1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(roufns,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
roufns1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(veg,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
veg1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
DO j = 1,ny1
DO i = 1,nx1
ix(i,j) = isx(i)
jy(i,j) = jsy(j)
xw(i,j) = &
MAX(0.,MIN(1.,(xs(isx(i)+1)-xs1(i))/(xs(isx(i)+1)-xs(isx(i)))))
yw(i,j) = &
MAX(0.,MIN(1.,(ys(jsy(j)+1)-ys1(j))/(ys(jsy(j)+1)-ys(jsy(j)))))
END DO
END DO
CALL intrp_soil(nx,ny,nx1,ny1,nstypin,nstyp1,xw,yw,ix,jy, &
tsfc,tsoil,wetsfc,wetdp,wetcanp, &
soiltyp,stypfrct,vegtyp, &
tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1, &
soiltyp1,stypfrct1,vegtyp1)
! CALL intrpxy3d(tsfc,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1, &
! wgtsx,isx,wgtsy,jsy,intrphopt, &
! tsfc1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
! CALL intrpxy3d(tsoil,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1, &
! wgtsx,isx,wgtsy,jsy,intrphopt, &
! tsoil1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
! CALL intrpxy3d(wetsfc,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1, &
! wgtsx,isx,wgtsy,jsy,intrphopt, &
! wetsfc1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
! CALL intrpxy3d(wetdp,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1, &
! wgtsx,isx,wgtsy,jsy,intrphopt, &
! wetdp1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
! CALL intrpxy3d(wetcanp,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1, &
! wgtsx,isx,wgtsy,jsy,intrphopt, &
! wetcanp1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
CALL intrpxy3d(snowdpth,nx,1,nx-1,ny,1,ny-1,1,1,1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
snowdpth1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
! DO i=1,nx1-1
! DO j=1,ny1-1
! vegtyp1 (i,j) = vegtyp (isx(i),jsy(j))
! END DO
! END DO
!
! DO iss=1,nstyp
! DO j=1,ny1-1
! DO i=1,nx1-1
! soiltyp1 (i,j,iss) = soiltyp (isx(i),jsy(j),iss)
! stypfrct1(i,j,iss) = stypfrct(isx(i),jsy(j),iss)
! END DO
! END DO
! END DO
!
!-----------------------------------------------------------------------
!
! Intepolate 3D arrays
!
!-----------------------------------------------------------------------
!
IF( z_or_p == 2) THEN ! This situation has not been tested since
! last major change to this program
DO k=1,nz
DO j=1,ny-1
DO i=1,nx-1
zp(i,j,k) = -ALOG( (pbar(i,j,k)+pprt(i,j,k))*1.0E-5 )
END DO
END DO
END DO
DO k=1,nz1
DO j=1,ny1-1
DO i=1,nx1-1
zp1(i,j,k) = -ALOG( plevels(k)*1.0E-3 )
END DO
END DO
END DO
END IF
IF( npoint_below_ground /= 0.AND.z_or_p == 1.AND. &
(bglopt == 4 .OR. bglopt == 5)) THEN
! Reconstruct horizontally homogeneous base state
redo_base_state = 1
ELSE
redo_base_state = 0
END IF
IF( redo_base_state /= 0 ) THEN
zpmin=0.5*(zp(2,2, 2)+zp(2,2,3))
zpmax=0.5*(zp(2,2,nz)+zp(2,2,nz-1))
DO i=2,nx-2
DO j=2,ny-2
zpmin=MIN(zpmin,0.5*(zp(i,j, 2)+zp(i,j, 3)))
zpmax=MAX(zpmax,0.5*(zp(i,j,nz)+zp(i,j,nz-1)))
END DO
END DO
zpmin1=MIN(zpmin1,0.5*(zp1(1,1,1)+zp1(1,1,2)))
DO i=2,nx1-2
DO j=2,ny1-2
zpmin1=MIN(zpmin1,0.5*(zp1(i,j,1)+zp1(i,j,2)))
END DO
END DO
snddelz = (zpmax-zpmin)/(lvlprof-2)
DO k=1,lvlprof
zsnd(k)= zpmin+(k-2)*snddelz
END DO
IF (bglopt == 4) THEN
DO k=2,lvlprof
temz1d3(k)=0.0
temz1d4(k)=0.0
temz1d5(k)=0.0
temz1d6(k)=0.0
temz1d7(k)=0.0
END DO
pbartop = 0.0
DO i=2,nx-2
DO j=2,ny-2
DO k=2,nz-1
temz1d1(k)=ptbar(i,j,k)
temz1d2(k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
CALL inte1d(temz1d1(2),temz1d2(2),nz-2, &
ptsnd(2),zsnd(2),lvlprof-1)
DO k=2,nz-1
temz1d1(k)=qvbar(i,j,k)
END DO
CALL inte1d(temz1d1(2),temz1d2(2),nz-2, &
qvsnd(2),zsnd(2),lvlprof-1)
DO k=2,nz-1
temz1d1(k)=(ubar(i,j,k)+ubar(i+1,j,k))*0.5
END DO
CALL inte1d(temz1d1(2),temz1d2(2),nz-2, &
ubsnd(2),zsnd(2),lvlprof-1)
DO k=2,nz-1
temz1d1(k)=(vbar(i,j,k)+vbar(i,j+1,k))*0.5
END DO
CALL inte1d(temz1d1(2),temz1d2(2),nz-2, &
vbsnd(2),zsnd(2),lvlprof-1)
DO k=2,nz-1
temz1d1(k)=ALOG(pbar(i,j,k))
END DO
CALL inte1d(temz1d1(2),temz1d2(2),nz-2, &
tem,zsnd(lvlprof-1),1)
pbartop = pbartop+tem
DO k=2,lvlprof
IF(zsnd(k) >= temz1d2(2)-1.0E-5) THEN
! Add only points above ground
temz1d3(k)=temz1d3(k)+1
temz1d4(k)=temz1d4(k)+ptsnd(k)
temz1d5(k)=temz1d5(k)+qvsnd(k)
temz1d6(k)=temz1d6(k)+ubsnd(k)
temz1d7(k)=temz1d7(k)+vbsnd(k)
END IF
END DO
END DO
END DO
!
! Obtain averages
!
DO k=2,lvlprof
ptsnd(k)=temz1d4(k)/temz1d3(k)
qvsnd(k)=temz1d5(k)/temz1d3(k)
ubsnd(k)=temz1d6(k)/temz1d3(k)
vbsnd(k)=temz1d7(k)/temz1d3(k)
END DO
pbartop=EXP(pbartop/temz1d3(lvlprof-1))
!
! Integrate hydrostatic equation down from top
!
DO k=2,lvlprof
temz1d1(k) = ptsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
END DO
psnd(lvlprof-1) = ( pbartop/p0 )**(rddcp)
DO k=lvlprof-2,2,-1
psnd(k)=psnd(k+1)-g*(zsnd(k)-zsnd(k+1)) &
/(cp*0.5*(temz1d1(k)+temz1d1(k+1)))
END DO
DO k=lvlprof,lvlprof
psnd(k)=psnd(k-1)-g*(zsnd(k)-zsnd(k-1)) &
/(cp*0.5*(temz1d1(k)+temz1d1(k-1)))
END DO
DO k = 2,lvlprof
psnd(k) = psnd(k)**(1/rddcp) * p0
END DO
IF( zpmin1 < zpmin ) THEN ! Extend 1d sounding to zpmin1 using
! constant temperature lapse rate
dtdz = -0.0068 ! Minus lapse rate
zsnd(1)=zpmin1
delz = zpmin1 - zpmin
tbark1 = ptsnd(2) * (psnd(2)/p0)**rddcp
tbark = tbark1+dtdz*delz
lnpbar = ALOG(psnd(2))-g*2.0*delz/(rd*(tbark1+tbark))
psnd(1) = EXP(lnpbar)
ptsnd(1)=tbark*(p0/psnd(1))**rddcp
qvsnd(1)=qvsnd(2)
ubsnd(1)=ubsnd(2)
vbsnd(1)=vbsnd(2)
ELSE ! Set for safty, won't be used.
psnd(1)=psnd(2)
ptsnd(1)=ptsnd(2)
qvsnd(1)=qvsnd(2)
ubsnd(1)=ubsnd(2)
vbsnd(1)=vbsnd(2)
END IF
ELSE ! bglopt=5
!wdt update
!
! Use extmnsnd (from ext2arps) to get average sounding
!
! tem1 - zp at zone center
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k) = 0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
END DO
END DO
! tem2 = log(pbar)
! tem3 = t
! tem4 = rhstar
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem2(i,j,k) = ALOG(pbar(i,j,k)) ! log pressure
tem3(i,j,k) = (ptprt(i,j,k)+ptbar(i,j,k)) &
*(((pprt(i,j,k)+pbar(i,j,k))/p0)**rddcp) ! temperature
qvsat=f_qvsat( pprt(i,j,k)+pbar(i,j,k), tem3(i,j,k) )
tem4(i,j,k)=SQRT(AMAX1(0.,(rhmax-(qv(i,j,k)/qvsat))))
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
xs_2d(i,j) = xs(i)
ys_2d(i,j) = ys(j)
END DO
END DO
i1mn = MAX(1,INT(xorig1/dx1-1.))
i1mx = MIN(nx-1,i1mn+nx1+1)
j1mn = MAX(1,INT(yorig1/dy1-1.))
j1mx = MIN(ny-1,j1mn+ny1+1)
!wdt update
CALL extmnsnd(nx1,ny1,nx,ny,nz,lvlprof, &
i1mn,i1mx,j1mn,j1mx,2,nz-1, &
xs1,ys1,xs_2d,ys_2d, &
tem1,tem2,tem3,tem4, &
u,v, & ! note: not vital that u & v be staggered
zsnd,temz1d1,psnd,temz1d2,ptsnd,temz1d3, &
qvsnd,temz1d4,ubsnd,vbsnd)
END IF
WRITE(6,'(/a/)')' Reconstructed horizontal average sounding:'
DO k = lvlprof,1,-1
WRITE(6,'(a,i4,4f10.2,f10.6)') 'k,z,p,pt,t,qv=', &
k,zsnd(k),psnd(k),ptsnd(k), &
ptsnd(k)*(psnd(k)/p0)**rddcp,qvsnd(k)
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Intepolate 3D arrays at u-point
!
!-----------------------------------------------------------------------
!
CALL avgz(zp1, 0 , nx1,ny1,nz1,1,nx1-1,1,ny1-1,1,nz1-1, tem11 )
CALL avgsu(tem11,nx1,ny1,nz1,1,ny1-1,1,nz1-1, tem21, tem11)
CALL avgz(zp, 0 , nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, tem1 )
CALL avgsu(tem1,nx,ny,nz,1,ny-1,1,nz-1, tem2, tem1)
CALL intrpxy3d(tem2,nx,1,nx,ny,1,ny-1,nz,1,nz-1, &
wgtux,iux,wgtsy,jsy,intrphopt, &
temx1y1zb,nx1,1,nx1, ny1,1,ny1-1, temx1yz)
DO i=1,nx1
DO j=1,ny1-1
k0 = nz-2
DO k=nz1-1,1,-1
z0 = tem21(i,j,k)
DO kk=k0,1,-1
IF(z0 >= temx1y1zb(i,j,kk)) THEN
k0 = kk
EXIT
END IF
END DO
1115 CONTINUE
IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
END DO
END DO
END DO
IF( intrpvopt == 1) THEN
DO i=1,nx1
DO j=1,ny1-1
DO k=1,nz1-1
wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/ &
(temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
END DO
END DO
END DO
ELSE
DO i=1,nx1
DO j=1,ny1-1
DO k=1,nz1-1
d=tem21(i,j,k)
a=temx1y1zb(i,j,kz(i,j,k)-1)
b=temx1y1zb(i,j,kz(i,j,k))
c=temx1y1zb(i,j,kz(i,j,k)+1)
wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
END DO
END DO
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
! IF( nfile == 1) THEN
IF(nfile == 1 .or. noutgrds > 0) THEN
IF( redo_base_state == 0 ) THEN ! get base-state from original arrays
CALL intrpxyz3d(ubar,nx,1,nx,ny,1,ny-1,nz,1,nz-1, &
wgtux,iux,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
ubar1,nx1,1,nx1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO i=1,nx1
DO j=1,ny1-1
DO k=1,nz1-1
IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
IF(bglopt == 2) THEN
ubar1(i,j,k)=temx1y1z(i,j,2)
ELSE IF(bglopt == 3) THEN
ubar1(i,j,k)=0.0
END IF
END IF
END DO
END DO
END DO
ELSE IF( redo_base_state == 1 ) THEN ! reconstruct base-state from 1d soundings
DO k=1,nz1-1
DO j=1,ny1-1
DO i=1,nx1
kk = MAX(1,MIN(lvlprof-1, &
INT((tem21(i,j,k)-zsnd(2))/snddelz)+2))
IF(tem21(i,j,k) < zsnd(2)) kk = 1
tem = (tem21(i,j,k)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
ubar1(i,j,k)= ubsnd(kk)+(ubsnd(kk+1)-ubsnd(kk))*tem
END DO
END DO
END DO
END IF
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
IF (ntagopt == 1) THEN ! adjust tem21 (heights) so that surface wind is
! same as input surface wind (by matching tem21
! to the input height, temx1y1zb, near the surface).
DO i=1,nx1
DO j=1,ny1-1
hght0 = tem21(i,j,2)
dhght = hght0 - temx1y1zb(i,j,2)
IF (dhght > 0) THEN
DO k=1,nz1-1
fac = MAX(0., tem21(i,j,k) - hght0)/aghght
IF (fac > 1) EXIT
tem21(i,j,k) = fac*fac*tem21(i,j,k) &
+ (1.-fac*fac)*(tem21(i,j,k)-dhght)
END DO
ENDIF
END DO
END DO
DO i=1,nx1
DO j=1,ny1-1
k0 = nz-2
DO k=nz1-1,1,-1
z0 = tem21(i,j,k)
DO kk=k0,1,-1
IF(z0 >= temx1y1zb(i,j,kk)) THEN
k0 = kk
EXIT
END IF
END DO
IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
END DO
END DO
END DO
IF( intrpvopt == 1) THEN
DO i=1,nx1
DO j=1,ny1-1
DO k=1,nz1-1
wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/ &
(temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
END DO
END DO
END DO
ELSE
DO i=1,nx1
DO j=1,ny1-1
DO k=1,nz1-1
d=tem21(i,j,k)
a=temx1y1zb(i,j,kz(i,j,k)-1)
b=temx1y1zb(i,j,kz(i,j,k))
c=temx1y1zb(i,j,kz(i,j,k)+1)
wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
END DO
END DO
END IF
END IF
CALL intrpxyz3d(u,nx,1,nx,ny,1,ny-1,nz,1,nz-1, &
wgtux,iux,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
u1,nx1,1,nx1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground
DO i=1,nx1
DO j=1,ny1-1
DO k=nz1-2,1,-1
IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
IF(bglopt == 2 .OR. bglopt == 4 .OR. bglopt == 5) THEN
u1 (i,j,k)=temx1y1z(i,j,2)
ELSE IF(bglopt == 3) THEN
u1 (i,j,k)=misvalue
END IF
END IF
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Intepolate 3D arrays at v-point
!
!-----------------------------------------------------------------------
!
CALL avgz(zp1, 0, nx1,ny1,nz1,1,nx1-1,1,ny1-1,1,nz1-1, tem11 )
CALL avgsv(tem11,nx1,ny1,nz1,1,nx1-1,1,nz1-1, tem21, tem11)
CALL avgz(zp, 0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, tem1 )
CALL avgsv(tem1,nx,ny,nz,1,nx-1,1,nz-1, tem2, tem1)
CALL intrpxy3d(tem2,nx,1,nx-1,ny,1,ny,nz,1,nz-1, &
wgtsx,isx,wgtvy,jvy,intrphopt, &
temx1y1zb,nx1,1,nx1-1, ny1,1,ny1, temx1yz)
DO i=1,nx1-1
DO j=1,ny1
k0=nz-2
DO k=nz1-1,1,-1
z0 = tem21(i,j,k)
DO kk=k0,1,-1
IF(z0 >= temx1y1zb(i,j,kk)) THEN
k0 = kk
EXIT
END IF
END DO
1215 CONTINUE
IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
END DO
END DO
END DO
IF( intrpvopt == 1) THEN
DO i=1,nx1-1
DO j=1,ny1
DO k=1,nz1-1
wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/ &
(temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
END DO
END DO
END DO
ELSE
DO i=1,nx1-1
DO j=1,ny1
DO k=1,nz1-1
d=tem21(i,j,k)
a=temx1y1zb(i,j,kz(i,j,k)-1)
b=temx1y1zb(i,j,kz(i,j,k))
c=temx1y1zb(i,j,kz(i,j,k)+1)
wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
END DO
END DO
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
! IF( nfile == 1) THEN
IF(nfile == 1 .or. noutgrds > 0) THEN
IF( redo_base_state == 0 ) THEN ! get base-state from original arrays
CALL intrpxyz3d(vbar,nx,1,nx-1,ny,1,ny,nz,1,nz-1, &
wgtsx,isx,wgtvy,jvy,wgtz,kz, &
intrphopt,intrpvopt, &
vbar1,nx1,1,nx1-1, ny1,1,ny1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO i=1,nx1-1
DO j=1,ny1
DO k=nz1-2,1,-1
IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
IF(bglopt == 2) THEN
vbar1(i,j,k)=temx1y1z(i,j,2)
ELSE IF(bglopt == 3) THEN
vbar1(i,j,k)=0.0
END IF
END IF
END DO
END DO
END DO
ELSE IF( redo_base_state == 1 ) THEN ! reconstruct base-state from 1d soundings
DO k=1,nz1-1
DO i=1,nx1-1
DO j=1,ny1
kk = MAX(1,MIN(lvlprof-1, &
INT((tem21(i,j,k)-zsnd(2))/snddelz)+2))
IF(tem21(i,j,k) < zsnd(2)) kk = 1
tem = (tem21(i,j,k)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
vbar1(i,j,k)= vbsnd(kk)+(vbsnd(kk+1)-vbsnd(kk))*tem
END DO
END DO
END DO
END IF
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
IF (ntagopt == 1) THEN ! adjust tem21 (heights) so that surface wind is
! same as input surface wind (by matching tem21
! to the input height, temx1y1zb, near the surface).
DO i=1,nx1-1
DO j=1,ny1
hght0 = tem21(i,j,2)
dhght = hght0 - temx1y1zb(i,j,2)
IF (dhght > 0) THEN
DO k=1,nz1-1
fac = MAX(0., tem21(i,j,k) - hght0)/aghght
IF (fac > 1) EXIT
tem21(i,j,k) = fac*fac*tem21(i,j,k) &
+ (1.-fac*fac)*(tem21(i,j,k)-dhght)
END DO
ENDIF
END DO
END DO
DO i=1,nx1-1
DO j=1,ny1
k0=nz-2
DO k=nz1-1,1,-1
z0 = tem21(i,j,k)
DO kk=k0,1,-1
IF(z0 >= temx1y1zb(i,j,kk)) THEN
k0 = kk
EXIT
END IF
END DO
IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
END DO
END DO
END DO
IF( intrpvopt == 1) THEN
DO i=1,nx1-1
DO j=1,ny1
DO k=1,nz1-1
wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/ &
(temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
END DO
END DO
END DO
ELSE
DO i=1,nx1-1
DO j=1,ny1
DO k=1,nz1-1
d=tem21(i,j,k)
a=temx1y1zb(i,j,kz(i,j,k)-1)
b=temx1y1zb(i,j,kz(i,j,k))
c=temx1y1zb(i,j,kz(i,j,k)+1)
wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
END DO
END DO
END IF
END IF
CALL intrpxyz3d(v,nx,1,nx-1,ny,1,ny,nz,1,nz-1, &
wgtsx,isx,wgtvy,jvy,wgtz,kz, &
intrphopt,intrpvopt, &
v1,nx1,1,nx1-1, ny1,1,ny1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground
DO i=1,nx1-1
DO j=1,ny1
DO k=nz1-2,1,-1
IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
IF(bglopt == 2 .OR. bglopt == 4 .OR. bglopt == 5) THEN
v1 (i,j,k)=temx1y1z(i,j,2)
ELSE IF(bglopt == 3) THEN
v1 (i,j,k)=misvalue
END IF
END IF
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Intepolate 3D arrays at w-point
!
!-----------------------------------------------------------------------
!
CALL intrpxy3d(zp,nx,1,nx-1,ny,1,ny-1,nz,1,nz, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
temx1y1zb,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!wdt update: use tem21 instead of zp1 below
tem21 = zp1
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
IF (ntagopt == 1) THEN ! adjust tem21 (heights) so that surface wind is
! same as input surface wind (by matching tem21
! to the input height, temx1y1zb, near the surface).
DO i=1,nx1-1
DO j=1,ny1-1
hght0 = tem21(i,j,2)
dhght = hght0 - temx1y1zb(i,j,2)
IF (dhght > 0) THEN
DO k=1,nz1-1
fac = MAX(0., tem21(i,j,k) - hght0)/aghght
IF (fac > 1) EXIT
tem21(i,j,k) = fac*fac*tem21(i,j,k) &
+ (1.-fac*fac)*(tem21(i,j,k)-dhght)
END DO
ENDIF
END DO
END DO
END IF
DO i=1,nx1-1
DO j=1,ny1-1
k0 = nz-1
DO k=nz1,1,-1
z0 = tem21(i,j,k)
DO kk=k0,1,-1
IF(z0 >= temx1y1zb(i,j,kk)) THEN
k0 = kk
EXIT
END IF
END DO
1315 CONTINUE
IF(intrpvopt == 1) kz(i,j,k)=MAX(2,MIN(k0,nz-1))
IF(intrpvopt == 2) kz(i,j,k)=MAX(3,MIN(k0,nz-1))
END DO
END DO
END DO
IF( intrpvopt == 1) THEN
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,nz1
wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/ &
(temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
END DO
END DO
END DO
ELSE
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,nz1
d=tem21(i,j,k)
a=temx1y1zb(i,j,kz(i,j,k)-1)
b=temx1y1zb(i,j,kz(i,j,k))
c=temx1y1zb(i,j,kz(i,j,k)+1)
wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
END DO
END DO
END IF
CALL intrpxyz3d(w,nx,1,nx-1,ny,1,ny-1,nz,1,nz, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
w1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1, &
temx1yz,temx1y1z)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground
DO i=1,nx1-1
DO j=1,ny1-1
DO k=nz1-1,1,-1
IF(zp1(i,j,k) < temx1y1zb(i,j,2)) THEN ! Blow input grid ground
IF(bglopt == 2 .OR. bglopt == 4 .OR. bglopt == 5) THEN
w1(i,j,k)=temx1y1z(i,j,2)
ELSE IF(bglopt == 3) THEN
w1(i,j,k)=misvalue
END IF
END IF
END DO
END DO
END DO
END IF
!-----------------------------------------------------------------------
!
! Interpolate 3D arrays at scalar point
!
!-----------------------------------------------------------------------
CALL avgz(zp1,0,nx1,ny1,nz1,1,nx1-1,1,ny1-1,1,nz1-1,tem21)
CALL avgz(zp,0,nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem2)
CALL intrpxy3d(tem2,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,intrphopt, &
temx1y1zb,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
DO i=1,nx1-1
DO j=1,ny1-1
htrnx1y1(i,j)=temx1y1zb(i,j,2) ! Old terrain on new grid at scalar level
! e.g., z of first grid level AGL.
DO k=1,nz1-1
IF(tem21(i,j,k) >= htrnx1y1(i,j)) GO TO 1350
END DO
1350 CONTINUE
ktrnx1y1(i,j)= k-1 ! The first scale-level below the first scale level
! above input grid ground
END DO
END DO
DO i=1,nx1-1
DO j=1,ny1-1
k0=nz-2
DO k=nz1-1,1,-1
z0 = tem21(i,j,k)
DO kk=k0,1,-1
IF(z0 >= temx1y1zb(i,j,kk) ) THEN
k0 = kk
EXIT
END IF
END DO
1415 CONTINUE
IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
END DO
END DO
END DO
IF( intrpvopt == 1) THEN
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,nz1-1
wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/ &
(temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
END DO
END DO
END DO
ELSE
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,nz1-1
d=tem21(i,j,k)
a=temx1y1zb(i,j,kz(i,j,k)-1)
b=temx1y1zb(i,j,kz(i,j,k))
c=temx1y1zb(i,j,kz(i,j,k)+1)
wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
END DO
END DO
END DO
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
! IF( nfile == 1) THEN
IF(nfile == 1 .or. noutgrds > 0) THEN
IF( redo_base_state == 0) THEN ! Get base-state from original arrays
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=ALOG(pbar(i,j,k))
END DO
END DO
END DO
CALL intrpxyz3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
pbar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO k=1,nz1-1
DO j=1,ny1-1
DO i=1,nx1-1
pbar1(i,j,k)=EXP(pbar1(i,j,k))
END DO
END DO
END DO
DO i=1,nx1-1
DO j=1,ny1-1
pbsfc(i,j)= EXP(temx1y1z(i,j,2))
END DO
END DO
CALL intrpxyz3d(ptbar,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
ptbar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO i=1,nx1-1
DO j=1,ny1-1
ptbsfc(i,j)= temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(qvbar,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
qvbar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO i=1,nx1-1
DO j=1,ny1-1
qvbsfc(i,j)= temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(rhobar,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
rhobar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
ELSE IF( redo_base_state == 1 ) THEN ! reconstruct base-state from 1d soundings
DO k=1,nz1-1
DO i=1,nx1-1
DO j=1,ny1-1
kk = MAX(1,MIN(lvlprof-1, &
INT((tem21(i,j,k)-zsnd(2))/snddelz)+2))
IF(tem21(i,j,k) < zsnd(2)) kk = 1
tem = (tem21(i,j,k)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
ptbar1(i,j,k)= ptsnd(kk)+(ptsnd(kk+1)-ptsnd(kk))*tem
qvbar1(i,j,k)= qvsnd(kk)+(qvsnd(kk+1)-qvsnd(kk))*tem
tema = ALOG(psnd (kk) )
pbar1 (i,j,k)=EXP(tema+(ALOG(psnd(kk+1))-tema)*tem)
rhobar1(i,j,k)=pbar1(i,j,k)/ &
(rd*ptbar1(i,j,k)*(pbar1(i,j,k)/p0)**rddcp)
END DO
END DO
END DO
DO i=1,nx1-1
DO j=1,ny1-1
kk= MAX(1,MIN(lvlprof-1, &
INT((htrnx1y1(i,j)-zsnd(2))/snddelz)+2))
IF(htrnx1y1(i,j) < zsnd(2)) k = 1
tem = (htrnx1y1(i,j)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
ptbsfc(i,j)= ptsnd(kk)+(ptsnd(kk+1)-ptsnd(kk))*tem
qvbsfc(i,j)= qvsnd(kk)+(qvsnd(kk+1)-qvsnd(kk))*tem
pbsfc (i,j)= EXP(ALOG(psnd (kk))+ &
(ALOG(psnd (kk+1))-ALOG(psnd (kk)))*tem)
END DO
END DO
END IF
END IF
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=ALOG(pbar(i,j,k)+pprt(i,j,k))
END DO
END DO
END DO
CALL intrpxyz3d(tem2,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
pprt1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO k=1,nz1-1
DO j=1,ny1-1
DO i=1,nx1-1
pprt1(i,j,k)=EXP(pprt1(i,j,k))-pbar1(i,j,k)
END DO
END DO
END DO
DO j=1,ny1-1
DO i=1,nx1-1
ppsfc(i,j)=EXP(temx1y1z(i,j,2))-pbsfc(i,j)
END DO
END DO
IF ( redo_base_state == 1 ) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
END DO
END DO
END DO
CALL intrpxyz3d(tem2,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
ptprt1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
ptpsfc(i,j)=temx1y1z(i,j,2)-ptbsfc(i,j)
END DO
END DO
DO k=1,nz1-1
DO j=1,ny1-1
DO i=1,nx1-1
ptprt1(i,j,k)=ptprt1(i,j,k)-ptbar1(i,j,k)
END DO
END DO
END DO
ELSE
CALL intrpxyz3d(ptprt,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
ptprt1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
ptpsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
END IF
CALL intrpxyz3d(qv,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
qv1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
qvsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(qc,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
qc1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
qcsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(qr,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
qr1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
qrsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(qi,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
qi1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
qisfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(qs,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
qs1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
qssfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(qh,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
qh1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
qhsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(tke,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
tke1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
tkesfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(kmh,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
kmh1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
kmhsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(kmv,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
kmv1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
kmvsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
CALL intrpxyz3d(radfrc,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1, &
wgtsx,isx,wgtsy,jsy,wgtz,kz, &
intrphopt,intrpvopt, &
radfrc1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1, &
temx1yz,temx1y1z)
DO j=1,ny1-1
DO i=1,nx1-1
radsfc(i,j)=temx1y1z(i,j,2)
END DO
END DO
!
! Reset extrapolated values below ground for certain options
!
IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground
IF( bglopt == 2 ) THEN ! Constant extension below ground
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,ktrnx1y1(i,j) ! only below ground
ptprt1(i,j,k) = ptpsfc(i,j)
pprt1 (i,j,k) = ppsfc (i,j)
qv1 (i,j,k) = qvsfc (i,j)
qc1 (i,j,k) = qcsfc (i,j)
qr1 (i,j,k) = qrsfc (i,j)
qi1 (i,j,k) = qisfc (i,j)
qs1 (i,j,k) = qssfc (i,j)
qh1 (i,j,k) = qhsfc (i,j)
tke1 (i,j,k) = tkesfc(i,j)
kmh1 (i,j,k) = kmhsfc(i,j)
kmv1 (i,j,k) = kmhsfc(i,j)
qvbar1(i,j,k) = qvbsfc(i,j)
ptbar1(i,j,k) = ptbsfc(i,j)
pbar1 (i,j,k) = pbsfc (i,j)
rhobar1(i,j,k)= pbsfc(i,j)/ &
(rd*ptbsfc(i,j)*(pbsfc(i,j)/p0)**rddcp)
radfrc1(i,j,k)= radsfc(i,j)
END DO
END DO
END DO
ELSE IF( bglopt == 3 ) THEN ! Set to a missing value
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,ktrnx1y1(i,j) ! only below ground
ptprt1(i,j,k) = misvalue
pprt1 (i,j,k) = misvalue
qvbar1(i,j,k) = 0.0
qv1 (i,j,k) = misvalue
qc1 (i,j,k) = misvalue
qr1 (i,j,k) = misvalue
qi1 (i,j,k) = misvalue
qs1 (i,j,k) = misvalue
qh1 (i,j,k) = misvalue
tke1 (i,j,k) = misvalue
kmh1 (i,j,k) = misvalue
kmv1 (i,j,k) = misvalue
ptbar1(i,j,k) = 0.0
pbar1 (i,j,k) = 0.0
rhobar1(i,j,k)= misvalue
radfrc1(i,j,k)= misvalue
END DO
END DO
END DO
!
ELSE IF( bglopt == 4 .OR. bglopt == 5) THEN
! Using a constant lapse rate for temperature
! and then the hydrostatic balance for pressure
DO i=1,nx1-1
DO j=1,ny1-1
DO k=1,ktrnx1y1(i,j) ! only below ground
qv1 (i,j,k) = qvsfc (i,j)
qc1 (i,j,k) = qcsfc (i,j)
qr1 (i,j,k) = qrsfc (i,j)
qi1 (i,j,k) = qisfc (i,j)
qs1 (i,j,k) = qssfc (i,j)
qh1 (i,j,k) = qhsfc (i,j)
tke1 (i,j,k) = tkesfc(i,j)
kmh1 (i,j,k) = kmhsfc(i,j)
kmv1 (i,j,k) = kmhsfc(i,j)
radfrc1(i,j,k)= radsfc(i,j)
! ptprt1(i,j,k) = ptpsfc(i,j)
! pprt1 (i,j,k) = ppsfc (i,j)
qvbar1(i,j,k) = qvbsfc(i,j)
END DO
END DO
END DO
IF( redo_base_state == 1 ) THEN ! There are points below ground
dtdz = -0.0068 ! Minus lapse rate
DO i=1,nx1-1
DO j=1,ny1-1
k=ktrnx1y1(i,j) ! First level below input grid ground
IF( k >= 1) THEN
delz = tem21(i,j,k)-htrnx1y1(i,j)
ttotk1 = (ptbsfc(i,j)+ptpsfc(i,j)) * &
((pbsfc (i,j)+ppsfc (i,j))/p0)**rddcp
ttotk = ttotk1+dtdz*delz
lnptot = ALOG(pbsfc (i,j)+ppsfc (i,j)) &
-g*2.0*delz/(rd*(ttotk1+ttotk))
pprt1(i,j,k) = EXP(lnptot)-pbar1(i,j,k)
ptprt1(i,j,k)= ttotk*(p0/(pbar1(i,j,k)+pprt1(i,j,k))) &
**rddcp-ptbar1(i,j,k)
END IF
DO k=ktrnx1y1(i,j)-1,1,-1 ! other below-ground levels
delz = tem21(i,j,k)-tem21(i,j,k+1)
ttotk1 = (ptbar1(i,j,k+1)+ptprt1(i,j,k+1)) * &
((pbar1 (i,j,k+1)+pprt1 (i,j,k+1))/p0)**rddcp
ttotk = ttotk1+dtdz*delz
lnptot = ALOG(pbar1(i,j,k+1)+pprt1(i,j,k+1)) &
-g*2.0*delz/(rd*(ttotk1+ttotk))
pprt1(i,j,k) = EXP(lnptot)-pbar1(i,j,k)
ptprt1(i,j,k)= ttotk*(p0/(pbar1(i,j,k)+pprt1(i,j,k))) &
**rddcp-ptbar1(i,j,k)
END DO
END DO
END DO
END IF
END IF
END IF ! Check on npoint_below_ground
800 CONTINUE
xgrdorg = xorig1
ygrdorg = yorig1
CALL xytoll(1,1,xgrdorg,ygrdorg,origlat,origlon)
CALL setorig( 2, origlat, origlon) ! set up the model origin on new grid
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
! IF( nfile == 1) THEN
IF(nfile == 1 .or. noutgrds > 0) THEN
!wdt update
DO i=1,nx1
x1_out(i) = x1(i)-xgrdorg
END DO
!wdt update
DO j=1,ny1
y1_out(j) = y1(j)-ygrdorg
END DO
WRITE (*,*) "ctrlat & ctrlon:",ctrlat1,ctrlon1
!wdt update: code block moved here
!wdt begin block
WRITE(6,'(/2x,a)')'For the output grid:'
WRITE(6,'(2x,a,f15.7)') &
'The latitude of the output grid center, ctrlat1=',ctrlat1
WRITE(6,'(2x,a,f15.7/)') &
'The longitude of the output grid center, ctrlon1=',ctrlon1
WRITE(6,'(2x,a,f15.7)') &
'The x coordinate of the output grid center, xctr1=',xctr1
WRITE(6,'(2x,a,f15.7)') &
'The y coordinate of the output grid center, yctr1=',yctr1
IF( ABS(xctr1-xctr1_old) > 0.0001*dx .OR. &
ABS(yctr1-yctr1_old) > 0.0001*dy ) THEN
WRITE(6,'(/1x,a/)') &
'################################################################'
WRITE(6,'(/2x,a/2x,a)') &
'Note that xctr1 and/or yctr1 had been changed by the program', &
'due to snap_to_grid option. Check earlier messages.'
WRITE(6,'(/1x,a/)') &
'################################################################'
END IF
WRITE(6,'(/2x,a/2x,a,2f15.7,a)') &
'The SW corner (i,j)=(2,2) of the grid is located at ', &
'(',xgrdorg,ygrdorg,') of the input grid.'
!wdt begin block
END IF
!-----------------------------------------------------------------------
!
! Copy grid information to the common block variables so that output
! routines have the correct values for the new grid.
!
!-----------------------------------------------------------------------
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
! save global setting before setting them to new grid values
ctrlat_sv = ctrlat
ctrlon_sv = ctrlon
latitud_sv = latitud
dx_sv = dx
dy_sv = dy
dz_sv = dz
dzmin_sv = dzmin
strhopt_sv = strhopt
zrefsfc_sv = zrefsfc
dlayer1_sv = dlayer1
dlayer2_sv = dlayer2
zflat_sv = zflat
strhtune_sv = strhtune
nstyp_sv = nstyp
ctrlat = ctrlat1
ctrlon = ctrlon1
latitud = ctrlat
dx = dx1
dy = dy1
dz = dz1
dzmin = dzmin1
strhopt = strhopt1
zrefsfc = zrefsfc1
dlayer1 = dlayer11
dlayer2 = dlayer21
zflat = zflat1
strhtune = strhtune1
nstyp = nstyp1
WRITE (cmnt(nocmnt),'(a,i4,a,i4,a,i4)') &
' nx =',nx1,', ny =',ny1,', nz =',nz1
!
!-----------------------------------------------------------------------
!
! Print out the max/min of output varaibles.
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/1x,a/)') &
'Min. and max. of input data interpolated to the new grid:'
!wdt update
CALL a3dmax0(x1_out,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
!wdt update
CALL a3dmax0(y1_out,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
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
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
CALL a3dmax0(ubar1,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
CALL a3dmax0(vbar1,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
CALL a3dmax0(ptbar1,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
CALL a3dmax0(pbar1,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
CALL a3dmax0(rhobar1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'rhobarmin=', amin,', rhobarmax=',amax
CALL a3dmax0(qvbar1,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
CALL a3dmax0(u1,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
CALL a3dmax0(v1,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
CALL a3dmax0(w1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'wmin = ', amin,', wmax =',amax
CALL a3dmax0(ptprt1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,', ptprtmax=',amax
CALL a3dmax0(pprt1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,', pprtmax =',amax
CALL a3dmax0(qv1,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
CALL a3dmax0(qc1,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
CALL a3dmax0(qr1,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
CALL a3dmax0(qi1,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
CALL a3dmax0(qs1,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
CALL a3dmax0(qh1,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
CALL a3dmax0(tke1,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
CALL a3dmax0(kmh1,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
CALL a3dmax0(kmv1,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
CALL a3dmax0(raing1,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
CALL a3dmax0(rainc1,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
CALL a3dmax0(prcrate1(1,1,1),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
CALL a3dmax0(prcrate1(1,1,2),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,', prcr2max=',amax
CALL a3dmax0(prcrate1(1,1,3),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,', prcr3max=',amax
CALL a3dmax0(prcrate1(1,1,4),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,', prcr4max=',amax
DO iss = 0, nstyp1
CALL a3dmax0(tsfc1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'tsfcmin = ', amin,', tsfcmax =',amax,' for soil type=',iss
CALL a3dmax0(tsoil1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'tsoilmin= ', amin,', tsoilmax=',amax,' for soil type=',iss
CALL a3dmax0(wetsfc1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'wetsmin = ', amin,', wetsmax =',amax,' for soil type=',iss
CALL a3dmax0(wetdp1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'wetdmin = ', amin,', wetdmax =',amax,' for soil type=',iss
CALL a3dmax0(wetcanp1(1,1,iss), &
1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6),a,i3)') &
'wetcmin = ', amin,', wetcmax =',amax,' for soil type=',iss
END DO
CALL a3dmax0(roufns1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'roufnmin= ', amin,', roufnmax =',amax
CALL a3dmax0(veg1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vegmin = ', amin,', vegmax =',amax
CALL a3dmax0(radfrc1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'radfnmin= ', amin,', radfnmax =',amax
CALL a3dmax0(radsw1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'radswmin= ', amin,', radswmax =',amax
CALL a3dmax0(rnflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'rnflxmin= ', amin,', rnflxmax =',amax
CALL a3dmax0(usflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'usflxmin= ', amin,', usflxmax =',amax
CALL a3dmax0(vsflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'vsflxmin= ', amin,', vsflxmax =',amax
CALL a3dmax0(ptsflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'ptflxmin= ', amin,', ptflxmax =',amax
CALL a3dmax0(qvsflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') 'qvflxmin= ', amin,', qvflxmax =',amax
tem11= 0.0 ! To be put in place of wbar
!
!-----------------------------------------------------------------------
!
! Data dump of the model grid and base state arrays:
!
! First find a unique name basdmpfn(1:lbasdmpf) for the grid and
! base state array dump file
!
! If grid/base state data has been written out once, skip
! the following writing block. Also no need to write out
! separate data for Savi3D dump. The same for GrADS dump.
!
!-----------------------------------------------------------------------
!
runname = new_runname
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
IF (noutgrds > 0) THEN
runname = name_grd(ng)
ENDIF
CALL gtlfnkey(runname, lfnkey)
IF(houtfmt /= 9 ) THEN
!wdt update
IF( nfile > 1 ) GO TO 300 ! If done already, skip this part.
CALL gtbasfn(runname(1:lfnkey),dirname,ldirnam,hdmpfmt, &
1,0,basdmpfn,lbasdmpf)
PRINT*
PRINT*,'Output grid/base state file is ', basdmpfn(1:lbasdmpf)
grdbas = 1 ! Dump out grd and base state arrays only
!wdt update
CALL setcornerll( nx1,ny1,x1_out,y1_out )
!wdt update
CALL dtadump(nx1,ny1,nz1,nstyp1,hdmpfmt,nchout, &
basdmpfn(1:lbasdmpf),grdbas,filcmprs, &
u1,v1,w1,ptprt1,pprt1,qv1,qc1,qr1,qi1,qs1,qh1, &
tke1,kmh1,kmv1, &
ubar1,vbar1,tem11,ptbar1,pbar1,rhobar1,qvbar1, &
x1_out,y1_out,z1,zp1,hterain1,j11,j21,j31, &
soiltyp1,stypfrct1,vegtyp1,lai1,roufns1,veg1, &
tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,snowdpth1, &
raing1,rainc1,prcrate1, &
radfrc1,radsw1,rnflx1, &
usflx1,vsflx1,ptsflx1,qvsflx1, &
tem21,wgtz,kz)
!wdt update - removed gboutcnt
300 CONTINUE
END IF
!
!-----------------------------------------------------------------------
!
! Then the time dependent fields:
!
!-----------------------------------------------------------------------
!
IF( .NOT. (houtfmt == 9 .AND. vroutcnt == 1) ) THEN
!
!-----------------------------------------------------------------------
!
! Reconstruct the file name using the specified directory name
!
!-----------------------------------------------------------------------
!
CALL gtdmpfn(runname(1:lfnkey),dirname, &
ldirnam,curtim,hdmpfmt,1,0, hdmpfn, ldmpf)
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
IF ( exbcdmp == 4 ) rayklow = -1
WRITE(6,'(a,a)') 'Writing t-dependent variable history dump ', &
hdmpfn(1:ldmpf)
grdbas = 0
!wdt update
CALL dtadump(nx1,ny1,nz1,nstyp1,hdmpfmt,nchout, &
hdmpfn(1:ldmpf),grdbas,filcmprs, &
u1,v1,w1,ptprt1,pprt1,qv1,qc1,qr1,qi1,qs1,qh1, &
tke1,kmh1,kmv1, &
ubar1,vbar1,tem11,ptbar1,pbar1,rhobar1,qvbar1, &
x1_out,y1_out,z1,zp1,hterain1,j11,j21,j31, &
soiltyp1,stypfrct1,vegtyp1,lai1,roufns1,veg1, &
tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,snowdpth1, &
raing1,rainc1,prcrate1, &
radfrc1,radsw1,rnflx1, &
usflx1,vsflx1,ptsflx1,qvsflx1, &
tem21,wgtz,kz)
!
!-----------------------------------------------------------------------
!
! Write out soil model variable file
!
!-----------------------------------------------------------------------
!
IF (sfcin == 1 .AND. soildmp > 0) THEN
CALL cvttsnd( curtim, timsnd, tmstrln )
soiloutfl = runname(1:lfnkey)//".soilvar."//timsnd(1:tmstrln)
lfn = lfnkey + 9 + tmstrln
IF( dirname /= ' ' ) THEN
temchar = soiloutfl
soiloutfl = dirname(1:ldirnam)//'/'//temchar
lfn = lfn + ldirnam + 1
END IF
CALL fnversn(soiloutfl, lfn)
PRINT *, 'Write soil initial data to ',soiloutfl(1:lfn)
CALL wrtsoil(nx1,ny1,nstyp1, soiloutfl(1:lfn), dx1,dy1, &
mapproj,trulat1,trulat2,trulon,sclfct, &
ctrlat1,ctrlon1, &
1,1,1,1,1,1, &
tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,snowdpth1,soiltyp1)
!wdt update
IF (soildmp == 1) CALL soilcntl(nx1,ny1, soiloutfl(1:lfn), &
1,1,1,1,1,1, x1_out,y1_out)
END IF ! sfcin.eq.1
IF(nfile == 1) THEN ! Need to do this only once
!
!-----------------------------------------------------------------------
!
! Write out terrain data
!
!-----------------------------------------------------------------------
!
IF (terndmp > 0) THEN
CALL getunit( nunit )
ternfn = 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 writtrn(nx1,ny1,ternfn(1:lternfn), dx1,dy1, &
mapproj,trulat1,trulat2,trulon,sclfct, &
ctrlat1,ctrlon1,hterain1)
!wdt update
IF (terndmp == 1) &
CALL trncntl(nx1,ny1,ternfn(1:lternfn), x1_out,y1_out)
END IF
!-----------------------------------------------------------------------
!
! Write out surface property data file: sfcoutfl .
!
!-----------------------------------------------------------------------
!
IF (landin == 1 .AND. sfcdmp > 0) THEN
sfcoutfl = runname(1:lfnkey)//".sfcdata"
lfn = lfnkey + 8
IF( dirname /= ' ' ) THEN
temchar = sfcoutfl
sfcoutfl = dirname(1:ldirnam)//'/'//temchar
lfn = lfn + ldirnam + 1
END IF
CALL fnversn(sfcoutfl, lfn)
PRINT *, 'Write surface property data in ',sfcoutfl(1:lfn)
CALL wrtsfcdt(nx1,ny1,nstyp1,sfcoutfl(1:lfn), dx1,dy1, &
mapproj,trulat1,trulat2,trulon,sclfct, &
ctrlat1,ctrlon1, &
1,1,1,1,1,0, &
soiltyp1,stypfrct1,vegtyp1,lai1,roufns1,veg1,veg1)
!wdt udpate
IF (sfcdmp == 1) CALL sfccntl(nx1,ny1, sfcoutfl(1:lfn), &
1,1,1,1,1,0, x1_out,y1_out)
END IF ! landin.eq.1
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
ctrlat = ctrlat_sv
ctrlon = ctrlon_sv
latitud = latitud_sv
dx = dx_sv
dy = dy_sv
dz = dz_sv
dzmin = dzmin_sv
strhopt = strhopt_sv
zrefsfc = zrefsfc_sv
dlayer1 = dlayer1_sv
dlayer2 = dlayer2_sv
zflat = zflat_sv
strhtune = strhtune_sv
nstyp = nstyp_sv
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
600 CONTINUE !wdt
IF (noutgrds > 0 .and. ng < noutgrds) GOTO 500
vroutcnt = 1 ! Variables have been written out at least once
IF (hinfmt == 9) GO TO 102
END DO
9001 CONTINUE
!wdt update: moved code block to above
WRITE (6,'(/a,g13.3,a)') "Current memory allocation=", &
current_memory_use,' (Mw).'
WRITE (6,'(a,g13.3,a/)') "Maxumum memory allocation=", &
max_memory_use,' (Mw).'
cpu_time = f_cputime() - cpu0
PRINT*,'Total CPU time used =', cpu_time,' (s).'
STOP
9002 CONTINUE
WRITE(6,'(1x,a,i2,/1x,a)') &
'Data read was unsuccessful. ireturn =', ireturn, &
'Job stopped in ARPSINTRP.'
STOP 9002
END PROGRAM arpsintrp
SUBROUTINE intrpx3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgtx,ix, & 1
intrphopt,aout,nx1,is1,ie1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in the first dimension
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: nx1,is1,ie1
REAL :: ain (nx ,ny,nz)
REAL :: aout(nx1,ny,nz)
REAL :: wgtx(nx1,3)
INTEGER :: ix(nx1)
INTEGER :: intrphopt
INTEGER :: i,j,k
IF(intrphopt == 1) THEN
DO k=ks ,ke
DO j=js ,je
DO i=is1,ie1
aout(i,j,k)= wgtx(i,1) *ain(ix(i) ,j,k) &
+(1.0-wgtx(i,1))*ain(ix(i)+1,j,k)
END DO
END DO
END DO
ELSE
DO k=ks ,ke
DO j=js ,je
DO i=is1,ie1
aout(i,j,k)=wgtx(i,1)*ain(ix(i)-1,j,k) &
+wgtx(i,2)*ain(ix(i) ,j,k) &
+wgtx(i,3)*ain(ix(i)+1,j,k)
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE intrpx3d
SUBROUTINE intrpy3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgty,jy, & 1
intrphopt,aout,ny1,js1,je1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in the second dimension
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: ny1,js1,je1
REAL :: ain (nx,ny ,nz)
REAL :: aout(nx,ny1,nz)
REAL :: wgty(ny1,3)
INTEGER :: jy(ny1)
INTEGER :: intrphopt
INTEGER :: i,j,k
IF(intrphopt == 1) THEN
DO k=ks ,ke
DO j=js1,je1
DO i=is ,ie
aout(i,j,k)= wgty(j,1) *ain(i,jy(j) ,k) &
+(1.0-wgty(j,1))*ain(i,jy(j)+1,k)
END DO
END DO
END DO
ELSE
DO k=ks ,ke
DO j=js1,je1
DO i=is ,ie
aout(i,j,k)=wgty(j,1)*ain(i,jy(j)-1,k) &
+wgty(j,2)*ain(i,jy(j) ,k) &
+wgty(j,3)*ain(i,jy(j)+1,k)
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE intrpy3d
SUBROUTINE intrpz3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgtz,kz, & 1
intrpvopt,aout,nz1,ks1,ke1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in the third dimension
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: nz1,ks1,ke1
REAL :: ain (nx,ny,nz)
REAL :: aout(nx,ny,nz1)
REAL :: wgtz(nx,ny,nz1,3)
INTEGER :: kz(nx,ny,nz1)
INTEGER :: intrpvopt
INTEGER :: i,j,k
IF(intrpvopt == 1) THEN
DO k=ks1,ke1
DO i=is ,ie
DO j=js ,je
aout(i,j,k)= wgtz(i,j,k,1) *ain(i,j,kz(i,j,k) ) &
+(1.0-wgtz(i,j,k,1))*ain(i,j,kz(i,j,k)+1)
END DO
END DO
END DO
ELSE
DO k=ks1,ke1
DO i=is ,ie
DO j=js ,je
aout(i,j,k)=wgtz(i,j,k,1)*ain(i,j,kz(i,j,k)-1) &
+wgtz(i,j,k,2)*ain(i,j,kz(i,j,k) ) &
+wgtz(i,j,k,3)*ain(i,j,kz(i,j,k)+1)
END DO
END DO
END DO
END IF
RETURN
END SUBROUTINE intrpz3d
SUBROUTINE intrpxy3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, & 20,2
wgtx,ix,wgty,jy,intrphopt, &
aout,nx1,is1,ie1, ny1,js1,je1, &
temx1yz)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in the first and second dimensions
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: nx1,is1,ie1, ny1,js1,je1
REAL :: ain (nx ,ny ,nz)
REAL :: aout(nx1,ny1,nz)
REAL :: wgtx(nx1,3),wgty(ny1,3)
INTEGER :: ix(nx1),jy(ny1)
INTEGER :: intrphopt
REAL :: temx1yz(nx1,ny,nz)
CALL intrpx3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, &
wgtx,ix,intrphopt, temx1yz,nx1,is1,ie1)
CALL intrpy3d(temx1yz,nx1,is1,ie1, ny,js,je, nz,ks,ke, &
wgty,jy,intrphopt, aout, ny1,js1,je1)
RETURN
END SUBROUTINE intrpxy3d
SUBROUTINE intrpxyz3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, & 22,2
wgtx,ix,wgty,jy,wgtz,kz, &
intrphopt,intrpvopt, &
aout,nx1,is1,ie1, ny1,js1,je1, nz1,ks1,ke1, &
temx1yz,temx1y1z)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Perform interpolation in all three dimensions
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue
! 4/1/1999.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
INTEGER :: nx1,is1,ie1, ny1,js1,je1, nz1,ks1,ke1
REAL :: ain (nx ,ny ,nz)
REAL :: aout(nx1,ny1,nz1)
REAL :: wgtx(nx1,3),wgty(ny1,3),wgtz(nx1,ny1,nz1,3)
INTEGER :: ix(nx1),jy(ny1),kz(nx1,ny1,nz1)
INTEGER :: intrphopt
INTEGER :: intrpvopt
REAL :: temx1yz (nx1,ny ,nz)
REAL :: temx1y1z(nx1,ny1,nz)
CALL intrpxy3d(ain, nx,is,ie, ny,js,je, nz,ks,ke, &
wgtx,ix,wgty,jy,intrphopt, &
temx1y1z, nx1,is1,ie1, ny1,js1,je1, temx1yz)
CALL intrpz3d (temx1y1z, nx1,is1,ie1, ny1,js1,je1, nz,ks,ke, &
wgtz,kz,intrpvopt, &
aout, nz1,ks1,ke1)
RETURN
END SUBROUTINE intrpxyz3d