!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM EXT2ARPS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
PROGRAM ext2arps,172
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Converts external forecast files coordinates and variables
! to those used in ARPS and writes the information in a
! history dump file.
!
! The actual reading of the external files is done by
! subroutine rdextfil. In many cases, it will only be necessary
! to make changes in extdims.inc and rdextfil to customize the
! processing for the user's specific data source.
!
! This program could become more general in the future, but
! for now we assume incoming data are supplied (or converted to)
! pressure(Pa), height(m), temperature(K), specific humidty (kg/kg),
! and u,v (m/s). Units should be converted from their original
! forms to these in rdextfil.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! Sept, 1994 after MAPS2ARPS
!
! MODIFICATION HISTORY:
! Oct, 1994 (KB)
! Changed input times to be more flexible.
!
! 17 Jan, 1995 (KB)
! Modified processing so that mean sounding would be for constant
! pressure surfaces rather than constant ARPS levels.
!
! 23 Jan, 1995 (KB)
! Added temporary array avgzs_ext, replacing the use of tem1 which
! would not work if number of ext levels is greater than ARPS.
!
! 9 August, 1995 (KB)
! Modified creation of mean soundings and storage of external
! pressure arrays to allow for incoming data on other than
! pressure vertical coordinates. Added switch for buoyancy
! balancing of vertical pressure gradient.
!
! 26 April, 1996 (KB)
! Version 2.0, including: Corrects divergence error
! using a function either linear in k-level, physical height or
! potential temperature. Also, moisture is interpolated and
! extrapolated using RHstar [ RHstar=sqrt(1. - RH) ], which should
! be more linear in height, rather than qv. Extrapolation above
! and below available model data are constant-in-height perturbation
! fields for all variables except assigned a standard lapse rate near
! the ground and dT/dz=0 at the top.
!
! 30 April 1996 (KB)
! Final preparation of version 2.0, added O'Brien option and
! order of interpolation option to input namelist.
!
! 8 November 1996 (KB)
! Version 2.1
! Changed the polynomial interpolation scheme for better
! accuracy and eliminated chance for a discontinuity at external
! grid cell boundaries. Linear or quadratic interpolation
! using local polynomial fit are supported (iorder=1 or 2).
!
! 3 March 1996 (KB)
! Pressure is now interpolated with polynomials of pressure
! rather than ln(p).
!
! 10 June 1997 (Fanyou Kong -- CMRP)
! Include surface 2D data analysis from external data sets, a new
! subroutine MKARPS2D is then added in file 'extlib23.f'.
!
! 16 July 1997 (Zonghui Huo)
! Added 'getcoamps.f' to read the COAMPS data on sigma-z levels.
!
! 06/16/1998 (Dennis Moon, SSESCO)
! Added support for RUC 2 on the AWIPS 211 grid. extdopt=7
!
! 06/16/1998 (Dennis Moon, SSESCO)
! Added support for NCEP global re-analysis on T62 Gaussian lat/lon
! grid. This does not work for ARPS grids which cross the 0 degree
! meridian. extdopt=8
!
! 29 Oct 1998 (R. Carpenter, G. Bassett)
! Added RUC-2 GEMPAK & ETA GEMPAK grid #104. Four digit years and
! extdname now used in constructing filenames of GEMPAK files.
!
! 15 Sep 1999 (K. Brewster)
! Added the microphysical variables to the dtadump calls so they
! will be used in those cases when available. Added check for
! positive definiteness to microphysical variables.
!
! 2000/03/23 (Donghai Wang)
! Added NCEP AVN GRIB global data, grid #3.
!
! 2000/07/28 (Ming Xue)
! Converted to F90 free format and implemented dynamic memory
! allocation.
!
! 2000/08/14 (Gene Bassett)
! Added multiple soil types, sfcdata variables and grid staggering
! for use with arps history format of external model data. Added options
! allowing ext2arps to produce results similar to arpsintrp (intropt,
! ext_lbc, ext_vbc).
!
! 2001/05/31 (Gene Bassett)
! Added exttrnopt, an option to use input grid terrain for output grid.
! Also corrected a bug introduced when updating the use of arps as an
! external model (2000/08) which caused only the first external file
! to have u&v rotated to arps map projection.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Input and output grid dimensions
!
!-----------------------------------------------------------------------
INTEGER :: nx ! Number of grid points in the x-dir of output arps grid
INTEGER :: ny ! Number of grid points in the y-dir of output arps grid
INTEGER :: nz ! Number of grid points in the z-dir of output arps grid
!
INTEGER :: nx_ext, ny_ext, nz_ext ! dimensions of external data grid
INTEGER :: nstyps
!
!-----------------------------------------------------------------------
!
! Space for mean sounding
!
!-----------------------------------------------------------------------
!
INTEGER :: lvlprof
REAL :: depthp
PARAMETER( lvlprof=601, depthp = 3.0E4 )
!
!-----------------------------------------------------------------------
!
! lvlprof:
!
! The ARPS interpolates unevenly-spaced data from the external
! data set to evenly-spaced data for use in defining the base
! state atmosphere. In this process, an intermediate sounding is
! generated at evenly-spaced altitudes, with the accuracy of the
! associated interpolation controlled by the parameter lvlprof.
! The larger lvlprof, the more accurate the interpolation
! (we recommend using lvlprof=200 for a model run with about 50
! points in the vertical direction). Using the intermediate
! sounding, the ARPS then generates a base state model sounding
! for the particular vertical grid resolution chosen (i.e.,
! the number of points in the vertical, nz, and the vertical grid
! spacing, dz).
!
! depthp:
!
! The depth of atmosphere over which the interpolated profiles
! will be defined. depthp should be greater than or equal to
! (nz-3)*dz, i.e., larger than the physical depth of the model
! domain. Otherwise, the code will extrapolate for gridpoints
! outside the domain, leading to possible inconsistencies. At all
! costs, any such extrapolation should be avoided.
!
!-----------------------------------------------------------------------
!
REAL :: psnd(lvlprof), &
zsnd(lvlprof), &
tsnd(lvlprof), &
ptsnd(lvlprof), &
rhssnd(lvlprof), &
qvsnd(lvlprof), &
rhosnd(lvlprof), &
usnd(lvlprof), &
vsnd(lvlprof), &
dumsnd(lvlprof), &
plsnd(lvlprof)
!
!-----------------------------------------------------------------------
!
! ARPS grid variables
!
!
! u x component of velocity at a given time level (m/s)
! v y component of velocity at a given time level (m/s)
! w Vertical component of Cartesian velocity at a given
! time level (m/s)
! ptprt Perturbation potential temperature at a given time
! level (K)
! pprt Perturbation pressure at a given time level (Pascal)
! qv Water vapor specific humidity at a given time level (kg/kg)
! qc Cloud water mixing ratio at a given time level (kg/kg)
! qr Rainwater mixing ratio at a given time level (kg/kg)
! qi Cloud ice mixing ratio at a given time level (kg/kg)
! qs Snow mixing ratio at a given time level (kg/kg)
! qh Hail mixing ratio at a given time level (kg/kg)
!
! ubar Base state x-velocity component (m/s)
! vbar Base state y-velocity component (m/s)
! wbar Base state vertical velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhobar Base state density (kg/m**3)
! qvbar Base state water vapor specific humidity (kg/kg)
!
! x x coordinate of grid points in physical/comp. space (m)
! y y coordinate of grid points in physical/comp. space (m)
! z z coordinate of grid points in computational space (m)
! zp Vertical coordinate of grid points in physical space (m)
! hterain Terrain height (m)
!
! j1 Coordinate transformation Jacobian -d(zp)/dx
! j2 Coordinate transformation Jacobian -d(zp)/dy
! j3 Coordinate transformation Jacobian d(zp)/dz
!
! tsfc Temperature at ground (K) (in top 1 cm layer)
! tsoil Deep soil temperature (K) (in deep 1 m layer)
! wetsfc Surface soil moisture in the top 1 cm layer
! wetdp Deep soil moisture in the deep 1 m layer
! wetcanp Canopy water amount
!
!-----------------------------------------------------------------------
!
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 :: xscl(:) ! The x-coordinate of scalar points.
REAL, allocatable :: yscl(:) ! The y-coordinate of scalar points.
REAL, allocatable :: zp(:,:,:) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, allocatable :: zs(:,:,:) ! The physical height of scalar points.
REAL, allocatable :: j1(:,:,:) ! Coordinate transformation Jacobian -d(zp)/dx.
REAL, allocatable :: j2(:,:,:) ! Coordinate transformation Jacobian -d(zp)/dy.
REAL, allocatable :: j3(:,:,:) ! Coordinate transformation Jacobian d(zp)/dz.
REAL, allocatable :: aj3z(:,:,:)! Coordinate transformation Jacobian defined
! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
REAL, allocatable :: hterain(:,:) ! The height of the terrain. (m)
!wdt update
REAL, allocatable :: htrn_ext(:,:) ! The height of the terrain. (m)
REAL, allocatable :: mapfct(:,:,:)! Map factor at scalar, u, and v points
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 :: pprt(:,:,:) ! Perturbation pressure from that
! of base state atmosphere (Pascal)
REAL, allocatable :: ptprt(:,:,:) ! Perturbation potential temperature
! from that of base state atmosphere (K)
!
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 :: pbar(:,:,:) ! Base state pressure (Pascal)
REAL, allocatable :: ptbar(:,:,:) ! Base state potential temperature (K)
REAL, allocatable :: qvbar(:,:,:) ! Base state water vapor specific humidity
! (kg/kg)
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 :: rhobar(:,:,:)! Base state density (kg/m3).
REAL, allocatable :: rhostr(:,:,:)! Base state density rhobar times j3.
REAL, allocatable :: wcont(:,:,:) ! Contravariant vertical velocity (m/s)
REAL, allocatable :: csndsq(:,:,:)! Speed of sound squared (m**2/s**2)
REAL, allocatable :: tsfc(:,:,:) ! Ground surface temperature (K)
REAL, allocatable :: tsoil(:,:,:) ! Deep soil temperature (K)
REAL, allocatable :: wetsfc(:,:,:) ! Surface soil moisture
REAL, allocatable :: wetdp(:,:,:) ! Deep soil moisture
REAL, allocatable :: wetcanp(:,:,:)! Water amount on canopy
REAL, allocatable :: snowdpth(:,:) ! Snow depth (m)
INTEGER, allocatable :: soiltyp (:,:,:)! Soil type
REAL, allocatable :: stypfrct(:,:,:)! Soil type fraction
INTEGER, allocatable :: vegtyp (:,:) ! Vegetation type
REAL, allocatable :: lai (:,:) ! Leaf Area Index
REAL, allocatable :: roufns (:,:) ! Surface roughness
REAL, allocatable :: veg (:,:) ! Vegetation fraction
REAL, allocatable :: tem1(:,:,:) ! Temporary work array.
REAL, allocatable :: tem2(:,:,:) ! Temporary work array.
REAL, allocatable :: tem3(:,:,:) ! Temporary work array.
REAL, allocatable :: tem4(:,:,:) ! Temporary work array.
REAL, allocatable :: tem5(:,:,:) ! Temporary work array.
REAL, allocatable :: tem6(:,:,:) ! Temporary work array.
REAL, allocatable :: tem7(:,:,:) ! Temporary work array.
REAL, allocatable :: tem8(:,:,:) ! Temporary work array.
REAL, allocatable :: xs2d(:,:)
REAL, allocatable :: ys2d(:,:)
REAL, allocatable :: xu2d(:,:)
REAL, allocatable :: yu2d(:,:)
REAL, allocatable :: xv2d(:,:)
REAL, allocatable :: yv2d(:,:)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: xw & yw
REAL, allocatable :: xw(:,:)
REAL, allocatable :: yw(:,:)
INTEGER, allocatable :: iscl(:,:) ! i index of scalar point
! in external array.
INTEGER, allocatable :: jscl(:,:) ! j index of scalar point
! in external array.
INTEGER, allocatable :: iu(:,:) ! i index of u-velocity point
! in external array.
INTEGER, allocatable :: ju(:,:) ! j index of u-velocity point
! in external array.
INTEGER, allocatable :: iv(:,:) ! i index of v-velocity point
! in external array.
INTEGER, allocatable :: jv(:,:) ! j index of v-velocity point
! in external array.
INTEGER :: iproj_arps ! ARPS map projection (tem storage)
REAL :: scale_arps ! ARPS map scale (tem storage)
REAL :: trlon_arps ! ARPS true-longitude (tem storage)
REAL :: latnot_arps(2) ! ARPS true-latitude(s) (tem storage)
REAL :: xorig_arps,yorig_arps ! ARPS map projection origin (tem storage)
!
!-----------------------------------------------------------------------
!
! NAMELIST parameter (In arps.input)
!
!-----------------------------------------------------------------------
!
INCLUDE 'ext2arps.inc'
!
!-----------------------------------------------------------------------
!
! External forecast variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext ! external data map projection
REAL :: scale_ext ! external data map scale factor
REAL :: trlon_ext ! external data true longitude
REAL :: latnot_ext(2) ! external data true latitude(s)
REAL :: x0_ext,y0_ext ! external data origin
REAL, allocatable :: x_ext(:) ! external data x-coordinate
REAL, allocatable :: y_ext(:) ! external data y-coordinate
REAL, allocatable :: xu_ext(:) ! external data u x-coordinate
REAL, allocatable :: yu_ext(:) ! external data u y-coordinate
REAL, allocatable :: xv_ext(:) ! external data v x-coordinate
REAL, allocatable :: yv_ext(:) ! external data v y-coordinate
REAL, allocatable :: lat_ext(:,:) ! external data latidude
REAL, allocatable :: lon_ext(:,:) ! external data longitude
REAL, allocatable :: latu_ext(:,:) ! external data latidude (x-stag)
REAL, allocatable :: lonu_ext(:,:) ! external data longitude (x-stag)
REAL, allocatable :: latv_ext(:,:) ! external data latidude (y-stag)
REAL, allocatable :: lonv_ext(:,:) ! external data longitude (y-stag)
REAL, allocatable :: zs_ext(:,:,:) ! external data physical height (m)
REAL, allocatable :: p_ext(:,:,:) ! Pressure (Pascals)
REAL, allocatable :: hgt_ext(:,:,:) ! Height (m)
REAL, allocatable :: zp2_ext(:,:,:) ! external w physical height (m)
REAL, allocatable :: zp_ext(:,:,:) ! Height (m) (on arps grid)
REAL, allocatable :: t_ext(:,:,:) ! Temperature (K)
REAL, allocatable :: u_ext(:,:,:) ! Eastward wind component
REAL, allocatable :: v_ext(:,:,:) ! Northward wind component
REAL, allocatable :: w_ext(:,:,:) ! Vertical wind component
REAL, allocatable :: qv_ext(:,:,:) ! Specific humidity (kg/kg)
REAL, allocatable :: qc_ext(:,:,:) ! Cloud H2O mixing ratio (kg/kg)
REAL, allocatable :: qr_ext(:,:,:) ! Rain H2O mixing ratio (kg/kg)
REAL, allocatable :: qi_ext(:,:,:) ! Ice H2O mixing ratio (kg/kg)
REAL, allocatable :: qs_ext(:,:,:) ! Snow H2O mixing ratio (kg/kg)
REAL, allocatable :: qh_ext(:,:,:) ! Hail H2O mixing ratio (kg/kg)
REAL, allocatable :: tsfc_ext (:,:,:) ! Temperature at surface (K)
REAL, allocatable :: tsoil_ext (:,:,:) ! Deep soil temperature (K)
REAL, allocatable :: wetsfc_ext (:,:,:) ! Surface soil moisture
REAL, allocatable :: wetdp_ext (:,:,:) ! Deep soil moisture
REAL, allocatable :: wetcanp_ext(:,:,:) ! Canopy water amount
REAL, allocatable :: snowdpth_ext(:,:) ! Snow depth (m)
INTEGER, allocatable :: soiltyp_ext (:,:,:)! Soil type
REAL, allocatable :: stypfrct_ext(:,:,:)! Soil type fraction
INTEGER, allocatable :: vegtyp_ext (:,:) ! Vegetation type
REAL, allocatable :: lai_ext (:,:) ! Leaf Area Index
REAL, allocatable :: roufns_ext (:,:) ! Surface roughness
REAL, allocatable :: veg_ext (:,:) ! Vegetation fraction
REAL, allocatable :: trn_ext (:,:) ! External terrain (m)
REAL, allocatable :: psfc_ext (:,:) ! Surface pressure (Pa)
REAL, allocatable :: t_2m_ext (:,:)
REAL, allocatable :: rh_2m_ext(:,:)
REAL, allocatable :: u_10m_ext(:,:)
REAL, allocatable :: v_10m_ext(:,:)
REAL, allocatable :: dxfld(:)
REAL, allocatable :: dyfld(:)
REAL, allocatable :: rdxfld(:)
REAL, allocatable :: rdyfld(:)
REAL, allocatable :: dxfldu(:)
REAL, allocatable :: dyfldu(:)
REAL, allocatable :: rdxfldu(:)
REAL, allocatable :: rdyfldu(:)
REAL, allocatable :: dxfldv(:)
REAL, allocatable :: dyfldv(:)
REAL, allocatable :: rdxfldv(:)
REAL, allocatable :: rdyfldv(:)
REAL, allocatable :: tem1_ext(:,:,:) ! Temporary work array
REAL, allocatable :: tem2_ext(:,:,:) ! Temporary work array
REAL, allocatable :: tem3_ext(:,:,:) ! Temporary work array
REAL, allocatable :: tem4_ext(:,:,:) ! Temporary work array
REAL, allocatable :: tem5_ext(:,:,:) ! Temporary work array
REAL, allocatable :: xa_ext(:,:)
REAL, allocatable :: ya_ext(:,:)
REAL, allocatable :: avgzs_ext(:,:,:)
!
!-----------------------------------------------------------------------
!
! include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.inc'
INCLUDE 'bndry.inc'
INCLUDE 'adjust.inc'
!
!-----------------------------------------------------------------------
!
! Non-dimensional smoothing coefficient
!
!-----------------------------------------------------------------------
!
REAL :: smfct1,smfct2
PARAMETER( smfct1=0.5, smfct2=-0.5 )
!
REAL :: rhmax
PARAMETER ( rhmax = 1.0)
!
!-----------------------------------------------------------------------
!
! Latitude and longitude for some diagnostic printing,
! e.g. to compare to an observed sounding
!
!-----------------------------------------------------------------------
!
REAL :: latdiag,londiag
PARAMETER (latdiag=34.5606,londiag=-103.0820)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: basdmpfn
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfinit
CHARACTER (LEN=9) :: julfname
!
INTEGER :: i,j,k,ksmth,istatus
INTEGER :: iyr,imo,iday,ihr,imin,isec,jldy
INTEGER :: ifhr,ifmin,ifsec,mfhr
INTEGER :: myr,initsec,iabssec,jabssec,kftime
INTEGER :: ifile,iprtopt,iniotfu,lbasdmpf,onvf,grdbas
INTEGER :: iextmn,iextmx,jextmn,jextmx
INTEGER :: idiag,jdiag
REAL :: latnot(2)
REAL :: amin,amax
REAL :: qvmin,qvmax,qvval
REAL :: csconst,pconst
REAL :: deltaz,tv_ext
REAL :: pres,temp,qvsat,rh,tvbar,qvprt,qtot
REAL :: xdiag,ydiag,dd,dmin,latd,lond
REAL :: ppasc,pmb,tc,tdc,theta,smix,e,bige,alge,dir,spd
!wdt update
CHARACTER (LEN=80) :: soiloutfl,temchar,ternfn
CHARACTER (LEN=80) :: timsnd
INTEGER :: lfn, tmstrln, lternfn
INTEGER :: tsfcout,tsoilout,wetsout,wetdout,wetcout,snowdout
INTEGER :: isnow,jsnow,ii,jj
REAL :: xumin,xumax,yvmin,yvmax
INTEGER :: strlen, ireturn
CHARACTER (LEN=3) :: FMT
CHARACTER (LEN=100) :: tmp_ch
INTEGER :: is,nstyp_ext
!wdt-williams 2002-01-17 GMB ntmerge
DOUBLE PRECISION :: ntmergeinv, mfac
INTEGER :: idist
!
!-----------------------------------------------------------------------
!
! 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...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
mgrid = 1
nestgrd = 0
nstyps = 4
nstyp_ext = 0
DO k=1,lvlprof
dumsnd(k) = 0.0
END DO
!
!-----------------------------------------------------------------------
!
! Call initpara to read in ARPS parameters of namelists
!
!-----------------------------------------------------------------------
!
CALL initpara
(nx,ny,nz,nstyps)
CALL initadas ! The adjust namelist block is used by ext2arps. Other
! adas parameters are not used by this program,
! but the name list blocks need to be read in sequence on
! Cray machines such as J90.
myr=MOD(year,100)
WRITE(julfinit,'(i2.2,i3.3,i2.2,i2.2)') myr,jday,hour,minute
CALL ctim2abss(year,month,day,hour,minute,second,initsec)
!
!-----------------------------------------------------------------------
!
! Read in additional namelists for external file specifications.
!
!-----------------------------------------------------------------------
!
extdopt = 0
extdfmt = 10
dir_extd = './'
extdname = 'may20'
nextdfil = 1
extdtime(1) = '1977-05-20.21:00:00+000:00:00'
iorder = 2
intropt = 1
nsmooth = 1
ext_lbc = 1
ext_vbc = 1
exttrnopt = 0
!wdt-williams 2002-01-17 GMB ntmerge
extntmrg = 1
READ (5,extdfile)
WRITE (6, '(/1x,a)') '&extdfile'
strlen = LEN( dir_extd )
CALL strlnth( dir_extd, strlen)
WRITE (6, '(3x,a,a,a)') &
'dir_extd = ''', dir_extd(1:strlen), ''','
strlen = LEN( extdname )
CALL strlnth( extdname, strlen)
WRITE (6, '(3x,a,a,a)') &
'extdname = ''', extdname(1:strlen), ''','
WRITE (6, '(3x,a,i4,a)') 'extdopt = ', extdopt, ','
WRITE (6, '(3x,a,i4,a)') 'extdfmt = ', extdfmt, ','
WRITE (6, '(3x,a,i4,a)') 'nextdfil = ', nextdfil, ','
DO i=1,nextdfil
strlen = LEN( extdtime(i) )
CALL strlnth( extdtime(i), strlen)
WRITE (6, '(3x,a,i2.2,a,a,a)') &
'extdtime(',i,') = ''', extdtime(i), ''','
END DO
WRITE (6, '(3x,a,i4,a)') 'iorder = ', iorder, ','
WRITE (6, '(3x,a,i4,a)') 'intropt = ', intropt, ','
WRITE (6, '(3x,a,i4,a)') 'nsmooth = ', nsmooth, ','
WRITE (6, '(3x,a,i4,a)') 'hydradj = ', hydradj, ','
WRITE (6, '(3x,a,i4,a)') 'wndadj = ', wndadj, ','
WRITE (6, '(3x,a,i4,a)') 'obropt = ', obropt, ','
WRITE (6, '(3x,a,f16.4,a)') 'obrzero = ', obrzero, ','
WRITE (6, '(3x,a,i4,a)') 'ext_lbc = ', ext_lbc, ','
WRITE (6, '(3x,a,i4,a)') 'ext_vbc = ', ext_vbc, ','
WRITE (6, '(3x,a,i4,a)') 'exttrnopt = ', exttrnopt,','
!wdt-williams 2002-01-17 GMB ntmerge
WRITE (6, '(3x,a,i4,a)') 'extntmrg = ', extntmrg, ','
WRITE (6, '(1x,a)') '&END'
! 2 CONTINUE
!-----------------------------------------------------------------------
!
! Allocate and initalize arrays based on dimension parameters
! read in from the input file
!
!-----------------------------------------------------------------------
select case ( extdopt )
case ( 0 ) ! ARPS grid
select case ( extdfmt )
case ( 1 )
FMT = 'bin'
case ( 2 )
FMT = 'asc'
case ( 3 )
FMT = 'hdf'
case ( 4 )
FMT = 'pak'
case ( 6 )
FMT = 'bn2'
case ( 7 )
FMT = 'net'
case ( 8 )
FMT = 'npk'
case ( 10 )
FMT = 'grb'
case default
PRINT*,'Invalid input history data format given by extdfmt.'
STOP
END select
tmp_ch=trim(dir_extd)//trim(extdname)//'.'//FMT//'grdbas'
PRINT*,'tmp_ch =', trim(tmp_ch)
CALL get_dims_from_data(extdfmt, &
trim(dir_extd)//trim(extdname)//'.'//FMT//'grdbas', &
nx_ext,ny_ext,nz_ext,nstyp_ext, ireturn)
nstyp_ext = min(nstyp_ext,nstyps)
IF( ireturn /= 0 ) THEN
PRINT*, &
'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
case ( 1 ) ! RUC (Hybrid-B, GRIB #87) from NMC is 81x62x25
nx_ext=81
ny_ext=62
nz_ext=25
case ( 2 ) ! ETA (GRIB #212, 40km) from NMC is 185x129x39
nx_ext = 185
ny_ext = 129
nz_ext = 39
case ( 3 ) ! OLAPS for 95 is 91x73x41
nx_ext=91
ny_ext=73
nz_ext=41
case ( 4 ) ! RUC in GEMPAK format from Rossby (e.g., those of 1995)
nx_ext=81
ny_ext=62
nz_ext=41
case ( 5 ) ! GEMPAK ETA data
PRINT*,'Need to find out the grid size and set in ext2arps.f.'
STOP
! nx_ext = ??
! ny_ext = ??
! nz_ext = ??
case ( 6 ) ! Special COAMPS files
PRINT*,'Please specify the size of COAMPS input data grid '
PRINT*,'by editing file ext2arps.f.'
STOP
! nx_ext = ??
! ny_ext = ??
! nz_ext = ??
case ( 7 ) ! Isobaric RUC2 GRIB file (AWIPS #211) from OSO
nx_ext=93
ny_ext=65
nz_ext=37
case ( 8 ) ! Global Reanalysis on T62 Gaussian grid
nx_ext = 192
ny_ext = 94
nz_ext = 28
case ( 9 ) ! RUC-2 (GEMPAK) from Doplight is 151x113x41
nx_ext=151
ny_ext=113
nz_ext=41
case ( 10 ) ! ETA (GEMPAK #104, 80km) from NMC is 147x110x39
nx_ext=147
ny_ext=110
nz_ext=39
case ( 11 ) ! Native RUC2 GRIB file (Grid #236) from OSO
nx_ext=151
ny_ext=113
nz_ext=40
case ( 12 ) ! Isobaric RUC2 GRIB file (Grid #236) from OSO
nx_ext=151
ny_ext=113
nz_ext=37
case ( 13 ) ! AVN (GRIB #3, 1x1 degree) from NCEP is 360x181x26
nx_ext = 360
ny_ext = 181
nz_ext = 26
case ( 14 ) ! Binary format lat/lon data
nx_ext = 26
ny_ext = 16
nz_ext = 26
case default
WRITE(6,'(/a/)') &
'extdopt option was invalid. Please specify a new option.'
STOP
END select
PRINT*,'nx_ext, ny_ext, nz_ext = ', nx_ext, ny_ext, nz_ext
allocate( allocate( allocate(
allocate(xscl(nx),stat=istatus)
allocate(yscl(ny),stat=istatus)
allocate(zp(nx,ny,nz),stat=istatus)
allocate(zs(nx,ny,nz),stat=istatus)
allocate(j1(nx,ny,nz),stat=istatus)
allocate(j2(nx,ny,nz),stat=istatus)
allocate(j3(nx,ny,nz),stat=istatus)
allocate(aj3z(nx,ny,nz),stat=istatus)
allocate(hterain(nx,ny),stat=istatus)
!wdt update
allocate(htrn_ext(nx,ny),stat=istatus)
allocate(mapfct(nx,ny,8),stat=istatus)
allocate( allocate( allocate( allocate(pprt(nx,ny,nz),stat=istatus)
allocate(ptprt(nx,ny,nz),stat=istatus)
allocate(qv(nx,ny,nz),stat=istatus)
allocate(qc(nx,ny,nz),stat=istatus)
allocate(qr(nx,ny,nz),stat=istatus)
allocate(qi(nx,ny,nz),stat=istatus)
allocate(qs(nx,ny,nz),stat=istatus)
allocate(qh(nx,ny,nz),stat=istatus)
!
allocate(pbar(nx,ny,nz),stat=istatus)
allocate(ptbar(nx,ny,nz),stat=istatus)
allocate(qvbar(nx,ny,nz),stat=istatus)
allocate(ubar(nx,ny,nz),stat=istatus)
allocate(vbar(nx,ny,nz),stat=istatus)
allocate(wbar(nx,ny,nz),stat=istatus)
allocate(rhobar(nx,ny,nz),stat=istatus)
allocate(rhostr(nx,ny,nz),stat=istatus)
allocate(wcont(nx,ny,nz),stat=istatus)
allocate(csndsq(nx,ny,nz),stat=istatus)
allocate(tsfc(nx,ny,0:nstyps),stat=istatus)
allocate(tsoil(nx,ny,0:nstyps),stat=istatus)
allocate(wetsfc(nx,ny,0:nstyps),stat=istatus)
allocate(wetdp(nx,ny,0:nstyps),stat=istatus)
allocate(wetcanp(nx,ny,0:nstyps),stat=istatus)
allocate(snowdpth(nx,ny),stat=istatus)
allocate(soiltyp (nx,ny,nstyps),stat=istatus)
allocate(stypfrct(nx,ny,nstyps),stat=istatus)
allocate(vegtyp (nx,ny),stat=istatus)
allocate(lai (nx,ny),stat=istatus)
allocate(roufns (nx,ny),stat=istatus)
allocate(veg (nx,ny),stat=istatus)
allocate(tem1(nx,ny,nz),stat=istatus)
allocate(tem2(nx,ny,nz),stat=istatus)
allocate(tem3(nx,ny,nz),stat=istatus)
allocate(tem4(nx,ny,nz),stat=istatus)
allocate(tem5(nx,ny,nz),stat=istatus)
allocate(tem6(nx,ny,nz),stat=istatus)
allocate(tem7(nx,ny,nz),stat=istatus)
allocate(tem8(nx,ny,nz),stat=istatus)
allocate(xs2d(nx,ny),stat=istatus)
allocate(ys2d(nx,ny),stat=istatus)
allocate(xu2d(nx,ny),stat=istatus)
allocate(yu2d(nx,ny),stat=istatus)
allocate(xv2d(nx,ny),stat=istatus)
allocate(yv2d(nx,ny),stat=istatus)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: xw & yw
allocate(xw(nx,ny),stat=istatus)
allocate(yw(nx,ny),stat=istatus)
allocate(iscl(nx,ny),stat=istatus)
allocate(jscl(nx,ny),stat=istatus)
allocate(iu(nx,ny),stat=istatus)
allocate(ju(nx,ny),stat=istatus)
allocate(iv(nx,ny),stat=istatus)
allocate(jv(nx,ny),stat=istatus)
x=0.0
y=0.0
z=0.0
xscl=0.0
yscl=0.0
zp=0.0
zs=0.0
j1=0.0
j2=0.0
j3=0.0
aj3z=0.0
hterain=0.0
!wdt update
htrn_ext=0.0
mapfct=0.0
u=0.0
v=0.0
w=0.0
pprt=0.0
ptprt=0.0
qv=0.0
qc=0.0
qr=0.0
qi=0.0
qs=0.0
qh=0.0
!
pbar=0.0
ptbar=0.0
qvbar=0.0
ubar=0.0
vbar=0.0
wbar=0.0
rhobar=0.0
rhostr=0.0
wcont=0.0
csndsq=0.0
tsfc=0.0
tsoil=0.0
wetsfc=0.0
wetdp=0.0
wetcanp=0.0
snowdpth=0.0
soiltyp =0
stypfrct=0.0
stypfrct(:,:,1) =1.0
vegtyp =0
lai =0.0
roufns =0.0
veg =0.0
tem1=0.0
tem2=0.0
tem3=0.0
tem4=0.0
tem5=0.0
tem6=0.0
tem7=0.0
tem8=0.0
xs2d=0.0
ys2d=0.0
xu2d=0.0
yu2d=0.0
xv2d=0.0
yv2d=0.0
iscl=0
jscl=0
iu=0
ju=0
iv=0
jv=0
!
!-----------------------------------------------------------------------
!
! Allocate and initialize external grid variables
!
!-----------------------------------------------------------------------
!
allocate(x_ext(nx_ext),stat=istatus)
allocate(y_ext(ny_ext),stat=istatus)
allocate(xu_ext(nx_ext),stat=istatus)
allocate(yu_ext(ny_ext),stat=istatus)
allocate(xv_ext(nx_ext),stat=istatus)
allocate(yv_ext(ny_ext),stat=istatus)
allocate(lat_ext(nx_ext,ny_ext),stat=istatus)
allocate(lon_ext(nx_ext,ny_ext),stat=istatus)
allocate(latu_ext(nx_ext,ny_ext),stat=istatus)
allocate(lonu_ext(nx_ext,ny_ext),stat=istatus)
allocate(latv_ext(nx_ext,ny_ext),stat=istatus)
allocate(lonv_ext(nx_ext,ny_ext),stat=istatus)
allocate(zs_ext(nx,ny,nz_ext),stat=istatus)
!
allocate(p_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(hgt_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(zp2_ext(nx,ny,nz_ext),stat=istatus)
allocate(zp_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(t_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(u_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(v_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(w_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(qv_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(qc_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(qr_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(qi_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(qs_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(qh_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(tsfc_ext (nx_ext,ny_ext,0:nstyps),stat=istatus)
allocate(tsoil_ext (nx_ext,ny_ext,0:nstyps),stat=istatus)
allocate(wetsfc_ext (nx_ext,ny_ext,0:nstyps),stat=istatus)
allocate(wetdp_ext (nx_ext,ny_ext,0:nstyps),stat=istatus)
allocate(wetcanp_ext(nx_ext,ny_ext,0:nstyps),stat=istatus)
allocate(snowdpth_ext(nx_ext,ny_ext),stat=istatus)
allocate(soiltyp_ext (nx_ext,ny_ext,nstyps),stat=istatus)
soiltyp_ext = 0
allocate(stypfrct_ext(nx_ext,ny_ext,nstyps),stat=istatus)
allocate(vegtyp_ext (nx_ext,ny_ext),stat=istatus)
allocate(lai_ext (nx_ext,ny_ext),stat=istatus)
allocate(roufns_ext (nx_ext,ny_ext),stat=istatus)
allocate(veg_ext (nx_ext,ny_ext),stat=istatus)
allocate(trn_ext (nx_ext,ny_ext),stat=istatus)
allocate(psfc_ext (nx_ext,ny_ext),stat=istatus)
allocate(t_2m_ext (nx_ext,ny_ext),stat=istatus)
allocate(rh_2m_ext(nx_ext,ny_ext),stat=istatus)
allocate(u_10m_ext(nx_ext,ny_ext),stat=istatus)
allocate(v_10m_ext(nx_ext,ny_ext),stat=istatus)
allocate(dxfld(nx_ext),stat=istatus)
allocate(dyfld(ny_ext),stat=istatus)
allocate(rdxfld(nx_ext),stat=istatus)
allocate(rdyfld(ny_ext),stat=istatus)
allocate(dxfldu(nx_ext),stat=istatus)
allocate(dyfldu(ny_ext),stat=istatus)
allocate(rdxfldu(nx_ext),stat=istatus)
allocate(rdyfldu(ny_ext),stat=istatus)
allocate(dxfldv(nx_ext),stat=istatus)
allocate(dyfldv(ny_ext),stat=istatus)
allocate(rdxfldv(nx_ext),stat=istatus)
allocate(rdyfldv(ny_ext),stat=istatus)
allocate(tem1_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(tem2_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(tem3_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(tem4_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(tem5_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
allocate(xa_ext(nx_ext,ny_ext),stat=istatus)
allocate(ya_ext(nx_ext,ny_ext),stat=istatus)
allocate(avgzs_ext(nx,ny,nz_ext),stat=istatus)
!
x_ext=0.0
y_ext=0.0
xu_ext=0.0
yu_ext=0.0
xv_ext=0.0
yv_ext=0.0
lat_ext=0.0
lon_ext=0.0
latu_ext=0.0
lonu_ext=0.0
latv_ext=0.0
lonv_ext=0.0
zs_ext=0.0
!
p_ext=0.0
hgt_ext=0.0
zp2_ext=0.0
zp_ext=0.0
t_ext=0.0
u_ext=0.0
v_ext=0.0
w_ext=0.0
qv_ext=0.0
qc_ext=0.0
qr_ext=0.0
qi_ext=0.0
qs_ext=0.0
qh_ext=0.0
tsfc_ext =-999.0
tsoil_ext =-999.0
wetsfc_ext =-999.0
wetdp_ext =-999.0
wetcanp_ext =-999.0
snowdpth_ext=-999.0
soiltyp = -999.0
stypfrct= -999.0
vegtyp = -999
lai = -999.0
roufns = -999.0
veg = -999.0
trn_ext =0.0
psfc_ext =0.0
t_2m_ext =0.0
rh_2m_ext=0.0
u_10m_ext=0.0
v_10m_ext=0.0
dxfld=0.0
dyfld=0.0
rdxfld=0.0
rdyfld=0.0
dxfldu=0.0
dyfldu=0.0
rdxfldu=0.0
rdyfldu=0.0
dxfldv=0.0
dyfldv=0.0
rdxfldv=0.0
rdyfldv=0.0
tem1_ext=0.0
tem2_ext=0.0
tem3_ext=0.0
tem4_ext=0.0
tem5_ext=0.0
xa_ext=0.0
ya_ext=0.0
avgzs_ext=0.0
!
!-----------------------------------------------------------------------
!
! Set up ARPS map projection
!
!-----------------------------------------------------------------------
!
latnot(1)=trulat1
latnot(2)=trulat2
CALL setmapr(mapproj,sclfct,latnot,trulon)
!
!-----------------------------------------------------------------------
!
! Set up ARPS grid
!
!-----------------------------------------------------------------------
!
IF (exttrnopt == 1) THEN ! don't really do grid here, just set map proj
ternopt = -1
hterain = 0
END IF
CALL inigrd(nx,ny,nz, &
x,y,z,zp,hterain,mapfct,j1,j2,j3,tem1,tem2,tem3)
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
aj3z(i,j,k)=0.5*(j3(i,j,k)+j3(i,j,k-1))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set up the height levels on which to define the base-state.
!
!-----------------------------------------------------------------------
!
deltaz = depthp/(lvlprof-1)
WRITE(6,'(/a,/a,f7.2,a/)') &
' The base state is formed from a mean sounding created', &
' with a ',deltaz,' meter interval.'
DO k=1,lvlprof
zsnd(k) = (k-1)*deltaz
END DO
!
!-----------------------------------------------------------------------
!
! Loop through the data times provided via NAMELIST.
!
!-----------------------------------------------------------------------
!
iniotfu = 21 ! FORTRAN unit number used for data output
DO ifile=1,nextdfil
!
!-----------------------------------------------------------------------
!
! Time conversions.
! Formats: extdtime='1994-05-06.18:00:00+000:00:00'
! julfname='941261800'
!
!-----------------------------------------------------------------------
!
READ(extdtime(ifile),'(a19,1x,a9)') extdinit,extdfcst
IF(extdfcst == ' ') extdfcst='000:00:00'
READ(extdinit, &
'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)',ERR=920,END=920) &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ(extdfcst, &
'(i3,1x,i2,1x,i2)',ERR=4,END=4) ifhr,ifmin,ifsec
4 CONTINUE
mfhr=MOD(ifhr,24)
jldy = jldy + ifhr/24
WRITE(julfname, &
'(i2.2,i3.3,i2.2,i2.2)') myr,jldy,ihr,mfhr
CALL ctim2abss(iyr,imo,iday,ihr,imin,isec,iabssec)
jabssec=(ifhr*3600) + (ifmin*60) + ifsec + iabssec
kftime=jabssec-initsec
curtim=FLOAT(kftime)
WRITE(6,'(a,a9,a,/19x,a,a19,a/a,a/a,i16,a,/,i26,a)') &
' Calling rdextfil, looking for ', &
extdfcst,' hour forecast ', &
'initialized at ',extdinit,' UTC', &
' Julian filename: ',julfname, &
' Which is ',jabssec,' abs seconds or ',kftime, &
' seconds from the ARPS initial time.'
!
!-----------------------------------------------------------------------
!
! Call getextd to reads and converts data to ARPS units
!
!-----------------------------------------------------------------------
!
select case ( extdopt )
case ( 0 ) ! ARPS grid
CALL getarps(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname,nstyps, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext,latu_ext,lonu_ext,latv_ext,lonv_ext, &
p_ext,hgt_ext,zp_ext,t_ext,qv_ext, &
u_ext,tem4_ext,v_ext,tem3_ext,w_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
soiltyp_ext,stypfrct_ext,vegtyp_ext, &
lai_ext,roufns_ext,veg_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext, &
tem1_ext,tem2_ext,istatus)
! tem3_ext holds uatv_ext and tem4_ext holds vatu_ext until
! u & v rotated to new map projection below (see calls to uvetomp)
trn_ext(:,:) = zp_ext(:,:,2)
case ( 1 )
!
!-----------------------------------------------------------------------
!
! Get NMC RUC GRIB #87
!
!-----------------------------------------------------------------------
!
CALL getnmcruc87(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext, &
istatus)
case ( 2 )
!
!-----------------------------------------------------------------------
!
! Get NMC ETA GRIB #212
!
!-----------------------------------------------------------------------
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. soiltyp_ext
CALL getnmceta212(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext, &
snowdpth_ext,trn_ext,psfc_ext,soiltyp_ext, &
istatus)
nstyp_ext = 1
stypfrct_ext(:,:,1) = 1.
stypfrct_ext(:,:,2:nstyps) = 0.
case ( 3 )
!
!-----------------------------------------------------------------------
!
! Get LAPS data. Need LAPS library (see makearps for ext2arps.laps)
!
!-----------------------------------------------------------------------
!
CALL getlaps(nx_ext,ny_ext,nz_ext,dir_extd, &
extdinit,extdfcst,julfname,iabssec, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
istatus)
case ( 4 )
!
!-----------------------------------------------------------------------
!
! Get RUC/MAPS data in GEMPAK format with a special setup for GEMPAK
! path and library (see makearps for ext2arps.gemruc).
!
!-----------------------------------------------------------------------
!
CALL getgemruc(nx_ext,ny_ext,nz_ext,dir_extd,extdname, &
extdinit,extdfcst,julfname,iabssec, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
istatus, tem1_ext)
case ( 5 )
CALL getgemeta(nx_ext,ny_ext,nz_ext,dir_extd,extdname, &
extdinit,extdfcst,julfname,iabssec, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
istatus, tem1_ext)
case ( 6 )
!
!-----------------------------------------------------------------------
!
! Get COAMPS data
!
!-----------------------------------------------------------------------
!
CALL getcoamps(nx_ext, ny_ext, nz_ext, dir_extd, &
extdinit,extdfcst, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext,istatus)
case ( 7 )
!
!-----------------------------------------------------------------------
!
! Get NMC RUC2 GRIB #211
!
!-----------------------------------------------------------------------
!
CALL getnmcruc211(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext, t_2m_ext, rh_2m_ext, &
u_10m_ext, v_10m_ext, istatus)
case ( 8 )
!
!-----------------------------------------------------------------------
!
! Get NCEP global re-analysis on T62 Gaussian lat/lon grid
!
!-----------------------------------------------------------------------
!
CALL getreanalt62(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext, &
wetcanp_ext, trn_ext,psfc_ext, istatus)
case ( 9 )
!
!-----------------------------------------------------------------------
!
! Get NMC RUC2 GEM
!
!-----------------------------------------------------------------------
!
CALL getgemruc2(nx_ext,ny_ext,nz_ext,dir_extd,extdname, &
extdinit,extdfcst,julfname,iabssec, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
istatus, tem1_ext)
case ( 10 )
!
!-----------------------------------------------------------------------
!
! Get NMC ETA GEMPAK #104
!
!-----------------------------------------------------------------------
!
CALL getgemeta2(nx_ext,ny_ext,nz_ext,dir_extd,extdname, &
extdinit,extdfcst,julfname,iabssec, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
istatus, tem1_ext)
!
!-----------------------------------------------------------------------
!
! Get NCEP RUC2 Native Coordinate GRIB (#236)
!
!-----------------------------------------------------------------------
!
case ( 11 )
CALL getnmcrucn236(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext,snowdpth_ext, &
istatus)
!
!-----------------------------------------------------------------------
!
! Get NCEP RUC2 Isobaric Coordinate GRIB (#236)
!
!-----------------------------------------------------------------------
!
case ( 12 )
CALL getnmcrucp236(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext, t_2m_ext, rh_2m_ext, &
u_10m_ext, v_10m_ext, snowdpth_ext, &
istatus)
case ( 13 )
!
!-----------------------------------------------------------------------
!
! Get NCEP AVN GRIB #3
!
!-----------------------------------------------------------------------
!
CALL getncepavn3(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext, &
wetsfc_ext,wetdp_ext,wetcanp_ext, &
snowdpth_ext,trn_ext,psfc_ext, &
istatus)
case ( 14 ) ! Binary format lat/lon data
CALL get_avn_bin(nx_ext,ny_ext,nz_ext, &
dir_extd,extdname,extdinit,extdfcst,julfname, &
iproj_ext,scale_ext,trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext,p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
snowdpth_ext,trn_ext,psfc_ext, &
istatus)
case default
WRITE(6,'(/a/)') &
'extdopt option was invalid. Please specify a new option.'
STOP
END select
IF(istatus /= 1) GO TO 999
PRINT*,' '
CALL a3dmax0(lat_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'lat_ext_min= ', amin,', lat_ext_max=',amax
CALL a3dmax0(lon_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'lon_ext_min= ', amin,', lon_ext_max=',amax
!
CALL a3dmax0(p_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'p_ext_min = ', amin,', p_ext_max =',amax
CALL a3dmax0(hgt_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'hgt_ext_min= ', amin,', hgt_ext_max=',amax
CALL a3dmax0(t_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
't_ext_min = ', amin,', t_ext_max =',amax
CALL a3dmax0(u_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'u_ext_min = ', amin,', u_ext_max =',amax
CALL a3dmax0(v_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'v_ext_min = ', amin,', v_ext_max =',amax
CALL a3dmax0(w_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'w_ext_min = ', amin,', w_ext_max =',amax
CALL a3dmax0(qv_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qv_ext_min = ', amin,', qv_ext_max =',amax
CALL a3dmax0(qc_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qc_ext_min = ', amin,', qc_ext_max =',amax
CALL a3dmax0(qr_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qr_ext_min = ', amin,', qr_ext_max =',amax
CALL a3dmax0(qi_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qi_ext_min = ', amin,', qi_ext_max =',amax
CALL a3dmax0(qs_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qs_ext_min = ', amin,', qs_ext_max =',amax
CALL a3dmax0(qh_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,nz_ext,1,nz_ext,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qh_ext_min = ', amin,', qh_ext_max =',amax
!wdt note: FIXME add soiltype, etc, plus add nstyp
CALL a3dmax0(tsfc_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'tsfc_ext_min= ', amin,', tsfc_ext_max=',amax
CALL a3dmax0(tsoil_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'tsoil_ext_min= ', amin,', tsoil_ext_max=',amax
CALL a3dmax0(wetsfc_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'wetsfc_ext_min= ', amin,', wetsfc_ext_max=',amax
CALL a3dmax0(wetdp_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'wetdp_ext_min= ', amin,', wetdp_ext_max=',amax
CALL a3dmax0(wetcanp_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'wetcanp_ext_min= ', amin,', wetcanp_ext_max=',amax
CALL a3dmax0(snowdpth_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'snowd_ext_min= ', amin,', snow_ext_max=',amax
CALL a3dmax0(trn_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'trn_ext_min= ', amin,', trn_ext_max=',amax
CALL a3dmax0(psfc_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'psfc_ext_min= ', amin,', psfc_ext_max=',amax
CALL a3dmax0(t_2m_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'T_2m_ext_min= ', amin,', T_2m_ext_max=',amax
CALL a3dmax0(rh_2m_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'RH_2m_ext_min= ', amin,', RH_2m_ext_max=',amax
CALL a3dmax0(u_10m_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'U_10m_ext_min= ', amin,', U_10m_ext_max=',amax
CALL a3dmax0(v_10m_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'V_10m_ext_min= ', amin,', V_10m_ext_max=',amax
PRINT*,' '
!
!-----------------------------------------------------------------------
!
! First time through the time loop, calculate grid
! transformation info.
!
!-----------------------------------------------------------------------
!
IF(ifile == 1) THEN
DO i=1,nx-1
xscl(i)=0.5*(x(i)+x(i+1))
END DO
xscl(nx)=2.*xscl(nx-1) - xscl(nx-2)
DO j=1,ny-1
yscl(j)=0.5*(y(j)+y(j+1))
END DO
yscl(ny)=2.*yscl(ny-1) - yscl(ny-2)
!
!-----------------------------------------------------------------------
!
! Find lat,lon locations of ARPS scalar, u and v grids.
!
!-----------------------------------------------------------------------
!
CALL xytoll(nx,ny,xscl,yscl,tem1,tem2)
CALL xytoll(nx,ny, x,yscl,tem3,tem4)
CALL xytoll(nx,ny,xscl, y,tem5,tem6)
CALL getmapr(iproj_arps,scale_arps,latnot_arps, &
trlon_arps,xorig_arps,yorig_arps)
!
!-----------------------------------------------------------------------
!
! Find x,y locations of external grid.
!
!-----------------------------------------------------------------------
!
CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
CALL setorig(1,x0_ext,y0_ext)
DO j=1,ny_ext
CALL lltoxy(1,1,lat_ext(1,j),lon_ext(1,j), &
x_ext(1),y_ext(j))
END DO
DO i=1,nx_ext
CALL lltoxy(1,1,lat_ext(i,1),lon_ext(i,1), &
x_ext(i),y_ext(1))
END DO
IF ( extdopt == 0 ) THEN
DO j=1,ny_ext
CALL lltoxy(1,1,latu_ext(1,j),lonu_ext(1,j), &
xu_ext(1),yu_ext(j))
END DO
DO i=1,nx_ext
CALL lltoxy(1,1,latu_ext(i,1),lonu_ext(i,1), &
xu_ext(i),yu_ext(1))
END DO
DO j=1,ny_ext
CALL lltoxy(1,1,latv_ext(1,j),lonv_ext(1,j), &
xv_ext(1),yv_ext(j))
END DO
DO i=1,nx_ext
CALL lltoxy(1,1,latv_ext(i,1),lonv_ext(i,1), &
xv_ext(i),yv_ext(1))
END DO
ENDIF
IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
DO i=1,nx_ext
IF(x_ext(i) < 0.0) x_ext(i)=x_ext(i)+360.
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Find x,y locations of ARPS scalar grid in terms of the
! external grid.
!
!-----------------------------------------------------------------------
!
CALL lltoxy(nx,ny,tem1,tem2,xs2d,ys2d)
IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
DO i=1,nx
DO j=1,ny
IF( xs2d(i,j) < 0. ) xs2d(i,j)= xs2d(i,j) + 360.
END DO
END DO
END IF
CALL setijloc(nx,ny,nx_ext,ny_ext,xs2d,ys2d, &
x_ext,y_ext,iscl,jscl)
!
!-----------------------------------------------------------------------
!
! Find x,y locations of ARPS u grid in terms of the
! external grid.
!
!-----------------------------------------------------------------------
!
CALL lltoxy(nx,ny,tem3,tem4,xu2d,yu2d)
IF( extdopt == 8 .OR. extdopt == 13 ) THEN
DO i=1,nx
DO j=1,ny
IF ( xu2d(i,j) < 0. ) xu2d(i,j)= xu2d(i,j) + 360.
END DO
END DO
END IF
IF ( extdopt == 0 ) THEN
CALL setijloc(nx,ny,nx_ext,ny_ext,xu2d,yu2d, &
xu_ext,yu_ext,iu,ju)
ELSE
CALL setijloc(nx,ny,nx_ext,ny_ext,xu2d,yu2d, &
x_ext,y_ext,iu,ju)
ENDIF
!
!-----------------------------------------------------------------------
!
! Find x,y locations of ARPS v grid in terms of the
! external grid.
!
!-----------------------------------------------------------------------
!
CALL lltoxy(nx,ny,tem5,tem6,xv2d,yv2d)
IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
DO i=1,nx
DO j=1,ny
IF ( xv2d(i,j) < 0. ) xv2d(i,j)= xv2d(i,j) + 360.
END DO
END DO
END IF
IF ( extdopt == 0 ) THEN
CALL setijloc(nx,ny,nx_ext,ny_ext,xv2d,yv2d, &
xv_ext,yv_ext,iv,jv)
ELSE
CALL setijloc(nx,ny,nx_ext,ny_ext,xv2d,yv2d, &
x_ext,y_ext,iv,jv)
ENDIF
xumin=xu2d(1,1)
xumax=xu2d(1,1)
DO j=1,ny-1
xumin=MIN(xumin,xu2d(1 ,j))
xumax=MAX(xumax,xu2d(nx,j))
END DO
yvmin=yv2d(1,1)
yvmax=yv2d(1,1)
DO i=1,nx-1
yvmin=MIN(yvmin,yv2d(i,1 ))
yvmax=MAX(yvmax,yv2d(i,ny))
END DO
IF(xumin < x_ext(1).OR.xumax > x_ext(nx_ext).OR. &
yvmin < y_ext(1).OR.yvmax > y_ext(ny_ext)) THEN
WRITE(6,'(a/a)') &
' ARPS domain extends outside available external data', &
' domain, ext2arps aborted.'
WRITE (6,*) "xu min & ext:",xumin,x_ext(1)
WRITE (6,*) "xu max & ext:",xumax,x_ext(nx_ext)
WRITE (6,*) "yv min & ext:",yvmin,y_ext(1)
WRITE (6,*) "yv max & ext:",yvmax,y_ext(ny_ext)
STOP
END IF
!
!-----------------------------------------------------------------------
!
! Test code, for diagnostic testing.
! Find x,y of Norman sounding in external grid.
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,latdiag,londiag,xdiag,ydiag)
IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
IF ( xdiag < 0. ) xdiag= xdiag + 360.
END IF
dmin=((xdiag-x_ext(1))*(xdiag-x_ext(1))+ &
(ydiag-y_ext(1))*(ydiag-y_ext(1)))
idiag=1
jdiag=1
!
DO j=1,ny_ext
DO i=1,nx_ext
dd=((xdiag-x_ext(i))*(xdiag-x_ext(i))+ &
(ydiag-y_ext(j))*(ydiag-y_ext(j)))
IF(dd < dmin) THEN
dmin=dd
idiag=i
jdiag=j
END IF
END DO
END DO
WRITE(6,'(a,f10.4,f10.4,/a,i5,i5,a,f10.4,f10.4)') &
' Nearest ext pt to diagnostic lat,lon: ', &
latdiag,londiag, &
' Diagnostic i,j: ', &
idiag,jdiag,' lat,lon= ', &
lat_ext(idiag,jdiag),lon_ext(idiag,jdiag)
WRITE(6,'(///a,/2x,a)') ' External sounding at idiag,jdiag', &
'k pres hgt temp theta dewp u v dir spd'
!
!-----------------------------------------------------------------------
!
! Convert units of external data and write as a sounding.
!
!-----------------------------------------------------------------------
!
DO k=nz_ext,1,-1
pmb=.01*p_ext(idiag,jdiag,k)
tc=t_ext(idiag,jdiag,k)-273.15
theta=t_ext(idiag,jdiag,k)*(p0/p_ext(idiag,jdiag,k))**rddcp
IF( qv_ext(idiag,jdiag,k) > 0.) THEN
smix=qv_ext(idiag,jdiag,k)/(1.-qv_ext(idiag,jdiag,k))
e=(pmb*smix)/(0.62197 + smix)
bige=e/( 1.001 + ( (pmb - 100.) / 900.) * 0.0034)
alge = ALOG(bige/6.112)
tdc = (alge * 243.5) / (17.67 - alge)
ELSE
tdc = tc-30.
END IF
CALL uvrotdd(1,1,lon_ext(idiag,jdiag), &
u_ext(idiag,jdiag,k),v_ext(idiag,jdiag,k), &
dir,spd)
WRITE(6,'(i4,f6.0,f9.0,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1)') &
k,pmb, &
hgt_ext(idiag,jdiag,k), &
tc,theta,tdc, &
u_ext(idiag,jdiag,k), &
v_ext(idiag,jdiag,k), &
dir,spd
END DO
!
!-----------------------------------------------------------------------
!
! Restore map projection to ARPS grid.
!
!-----------------------------------------------------------------------
!
CALL setmapr(iproj_arps,scale_arps,latnot_arps, &
trlon_arps)
CALL setorig(1,xorig_arps,yorig_arps)
!
!-----------------------------------------------------------------------
!
! Find the min and max of the external indices that will be used
! to be sure interpolation is within external dataset.
!
!-----------------------------------------------------------------------
!
iextmn=iscl(1,1)
jextmn=jscl(1,1)
iextmx=iscl(1,1)
jextmx=jscl(1,1)
DO j=1,ny
DO i=1,nx
iextmn=MIN(iextmn,iscl(i,j))
jextmn=MIN(jextmn,jscl(i,j))
iextmx=MAX(iextmx,iscl(i,j))
jextmx=MAX(jextmx,jscl(i,j))
iextmn=MIN(iextmn,iu(i,j))
jextmn=MIN(jextmn,ju(i,j))
iextmx=MAX(iextmx,iu(i,j))
jextmx=MAX(jextmx,ju(i,j))
iextmn=MIN(iextmn,iv(i,j))
jextmn=MIN(jextmn,jv(i,j))
iextmx=MAX(iextmx,iv(i,j))
jextmx=MAX(jextmx,jv(i,j))
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! External data comes in oriented to true north.
! Change that to u,v in the new coordinate system.
! (tem3_ext & tem4_ext from above hold uatv_ext & vatu_ext for use here)
!
!-----------------------------------------------------------------------
!
IF ( extdopt == 0 ) THEN
DO k=1,nz_ext
CALL uvetomp(nx_ext,ny_ext,u_ext(1,1,k),tem4_ext(1,1,k), &
lonu_ext,tem1_ext(1,1,k),tem2_ext(1,1,k))
END DO
u_ext = tem1_ext
DO k=1,nz_ext
CALL uvetomp(nx_ext,ny_ext,tem3_ext(1,1,k),v_ext(1,1,k), &
lonv_ext,tem1_ext(1,1,k),tem2_ext(1,1,k))
END DO
v_ext = tem2_ext
ELSE
DO k=1,nz_ext
CALL uvetomp(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k), &
lon_ext,tem1_ext(1,1,k),tem2_ext(1,1,k))
END DO
u_ext = tem1_ext
v_ext = tem2_ext
ENDIF
!
!-----------------------------------------------------------------------
!
! Process height data
!
! Interpolate the external heights horizontally to the
! ARPS x,y. This is for scalar pts.
!
!-----------------------------------------------------------------------
!
CALL setdxdy(nx_ext,ny_ext, &
1,nx_ext,1,ny_ext, &
x_ext,y_ext,dxfld,dyfld,rdxfld,rdyfld)
IF ( extdopt == 0 ) THEN
CALL setdxdy(nx_ext,ny_ext, &
1,nx_ext,1,ny_ext, &
xu_ext,yu_ext,dxfldu,dyfldu,rdxfldu,rdyfldu)
CALL setdxdy(nx_ext,ny_ext, &
1,nx_ext,1,ny_ext, &
xv_ext,yv_ext,dxfldv,dyfldv,rdxfldv,rdyfldv)
DO k=1,nz_ext
CALL fldint2d(nx,ny,nx_ext,ny_ext, &
1,nx,1,ny, &
1,nx_ext,1,ny_ext, &
iorder,xs2d,ys2d,zp_ext(1,1,k), &
x_ext,y_ext,iscl,jscl, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
zp2_ext(1,1,k))
END DO
ENDIF
DO k=1,nz_ext
CALL fldint2d(nx,ny,nx_ext,ny_ext, &
1,nx,1,ny, &
1,nx_ext,1,ny_ext, &
iorder,xs2d,ys2d,hgt_ext(1,1,k), &
x_ext,y_ext,iscl,jscl, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
zs_ext(1,1,k))
END DO
!
!-----------------------------------------------------------------------
!
! Calculate x and y coordinates of external grid in ARPS coordinate
! system. Store them in xa_ext and ya_ext.
!
!-----------------------------------------------------------------------
!
CALL lltoxy(nx_ext,ny_ext,lat_ext,lon_ext, &
xa_ext,ya_ext)
IF(ifile == 1) THEN
!
!-----------------------------------------------------------------------
!
! Interpolate external terrain to arps grid (if exttrnopt=1).
!
!-----------------------------------------------------------------------
!
!wdt update exttrnopt
!wdt-williams 2002-01-17 GMB ntmerge
IF (exttrnopt > 0) THEN ! set arps terrain to match external terrain
CALL mkarps2d (nx_ext,ny_ext,nx,ny, &
iorder,iscl,jscl,x_ext,y_ext, &
xs2d,ys2d,trn_ext,htrn_ext, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext(1,1,1),tem1_ext(1,1,2), &
tem1_ext(1,1,3),tem1_ext(1,1,4))
IF (exttrnopt == 1) THEN
hterain = htrn_ext
!wdt-williams 2002-01-17 GMB ntmerge
ELSE
ntmergeinv = 1.d0/extntmrg
DO j = 1,ny
DO i = 1,nx
idist = max(0,min(extntmrg,i-2,nx-2-i,j-2,ny-2-j))
mfac = idist*ntmergeinv
hterain(i,j) = (1.d0-mfac)*htrn_ext(i,j) &
+ mfac*hterain(i,j)
END DO
END DO
END IF
ternopt = -1
CALL inigrd(nx,ny,nz, &
x,y,z,zp,hterain,mapfct,j1,j2,j3,tem1,tem2,tem3)
DO k=2,nz-1
DO j=1,ny-1
DO i=1,nx-1
aj3z(i,j,k)=0.5*(j3(i,j,k)+j3(i,j,k-1))
END DO
END DO
END DO
END IF
!wdt update 2002-01-17
!-----------------------------------------------------------------------
!
! Write out terrain data
!
!-----------------------------------------------------------------------
IF (terndmp > 0) THEN
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(nx,ny,ternfn(1:lternfn), dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct, &
ctrlat,ctrlon,hterain)
!FIXME: add grads control file
!IF (terndmp == 1) &
! CALL trncntl(nx,ny,ternfn(1:lternfn), x1_out,y1_out, &
! tem11,tem21)
END IF
!
!-----------------------------------------------------------------------
!
! Set up z grid at scalar vertical levels.
!
!-----------------------------------------------------------------------
!
onvf = 0
CALL avgz(zp, onvf, nx,ny,nz, 1,nx, 1,ny, 1,nz-1, zs)
!
DO j=1,ny-1
DO i=1,nx-1
zs(i,j,nz)=2.*zs(i,j,nz-1) - zs(i,j,nz-2)
END DO
END DO
CALL edgfill(zs, 1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
CALL edgfill(zp, 1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
END IF
!
!-----------------------------------------------------------------------
!
! Data transformations
! Calculate log of pressure for external grid.
! Change qv to RHstar RHStar=sqrt(1.-relative humidity)
!
!-----------------------------------------------------------------------
!
qvmin=999.
qvmax=-999.
DO k=1,nz_ext
DO j=1,ny_ext
DO i=1,nx_ext
qvmax=AMAX1(qvmax,qv_ext(i,j,k))
qvmin=AMIN1(qvmin,qv_ext(i,j,k))
tem5_ext(i,j,k)=ALOG(p_ext(i,j,k))
qvsat=f_qvsat( p_ext(i,j,k), t_ext(i,j,k) )
qv_ext(i,j,k)=SQRT(AMAX1(0.,(rhmax-(qv_ext(i,j,k)/qvsat))))
! u_ext(i,j,k)=tem1_ext(i,j,k)
! v_ext(i,j,k)=tem2_ext(i,j,k)
END DO
END DO
END DO
PRINT *, ' qv_ext min = ',qvmin
PRINT *, ' qv_ext max = ',qvmax
!
!-----------------------------------------------------------------------
!
! Calculate base-state sounding (vertical profile)
!
!-----------------------------------------------------------------------
!
CALL extmnsnd(nx,ny,nx_ext,ny_ext,nz_ext,lvlprof, &
iextmn,iextmx,jextmn,jextmx,1,nz_ext, &
xscl,yscl,xa_ext,ya_ext, &
hgt_ext,tem5_ext,t_ext, &
qv_ext,u_ext,v_ext, &
zsnd,plsnd,psnd,tsnd,ptsnd,rhssnd,qvsnd, &
rhosnd,usnd,vsnd)
!
!-----------------------------------------------------------------------
!
! Process Pressure data
!
!-----------------------------------------------------------------------
!
iprtopt=1
CALL mkarpsvlz(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,p_ext, &
zsnd,plsnd,pbar,pprt, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
DO k=1,nz
CALL smooth9p(pprt(1,1,k), nx,ny, 1,nx,1,ny, tem1)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Process potential temperature data
!
!-----------------------------------------------------------------------
!
DO k=1,nz_ext
DO j=1,ny_ext
DO i=1,nx_ext
tem5_ext(i,j,k)=t_ext(i,j,k)*((p0/p_ext(i,j,k))**rddcp)
END DO
END DO
END DO
iprtopt=1
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,tem5_ext, &
zsnd,ptsnd,ptbar,ptprt, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs,ptprt,tem1,ptprt)
! CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
! : smfct2,zs,ptprt,tem1,ptprt)
END DO
!
!-----------------------------------------------------------------------
!
! Process qv data.
!
!-----------------------------------------------------------------------
!
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,qv_ext, &
zsnd,rhssnd,qvbar,qv, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs, qv,tem1, qv)
! CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
! : smfct2,zs, qv,tem1, qv)
END DO
!
!-----------------------------------------------------------------------
!
! Convert rhstar back to qv for writing, further calculations
!
!-----------------------------------------------------------------------
!
qvmax=-999.
qvmin=999.
DO k=1,nz_ext
DO j=1,ny_ext
DO i=1,nx_ext
qvsat=f_qvsat( p_ext(i,j,k), t_ext(i,j,k) )
rh=AMAX1(0.01,rhmax-(qv_ext(i,j,k)*qv_ext(i,j,k)))
qv_ext(i,j,k)=rh*qvsat
qvmax=AMAX1(qvmax,qv_ext(i,j,k))
qvmin=AMIN1(qvmin,qv_ext(i,j,k))
END DO
END DO
END DO
PRINT *, ' qv_ext min = ',qvmin
PRINT *, ' qv_ext max = ',qvmax
!
qvmax=-999.
qvmin=999.
DO k=1,nz
DO j=1,ny
DO i=1,nx
pres = pbar(i,j,k)+pprt(i,j,k)
temp = (ptbar(i,j,k)+ptprt(i,j,k))*((pres/p0)**rddcp)
qvsat=f_qvsat( pres, temp )
IF(qvsat > 1.) THEN
PRINT *, ' qvsat trouble at: ',i,j,k
PRINT *, ' temp,press,qvsat: ',temp,pres,qvsat
END IF
rh=AMAX1(0.01,rhmax-(qv(i,j,k)*qv(i,j,k)))
rh=AMIN1(rh,rhmax)
qv(i,j,k)=rh*qvsat
IF(qv(i,j,k) > 0.5) THEN
PRINT *, ' qvsat trouble at: ',i,j,k
PRINT *, ' temp,press,qvsat: ',temp,pres,qv(i,j,k)
PRINT *, ' rh= ',rh
END IF
qvmax=AMAX1(qvmax,qv(i,j,k))
qvmin=AMIN1(qvmin,qv(i,j,k))
pres = pbar(i,j,k)
temp = ptbar(i,j,k)*((pres/p0)**rddcp)
qvsat=f_qvsat( pres, temp )
rh=AMAX1(0.01,rhmax-(qvbar(i,j,k)*qvbar(i,j,k)))
rh=AMIN1(rh,rhmax)
qvbar(i,j,k)=rh*qvsat
END DO
END DO
END DO
PRINT *, ' qv min = ',qvmin
PRINT *, ' qv max = ',qvmax
!
!-----------------------------------------------------------------------
!
! Process qc data.
! If first data point is less than or equal to -1, then
! no qc info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
IF(qc_ext(1,1,1) > -1.0) THEN
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,qc_ext, &
zsnd,dumsnd,tem1,qc, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs, qc,tem1, qc)
END DO
ELSE
DO k=1,nz
DO j=1,ny
DO i=1,nx
qc(i,j,k)=0.
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Process qr data.
! If first data point is less than or equal to -1, then
! no qr info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
IF(qr_ext(1,1,1) > -1.0) THEN
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,qr_ext, &
zsnd,dumsnd,tem1,qr, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs, qr,tem1, qr)
END DO
DO k=1,nz
DO j=1,ny
DO i=1,nx
qr(i,j,k)=MAX(qr(i,j,k),0.)
END DO
END DO
END DO
ELSE
DO k=1,nz
DO j=1,ny
DO i=1,nx
qr(i,j,k)=0.
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Process qi data.
! If first data point is less than or equal to -1, then
! no qi info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
iceout=0
IF(qi_ext(1,1,1) > -1.0) THEN
iceout=1
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,qi_ext, &
zsnd,dumsnd,tem1,qi, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs, qi,tem1, qi)
END DO
DO k=1,nz
DO j=1,ny
DO i=1,nx
qi(i,j,k)=MAX(qi(i,j,k),0.)
END DO
END DO
END DO
ELSE
DO k=1,nz
DO j=1,ny
DO i=1,nx
qi(i,j,k)=0.
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Process qs data.
! If first data point is less than or equal to -1, then
! no qs info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
IF(qs_ext(1,1,1) > -1.0) THEN
iceout=1
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,qs_ext, &
zsnd,dumsnd,tem1,qs, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs, qs,tem1, qs)
END DO
DO k=1,nz
DO j=1,ny
DO i=1,nx
qs(i,j,k)=MAX(qs(i,j,k),0.)
END DO
END DO
END DO
ELSE
DO k=1,nz
DO j=1,ny
DO i=1,nx
qs(i,j,k)=0.
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Process qh data.
! If first data point is less than or equal to -1, then
! no qh info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
IF(qh_ext(1,1,1) > -1.0) THEN
iceout=1
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,qh_ext, &
zsnd,dumsnd,tem1,qh, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs, qh,tem1, qh)
END DO
DO k=1,nz
DO j=1,ny
DO i=1,nx
qh(i,j,k)=MAX(qh(i,j,k),0.)
END DO
END DO
END DO
ELSE
DO k=1,nz
DO j=1,ny
DO i=1,nx
qh(i,j,k)=0.
END DO
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Process density which has been stuffed into tem2_ext array.
! We really only need rhobar, so set iprtopt=1 and pass
! a temporary array in place of rhoprt.
!
!-----------------------------------------------------------------------
!
DO k=1,nz_ext
DO j=1,ny_ext
DO i=1,nx_ext
tv_ext=t_ext(i,j,k) / &
((1.0+qv_ext(i,j,k))* &
(1.0-(qv_ext(i,j,k)/(rddrv+qv_ext(i,j,k)))))
tem5_ext(i,j,k)=p_ext(i,j,k)/(rd*tv_ext)
END DO
END DO
END DO
!
iprtopt=1
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
hgt_ext,zs_ext,xs2d,ys2d,zs,tem5_ext, &
zsnd,rhosnd,rhobar,tem1, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zs,rhobar,tem1,rhobar)
! CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
! : smfct2,zs,rhobar,tem1,rhobar)
END DO
!
!-----------------------------------------------------------------------
!
! Calculate and store the sound wave speed squared in csndsq.
!
!-----------------------------------------------------------------------
!
IF(csopt == 1) THEN ! Original definition of sound speed.
DO k=1,nz
DO j=1,ny
DO i=1,nx
csndsq(i,j,k)= cpdcv*pbar(i,j,k)/rhobar(i,j,k)
END DO
END DO
END DO
ELSE IF(csopt == 2) THEN ! Original sound speed multiplied
! by a factor
csconst = csfactr**2*cpdcv
DO k=1,nz
DO j=1,ny
DO i=1,nx
csndsq(i,j,k)= csconst * pbar(i,j,k)/rhobar(i,j,k)
END DO
END DO
END DO
ELSE ! Specified constant sound speed.
DO k=1,nz
DO j=1,ny
DO i=1,nx
csndsq(i,j,k)= csound * csound
END DO
END DO
END DO
END IF
!
IF(hydradj == 1 .OR. hydradj == 2) THEN
pconst=0.5*g*dz
!
!-----------------------------------------------------------------------
!
! Create thermal buoyancy at each scalar point,
! which is stored in tem2
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
qvprt=qv(i,j,k)-qvbar(i,j,k)
qtot=qc(i,j,k)+qr(i,j,k)+qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
tem2(i,j,k)=j3(i,j,k)*rhobar(i,j,k)* &
g*( (ptprt(i,j,k)/ptbar(i,j,k)) + &
(qvprt/(rddrv+qvbar(i,j,k))) - &
((qvprt+qtot)/(1.+qvbar(i,j,k))) )
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Average thermal buoyancy to w points
! which is stored in tem3
!
!-----------------------------------------------------------------------
!
CALL avgsw(tem2,nx,ny,nz,1,nx,1,ny, tem3)
IF(hydradj == 1) THEN
DO i=1,nx
DO j=1,ny
tem1(i,j,1)=pprt(i,j,1)
DO k=2,nz-2
tem1(i,j,k)= &
( (1. - (pconst*j3(i,j,k-1)/csndsq(i,j,k-1)) )* &
tem1(i,j,k-1) + &
dz*tem3(i,j,k) ) / &
(1. + (pconst*j3(i,j,k)/csndsq(i,j,k)) )
IF(j == 26 .AND. MOD(i,10) == 0) THEN
IF(k == nz-2) PRINT *, ' Point i= ',i, ' j=26'
PRINT *, ' k,pprt,tem1,diff =', &
k,pprt(i,j,k),tem1(i,j,k), &
(tem1(i,j,k)-pprt(i,j,k))
END IF
pprt(i,j,k)=tem1(i,j,k)
END DO
pprt(i,j,nz-1)=pprt(i,j,nz-2)
END DO
END DO
ELSE IF(hydradj == 2) THEN
DO i=1,nx
DO j=1,ny
tem1(i,j,nz-1)=pprt(i,j,nz-1)
DO k=nz-2,2,-1
tem1(i,j,k)= &
( (1.+ (pconst*j3(i,j,k+1)/csndsq(i,j,k+1)) )* &
tem1(i,j,k+1) - &
dz*tem3(i,j,k+1) ) / &
(1.- (pconst*j3(i,j,k )/csndsq(i,j,k )) )
IF(j == 26 .AND. MOD(i,10) == 0) THEN
IF(k == nz-2) PRINT *, ' Point i= ',i, ' j=26'
PRINT *, ' k,pprt,tem1,diff =', &
k,pprt(i,j,k),tem1(i,j,k), &
(tem1(i,j,k)-pprt(i,j,k))
END IF
pprt(i,j,k)=tem1(i,j,k)
END DO
pprt(i,j,1)=pprt(i,j,2)
END DO
END DO
END IF
ELSE IF (hydradj == 3) THEN
!
!-----------------------------------------------------------------------
!
! Calculate total pressure, store in tem1.
! Calculate virtual temperature, store in tem2.
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k) = pbar(i,j,k)+pprt(i,j,k)
temp = (ptbar(i,j,k)+ptprt(i,j,k))* &
((tem1(i,j,k)/p0)**rddcp)
tem2(i,j,k) = temp*(1.0+rvdrd*qv(i,j,k))/(1.0+qv(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Integrate hydrostatic equation, begining at k=2
!
!-----------------------------------------------------------------------
!
pconst=-g/rd
DO k=3,nz-1
DO j=1,ny
DO i=1,nx
tvbar=0.5*(tem2(i,j,k)+tem2(i,j,k-1))
tem1(i,j,k)=tem1(i,j,k-1)* &
EXP(pconst*(zs(i,j,k)-zs(i,j,k-1))/tvbar)
pprt(i,j,k)=tem1(i,j,k)-pbar(i,j,k)
END DO
END DO
END DO
DO j=1,ny
DO i=1,nx
tvbar=0.5*(tem2(i,j,2)+tem2(i,j,1))
tem1(i,j,1)=tem1(i,j,2)* &
EXP(pconst*(zs(i,j,1)-zs(i,j,2))/tvbar)
pprt(i,j,1)=tem1(i,j,1)-pbar(i,j,1)
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Process wind data
!
!-----------------------------------------------------------------------
!
CALL avgsu(zs_ext,nx,ny,nz_ext, 1,ny, 1,nz_ext, avgzs_ext, tem3_ext)
CALL avgsu(zs,nx,ny,nz, 1,ny, 1,nz, tem2, tem3)
DO k=1,nz_ext
DO j=1,ny
avgzs_ext(1,j,k)=zs_ext(1,j,k)
avgzs_ext(nx,j,k)=zs_ext(nx-1,j,k)
END DO
END DO
DO k=1,nz
DO j=1,ny
tem2(1,j,k)=zs(1,j,k)
tem2(nx,j,k)=zs(nx-1,j,k)
END DO
END DO
IF ( extdopt == 0 ) THEN ! u is staggered as in arps
CALL avgsu(hgt_ext,nx_ext,ny_ext,nz_ext,1,ny_ext,1,nz_ext,tem5_ext, &
tem3_ext)
tem5_ext(1,1:ny_ext,1:nz_ext) = hgt_ext(1,1:ny_ext,1:nz_ext)
tem5_ext(nx,1:ny_ext,1:nz_ext) = hgt_ext(nx-1,1:ny_ext,1:nz_ext)
!
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iu,ju,xu_ext,yu_ext, &
tem5_ext,avgzs_ext,xu2d,yu2d,tem2,u_ext, &
zsnd,usnd,ubar,u, &
dxfldu,dyfldu,rdxfldu,rdyfldu, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
ELSE ! u is on a scalar point
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iu,ju,x_ext,y_ext, &
hgt_ext,avgzs_ext,xu2d,yu2d,tem2,u_ext, &
zsnd,usnd,ubar,u, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
ENDIF
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,tem2, u,tem1, u)
! call smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
! : smfct2,tem2, u,tem1, u)
END DO
CALL a3dmax0(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'umin final= ', amin,', umax final=',amax
!
!-----------------------------------------------------------------------
!
! Process v component
!
!-----------------------------------------------------------------------
!
CALL avgsv(zs_ext,nx,ny,nz_ext, 1,nx, 1,nz_ext, avgzs_ext, tem3_ext)
CALL avgsv(zs,nx,ny,nz, 1,nx, 1,nz, tem2, tem3)
DO k=1,nz_ext
DO i=1,nx
avgzs_ext(i,1,k)=zs_ext(i,1,k)
avgzs_ext(i,ny,k)=zs_ext(i,ny-1,k)
END DO
END DO
DO k=1,nz
DO i=1,nx
tem2(i,1,k)=zs(i,1,k)
tem2(i,ny,k)=zs(i,ny-1,k)
END DO
END DO
IF ( extdopt == 0 ) THEN ! v is staggered as in arps
CALL avgsv(hgt_ext,nx_ext,ny_ext,nz_ext,1,nx_ext,1,nz_ext,tem5_ext, &
tem3_ext)
tem5_ext(1:nx_ext,1,1:nz_ext) = hgt_ext(1:nx_ext,1,1:nz_ext)
tem5_ext(1:nx_ext,ny,1:nz_ext) = hgt_ext(1:nx_ext,ny-1,1:nz_ext)
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iv,jv,xv_ext,yv_ext, &
tem5_ext,avgzs_ext,xv2d,yv2d,tem2,v_ext, &
zsnd,vsnd,vbar,v, &
dxfldv,dyfldv,rdxfldv,rdyfldv, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
ELSE ! v is on a scalar point
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iv,jv,x_ext,y_ext, &
hgt_ext,avgzs_ext,xv2d,yv2d,tem2,v_ext, &
zsnd,vsnd,vbar,v, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
ENDIF
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,tem2, v,tem1, v)
END DO
CALL a3dmax0(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'vmin final= ', amin,', vmax final=',amax
!
!-----------------------------------------------------------------------
!
! Process 2D surface fields if required (sfcout = 1)
!
!-----------------------------------------------------------------------
!
! determine the number of soil types supplied
IF (nstyp_ext <= 0) THEN
DO is=1,nstyps
IF (tsfc_ext(1,1,is) > 0) nstyp_ext = is
ENDDO
ENDIF
! GMB 2001-03-13 wdt
! This should never be the case (or else ext sfc arrays allocated to small!):
nstyp_ext = min(nstyp_ext,nstyps)
tsfcout = 0
tsoilout = 0
wetsout = 0
wetdout = 0
wetcout = 0
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
!wdt begin block:
!commented out mkarps2d & smooth9's, add call to intrp_soil
DO is=0,nstyp_ext
CALL a3dmax0(tsfc_ext(1,1,is),1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext, &
1,1,1,1,amax,amin)
WRITE(6,'(1x,a,i3,2(a,e13.6))') 'is=',is, &
'tsfc_ext_min= ', amin,', tsfc_ext_max=',amax
IF ( tsfc_ext(1,1,is) > 0.0 ) THEN
tsfcout = 1
! CALL mkarps2d (nx_ext,ny_ext,nx,ny, &
! iorder,iscl,jscl,x_ext,y_ext, &
! xs2d,ys2d,tsfc_ext(1,1,is),tsfc(1,1,is), &
! dxfld,dyfld,rdxfld,rdyfld, &
! tem1_ext(1,1,1),tem1_ext(1,1,2), &
! tem1_ext(1,1,3),tem1_ext(1,1,4))
! DO ksmth=1,nsmooth
! CALL smooth9p(tsfc(1,1,is), nx,ny, 1,nx,1,ny, tem1)
! END DO
END IF
IF ( tsoil_ext(1,1,is) > 0.0 ) THEN
tsoilout = 1
! CALL mkarps2d (nx_ext,ny_ext,nx,ny, &
! iorder,iscl,jscl,x_ext,y_ext, &
! xs2d,ys2d,tsoil_ext(1,1,is),tsoil(1,1,is), &
! dxfld,dyfld,rdxfld,rdyfld, &
! tem1_ext(1,1,1),tem1_ext(1,1,2), &
! tem1_ext(1,1,3),tem1_ext(1,1,4))
! DO ksmth=1,nsmooth
! CALL smooth9p(tsoil(1,1,is), nx,ny, 1,nx,1,ny, tem1)
! END DO
END IF
IF ( wetsfc_ext(1,1,is) >= 0.0 ) THEN
wetsout = 1
! CALL mkarps2d (nx_ext,ny_ext,nx,ny, &
! iorder,iscl,jscl,x_ext,y_ext, &
! xs2d,ys2d,wetsfc_ext(1,1,is),wetsfc(1,1,is), &
! dxfld,dyfld,rdxfld,rdyfld, &
! tem1_ext(1,1,1),tem1_ext(1,1,2), &
! tem1_ext(1,1,3),tem1_ext(1,1,4))
! DO ksmth=1,nsmooth
! CALL smooth9p(wetsfc(1,1,is), nx,ny, 1,nx,1,ny, tem1)
! END DO
END IF
IF ( wetdp_ext(1,1,is) >= 0.0 ) THEN
wetdout = 1
! CALL mkarps2d (nx_ext,ny_ext,nx,ny, &
! iorder,iscl,jscl,x_ext,y_ext, &
! xs2d,ys2d,wetdp_ext(1,1,is),wetdp(1,1,is), &
! dxfld,dyfld,rdxfld,rdyfld, &
! tem1_ext(1,1,1),tem1_ext(1,1,2), &
! tem1_ext(1,1,3),tem1_ext(1,1,4))
! DO ksmth=1,nsmooth
! CALL smooth9p(wetdp(1,1,is), nx,ny, 1,nx,1,ny, tem1)
! END DO
END IF
IF ( wetcanp_ext(1,1,is) >= 0.0 ) THEN
wetcout = 1
! CALL mkarps2d (nx_ext,ny_ext,nx,ny, &
! iorder,iscl,jscl,x_ext,y_ext, &
! xs2d,ys2d,wetcanp_ext(1,1,is),wetcanp(1,1,is), &
! dxfld,dyfld,rdxfld,rdyfld, &
! tem1_ext(1,1,1),tem1_ext(1,1,2), &
! tem1_ext(1,1,3),tem1_ext(1,1,4))
! DO ksmth=1,nsmooth
! CALL smooth9p(wetcanp(1,1,is), nx,ny, 1,nx,1,ny, tem1)
! END DO
END IF
END DO ! is
!wdt end
CALL fix_soil_nstyp(nx,ny,nstyp_ext,nstyp,tsfc,tsoil, &
wetsfc,wetdp,wetcanp)
IF ( snowdpth_ext(1,1) >= 0 ) THEN
snowdout = 1
!wdt
! CALL mkarps2d (nx_ext,ny_ext,nx,ny, &
! iorder,iscl,jscl,x_ext,y_ext, &
! xs2d,ys2d,snowdpth_ext,snowdpth, &
! dxfld,dyfld,rdxfld,rdyfld, &
! tem1_ext(1,1,1),tem1_ext(1,1,2), &
! tem1_ext(1,1,3),tem1_ext(1,1,4))
DO j=1,ny
DO i=1,nx
dmin=((xs2d(i,j)-x_ext(1))*(xs2d(i,j)-x_ext(1)) + &
(ys2d(i,j)-y_ext(1))*(ys2d(i,j)-y_ext(1)))
isnow=1
jsnow=1
DO jj=jscl(i,j),MAX(jscl(i,j)+1,ny_ext)
DO ii=iscl(i,j),MAX(iscl(i,j)+1,nx_ext)
dd=((xs2d(i,j)-x_ext(ii))*(xs2d(i,j)-x_ext(ii)) + &
(ys2d(i,j)-y_ext(jj))*(ys2d(i,j)-y_ext(jj)))
IF(dd < dmin) THEN
dmin=dd
isnow=ii
jsnow=jj
END IF
END DO
END DO
snowdpth(i,j) = snowdpth_ext(isnow,jsnow)
END DO
END DO
END IF
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
!wdt begin block:
xw = 0.
yw = 0.
DO j = 1,ny-1
DO i = 1,nx-1
xw(i,j) = &
MAX(0.,MIN(1.,(x_ext(iscl(i,j)+1)-xs2d(i,j)) &
/(x_ext(iscl(i,j)+1)-x_ext(iscl(i,j)))))
yw(i,j) = &
MAX(0.,MIN(1.,(y_ext(jscl(i,j)+1)-ys2d(i,j)) &
/(y_ext(jscl(i,j)+1)-y_ext(jscl(i,j)))))
END DO
END DO
!tmp wdt - for reference
!CALL intrp_soil(nx,ny,nx1,ny1,nstyp,nstyp,xw,yw,ix,jy, &
! tsfc,tsoil,wetsfc,wetdp,wetcanp, &
! soiltyp,stypfrct,vegtyp, &
! tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1, &
! soiltyp1,stypfrct1,vegtyp1)
CALL intrp_soil(nx_ext,ny_ext,nx,ny,nstyp_ext,nstyp,xw,yw,iscl,jscl, &
tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
soiltyp_ext,stypfrct_ext,vegtyp_ext, &
tsfc,tsoil,wetsfc,wetdp,wetcanp, &
soiltyp,stypfrct,vegtyp)
!wdt FIXME other data types ...; make soiltyp consistant with
! soil temps, etc.; move this section before soil temp/moist section
! INTEGER, allocatable :: soiltyp (:,:,:)! Soil type
! REAL, allocatable :: stypfrct(:,:,:)! Soil type fraction
! INTEGER, allocatable :: vegtyp (:,:) ! Vegetation type
! REAL, allocatable :: lai (:,:) ! Leaf Area Index
! REAL, allocatable :: roufns (:,:) ! Surface roughness
! REAL, allocatable :: veg (:,:) ! Vegetation fraction
!
! CALL fix_stypfrct_nstyp(nx,ny,nstyp_ext,nstyp,stypfrct)
!wdt FIXME
! Add capability to write out sfcdata file here? ...
!wdt end
IF ( tsfcout == 1 .OR. tsoilout == 1 .OR. &
wetsout == 1 .OR. wetdout == 1 .OR. &
wetcout == 1 .OR. snowdout == 1 ) THEN
CALL cvttsnd( curtim, timsnd, tmstrln )
soiloutfl = runname(1:lfnkey)//".soilvar."//timsnd(1:tmstrln)
lfn = lfnkey + 9 + tmstrln
IF( dirname /= ' ' ) THEN
temchar = soiloutfl
soiloutfl = dirname(1:ldirnam)//'/'//temchar
lfn = lfn + ldirnam + 1
END IF
IF (soildmp > 0) THEN
CALL fnversn(soiloutfl, lfn)
PRINT *, 'Write soil initial data to ',soiloutfl(1:lfn)
CALL wrtsoil(nx,ny,nstyp, soiloutfl,dx,dy, &
mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, &
tsfcout,tsoilout,wetsout, &
wetdout,wetcout,snowdout, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,soiltyp)
IF (soildmp == 1) CALL soilcntl(nx,ny, soiloutfl, &
tsfcout,tsoilout,wetsout, &
wetdout,wetcout,snowdout, &
x,y)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Test code, for diagnostic testing.
! Find x,y of diagnostic sounding location in ARPS grid.
!
!-----------------------------------------------------------------------
!
CALL lltoxy(1,1,latdiag,londiag,xdiag,ydiag)
!
dmin=((xdiag-xscl(1))*(xdiag-xscl(1))+ &
(ydiag-yscl(1))*(ydiag-yscl(1)))
idiag=1
jdiag=1
!
DO j=2,ny-2
DO i=2,nx-2
dd=((xdiag-xscl(i))*(xdiag-xscl(i))+ &
(ydiag-yscl(j))*(ydiag-yscl(j)))
IF(dd < dmin) THEN
dmin=dd
idiag=i
jdiag=j
END IF
END DO
END DO
CALL xytoll(1,1,xscl(idiag),yscl(jdiag), &
latd,lond)
WRITE(6,'(a,f10.4,f10.4,/a,i5,i5,a,f10.4,f10.4)') &
' Nearest ARPS pt to diagnostic lat,lon: ', &
latdiag,londiag, &
' Diagnostic i,j: ', &
idiag,jdiag,' lat,lon= ',latd,lond
WRITE(6,'(///a,/2x,a)') &
' ARPS extracted sounding at idiag,jdiag', &
'k pres hgt temp theta dewp u v dir spd'
!
!-----------------------------------------------------------------------
!
! Convert units of ARPS data and write as a sounding.
!
!-----------------------------------------------------------------------
!
DO k=nz-2,1,-1
ppasc=pbar(idiag,jdiag,k)+pprt(idiag,jdiag,k)
pmb=.01*(pbar(idiag,jdiag,k)+pprt(idiag,jdiag,k))
theta=ptbar(idiag,jdiag,k)+ptprt(idiag,jdiag,k)
tc=(theta*((ppasc/p0)**rddcp))-273.15
IF( qv(idiag,jdiag,k) > 0.) THEN
smix=qv(idiag,jdiag,k)/(1.-qv(idiag,jdiag,k))
e=(pmb*smix)/(0.62197 + smix)
bige=e/( 1.001 + ( (pmb - 100.) / 900.) * 0.0034)
alge = ALOG(bige/6.112)
tdc = (alge * 243.5) / (17.67 - alge)
ELSE
tdc = tc-30.
END IF
CALL uvrotdd(1,1,londiag, &
u(idiag,jdiag,k), &
v(idiag,jdiag,k), &
dir,spd)
WRITE(6,'(i4,f6.0,f9.0,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1)') &
k,pmb, &
zs(idiag,jdiag,k), &
tc,theta,tdc, &
u(idiag,jdiag,k), &
v(idiag,jdiag,k), &
dir,spd
END DO
!
!-----------------------------------------------------------------------
!
! Find vertical velocity and make any u-v adjustments
! according to wndadj option.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
rhostr(i,j,k)=ABS(j3(i,j,k))*rhobar(i,j,k)
END DO
END DO
END DO
IF(wndadj == 0) THEN
IF ( extdopt == 0 ) THEN ! ARPS history data
iprtopt=0
CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof, &
iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext, &
zp_ext,zp2_ext,xs2d,ys2d,zp,w_ext, &
zsnd,dumsnd,wbar,w, &
dxfld,dyfld,rdxfld,rdyfld, &
tem1_ext,tem2_ext,tem3_ext, &
tem4_ext)
DO ksmth=1,nsmooth
CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz, &
smfct1,zp,w,tem1,w)
END DO
ELSE
w = 0.
ENDIF
ELSE
CALL adjuvw( nx,ny,nz, u,v,w,wcont,ptprt,ptbar, &
zp,j1,j2,j3,aj3z,mapfct,rhostr,tem1, &
wndadj,obropt,obrzero,0, &
tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
END IF
!
!-----------------------------------------------------------------------
!
! Enforce vertical w boundary conditions
!
!-----------------------------------------------------------------------
!
IF (ext_lbc == 1 .or. ext_vbc == 1) &
CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
!
IF (ext_vbc == 1) CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v, &
rhostr,tem1,tem2,tem3,j1,j2,j3)
!
!
!-----------------------------------------------------------------------
!
! Assign zero-gradient horizontal boundary conditions
! to the horizontal & vertical winds.
!
!-----------------------------------------------------------------------
!
IF (ext_lbc == 1) THEN
CALL lbcw(nx,ny,nz,1.0, w,wcont,tem1,tem2,tem3,tem4,3,3,3,3)
CALL bcu(nx,ny,nz,1.0, u, tem1,tem2,tem3,tem4, 3,3,3,3,1,1)
CALL bcv(nx,ny,nz,1.0, v, tem1,tem2,tem3,tem4, 3,3,3,3,1,1)
ENDIF
!
!-----------------------------------------------------------------------
!
! Zero out uninitialized fields
!
!-----------------------------------------------------------------------
!
tem1=0.
!
!-----------------------------------------------------------------------
!
! Print out the max/min of output variables.
!
!-----------------------------------------------------------------------
!
WRITE(6,'(/1x,a/)') &
'Min and max of External data interpolated to the ARPS grid:'
CALL a3dmax0(x,1,nx,1,nx,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(/1x,2(a,e13.6))') &
'xmin = ', amin,', xmax =',amax
CALL a3dmax0(y,1,ny,1,ny,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'ymin = ', amin,', ymax =',amax
CALL a3dmax0(z,1,nz,1,nz,1,1,1,1,1,1,1,1, amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'zmin = ', amin,', zmax =',amax
CALL a3dmax0(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'zpmin = ', amin,', zpmax =',amax
CALL a3dmax0(rhobar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'rhobarmin=', amin,', rhobarmax=',amax
CALL a3dmax0(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'umin = ', amin,', umax =',amax
CALL a3dmax0(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'vmin = ', amin,', vmax =',amax
CALL a3dmax0(w,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'wmin = ', amin,', wmax =',amax
CALL a3dmax0(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'ptbarmin= ', amin,', ptbarmax=',amax
CALL a3dmax0(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'ptprtmin= ', amin,', ptprtmax=',amax
CALL a3dmax0(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'pbarmin= ', amin,', pbarmax =',amax
CALL a3dmax0(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'pprtmin = ', amin,', pprtmax =',amax
CALL a3dmax0(qv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1, &
amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'qvmin = ', amin,', qvmax =',amax
IF(sfcout == 1) THEN
CALL a3dmax0(tsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'tsfcmin = ', amin,', tsfcmax =',amax
CALL a3dmax0(tsoil,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'tsoilmin= ', amin,', tsoilmax=',amax
CALL a3dmax0(wetsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'wetsfcmin = ', amin,', wetsfcmax =',amax
CALL a3dmax0(wetdp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'wetdpmin= ', amin,', wetdpmax=',amax
CALL a3dmax0(wetcanp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
WRITE(6,'(1x,2(a,e13.6))') &
'wetcanpmin = ', amin,', wetcanpmax =',amax
END IF
!
!-----------------------------------------------------------------------
!
! Print the mean sounding that was used in setting the
! mean ARPS variables.
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a/a)') &
' Mean sounding for ARPS base state variables', &
' k p(mb) z(m) pt(mb) T(C) qv(kg/kg) ', &
' RH % u(m/s) v(m/s)'
DO k=lvlprof,1,-50
pres = psnd(k)
temp = ptsnd(k)*((pres/p0)**rddcp)
rh=AMAX1(0.01,rhmax-(rhssnd(k)*rhssnd(k)))
qvsat=f_qvsat( pres, temp )
qvval=rh*qvsat
WRITE(6,'(i4,f9.1,f9.1,f9.1,f9.1,f9.5,f9.1,f9.1,f9.1)') &
k,(0.01*psnd(k)),zsnd(k),ptsnd(k),(temp-273.15), &
qvval,(100.*rh),usnd(k),vsnd(k)
END DO
!
!-----------------------------------------------------------------------
!
! Data dump of the model grid and base state arrays:
!
! First find a unique name basdmpfn(1:lbasdmpf) for the grid and
! base state array dump file
!
!-----------------------------------------------------------------------
!
tem1=0.
ldirnam=LEN(dirname)
CALL strlnth( dirname, ldirnam)
IF ( hdmpfmt == 5 ) THEN
WRITE (6,'(a/a)') &
'Program ext2arps does not support Savi3D format.', &
'Reset to binary format, hdmpfmt=1.'
hdmpfmt = 1
END IF
IF ( hdmpfmt == 9 ) GO TO 450
CALL gtbasfn(runname(1:lfnkey),dirname,ldirnam,hdmpfmt, &
mgrid,nestgrd,basdmpfn,lbasdmpf)
PRINT*,'Writing base state history dump ',basdmpfn(1:lbasdmpf)
grdbas =1
mstout =1
rainout=0
prcout =0
trbout =0
IF (nstyp < 2) landout = 0
CALL dtadump(nx,ny,nz,nstyp, &
hdmpfmt,iniotfu,basdmpfn(1:lbasdmpf),grdbas,filcmprs, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
tem1,tem1,tem1, &
ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain,j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
tem1,tem1,tem1, &
tem1,tem1,tem1, &
tem1,tem1,tem1,tem1, &
tem2,tem3,tem4)
450 CONTINUE
PRINT*,'curtim = ', curtim
CALL gtdmpfn(runname(1:lfnkey),dirname,ldirnam,curtim,hdmpfmt, &
mgrid,nestgrd, hdmpfn, ldmpf)
WRITE(6,'(1x,a,a)') 'History data dump in file ',hdmpfn(1:ldmpf)
grdbas=0
CALL dtadump(nx,ny,nz,nstyp, &
hdmpfmt,iniotfu,hdmpfn(1:ldmpf),grdbas,filcmprs, &
u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
tem1,tem1,tem1, &
ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar, &
x,y,z,zp,hterain,j1,j2,j3, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth, &
tem1,tem1,tem1, &
tem1,tem1,tem1, &
tem1,tem1,tem1,tem1, &
tem2,tem3,tem4)
END DO
!
!-----------------------------------------------------------------------
!
! Friendly exit message.
!
!-----------------------------------------------------------------------
!
WRITE(6,'(a/a,i4,a)') &
' Normal succesful completion of EXT2ARPS.', &
' Processed',nextdfil,' file(s)'
STOP
!
!-----------------------------------------------------------------------
!
! Problem doing time conversions.
!
!-----------------------------------------------------------------------
!
920 CONTINUE
WRITE(6,'(a,/,a,i4,a,i4/,a,a19)') &
' Aborting, error in time format for external file', &
' File number:',ifile,' of',nextdfil, &
' External file time provided:',extdtime(ifile)
STOP
!
!-----------------------------------------------------------------------
!
! Error status returned from rdextfil
!
!-----------------------------------------------------------------------
!
999 CONTINUE
WRITE(6,'(a,i6)') &
' Aborting, error reading external file. istatus=', &
istatus
STOP
END PROGRAM ext2arps