PROGRAM arps2ncdf,70
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARSP2NCDF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Converts ARPS history files to a netCDF format file.
!
! Reads in a history file produced by ARPS in any ARPS format.
! Converts variables and interpolates to a fixed set of pressure
! levels. Sets up variables needed for LDADS/AWIPS D2d display
! of data.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster
! 02/19/1999
!
! MODIFICATION HISTORY:
! 08/13/1999
! Fixed MSLP calculation (J. Levit, K. Brewster)
!
!
!-----------------------------------------------------------------------
!
! DATA ARRAYS READ IN:
!
! 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 z coordinate of grid points in physical space (m)
!
! w vertical component of velocity in Cartesian
! coordinates (m/s).
!
! ptprt perturbation potential temperature (K)
! pprt perturbation pressure (Pascal)
! uprt perturbation x velocity component (m/s)
! vprt perturbation y velocity component (m/s)
! wprt perturbation z velocity component (m/s)
!
! qv water vapor mixing ratio (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
!
! ubar Base state x velocity component (m/s)
! vbar Base state y velocity component (m/s)
! wbar Base state z 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 mixing ratio (kg/kg)
!
! soiltyp Soil type
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsfc Temperature at surface (K)
! tsoil Deep soil temperature (K)
! wetsfc Surface soil moisture
! wetdp Deep soil moisture
! wetcanp Canopy water amount
!
! rain Total rain (raing + rainc)
! raing Grid supersaturation rain
! rainc Cumulus convective rain
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Grid dimensions.
INTEGER :: nstyps ! Maximum number of soil types.
INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,nf,lenfil
PARAMETER (nhisfile_max=200)
CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.inc'
INCLUDE 'arps2ncdf.inc'
INCLUDE 'netcdf.inc'
!
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x(:) ! The x-coord. of the physical and
! computational grid.
! Defined at u-point.
REAL, ALLOCATABLE :: y(:) ! The y-coord. of the physical and
! computational grid.
! Defined at v-point.
REAL, ALLOCATABLE :: z(:) ! The z-coord. of the computational
! grid. Defined at w-point.
REAL, ALLOCATABLE :: zp(:,:,:) ! The height of the terrain.
REAL, ALLOCATABLE :: uprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: wprt (:,:,:) ! Perturbation w-velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal)
REAL, ALLOCATABLE :: qvprt (:,:,:) ! Perturbation 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 :: qv (:,:,:) ! Water vapor specific humidity (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 )
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
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! 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 :: tsfc (:,:,:) ! Temperature at surface (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 :: rain (:,:) ! Total rainfall
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 :: e_mb (:,:,:) ! vapor pressure in mb
REAL, ALLOCATABLE :: mix (:,:,:) ! total mixing ratio (kg/kg)
REAL, ALLOCATABLE :: esat_mb(:,:,:) ! saturation vapor pressure in mb
REAL, ALLOCATABLE :: rh (:,:,:) ! Relative humidity in %
REAL, ALLOCATABLE :: t_dew (:,:,:) ! dewpoint temp. in degrees K
!
!-----------------------------------------------------------------------
!
! 2-D stability index arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: lcl(:,:) ! lifting condensation level
REAL, ALLOCATABLE :: lfc(:,:) ! level of free convection
REAL, ALLOCATABLE :: el(:,:) ! equilibrium level
REAL, ALLOCATABLE :: twdf(:,:) ! max. wet bulb pot. temp. difference
REAL, ALLOCATABLE :: li(:,:) ! lifted index
REAL, ALLOCATABLE :: pbe(:,:) ! CAPE
REAL, ALLOCATABLE :: mbe(:,:) ! Moist CAPE
REAL, ALLOCATABLE :: nbe(:,:) ! CIN
REAL, ALLOCATABLE :: tcap(:,:) ! Cap Strength
!
REAL, ALLOCATABLE :: heli(:,:) ! helicity
REAL, ALLOCATABLE :: brn(:,:) ! Bulk Richardson Number (Weisman and Klemp)
REAL, ALLOCATABLE :: brnu(:,:) ! Shear parameter of BRN, "U"
REAL, ALLOCATABLE :: srlfl(:,:) ! storm-relative low-level flow (0-2km AGL)
REAL, ALLOCATABLE :: srmfl(:,:) ! storm-relative mid-level flow (2-9km AGL)
REAL, ALLOCATABLE :: shr37(:,:) ! 7km - 3km wind shear
REAL, ALLOCATABLE :: ustrm(:,:) ! Estimated storm motion (modified Bob Johns)
REAL, ALLOCATABLE :: vstrm(:,:) ! Estimated storm motion (modified Bob Johns)
REAL, ALLOCATABLE :: blcon(:,:) ! Boundary layer convergence
!
!-----------------------------------------------------------------------
!
! Work Arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
REAL, ALLOCATABLE :: tem4(:,:,:)
REAL, ALLOCATABLE :: tem2d1(:,:)
REAL, ALLOCATABLE :: tem2d2(:,:)
REAL, ALLOCATABLE :: tem2d3(:,:)
REAL, ALLOCATABLE :: tem2d4(:,:)
REAL, ALLOCATABLE :: wrk1(:),wrk2(:),wrk3(:),wrk4(:),wrk5(:),wrk6(:)
REAL, ALLOCATABLE :: wrk7(:),wrk8(:),wrk9(:),wrk10(:),wrk11(:),wrk12(:)
REAL, ALLOCATABLE :: xs(:)
REAL, ALLOCATABLE :: ys(:)
REAL, ALLOCATABLE :: zps(:,:,:)
REAL, ALLOCATABLE :: zps_km(:,:,:)
REAL, ALLOCATABLE :: mf2d(:,:)
REAL, ALLOCATABLE :: coriolis(:,:)
REAL, ALLOCATABLE :: spacing(:,:)
!
!-----------------------------------------------------------------------
!
! Misc ARPS variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nch
REAL :: time
REAL :: latnot(2)
LOGICAL :: init
!
!-----------------------------------------------------------------------
!
! netCDF dimensions
!
! nprlvl is difined in arps2cdf.inc
!
!-----------------------------------------------------------------------
!
INTEGER :: nxcdf,nycdf,mxtime
PARAMETER(mxtime=40)
INTEGER :: n3dvar,n2dvar,nvartot,nstvar
PARAMETER (n3dvar=6, &
n2dvar=28, &
nvartot=(n3dvar+n2dvar), &
nstvar=3)
!
!-----------------------------------------------------------------------
!
! netCDF variables
!
!-----------------------------------------------------------------------
!
INTEGER :: istatus,ncid
INTEGER :: dimidrec,dimidtot,dimidnt,dimidnav
INTEGER :: dimidnx,dimidny,dimid1,dimidnz,dimchrlvl,dimnamlen
INTEGER :: vtmrefid,valtimid,reftimid,origid,modelid
!
INTEGER :: vecistr(2)
INTEGER :: planeistr(4)
INTEGER :: volistart(4)
DATA vecistr /1,1/
DATA planeistr /1,1,1,1/
DATA volistart /1,1,1,1/
!
INTEGER :: planedim(4)
INTEGER :: voldim(4)
INTEGER :: lvldim(2)
INTEGER :: invdim(2)
INTEGER :: planelen(4)
INTEGER :: lvllen(2)
INTEGER :: sfclen(2)
INTEGER :: vollen(4)
INTEGER :: invlen(2)
!
REAL :: vrange(2)
!
!-----------------------------------------------------------------------
!
! netCDF output variables
!
!-----------------------------------------------------------------------
!
INTEGER :: v3did(n3dvar)
INTEGER :: v3dlvlid(n3dvar)
INTEGER :: v2dinvid(n2dvar)
INTEGER :: v3dinvid(n3dvar)
INTEGER :: v2did(n2dvar)
INTEGER :: vstid(nstvar)
INTEGER :: v2dlvlid(n2dvar)
INTEGER :: vtmref(mxtime)
REAL*8 valtime(mxtime)
REAL*8 reftime(mxtime)
INTEGER :: ntime
!
CHARACTER (LEN=4) :: v3dnam(n3dvar)
CHARACTER (LEN=4) :: v2dnam(n2dvar)
CHARACTER (LEN=14) :: vstnam(nstvar)
CHARACTER (LEN=10) :: v3dlvl(nprlvl)
CHARACTER (LEN=10) :: v2dlvl
CHARACTER (LEN=30) :: v3dlngnam(n3dvar)
CHARACTER (LEN=30) :: v2dlngnam(n2dvar)
CHARACTER (LEN=30) :: vstlngnam(nstvar)
CHARACTER (LEN=30) :: v3dunits(n3dvar)
CHARACTER (LEN=30) :: v2dunits(n2dvar)
CHARACTER (LEN=30) :: vstunits(nstvar)
CHARACTER (LEN=nprlvl) :: inventory
CHARACTER (LEN=1) :: inventory_sfc
REAL :: v3dmin(n3dvar)
REAL :: v2dmin(n2dvar)
REAL :: v3dmax(n3dvar)
REAL :: v2dmax(n2dvar)
!
CHARACTER (LEN=10) :: v2dlvlnam
CHARACTER (LEN=10) :: v3dlvlnam
CHARACTER (LEN=13) :: v2dinv
CHARACTER (LEN=13) :: v3dinv
!
DATA vtmref /mxtime*0/
!
!-----------------------------------------------------------------------
!
! Variable indexes and descriptors
!
!-----------------------------------------------------------------------
!
INTEGER :: idxgh,idxrh,idxt,idxuw,idxvw,idxww
PARAMETER (idxgh=1,idxrh=2,idxt=3, &
idxuw=4,idxvw=5,idxww=6)
!
DATA v3dnam /'gh','rh','t','uw','vw','ww'/
DATA v3dlngnam /'Geopotential Height', &
'Relative Humidity', &
'Temperature', &
'u Wind Component', &
'v Wind Component', &
'w Wind Component'/
DATA v3dunits /'geopotential meters', &
'percent', &
'degrees K', &
'meters/second', &
'meters/second', &
'meters/second'/
DATA v3dmin / 0., 0.,150.,-300.,-300.,-100./
DATA v3dmax /40000.,110.,400., 300., 300., 100./
!
INTEGER :: idxmslp,idxp,idxlgsp,idxcp,idxtp
INTEGER :: idxsfu,idxsfv,idxsft,idxdpt
INTEGER :: idxsh,idxept,idxsli,idxcape,idxcin
INTEGER :: idxlcl,idxlfc,idxel,idxtwdf,idxtcap
INTEGER :: idxshr37,idxustrm,idxvstrm
INTEGER :: idxsrlfl,idxsrmfl,idxheli,idxbrn,idxbrnu,idxblcon
PARAMETER (idxmslp=1, &
idxp =2, &
idxlgsp=3, &
idxcp =4, &
idxtp =5, &
idxsfu =6, &
idxsfv =7, &
idxsft =8, &
idxdpt =9, &
idxsh =10, &
idxept =11, &
idxsli =12, &
idxcape=13, &
idxcin =14, &
idxlcl =15, &
idxlfc =16, &
idxel =17, &
idxtwdf =18, &
idxtcap =19, &
idxshr37=20, &
idxustrm=21, &
idxvstrm=22, &
idxsrlfl=23, &
idxsrmfl=24, &
idxheli =25, &
idxbrn =26, &
idxbrnu =27, &
idxblcon=28)
DATA v2dnam /'mslp','p','lgsp','cp','tp','sfu','sfv', &
'sft','dpt','sh','ept','sli','cape','cin','lcl', &
'lfc','el','twdf','tcap','shr37','ustrm','vstrm', &
'srlfl','srmfl','heli','brn','brnu','blcon'/
DATA v2dlngnam /'Mean Sea Level Pressure', &
'Surface Pressure', &
'Grid Scale Precipitation', &
'Convective Precipitation', &
'Total Precipitation', &
'Surface u Wind Component', &
'Surface v Wind Component', &
'Surface Temperature', &
'Surface Dew Point', &
'Specific Humidity', &
'Theta-e', &
'Surface Lifted Index', &
'CAPE', &
'CIN', &
'Lifting Condensation Lvl', &
'Level of Free Convection', &
'Equilibrium Level', &
'Max Wet-Bulb Temp Diff', &
'Cap Strength', &
'3-7km Shear', &
'Est Storm Motion u comp', &
'Est Storm Motion v comp', &
'Storm Rel Low-Level Flow', &
'Storm Rel Mid-Level Flow', &
'Storm Rel Helicity', &
'Bulk Richardson Number', &
'BRN Shear u', &
'Bound-Layer Conv x 1000'/
DATA v2dunits /'Pascals','Pascals','mm','mm','mm', &
'meters/second','meters/second', &
'degrees K','degrees K','kg/kg', &
'degrees K','degrees C','J/kg','J/kg', &
'meters','meters','meters','degrees C','degrees C', &
'1/second','meters/second','meters/second', &
'meters/second','meters/second','m2/s2', &
'unitless','meters/second','1/second'/
DATA v2dmin / 80000.,50000.,0.,0.,0., &
-100.,-100., &
-100.,-100.,0., &
200.,-100.,0.,0., &
0.,0.,0.,-100.,0.,0.,-100.,-100., &
-100.,-100.,-900.,0.,0.,-10./
DATA v2dmax / 110000.,120000.,1000.,1000.,1000., &
100.,100., &
100.,100.,1., &
500.,100.,10000.,10000., &
40000.,40000.,40000.,100.,100.,1.,100.,100., &
100.,100.,900.,900.,100.,10./
!
INTEGER :: idxtop,idxcor,idxspa
PARAMETER (idxtop=1,idxcor=2,idxspa=3)
DATA vstnam /'staticTopo', &
'staticCoriolis', &
'staticSpacing'/
DATA vstlngnam /'Topography', &
'Coriolis parameter', &
'Grid spacing'/
DATA vstunits /'meters', &
'/second', &
'meters'/
!
CHARACTER (LEN=20) :: cproj
!
! The following are the projections supported by LDAD/AWIPS
! in their ordering/naming convention
!
CHARACTER (LEN=20) :: chproj(17)
DATA chproj /'STEREOGRAPHIC', &
'ORTHOGRAPHIC', &
'LAMBERT_CONFORMAL', &
'AZIMUTHAL_EQUAL_AREA', &
'GNOMONIC', &
'AZIMUTHAL_EQUIDISTANT', &
'SATELLITE_VIEW', &
'CYLINDRICAL_EQUIDISTANT', &
'MERCATOR', &
'MOLLWEIDE', &
'NULL', &
'NULL', &
'NULL', &
'NULL', &
'NULL', &
'AZIMUTH_RANGE', &
'RADAR_AZ_RAN'/
!
! kproj is the index in the LDAD/AWIPS map projection list
! for the mapproj number in the ARPS file
!
INTEGER :: kproj(3)
DATA kproj /1,3,9/
!
REAL, ALLOCATABLE :: static2d(:,:,:)
REAL, ALLOCATABLE :: cdf2d(:,:,:,:)
REAL, ALLOCATABLE :: cdf3d(:,:,:,:)
!
!-----------------------------------------------------------------------
!
! Constants
!
!-----------------------------------------------------------------------
!
INTEGER :: itim1970
PARAMETER (itim1970=315619200)
CHARACTER (LEN=8) :: cdldate
CHARACTER (LEN=40) :: depictor
CHARACTER (LEN=132) :: origin,model
PARAMETER (cdldate='19990301', &
depictor='TEST_ARPS', &
origin='HubCAPS', &
model='ARPS 4.4.1')
!
!-----------------------------------------------------------------------
!
! Namelists
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=132) :: cdfname
INTEGER :: ipktyp,nbits
NAMELIST /ncdf_options/ cdfname,ipktyp,nbits
!
!-----------------------------------------------------------------------
!
! External functions
!
!-----------------------------------------------------------------------
!
REAL :: wmr2td,oe,dpt
EXTERNAL wmr2td,oe,dpt
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=10) :: nzstr
REAL :: latll,lonll,latur,lonur,dxkm,dykm,latdxdy,londxdy,rotate
INTEGER :: ivar,ls,itime,itimref,iproj
!
INTEGER :: i,j,k,klev,ifile,ireturn,imid,jmid
REAL :: rlnpcdf,r2d
REAL :: ctrx,ctry,swx,swy
REAL :: gamma,ex1,ex2,rln700,p00,t00
REAL :: prmb,tdew
REAL, ALLOCATABLE :: p_mb(:,:,:)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Initializations
!
!-----------------------------------------------------------------------
!
!
WRITE(6,'(6(/a))') &
'###############################################################',&
'# #',&
'# Welcome to ARPS2NCDF, a program that reads in history files #',&
'# generated by ARPS and produces a netCDF format file. #',&
'# #',&
'###############################################################'
!
!-----------------------------------------------------------------------
!
! Get the names of the input data files.
!
!-----------------------------------------------------------------------
!
CALL get_input_file_names
(hinfmt,grdbasfn,hisfile,nhisfile)
lengbf = len_trim(grdbasfn)
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nstyps, ireturn)
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz
nxcdf=(nx-3)
nycdf=(ny-3)
planelen(1)=nxcdf
planelen(2)=nycdf
planelen(3)=1
planelen(4)=1
lvllen(1)=10
lvllen(2)=nprlvl
sfclen(1)=10
sfclen(2)=1
vollen(1)=nxcdf
vollen(2)=nycdf
vollen(3)=nprlvl
vollen(4)=1
!
!-----------------------------------------------------------------------
!
! Get output options and the name of the grid/base data set.
!
!-----------------------------------------------------------------------
!
ipktyp=1
nbits=16
READ(5,ncdf_options,END=900)
ntime=nhisfile
lengbf=LEN_trim(grdbasfn)
WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)
ALLOCATE(x (nx))
ALLOCATE(y (ny))
ALLOCATE(z (nz))
ALLOCATE(zp (nx,ny,nz))
ALLOCATE(uprt (nx,ny,nz))
ALLOCATE(vprt (nx,ny,nz))
ALLOCATE(wprt (nx,ny,nz))
ALLOCATE(ptprt (nx,ny,nz))
ALLOCATE(pprt (nx,ny,nz))
ALLOCATE(qvprt (nx,ny,nz))
ALLOCATE(qv (nx,ny,nz))
ALLOCATE(qc (nx,ny,nz))
ALLOCATE(qr (nx,ny,nz))
ALLOCATE(qi (nx,ny,nz))
ALLOCATE(qs (nx,ny,nz))
ALLOCATE(qh (nx,ny,nz))
ALLOCATE(tke (nx,ny,nz))
ALLOCATE(kmh (nx,ny,nz))
ALLOCATE(kmv (nx,ny,nz))
ALLOCATE(ubar (nx,ny,nz))
ALLOCATE(vbar (nx,ny,nz))
ALLOCATE(wbar (nx,ny,nz))
ALLOCATE(ptbar (nx,ny,nz))
ALLOCATE(pbar (nx,ny,nz))
ALLOCATE(rhobar (nx,ny,nz))
ALLOCATE(qvbar (nx,ny,nz))
ALLOCATE(soiltyp (nx,ny,nstyps))
ALLOCATE(stypfrct(nx,ny,nstyps))
ALLOCATE(vegtyp (nx,ny))
ALLOCATE(lai (nx,ny))
ALLOCATE(roufns (nx,ny))
ALLOCATE(veg (nx,ny))
ALLOCATE(tsfc (nx,ny,0:nstyps))
ALLOCATE(tsoil (nx,ny,0:nstyps))
ALLOCATE(wetsfc (nx,ny,0:nstyps))
ALLOCATE(wetdp (nx,ny,0:nstyps))
ALLOCATE(wetcanp(nx,ny,0:nstyps))
ALLOCATE(snowdpth(nx,ny))
ALLOCATE(rain (nx,ny))
ALLOCATE(raing(nx,ny))
ALLOCATE(rainc(nx,ny))
ALLOCATE(prcrate(nx,ny,4))
ALLOCATE(radfrc(nx,ny,nz))
ALLOCATE(radsw (nx,ny))
ALLOCATE(rnflx (nx,ny))
ALLOCATE(usflx (nx,ny))
ALLOCATE(vsflx (nx,ny))
ALLOCATE(ptsflx(nx,ny))
ALLOCATE(qvsflx(nx,ny))
ALLOCATE(e_mb (nx,ny,nz))
ALLOCATE(mix (nx,ny,nz))
ALLOCATE(esat_mb(nx,ny,nz))
ALLOCATE(rh (nx,ny,nz))
ALLOCATE(t_dew (nx,ny,nz))
ALLOCATE(lcl(nx,ny))
ALLOCATE(lfc(nx,ny))
ALLOCATE(el(nx,ny))
ALLOCATE(twdf(nx,ny))
ALLOCATE(li(nx,ny))
ALLOCATE(pbe(nx,ny))
ALLOCATE(mbe(nx,ny))
ALLOCATE(nbe(nx,ny))
ALLOCATE(tcap(nx,ny))
ALLOCATE(tem1(nx,ny,nz))
ALLOCATE(tem2(nx,ny,nz))
ALLOCATE(tem3(nx,ny,nz))
ALLOCATE(tem4(nx,ny,nz))
ALLOCATE(tem2d1(nx,ny))
ALLOCATE(tem2d2(nx,ny))
ALLOCATE(tem2d3(nx,ny))
ALLOCATE(tem2d4(nx,ny))
ALLOCATE(wrk1(nz),wrk2(nz),wrk3(nz),wrk4(nz),wrk5(nz),wrk6(nz))
ALLOCATE(wrk7(nz),wrk8(nz),wrk9(nz),wrk10(nz),wrk11(nz),wrk12(nz))
ALLOCATE(xs(nx))
ALLOCATE(ys(ny))
ALLOCATE(zps(nx,ny,nz))
ALLOCATE(zps_km(nx,ny,nz))
ALLOCATE(p_mb(nx,ny,nz))
ALLOCATE(heli (nx,ny))
ALLOCATE(brn (nx,ny))
ALLOCATE(brnu (nx,ny))
ALLOCATE(srlfl(nx,ny))
ALLOCATE(srmfl(nx,ny))
ALLOCATE(shr37(nx,ny))
ALLOCATE(ustrm(nx,ny))
ALLOCATE(vstrm(nx,ny))
ALLOCATE(blcon(nx,ny))
ALLOCATE(mf2d(nx,ny))
ALLOCATE(coriolis(nx,ny))
ALLOCATE(spacing(nx,ny))
ALLOCATE(static2d(nxcdf,nycdf,nstvar))
ALLOCATE(cdf2d(nxcdf,nycdf,1,n2dvar))
ALLOCATE(cdf3d(nxcdf,nycdf,nprlvl,n3dvar))
x =0.0
y =0.0
z =0.0
zp =0.0
uprt =0.0
vprt =0.0
wprt =0.0
ptprt =0.0
pprt =0.0
qvprt =0.0
qc =0.0
qr =0.0
qi =0.0
qs =0.0
qh =0.0
tke =0.0
kmh =0.0
kmv =0.0
ubar =0.0
vbar =0.0
wbar =0.0
ptbar =0.0
pbar =0.0
rhobar =0.0
qvbar =0.0
qv =0.0
soiltyp =0.0
stypfrct=0.0
vegtyp =0.0
lai =0.0
roufns =0.0
veg =0.0
tsfc =0.0
tsoil =0.0
wetsfc =0.0
wetdp =0.0
wetcanp=0.0
snowdpth=0.0
rain =0.0
raing=0.0
rainc=0.0
prcrate=0.0
radfrc=0.0
radsw =0.0
rnflx =0.0
usflx =0.0
vsflx =0.0
ptsflx=0.0
qvsflx=0.0
e_mb =0.0
mix =0.0
esat_mb =0.0
rh =0.0
t_dew =0.0
lcl =0.0
lfc =0.0
el =0.0
twdf=0.0
li =0.0
pbe =0.0
mbe =0.0
nbe =0.0
tcap=0.0
tem1 =0.0
tem2 =0.0
tem3 =0.0
tem4 =0.0
tem2d1=0.0
tem2d2=0.0
tem2d3=0.0
tem2d4=0.0
wrk1 = 0.0
wrk2 = 0.0
wrk3 = 0.0
wrk4 = 0.0
wrk5 = 0.0
wrk6 = 0.0
wrk7 = 0.0
wrk8 = 0.0
wrk9 = 0.0
wrk10= 0.0
wrk11= 0.0
wrk12= 0.0
xs=0.0
ys=0.0
zps=0.0
zps_km=0.0
p_mb =0.0
heli =0.0
brn =0.0
brnu =0.0
srlfl=0.0
srmfl=0.0
shr37=0.0
ustrm=0.0
vstrm=0.0
blcon=0.0
mf2d=0.0
coriolis=0.0
spacing=0.0
init=.false.
DO ifile=1,nhisfile
lenfil = LEN(hisfile(ifile))
CALL strlnth( hisfile(ifile), lenfil)
WRITE(6,'(/a,/1x,a)')' The data set name is ', &
hisfile(ifile)(1:lenfil)
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
CALL dtaread(nx,ny,nz,nstyps,hinfmt,nch,grdbasfn(1:lengbf), &
lengbf, &
hisfile(ifile)(1:lenfil),lenfil,time, &
x,y,z,zp, uprt ,vprt ,wprt ,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)
!
!-----------------------------------------------------------------------
!
! ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
IF( ireturn == 0 ) THEN ! successful read
IF( .NOT.init ) THEN
!
!-----------------------------------------------------------------------
!
! Establish coordinate for scalar fields.
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
xs(i)=0.5*(x(i)+x(i+1))
END DO
xs(nx)=2.*xs(nx-1)-xs(nx-2)
DO j=1,ny-1
ys(j)=0.5*(y(j)+y(j+1))
END DO
ys(ny)=2.*ys(ny-1)-ys(ny-2)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
zps(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
END DO
END DO
END DO
! nf_noerr = 0 ! correct??
!
!-----------------------------------------------------------------------
!
! Open netCDF file
!
!-----------------------------------------------------------------------
!
istatus=nf_create(cdfname,nf_clobber,ncid)
IF (istatus /= nf_noerr) CALL handle_err('nf_create',istatus)
!
!-----------------------------------------------------------------------
!
! Define dimensions
!
!-----------------------------------------------------------------------
!
IF(nprlvl < 10) THEN
WRITE(nzstr,'(a,i1)') 'levels_',nprlvl
ELSE IF(nprlvl < 100) THEN
WRITE(nzstr,'(a,i2)') 'levels_',nprlvl
ELSE
WRITE(nzstr,'(a,i3)') 'levels_',nprlvl
END IF
!
istatus=nf_def_dim(ncid,'record',nf_unlimited,dimidrec)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim record',istatus)
istatus=nf_def_dim(ncid,'n_valtimes',ntime,dimidnt)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim ntime',istatus)
istatus=nf_def_dim(ncid,'data_variables',nvartot,dimidtot)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim nvartot',istatus)
istatus=nf_def_dim(ncid,'charsPerLevel',10,dimchrlvl)
IF (istatus /= nf_noerr) &
CALL handle_err('nf_def_dim charsPerLevel',istatus)
istatus=nf_def_dim(ncid,'namelen',132,dimnamlen)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim namelen',istatus)
istatus=nf_def_dim(ncid,'x',nxcdf,dimidnx)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim nx',istatus)
istatus=nf_def_dim(ncid,'y',nycdf,dimidny)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim ny',istatus)
istatus=nf_def_dim(ncid,'levels_1',1,dimid1)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim 1',istatus)
istatus=nf_def_dim(ncid,nzstr,nprlvl,dimidnz)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim levels',istatus)
istatus=nf_def_dim(ncid,'nav',1,dimidnav)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim nav',istatus)
!
planedim(1)=dimidnx
planedim(2)=dimidny
planedim(3)=dimid1
planedim(4)=dimidrec
voldim(1)=dimidnx
voldim(2)=dimidny
voldim(3)=dimidnz
voldim(4)=dimidrec
lvldim(1)=dimchrlvl
lvldim(2)=dimidnz
invdim(1)=dimidnz
invdim(2)=dimidrec
!
!-----------------------------------------------------------------------
!
! Define 3-d variables and attributes
!
!-----------------------------------------------------------------------
!
DO ivar=1,n3dvar
vrange(1)=v3dmin(ivar)
vrange(2)=v3dmax(ivar)
istatus=nf_def_var(ncid,v3dnam(ivar),nf_float,4, &
voldim,v3did(ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var 3d',istatus)
ls=LEN(v3dlngnam(ivar))
CALL strlnth(v3dlngnam(ivar),ls)
istatus=nf_put_att_text(ncid,v3did(ivar),'long_name', &
ls,v3dlngnam(ivar)(1:ls))
ls=LEN(v3dunits(ivar))
CALL strlnth(v3dunits(ivar),ls)
istatus=nf_put_att_text(ncid,v3did(ivar),'units', &
ls,v3dunits(ivar)(1:ls))
istatus=nf_put_att_real(ncid,v3did(ivar),'valid_range', &
nf_float,2,vrange)
istatus=nf_put_att_real(ncid,v3did(ivar),'_FillValue', &
nf_float,1,-99999.)
istatus=nf_put_att_int(ncid,v3did(ivar),'_n3D', &
nf_int,1,nprlvl)
istatus=nf_put_att_text(ncid,v3did(ivar),'levels', &
2,'MB')
!
ls=LEN(v3dnam(ivar))
CALL strlnth(v3dnam(ivar),ls)
WRITE(v3dlvlnam,'(a,a)') v3dnam(ivar)(1:ls),'Levels'
istatus=nf_def_var(ncid,v3dlvlnam,nf_char,2, &
lvldim,v3dlvlid(ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var lvl 3d',istatus)
!
WRITE(v3dinv,'(a,a)') v3dnam(ivar)(1:ls),'Inventory'
istatus=nf_def_var(ncid,v3dinv,nf_char,2, &
invdim,v3dinvid(ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var inv 3d',istatus)
!
END DO
!
!-----------------------------------------------------------------------
!
! Define 2-d variables and attributes
!
!-----------------------------------------------------------------------
!
invdim(1)=dimid1
invdim(2)=dimidrec
lvldim(2)=dimid1
DO ivar=1,n2dvar
vrange(1)=v2dmin(ivar)
vrange(2)=v2dmax(ivar)
istatus=nf_def_var(ncid,v2dnam(ivar),nf_float,4, &
planedim,v2did(ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var 2d',istatus)
ls=LEN(v2dlngnam(ivar))
CALL strlnth(v2dlngnam(ivar),ls)
istatus=nf_put_att_text(ncid,v2did(ivar),'long_name', &
ls,v2dlngnam(ivar)(1:ls))
ls=LEN(v2dunits(ivar))
CALL strlnth(v2dunits(ivar),ls)
istatus=nf_put_att_text(ncid,v2did(ivar),'units', &
ls,v2dunits(ivar)(1:ls))
istatus=nf_put_att_real(ncid,v2did(ivar),'valid_range', &
nf_float,2,vrange)
istatus=nf_put_att_real(ncid,v2did(ivar),'_FillValue', &
nf_float,1,-99999.)
istatus=nf_put_att_int(ncid,v2did(ivar),'_n3D', &
nf_int,1,0)
IF(ivar == idxmslp) THEN
istatus=nf_put_att_text(ncid,v2did(ivar),'levels', &
3,'MSL')
ELSE
istatus=nf_put_att_text(ncid,v2did(ivar),'levels', &
3,'SFC')
END IF
!
!
ls=LEN(v2dnam(ivar))
CALL strlnth(v2dnam(ivar),ls)
WRITE(v2dlvlnam,'(a,a)') v2dnam(ivar)(1:ls),'Levels'
istatus=nf_def_var(ncid,v2dlvlnam,nf_char,2, &
lvldim,v2dlvlid(ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var lvl 2d',istatus)
!
WRITE(v2dinv,'(a,a)') v2dnam(ivar)(1:ls),'Inventory'
istatus=nf_def_var(ncid,v2dinv,nf_char,2, &
invdim,v2dinvid(ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var inv 2d',istatus)
END DO
!
!-----------------------------------------------------------------------
!
! Define 1-d variables
!
!-----------------------------------------------------------------------
!
istatus=nf_def_var(ncid,'valtimeMINUSreftime',nf_int,1, &
dimidnt,vtmrefid)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var vtmref',istatus)
istatus=nf_put_att_text(ncid,vtmrefid,'units', &
7,'seconds')
istatus=nf_def_var(ncid,'reftime',nf_double,1, &
dimidnt,reftimid)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var reftim',istatus)
istatus=nf_put_att_text(ncid,reftimid,'units', &
33,'seconds since 1970-1-1 00:00:00.0')
istatus=nf_def_var(ncid,'valtime',nf_double,1, &
dimidnt,valtimid)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var valtime',istatus)
istatus=nf_put_att_text(ncid,vtmrefid,'units', &
33,'seconds since 1970-1-1 00:00:00.0')
istatus=nf_def_var(ncid,'origin',nf_char,1, &
dimnamlen,origid)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var orig',istatus)
istatus=nf_def_var(ncid,'model',nf_char,1, &
dimnamlen,modelid)
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var model',istatus)
!
! Define static variables
!
DO ivar=1,nstvar
istatus=nf_def_var(ncid,vstnam(ivar),nf_float,2, &
planedim,vstid(ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_def_var static',istatus)
ls=LEN(vstlngnam(ivar))
CALL strlnth(vstlngnam(ivar),ls)
istatus=nf_put_att_text(ncid,vstid(ivar),'long_name', &
ls,vstlngnam(ivar)(1:ls))
ls=LEN(vstunits(ivar))
CALL strlnth(vstunits(ivar),ls)
istatus=nf_put_att_text(ncid,vstid(ivar),'units', &
ls,vstunits(ivar)(1:ls))
END DO
!
!-----------------------------------------------------------------------
!
! Global Attributes
!
!-----------------------------------------------------------------------
!
latnot(1)=trulat1
latnot(2)=trulat2
CALL setmapr( mapproj, 1.0 , latnot , trulon)
CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
dx=x(3)-x(2)
dy=y(3)-y(2)
dxkm=0.001*dx
dykm=0.001*dy
swx = ctrx - (FLOAT(nx-3)/2.) * dx
swy = ctry - (FLOAT(ny-3)/2.) * dy
CALL setorig( 1, swx, swy)
!
CALL xytoll(1,1,xs(2),ys(2),latll,lonll)
CALL xytoll(1,1,xs(nx-2),ys(ny-2),latur,lonur)
!
IF(mapproj > 0) THEN
iproj=kproj(mapproj)
cproj=chproj(iproj)
ELSE
iproj=0
cproj='NONE'
END IF
!
IF(mapproj == 2) THEN
latdxdy=0.5*(trulat1+trulat2)
ELSE
latdxdy=trulat1
END IF
londxdy=ctrlon
CALL xytoll(nx,ny,xs,ys,coriolis,tem1)
CALL xytomf(nx,ny,xs,ys,mf2d)
r2d = ATAN(1.0)*4.0/180.0
DO j=1,ny-1
DO i=1,nx-1
coriolis(i,j) = 2.*omega*SIN( r2d * coriolis(i,j) )
spacing(i,j) = dx*mf2d(i,j)
END DO
END DO
!
WRITE(6,'(a,a)') ' cdl date ',cdldate
!
istatus=nf_put_att_text(ncid,nf_global,'cdlDate', &
8,cdldate)
ls=LEN(depictor)
CALL strlnth(depictor,ls)
istatus=nf_put_att_text(ncid,nf_global,'depictorName', &
ls,depictor)
istatus=nf_put_att_int(ncid,nf_global,'projIndex', &
nf_int,1,iproj)
ls=LEN(cproj)
CALL strlnth(cproj,ls)
istatus=nf_put_att_text(ncid,nf_global,'projName', &
ls,cproj)
istatus=nf_put_att_real(ncid,nf_global,'centralLat', &
nf_float,1,ctrlat)
istatus=nf_put_att_real(ncid,nf_global,'centralLon', &
nf_float,1,ctrlon)
rotate=ctrlon-trulon
istatus=nf_put_att_real(ncid,nf_global,'rotation', &
nf_float,1,rotate)
istatus=nf_put_att_real(ncid,nf_global,'lat00', &
nf_float,1,latll)
istatus=nf_put_att_real(ncid,nf_global,'lon00', &
nf_float,1,lonll)
istatus=nf_put_att_real(ncid,nf_global,'latNxNy', &
nf_float,1,latur)
istatus=nf_put_att_real(ncid,nf_global,'lonNxNy', &
nf_float,1,lonur)
istatus=nf_put_att_real(ncid,nf_global,'dxKm', &
nf_float,1,dxkm)
istatus=nf_put_att_real(ncid,nf_global,'dyKm', &
nf_float,1,dykm)
istatus=nf_put_att_real(ncid,nf_global,'latDxDy', &
nf_float,1,latdxdy)
istatus=nf_put_att_real(ncid,nf_global,'lonDxDy', &
nf_float,1,londxdy)
!
!-----------------------------------------------------------------------
!
! End define mode
!
!-----------------------------------------------------------------------
!
istatus=nf_enddef(ncid)
IF (istatus /= nf_noerr) CALL handle_err('nf_enddef',istatus)
!
!-----------------------------------------------------------------------
!
! Set inventory arrays.
!
!-----------------------------------------------------------------------
!
DO k=1,nprlvl
inventory(k:k)='1'
IF(iprlvl(k) > 999) THEN
WRITE(v3dlvl(k),'(a,i4)') 'MB ',iprlvl(k)
ELSE IF(iprlvl(k) > 99) THEN
WRITE(v3dlvl(k),'(a,i3)') 'MB ',iprlvl(k)
ELSE
WRITE(v3dlvl(k),'(a,i2)') 'MB ',iprlvl(k)
END IF
END DO
!
inventory_sfc='1'
!
CALL ctim2abss(year,month,day,hour,minute,second,itimref)
itimref=itimref-itim1970
PRINT *, ' reference time: ',itimref
DO i=1,ntime
reftime(i)=FLOAT(itimref)
END DO
!
!-----------------------------------------------------------------------
!
! End of first-read intializations
!
!-----------------------------------------------------------------------
!
init=.true.
END IF
!
!-----------------------------------------------------------------------
!
! Begin ARPS data conversions
!
!-----------------------------------------------------------------------
!
itime=itimref+nint(time)
valtime(ifile)=FLOAT(itime)
vtmref(ifile)=nint(time)
DO i=1,nx
DO j=1,ny
rain(i,j)=raing(i,j)+rainc(i,j)
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
pprt(i,j,k)=pprt(i,j,k)+pbar(i,j,k)
ptprt(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
qvprt(i,j,k)=qvprt(i,j,k)+qvbar(i,j,k)
tem1(i,j,k)=0.5*(uprt(i,j,k)+ubar(i,j,k)+ &
uprt(i+1,j,k)+ubar(i+1,j,k))
tem2(i,j,k)=0.5*(vprt(i,j,k)+vbar(i,j,k)+ &
vprt(i,j+1,k)+vbar(i,j+1,k))
tem3(i,j,k)=0.5*(wprt(i,j,k)+wbar(i,j,k)+ &
wprt(i,j,k+1)+wbar(i,j,k+1))
qi(i,j,k)=qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Swap wind data back into wind arrays
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
uprt(i,j,k)=tem1(i,j,k)
vprt(i,j,k)=tem2(i,j,k)
wprt(i,j,k)=tem3(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Put temperature into tem2 and -ln(p) into tem3
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=ptprt(i,j,k)* &
(pprt(i,j,k)/p0)**rddcp
tem3(i,j,k)=-ALOG(pprt(i,j,k))
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
zps_km(i,j,k)=zps(i,j,k)/1000.0
p_mb(i,j,k)=0.01*pprt(i,j,k)
mix(i,j,k)=1000.0*(qvprt(i,j,k)/(1.-qvprt(i,j,k)))
e_mb(i,j,k)=qvprt(i,j,k)*p_mb(i,j,k)/(qvprt(i,j,k)+ &
287.0/461.0)
esat_mb(i,j,k)=6.1078*EXP((2.500780E6/461.0)*((1.0/273.15)- &
(1.0/tem2(i,j,k))))
t_dew(i,j,k)=wmr2td(p_mb(i,j,k),mix(i,j,k)) &
+ 273.15
rh(i,j,k)=100.0*(e_mb(i,j,k)/esat_mb(i,j,k))
rh(i,j,k)=MAX(0.,rh(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate stability indices.
! Use level k=2 as the "surface".
!
!-----------------------------------------------------------------------
!
imid=(nx/2)+1
jmid=(ny/2)+1
CALL arps_be(nx,ny,nz, &
pprt,zps_km,tem2,qvprt, &
lcl,lfc,el,twdf,li,pbe,mbe,nbe,tcap, &
wrk1,wrk2,wrk3,wrk4,wrk5,wrk6, &
wrk7,wrk8,wrk9,wrk10,wrk11,wrk12,tem2d1)
PRINT *, ' Sample stability values: '
PRINT *, ' lcl,lfc : ',lcl(imid,jmid),lfc(imid,jmid)
PRINT *, ' el, twdf: ',el(imid,jmid),twdf(imid,jmid)
PRINT *, ' li, pbe : ',li(imid,jmid),pbe(imid,jmid)
PRINT *, ' mbe, nbe: ',mbe(imid,jmid),nbe(imid,jmid)
PRINT *, ' tcap : ',tcap(imid,jmid)
CALL calcshr(nx,ny,nz,x,y,zp,mf2d, &
pprt,tem2,uprt,vprt,pbe, &
shr37,ustrm,vstrm,srlfl,srmfl,heli,brn,brnu,blcon, &
tem2d1,tem2d2,tem2d3,tem2d4,tem4)
PRINT *, ' Sample shear values: '
PRINT *, ' shr37,ustrm,vstrm: ', &
shr37(imid,jmid),ustrm(imid,jmid),vstrm(imid,jmid)
PRINT *, ' srlfl,srmfl,heli: ', &
srlfl(imid,jmid),srmfl(imid,jmid),heli(imid,jmid)
PRINT *, ' brn,brnu,blcon: ', &
brn(imid,jmid),brnu(imid,jmid),blcon(imid,jmid)
!
! Store k=2 theta-e and dewpt in tem1,
! level 1 and 2 respectively.
!
DO j=1,ny-1
DO i=1,nx-1
prmb=0.01*pprt(i,j,2)
tdew=wmr2td(prmb,(1000.*qvprt(i,j,2)))
tem1(i,j,1)=oe((tem2(i,j,2)-273.15),tdew,prmb) + 273.15
tem1(i,j,2)=tdew+273.15
tem1(i,j,3)=prmb
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Output near-sfc data.
! Simularity theory or something could be applied here
! to make these data be valid at sfc instrument height,
! for now we output level 2.
!
!-----------------------------------------------------------------------
!
PRINT *, ' Storing near-sfc pressure'
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,3), &
cdf2d(1,1,1,idxp) )
PRINT *,' Storing grid-scale rainfall'
CALL cdf_fill(nx,ny,nxcdf,nycdf, raing, &
cdf2d(1,1,1,idxlgsp) )
PRINT *,' Storing convective rainfall'
CALL cdf_fill(nx,ny,nxcdf,nycdf, rainc, &
cdf2d(1,1,1,idxcp) )
PRINT *,' Storing total accumulated rainfall'
CALL cdf_fill(nx,ny,nxcdf,nycdf, rain, &
cdf2d(1,1,1,idxtp) )
PRINT *, ' Storing near-sfc u wind'
CALL cdf_fill(nx,ny,nxcdf,nycdf, uprt(1,1,2), &
cdf2d(1,1,1,idxsfu) )
PRINT *, ' Storing near-sfc v wind'
CALL cdf_fill(nx,ny,nxcdf,nycdf, vprt(1,1,2), &
cdf2d(1,1,1,idxsfv) )
PRINT *, ' Storing near-sfc temperature'
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem2(1,1,2), &
cdf2d(1,1,1,idxsft) )
PRINT *, ' Storing near-sfc dew point temperature'
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,2), &
cdf2d(1,1,1,idxdpt) )
PRINT *,' Storing near-sfc specific humidity'
CALL cdf_fill(nx,ny,nxcdf,nycdf, qvprt(1,1,2), &
cdf2d(1,1,1,idxsh) )
PRINT *, ' Storing near-sfc theta-e'
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1), &
cdf2d(1,1,1,idxept) )
PRINT *, ' Storing near-sfc LI'
CALL cdf_fill(nx,ny,nxcdf,nycdf, li, &
cdf2d(1,1,1,idxsli) )
PRINT *, ' Storing near-sfc CAPE'
CALL cdf_fill(nx,ny,nxcdf,nycdf, pbe, &
cdf2d(1,1,1,idxcape) )
PRINT *, ' Storing near-sfc CIN'
CALL cdf_fill(nx,ny,nxcdf,nycdf, nbe, &
cdf2d(1,1,1,idxcin) )
PRINT *, ' Storing near-sfc LCL'
CALL cdf_fill(nx,ny,nxcdf,nycdf, lcl, &
cdf2d(1,1,1,idxlcl) )
PRINT *, ' Storing near-sfc LFC'
CALL cdf_fill(nx,ny,nxcdf,nycdf, lfc, &
cdf2d(1,1,1,idxlfc) )
PRINT *, ' Storing near-sfc EL'
CALL cdf_fill(nx,ny,nxcdf,nycdf, el, &
cdf2d(1,1,1,idxel) )
PRINT *, ' Storing wet bulb temp diff'
CALL cdf_fill(nx,ny,nxcdf,nycdf, twdf, &
cdf2d(1,1,1,idxtwdf) )
PRINT *, ' Storing Cap Strength'
CALL cdf_fill(nx,ny,nxcdf,nycdf, tcap, &
cdf2d(1,1,1,idxtcap) )
PRINT *, ' Storing 3-7km Shear'
CALL cdf_fill(nx,ny,nxcdf,nycdf, shr37, &
cdf2d(1,1,1,idxshr37) )
PRINT *, ' Storing u-storm'
CALL cdf_fill(nx,ny,nxcdf,nycdf, ustrm, &
cdf2d(1,1,1,idxustrm) )
PRINT *, ' Storing v-storm'
CALL cdf_fill(nx,ny,nxcdf,nycdf, vstrm, &
cdf2d(1,1,1,idxvstrm) )
PRINT *, ' Storing low-level SR flow'
CALL cdf_fill(nx,ny,nxcdf,nycdf, srlfl, &
cdf2d(1,1,1,idxsrlfl) )
PRINT *, ' Storing mid-level SR flow'
CALL cdf_fill(nx,ny,nxcdf,nycdf, srmfl, &
cdf2d(1,1,1,idxsrmfl) )
PRINT *, ' Storing helicity'
CALL cdf_fill(nx,ny,nxcdf,nycdf, heli, &
cdf2d(1,1,1,idxheli) )
PRINT *, ' Storing Bulk Richardson Number'
CALL cdf_fill(nx,ny,nxcdf,nycdf, brn, &
cdf2d(1,1,1,idxbrn) )
PRINT *, ' Storing BRN Shear'
CALL cdf_fill(nx,ny,nxcdf,nycdf, brnu, &
cdf2d(1,1,1,idxbrnu) )
PRINT *, ' Storing Boundary Layer Convergence'
CALL cdf_fill(nx,ny,nxcdf,nycdf, blcon, &
cdf2d(1,1,1,idxblcon) )
!
!-----------------------------------------------------------------------
!
! Calculate constant pressure level data
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Put temperature into tem2 and -ln(p) into tem3.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=ptprt(i,j,k)* &
(pprt(i,j,k)/p0)**rddcp
tem3(i,j,k)=-ALOG(pprt(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate temperature (K) at ARPS grid points
! and 700mb pressure level
!
!-----------------------------------------------------------------------
!
rln700=-ALOG(70000.0)
CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, &
tem2,tem3,rln700,tem1)
!
!-----------------------------------------------------------------------
!
! Calculate sea level pressure (mb)
! Reduction method: Benjamin and Miller: 1990, MWR, vol.118, No.10,
! Page: 2100-2101
!
!-----------------------------------------------------------------------
!
! gamma=.0065 ! std lapse rate per meter
! ex2=5.2558774
! DO 355 j=1,ny-1
! DO 355 i=1,nx-1
! p00 = 0.01*(pprt(i,j,2))
! tem1(i,j,1)=p00*
! : (1.0+gamma*zps(i,j,2)/tem1(i,j,1))**ex2
! 355 CONTINUE
gamma=.0065 ! std lapse rate per meter
ex1=0.1903643 ! R*gamma/g
ex2=5.2558774
DO j=1,ny-1
DO i=1,nx-1
p00 = 0.01*(pprt(i,j,2))
IF (p00 <= 700.0) THEN
t00 = tem2(i,j,2)
ELSE
t00 = tem1(i,j,1)*(p00/700.0)**ex1
END IF
tem1(i,j,1)=p00*(1.0+gamma*zps(i,j,2)/t00)**ex2
END DO
END DO
PRINT *, ' Storing MSL Pressure '
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1, &
cdf2d(1,1,1,idxmslp) )
!
!-----------------------------------------------------------------------
!
! Store terrain data.
!
!-----------------------------------------------------------------------
!
PRINT *, ' Storing terrain data.'
CALL cdf_fill(nx,ny,nxcdf,nycdf, zp(1,1,2), &
static2d(1,1,idxtop) )
PRINT *, ' Storing Coriolis data'
CALL cdf_fill(nx,ny,nxcdf,nycdf, coriolis, &
static2d(1,1,idxcor) )
PRINT *, ' Storing spacing data.'
CALL cdf_fill(nx,ny,nxcdf,nycdf, spacing, &
static2d(1,1,idxspa) )
!
!-----------------------------------------------------------------------
!
! Calculate atmospheric state variables.
!
!-----------------------------------------------------------------------
!
DO klev=1,nprlvl
rlnpcdf=-ALOG(100.*FLOAT(iprlvl(klev)))
PRINT *, ' Storing netCDF data at pr= ',iprlvl(klev)
!
! Store gh
!
CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, &
zps,tem3,rlnpcdf,tem1)
CALL extrph(nx,ny,nz,zps,tem2,pprt, &
iprlvl(klev),tem1)
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1, &
cdf3d(1,1,klev,idxgh) )
!
! Store temperature (tem2)
!
CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, &
tem2,tem3,rlnpcdf,tem1)
CALL extrpt(nx,ny,nz,tem2,pprt,zps, &
iprlvl(klev),tem1)
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1, &
cdf3d(1,1,klev,idxt) )
!
! Store relative humidity
!
CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, &
rh,tem3,rlnpcdf,tem1)
CALL extrpq(nx,ny,nz,rh,pprt,iprlvl(klev),tem1)
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1, &
cdf3d(1,1,klev,idxrh) )
!
! Store u and v wind components
!
CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, &
uprt,tem3,rlnpcdf,tem1(1,1,1))
CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, &
vprt,tem3,rlnpcdf,tem1(1,1,2))
CALL extrpuv(nx,ny,nz,uprt,vprt,pprt,zps, &
iprlvl(klev),tem1(1,1,1),tem1(1,1,2))
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1), &
cdf3d(1,1,klev,idxuw) )
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,2), &
cdf3d(1,1,klev,idxvw) )
!
! Store vertical velocity
!
CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, &
wprt,tem3,rlnpcdf,tem1)
CALL extrpq(nx,ny,nz,wprt,pprt,iprlvl(klev),tem1)
CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1), &
cdf3d(1,1,klev,idxww) )
END DO
END IF ! Good return from read
!
invlen(1)=nprlvl
invlen(2)=1
vecistr(2)=ifile
volistart(4)=ifile
DO ivar=1,n3dvar
WRITE(6,'(a,i4,2x,a)') ' Writing 3d ivar:',ivar, &
v3dlngnam(ivar)
!
IF(ifile == 1) THEN
istatus=nf_put_vara_text(ncid,v3dlvlid(ivar), &
planeistr,lvllen,v3dlvl)
IF (istatus /= nf_noerr) &
CALL handle_err('nf_put_vara_text 3d lvl',istatus)
END IF
istatus=nf_put_vara_text(ncid,v3dinvid(ivar), &
vecistr,invlen,inventory)
IF (istatus /= nf_noerr) &
CALL handle_err('nf_put_vara_text 3d inv',istatus)
istatus=nf_put_vara_real(ncid,v3did(ivar), &
volistart,vollen,cdf3d(1,1,1,ivar))
IF (istatus /= nf_noerr) &
CALL handle_err('nf_put_vara_real 3d var',istatus)
END DO
!
invlen(1)=1
invlen(2)=1
planeistr(4)=ifile
DO ivar=1,n2dvar
WRITE(6,'(a,i4,2x,a)') ' Writing 2d ivar:',ivar, &
v2dlngnam(ivar)
IF(ifile == 1) THEN
IF(ivar == idxmslp) THEN
v2dlvl='MSL'
ELSE
v2dlvl='SFC'
END IF
istatus=nf_put_vara_text(ncid,v2dlvlid(ivar), &
vecistr,sfclen,v2dlvl)
IF (istatus /= nf_noerr) &
CALL handle_err('nf_put_vara_text 2d lvl',istatus)
END IF
istatus=nf_put_vara_text(ncid,v2dinvid(ivar), &
vecistr,invlen,inventory_sfc)
IF (istatus /= nf_noerr) &
CALL handle_err('nf_put_vara_text 2d inv',istatus)
istatus=nf_put_vara_real(ncid,v2did(ivar), &
planeistr,planelen,cdf2d(1,1,1,ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_put_vara_real 2d',istatus)
END DO
END DO
!
! Write static variables
!
DO ivar=1,nstvar
WRITE(6,'(a,i4,2x,a)') ' Writing st ivar:',ivar, &
vstlngnam(ivar)
istatus=nf_put_vara_real(ncid,vstid(ivar), &
planeistr,planelen,static2d(1,1,ivar))
IF (istatus /= nf_noerr) CALL handle_err('nf_put_vara_real static',istatus)
END DO
!
! Write 1-d variables
!
istatus=nf_put_vara_int(ncid,vtmrefid, &
1,ntime,vtmref)
istatus=nf_put_vara_double(ncid,reftimid, &
1,ntime,reftime)
istatus=nf_put_vara_double(ncid,valtimid, &
1,ntime,valtime)
ls=LEN(origin)
CALL strlnth(origin,ls)
istatus=nf_put_vara_text(ncid,origid, &
1,ls,origin)
ls=LEN(model)
CALL strlnth(model,ls)
istatus=nf_put_vara_text(ncid,modelid, &
1,ls,model)
!
! Close file
!
istatus=nf_close(ncid)
IF (istatus /= nf_noerr) CALL handle_err('nf_close',istatus)
STOP
!
900 CONTINUE
WRITE(6,'(a)') ' Error reading input data'
950 CONTINUE
WRITE(6,'(a)') ' Error setting up netCDF file'
STOP
END PROGRAM arps2ncdf
!
SUBROUTINE extrph(nx,ny,nz,zps,t,pr,iprlvl,hgtcdf) 2
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EXTRPH ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Extrapolate height by using a std atmos lapse rate
! below the last physical level. Above the domain,
! assume a constant temperature above 300 mb, otherwise
! use the std atmos lapse rate.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster, from arps2gem
! 2/19/99
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: zps(nx,ny,nz)
REAL :: t(nx,ny,nz)
REAL :: pr(nx,ny,nz)
INTEGER :: iprlvl
REAL :: hgtcdf(nx,ny)
!
INCLUDE 'phycst.inc'
!
REAL :: gamma,rddg,const
PARAMETER ( gamma = 0.0065, & ! degrees/m lapse rate
rddg = (rd/g), &
const = (rd*gamma/g) )
!
INTEGER :: i,j
REAL :: prcdf
!
prcdf=100.*FLOAT(iprlvl)
DO j=1,ny-1
DO i=1,nx-1
IF(prcdf < pr(i,j,nz-1)) THEN
IF(pr(i,j,nz-1) <= 30000.) THEN
hgtcdf(i,j)=zps(i,j,nz-1) + &
rddg*t(i,j,nz-1)*ALOG(pr(i,j,nz-1)/prcdf)
ELSE
hgtcdf(i,j)=zps(i,j,nz-1) + (t(i,j,nz-1)/gamma)* &
(1.-(prcdf/pr(i,j,nz-1))**const)
END IF
ELSE IF(prcdf >= pr(i,j,2)) THEN
hgtcdf(i,j)=zps(i,j,2) + (t(i,j,2)/gamma)* &
(1.-(prcdf/pr(i,j,2))**const)
END IF
END DO
END DO
RETURN
END SUBROUTINE extrph
!
SUBROUTINE extrpt(nx,ny,nz,t,pr,zps,iprlvl,tcdf) 3
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EXTRPT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Extrapolate temperature by using a std atmos lapse rate
! below the last physical level. Above the domain,
! assume a constant temperature above 300 mb, otherwise
! use the std atmos lapse rate.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster, from arps2gem
! 2/19/99
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz
REAL :: t(nx,ny,nz)
REAL :: pr(nx,ny,nz)
REAL :: zps(nx,ny,nz)
INTEGER :: iprlvl
REAL :: tcdf(nx,ny)
!
INCLUDE 'phycst.inc'
!
REAL :: gamma,const
PARAMETER ( gamma = 0.0065, & ! degrees/m lapse rate
const = (rd*gamma/g) )
!
INTEGER :: i,j
REAL :: prcdf
!
prcdf=100.*FLOAT(iprlvl)
DO j=1,ny-1
DO i=1,nx-1
IF(prcdf <= pr(i,j,nz-1)) THEN
IF(pr(i,j,nz-1) <= 30000.) THEN
tcdf(i,j)=t(i,j,nz-1)
ELSE
tcdf(i,j)=t(i,j,nz-1)* &
((prcdf/pr(i,j,nz-1))**const)
END IF
ELSE IF(prcdf >= pr(i,j,2)) THEN
tcdf(i,j)=t(i,j,2)* &
((prcdf/pr(i,j,2))**const)
END IF
END DO
END DO
RETURN
END SUBROUTINE extrpt
!
SUBROUTINE extrpq(nx,ny,nz,q,pr,iprlvl,qcdf) 8
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EXTRPQ ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Assign a zero-gradient value to scalars below ground.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster, from arps2gem
! 2/19/99
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: q(nx,ny,nz)
REAL :: pr(nx,ny,nz)
INTEGER :: iprlvl
REAL :: qcdf(nx,ny)
!
INTEGER :: i,j
REAL :: prcdf
!
prcdf=100.*FLOAT(iprlvl)
DO j=1,ny-1
DO i=1,nx-1
IF(prcdf < pr(i,j,nz-1)) THEN
qcdf(i,j)=q(i,j,nz-1)
ELSE IF(prcdf > pr(i,j,2)) THEN
qcdf(i,j)=q(i,j,2)
END IF
END DO
END DO
RETURN
END SUBROUTINE extrpq
SUBROUTINE extrpuv(nx,ny,nz,us,vs,pr,zps, & 2
iprlvl,ucdf,vcdf)
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE EXTRPUV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Assign a constant value to sub-terrainian winds
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster, from arps2gem
! 2/19/99
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: us(nx,ny,nz)
REAL :: vs(nx,ny,nz)
REAL :: pr(nx,ny,nz)
REAL :: zps(nx,ny,nz)
INTEGER :: iprlvl
REAL :: ucdf(nx,ny)
REAL :: vcdf(nx,ny)
!
INTEGER :: i,j
REAL :: prcdf
!
prcdf=100.*FLOAT(iprlvl)
DO j=1,ny-1
DO i=1,nx-1
IF(prcdf < pr(i,j,nz-1)) THEN
ucdf(i,j)=us(i,j,nz-1)
vcdf(i,j)=vs(i,j,nz-1)
ELSE IF(prcdf > pr(i,j,2)) THEN
ucdf(i,j)=us(i,j,2)
vcdf(i,j)=vs(i,j,2)
END IF
END DO
END DO
RETURN
END SUBROUTINE extrpuv
!
SUBROUTINE handle_err(string,istatus) 32
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE HANDLE_ERR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Report netCDF error in plain English. Exit.
!
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster
! 2/19/99
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
CHARACTER (LEN=*) :: string
CHARACTER (LEN=80) :: nf_strerror
INTEGER :: istatus
WRITE(6,'(a,a,/a)') 'netCDF error: ',string, &
nf_strerror(istatus)
STOP
END SUBROUTINE handle_err
!
SUBROUTINE cdf_fill(nx,ny,nxcdf,nycdf, &
arpsvar,cdfvar)
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CDF_FILL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Fill netCDF output array which is the physical domain of the
! ARPS scalar grid.
!
!-----------------------------------------------------------------------
!
!
! AUTHOR: Keith Brewster
! 2/19/99
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny
INTEGER :: nxcdf,nycdf
REAL :: arpsvar(nx,ny)
REAL :: cdfvar(nxcdf,nycdf)
INTEGER :: i,j
!
DO j=2,ny-2
DO i=2,nx-2
cdfvar(i-1,j-1)=arpsvar(i,j)
END DO
END DO
RETURN
END SUBROUTINE cdf_fill