PROGRAM arps2ncdf,106 ! !################################################################## !################################################################## !###### ###### !###### 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) ! ! 1 June 2002 Eric Kemp ! Soil variable updates. ! ! 11/22/2003 (K. Brewster) ! Added logic to allow appending new time levels to existing file. ! Includes additions to input file. Corrected itim1970 constant. ! ! 12/02/2003 (T. Oram, NWS SFMG) ! Modifications for AWIPS compatibility plus comments to help ! myself understand what is going on ! a. This program works with the arps.cdl netCDF cdl file that is used ! for the AWIPS localization process. The arps.cdl file is new. ! b. Changed the tracking of the current and maximum number of files that ! can be appended to a netCDF file using the record and n_valtimes ! dimensions in the netCDF file. ! c. modified arps2ncdf.input for new variables and modified write of ! 1-d variables (valtime, reftime, valtimeMINUSreftime, etc) to only ! be written once for netcdf append files. This assumes that all forecast ! times are known when the first history file is read and converted to ! netCDF. Plan to move the number of pressure levels and the pressure level ! array into arps2ncdf.input for run-time determination rather than ! netcdf.inc for compile-time determination. ! d. The left most dimension of the inventory variables also was ! changed to n_valtimes rather than record. AWIPS chokes ! without a full inventory array for all valid times. The new dimension of ! the inventory array matches the AWIPS convention. ! e. added surface level to the 3-dimensional array to match the awips convention ! f. added cw, cice, and tke to 3-dimensional array (tdo, 02/02/2004) ! g. removed 2-d fields not included in the dataFieldTable.txt file. ! also changed reftime to remove the record dimension to match awips (tdo, 02/05/2004) ! ! 02/22/2005 (Y. Wang) ! Moved the pressure level specification to arps2ncdf.input for run-time ! determination as planed by T. Oram. ! ! Cleaned up documents and formats. ! !----------------------------------------------------------------------- ! ! 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) ! zpsoil z coordinate of soil model (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 ! ! tsoil Soil Temperature (K) ! qsoil 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 :: nzsoil ! Soil levels INTEGER :: nstyps ! Maximum number of soil types. INTEGER, PARAMETER :: nhisfile_max = 200 INTEGER :: hinfmt,nhisfile,lengbf,lenfil CHARACTER (LEN=256) :: grdbasfn,hisfile(nhisfile_max) ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INCLUDE 'phycst.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 :: zpsoil(:,:,:) ! The height of the soil model. 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 :: tsoil (:,:,:,:) ! Soil Temperature (K) REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! 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 :: radswnet(:,:) ! Net shortwave radiation REAL, ALLOCATABLE :: radlwin(:,:) ! Incoming longwave radiation 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(:,:,:) ! J.Case, ENSCO Inc. (9/29/2004) ! Additional variables. ! Note: radar_comp is maximum column reflecitivity and ! echo_top is maximum height of 18 dBZ echo in a column. REAL, ALLOCATABLE :: radar_3d(:,:,:),radar_comp(:,:),echo_top(:,:) REAL, ALLOCATABLE :: mf2d(:,:) REAL, ALLOCATABLE :: coriolis(:,:) REAL, ALLOCATABLE :: spacing(:,:) ! !----------------------------------------------------------------------- ! ! Misc ARPS variables ! !----------------------------------------------------------------------- ! INTEGER :: nch INTEGER :: first18 ! J.Case, ENSCO Inc. (9/29/2004) REAL :: time REAL :: latnot(2) LOGICAL :: init ! !----------------------------------------------------------------------- ! ! netCDF dimensions ! ! nprlvl defined in arps2cdf.inc is the number of pressure levels ! nxcdf number of points in x direction in netcdf file ! nycdf number of points in y direction in netcdf file ! nrecord number of records in netcdf file ! ntotaltime number of valid times that are to be or can be written into file ! ncdflvls number of levels in 3-d grids in file; nprlvl + 1 (surface) ! !----------------------------------------------------------------------- ! INTEGER, PARAMETER :: mxtime = 40 INTEGER :: nxcdf,nycdf,nrecord,nvaltimes,ntotaltimes INTEGER, PARAMETER :: nprlvl_max = 128 INTEGER :: iprlvl(nprlvl_max) INTEGER :: nprlvl, ncdflvls ! = nprlvl+1 INTEGER, PARAMETER :: n3dvar = 10, n2dvar = 16, nstvar = 3 INTEGER, PARAMETER :: nvartot = (n3dvar+n2dvar) ! J.Case, ENSCO Inc. (9/29/2004) -- Added 1 new 3D ! variable (rr) and 2 new 2D variables (cxr, mret) ! !----------------------------------------------------------------------- ! ! netCDF variables ! !----------------------------------------------------------------------- ! INTEGER :: istatus,ncid INTEGER :: dimidrec,dimidtot,dimidnt,dimidnav INTEGER :: dimidnx,dimidny,dimid1,dimidnz,dimchrlvl,dimnamlen INTEGER :: vtmrefid,valtimid,reftimid,origid,modelid,ntwrtid 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) INTEGER :: ntwrt 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) CHARACTER (LEN=4) :: v3dnam(n3dvar) CHARACTER (LEN=4) :: v2dnam(n2dvar) CHARACTER (LEN=14) :: vstnam(nstvar) CHARACTER (LEN=10) :: v3dlvl(nprlvl_max) 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_max) :: inventory CHARACTER (LEN=1) :: inventory_sfc CHARACTER (LEN=1) :: inventory_blank 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,idxcw,idxci,idxtk ! J.Case, ENSCO Inc. (9/29/2004) INTEGER :: idxrr PARAMETER (idxgh=1,idxrh=2,idxt=3, & idxuw=4,idxvw=5,idxww=6, & idxcw=7,idxci=8,idxtk=9,idxrr=10) ! ! J.Case, cont. (9/29/2004) -- Added radar reflectivity to 3D variables (rr). DATA v3dnam /'gh','rh','t','uw','vw','ww','cw','cice','tke','rr'/ DATA v3dlngnam /'Geopotential Height', & 'Relative Humidity', & 'Temperature', & 'u Wind Component', & 'v Wind Component', & 'w Wind Component', & 'Cloud Water', & 'Cloud Ice', & 'Turbulent Kinetic Energy', & 'Radar Reflectivity'/ DATA v3dunits /'geopotential meters', & 'percent', & 'degrees K', & 'meters/second', & 'meters/second', & 'meters/second', & 'kg/kg', & 'kg/kg', & '(meters/second)**2', & 'dBZ'/ DATA v3dmin / 0., 0.,150.,-300.,-300.,-100.,0.,0.,0.,-200./ DATA v3dmax /40000.,110.,400., 300., 300., 100.,1.,1.,1000.,200./ INTEGER :: idxmslp,idxp,idxlgsp,idxcp,idxtp INTEGER :: idxdpt INTEGER :: idxsh,idxept,idxsli,idxcape,idxcin INTEGER :: idxustrm,idxvstrm INTEGER :: idxheli ! J.Case, ENSCO Inc. (9/29/2004) -- new 2D variables cxr and mret INTEGER :: idxcxr,idxmret PARAMETER (idxmslp=1, & idxp =2, & idxlgsp=3, & idxcp =4, & idxtp =5, & idxdpt =6, & idxsh =7, & idxept =8, & idxsli =9, & idxcape=10, & idxcin =11, & idxustrm=12, & idxvstrm=13, & idxheli =14, & idxcxr =15, & idxmret =16) ! idxsfu =6, & ! idxsfv =7, & ! idxsft =8, & ! idxlcl =15, & ! idxlfc =16, & ! idxel =17, & ! idxtwdf =18, & ! idxtcap =19, & ! idxshr37=20, & ! idxsrlfl=23, & ! idxsrmfl=24, & ! idxbrn =26, & ! idxbrnu =27, & ! idxblcon=28) ! these 2d variable names must match cdl names in the AWIPS ! dataFieldTable.txt file ! DATA v2dnam /'mslp','p','lgsp','cp','tp','sfu','sfv', & ! 'sft','dpt','sh','ept','sli','cape','cin','lcl', & ! 'lfc','el','twdf','tcap','sh37','ustm','vstm', & ! 'srlf','srmf','heli','brn','brnu','blcn'/ !JLC (8/23/2004) -- added cxr,mret to 2D variables. DATA v2dnam /'mslp','p','lgsp','cp','tp', & 'dpt','sh','ept','sli','cape','cin', & 'ustm','vstm', & 'heli','cxr','mret'/ DATA v2dlngnam /'Mean Sea Level Pressure', & 'Surface Pressure', & 'Grid Scale Precipitation', & 'Convective Precipitation', & 'Total Precipitation', & 'Surface Dew Point', & 'Specific Humidity', & 'Theta-e', & 'Surface Lifted Index', & 'CAPE', & 'CIN', & 'Est Storm Motion u comp', & 'Est Storm Motion v comp', & 'Storm Rel Helicity', & 'Column Max Reflectivity', & 'Maximum Radar Echo Tops'/ ! 'Surface u Wind Component', & ! 'Surface v Wind Component', & ! 'Surface Temperature', & ! 'Lifting Condensation Lvl', & ! 'Level of Free Convection', & ! 'Equilibrium Level', & ! 'Max Wet-Bulb Temp Diff', & ! 'Cap Strength', & ! '3-7km Shear', & ! 'Storm Rel Low-Level Flow', & ! 'Storm Rel Mid-Level Flow', & ! 'Bulk Richardson Number', & ! 'BRN Shear u', & ! 'Bound-Layer Conv x 1000'/ ! J.Case, cont. (9/29/2004) DATA v2dunits /'Pascals','Pascals','mm','mm','mm', & 'degrees K','kg/kg', & 'degrees K','degrees C','J/kg','J/kg', & 'meters/second','meters/second', & 'm2/s2','dBZ','meters'/ DATA v2dmin / 80000.,50000.,0.,0.,0., & -100.,0., & 200.,-100.,0.,0., & -100.,-100., & -900.,0.,0./ DATA v2dmax / 110000.,120000.,1000.,1000.,1000., & 100.,1., & 500.,100.,10000.,10000., & 100.,100., & 900.,100.,20000./ 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=25) :: 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, PARAMETER :: itim1970=315622800 CHARACTER(LEN=8), PARAMETER :: cdldate = '20031202' CHARACTER(LEN=40), PARAMETER :: depictor = 'ARPS' CHARACTER(LEN=132), PARAMETER :: origin = 'CAPS' CHARACTER(LEN=132), PARAMETER :: model = 'ARPS 5.3' ! !----------------------------------------------------------------------- ! ! Namelists ! !----------------------------------------------------------------------- ! CHARACTER (LEN=256) :: cdfname INTEGER :: ipktyp,nbits,ncdapnd,ncdntime,ncdtinterval NAMELIST /ncdf_options/ nprlvl,iprlvl,cdfname,ipktyp,nbits,ncdapnd, & ncdntime,ncdtinterval ! !----------------------------------------------------------------------- ! ! 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,idxtime, ntime INTEGER :: i,j,k,klev,ifile,ireturn,imid,jmid,ncmode REAL :: rlnpcdf,r2d REAL :: ctrx,ctry,swx,swy REAL :: gamma,ex1,ex2,rln700,p00,t00 REAL :: prmb,tdew LOGICAL :: adddefs,ncexist REAL, ALLOCATABLE :: p_mb(:,:,:) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !----------------------------------------------------------------------- ! ! Initializations ! !----------------------------------------------------------------------- ! ! WRITE(6,'(9(/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. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! get_input_filenames reads the history file input options ! from the namelist then open the file and read some of ! the grid dimension information ! !----------------------------------------------------------------------- ! 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,nzsoil,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 ! !----------------------------------------------------------------------- ! ! Read the netcdf output options from the netcdf namelist options ! ncdapnd option to append to netcdf file ! ncdntime the number of valid times in the file ! the name of the grid/base data set. ! !----------------------------------------------------------------------- ! ncdapnd=0 ipktyp=1 nbits=16 READ(5,ncdf_options,END=900) lengbf=LEN_trim(grdbasfn) WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf) !----------------------------------------------------------------------- ! output grid is 3 gridpoints less than history file grid ! in both x and y directions !----------------------------------------------------------------------- ncdflvls = nprlvl + 1 ! include surface layer nxcdf=(nx-3) nycdf=(ny-3) planelen(1) = nxcdf planelen(2) = nycdf planelen(3) = 1 planelen(4) = 1 lvllen(1) = 10 lvllen(2) = ncdflvls sfclen(1) = 10 sfclen(2) = 1 vollen(1) = nxcdf vollen(2) = nycdf vollen(3) = ncdflvls vollen(4) = 1 ncdntime = max(nhisfile,ncdntime) !----------------------------------------------------------------------- ! ! now that you've read all the dimensions at execute time ! allocate the memory for the arrays ! !----------------------------------------------------------------------- ALLOCATE(x (nx)) ALLOCATE(y (ny)) ALLOCATE(z (nz)) ALLOCATE(zp (nx,ny,nz)) ALLOCATE(zpsoil (nx,ny,nzsoil)) 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(tsoil (nx,ny,nzsoil,0:nstyps)) ALLOCATE(qsoil (nx,ny,nzsoil,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(radswnet(nx,ny)) ALLOCATE(radlwin(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)) ! J.Case, ENSCO Inc. (9/29/2004) ALLOCATE(radar_3d(nx,ny,nz)) ALLOCATE(radar_comp(nx,ny)) ALLOCATE(echo_top(nx,ny)) 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,ncdflvls,n3dvar)) !----------------------------------------------------------------------- ! initialize a bunch of stuff !----------------------------------------------------------------------- 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 tsoil =0.0 qsoil =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 radswnet =0.0 radlwin =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 ! J.Case, ENSCO Inc. (9/29/2004) radar_3d=0.0 radar_comp=0.0 echo_top=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. nrecord=0 nvaltimes=0 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,nzsoil,nstyps,hinfmt,nch,grdbasfn(1:lengbf), & lengbf, & hisfile(ifile)(1:lenfil),lenfil,time, & x,y,z,zp,zpsoil,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, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn, tem1,tem2,tem3) ! !----------------------------------------------------------------------- ! ! ireturn = 0 for a successful read ! !----------------------------------------------------------------------- ! IF( ireturn == 0 ) THEN ! successful read ! !=============================================================== ! ! J.Case, ENSCO, Inc. 9/29/2004 -- added subroutine call for ! derived radar reflectivity products. ! !=============================================================== ! WRITE(*,*) ' ' WRITE (*,*) 'Deriving radar reflectivity.' CALL REFLEC (nx,ny,nz, rhobar, qr, qs, qh, radar_3d) DO j=1,ny-1 DO i=1,nx-1 radar_comp(i,j) = 0.0 echo_top(i,j) = 0.0 first18 = 1 DO k=1,nz-1 IF (radar_3d(i,j,k) > radar_comp(i,j)) THEN radar_comp(i,j) = radar_3d(i,j,k) ENDIF END DO DO k=nz-1,1,-1 IF (radar_3d(i,j,k) >= 18.0 .AND. first18 == 1) THEN echo_top(i,j) = zp(i,j,k) first18 = 0 ENDIF END DO END DO END DO 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 ! !----------------------------------------------------------------------- ! ! 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 = ACOS(-1.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 !----------------------------------------------------------------------- ! ! Open netCDF file ! !----------------------------------------------------------------------- ! INQUIRE(FILE=cdfname,EXIST=ncexist) !----------------------------------------------------------------------- ! if this is an append and the netCDF files exists ! then read the record and n_valtimes dimensions from the ! netcdf file !----------------------------------------------------------------------- IF(ncdapnd > 0 .AND. ncexist) THEN ncmode=NF_WRITE adddefs=.false. istatus=NF_OPEN(cdfname,ncmode,ncid) IF (istatus /= NF_NOERR) CALL handle_err('nf_open',istatus) ! i think we want to replace with a read of the ! record dimension to determine the record to write to ! and the n_valtimes dimension to determine the maximum ! number of forecast times that can be included in the file (tdo) istatus=NF_INQ_DIMID(ncid,'record',dimidrec) IF (istatus /= NF_NOERR) CALL handle_err('nf_inq_dimid record',istatus) istatus=NF_INQ_DIMLEN(ncid,dimidrec,nrecord) istatus=NF_INQ_DIMID(ncid,'n_valtimes',dimidnt) IF (istatus /= NF_NOERR) CALL handle_err('nf_inq_dimid ntime',istatus) istatus=NF_INQ_DIMLEN(ncid,dimidnt,nvaltimes) WRITE(6,'(2(a,i5))') 'current record =',nrecord, & ',number of valid times =',nvaltimes, & ',number of history files =',nhisfile ntotaltimes=nrecord+nhisfile IF (ntotaltimes > nvaltimes) THEN WRITE(6,'(a)') 'Trying to write too many hours into netcdf file' STOP ENDIF ntwrt=nrecord ! i think this is the simplest change to preserve Keith Brewsters ! append capability while living within the AWIPS netCDF conventions ! istatus=nf_inq_varid(ncid,'ntwrite',ntwrtid) ! IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus) ! istatus=nf_get_var_int(ncid,ntwrtid,ntwrt) ! IF (istatus /= nf_noerr) CALL handle_err('nf_get_var_int',istatus) ELSE ncmode=NF_CLOBBER adddefs=.true. istatus=NF_CREATE(cdfname,ncmode,ncid) ntwrt=0 IF (istatus /= NF_NOERR) CALL handle_err('nf_create',istatus) END IF ! !----------------------------------------------------------------------- ! ! Define dimensions ! !----------------------------------------------------------------------- ! IF(ncdflvls < 10) THEN WRITE(nzstr,'(a,i1)') 'levels_',ncdflvls ELSE IF(ncdflvls < 100) THEN WRITE(nzstr,'(a,i2)') 'levels_',ncdflvls ELSE WRITE(nzstr,'(a,i3)') 'levels_',ncdflvls END IF IF( adddefs ) THEN 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',ncdntime,dimidnt) IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim ncdntime',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,ncdflvls,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) ! ! The following arrays are used to create the dimensions of the variables ! that are defined next. The arrays consist of dimension ids (from above) ! When listed by the ncdump utility, the 1st dimension (i.e. arraydim(1)) is ! the right most dimension listed while the last dimension is the left most ! dimension. The order of dimensions here is set to match the AWIPS convention ! 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)=dimidnt ! !----------------------------------------------------------------------- ! ! 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,ncdflvls) istatus=nf_put_att_text(ncid,v3did(ivar),'levels', & 10,'MB and SFC') 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) ! ! create the inventory variable name by appending Inventory to variable ! then define the inventory variable ! 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)=dimidnt 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) ! ! create the inventory variable name by appending Inventory to variable ! then define the inventory variable ! 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,0, & 0,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) ! commented out the following definition of ntwrt variable ! probably should be a global attribute (tdo) ! istatus=nf_def_var(ncid,'ntwrite',nf_int,0,0,ntwrtid) ! IF (istatus /= nf_noerr) CALL handle_err('nf_def_var ntwrite',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 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) ELSE ! get variable ids from existing dataset ! ! Get dimension ids ! The inquery for record and n_valtimes is redundant now, but ok so won't touch ! istatus=nf_inq_dimid(ncid,'record',dimidrec) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid record',istatus) istatus=nf_inq_dimid(ncid,'n_valtimes',dimidnt) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid ntime',istatus) istatus=nf_inq_dimid(ncid,'data_variables',dimidtot) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid nvartot',istatus) istatus=nf_inq_dimid(ncid,'charsPerLevel',dimchrlvl) IF (istatus /= nf_noerr) & CALL handle_err('nf_inq_dimid charsPerLevel',istatus) istatus=nf_inq_dimid(ncid,'namelen',dimnamlen) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid namelen',istatus) istatus=nf_inq_dimid(ncid,'x',dimidnx) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid nx',istatus) istatus=nf_inq_dimid(ncid,'y',dimidny) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid ny',istatus) istatus=nf_inq_dimid(ncid,'levels_1',dimid1) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid 1',istatus) istatus=nf_inq_dimid(ncid,nzstr,dimidnz) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid levels',istatus) istatus=nf_inq_dimid(ncid,'nav',dimidnav) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid nav',istatus) ! ! Get variable ids ! DO ivar=1,n3dvar istatus=nf_inq_varid(ncid,v3dnam(ivar),v3did(ivar)) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 3d',istatus) ls=LEN(v3dnam(ivar)) CALL strlnth(v3dnam(ivar),ls) WRITE(v3dinv,'(a,a)') v3dnam(ivar)(1:ls),'Inventory' istatus=nf_inq_varid(ncid,v3dinv,v3dinvid(ivar)) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 3dinv',istatus) END DO DO ivar=1,n2dvar istatus=nf_inq_varid(ncid,v2dnam(ivar),v2did(ivar)) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 2d',istatus) ls=LEN(v2dnam(ivar)) CALL strlnth(v2dnam(ivar),ls) WRITE(v2dinv,'(a,a)') v2dnam(ivar)(1:ls),'Inventory' istatus=nf_inq_varid(ncid,v2dinv,v2dinvid(ivar)) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 2dinv',istatus) END DO istatus=nf_inq_varid(ncid,'valtimeMINUSreftime',vtmrefid) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus) ! we're going to have trouble here need to replace with a read of the ! number of the number of records in the dimensions istatus=nf_get_vara_int(ncid,vtmrefid,1,ntwrt,vtmref) IF (istatus /= nf_noerr) CALL handle_err('nf_get_vara_int',istatus) istatus=nf_inq_varid(ncid,'reftime',reftimid) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus) istatus=nf_get_vara_double(ncid,reftimid,1,1,reftime) IF (istatus /= nf_noerr) CALL handle_err('nf_get_vara_double',istatus) istatus=nf_inq_varid(ncid,'valtime',valtimid) IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus) istatus=nf_get_vara_double(ncid,valtimid,1,ntwrt,valtime) IF (istatus /= nf_noerr) CALL handle_err('nf_get_vara_double',istatus) END IF ! !----------------------------------------------------------------------- ! ! Set inventory arrays. ! inventory(k:k) appends the right number of 1's to the string ! !----------------------------------------------------------------------- ! 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 ! add surface to 3d level definition (tdo) k=ncdflvls inventory(k:k)='1' WRITE(v3dlvl(k),'(a)') 'SFC' ! ! added a blank inventory (tdo) ! inventory_sfc='1' inventory_blank='' CALL ctim2abss(year,month,day,hour,minute,second,itimref) itimref=itimref-itim1970 ! !----------------------------------------------------------------------- ! ! End of first-read intializations ! !----------------------------------------------------------------------- ! init=.true. END IF ! !----------------------------------------------------------------------- ! ! Begin ARPS data conversions ! !----------------------------------------------------------------------- ! itime=itimref+nint(time) idxtime=ifile+ntwrt valtime(idxtime)=FLOAT(itime) vtmref(idxtime)=nint(time) reftime(idxtime)=FLOAT(itimref) 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 ! !----------------------------------------------------------------------- ! uprt=0. vprt=0. wprt=0. 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. This is approximately 10m in ! our configuration ! !----------------------------------------------------------------------- ! ! addition of near surface gph (tdo) PRINT *, ' Storing near-sfc gh' CALL cdf_fill(nx,ny,nxcdf,nycdf, zps(1,1,2), & cdf3d(1,1,ncdflvls,idxgh) ) ! addition of near surface rh (tdo) PRINT *, ' Storing near-sfc rh' CALL cdf_fill(nx,ny,nxcdf,nycdf, rh(1,1,2), & cdf3d(1,1,ncdflvls,idxrh) ) 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) ) ! delete sfu from netcdf file !PRINT *, ' Storing near-sfc u wind' !CALL cdf_fill(nx,ny,nxcdf,nycdf, uprt(1,1,2), & ! cdf2d(1,1,1,idxsfu) ) ! add write of u wind into 3d array (tdo) PRINT *, ' Storing near-sfc u wind into 3 d array' CALL cdf_fill(nx,ny,nxcdf,nycdf, uprt(1,1,2), & cdf3d(1,1,ncdflvls,idxuw) ) ! delete sfv from netcdf file ! PRINT *, ' Storing near-sfc v wind' ! CALL cdf_fill(nx,ny,nxcdf,nycdf, vprt(1,1,2), & ! cdf2d(1,1,1,idxsfv) ) ! add write of v wind into 3d array (tdo) PRINT *, ' Storing near-sfc v wind into 3 d array' CALL cdf_fill(nx,ny,nxcdf,nycdf, vprt(1,1,2), & cdf3d(1,1,ncdflvls,idxvw) ) ! add write of w wind into 3d array (tdo) PRINT *, ' Storing near-sfc w wind into 3 d array' CALL cdf_fill(nx,ny,nxcdf,nycdf, wprt(1,1,2), & cdf3d(1,1,ncdflvls,idxww) ) ! delete sft from netcdf file ! PRINT *, ' Storing near-sfc temperature' ! CALL cdf_fill(nx,ny,nxcdf,nycdf, tem2(1,1,2), & ! cdf2d(1,1,1,idxsft) ) ! add the storage of surface data into the 3d array (tdo) PRINT *, ' Storing near-sfc temperature into 3d array' CALL cdf_fill(nx,ny,nxcdf,nycdf, tem2(1,1,2), & cdf3d(1,1,ncdflvls,idxt) ) 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) ) ! ! J.Case, ENSOC Inc. (9/29/2004) -- Added 2D radar variables to output. PRINT *, ' Storing Column Max Reflectivity' CALL cdf_fill(nx,ny,nxcdf,nycdf, radar_comp, & cdf2d(1,1,1,idxcxr) ) PRINT *, ' Storing Maximum Radar Echo Tops' CALL cdf_fill(nx,ny,nxcdf,nycdf, echo_top, & cdf2d(1,1,1,idxmret) ) !----------------------------------------------------------------------- ! ! 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) ) ! ! Store cloud water ! CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, & qc,tem3,rlnpcdf,tem1) CALL extrpq(nx,ny,nz,qc,pprt,iprlvl(klev),tem1) CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1), & cdf3d(1,1,klev,idxcw) ) ! ! Store cloud ice ! CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, & qi,tem3,rlnpcdf,tem1) CALL extrpq(nx,ny,nz,qi,pprt,iprlvl(klev),tem1) CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1), & cdf3d(1,1,klev,idxci) ) ! ! Store turbulent kinetic energy ! CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, & tke,tem3,rlnpcdf,tem1) CALL extrpq(nx,ny,nz,tke,pprt,iprlvl(klev),tem1) CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1), & cdf3d(1,1,klev,idxtk) ) ! ! J.Case, ENSCO Inc. (9/29/2004) ! Store 3D radar reflectivity ! CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, & radar_3d,tem3,rlnpcdf,tem1) CALL extrpq(nx,ny,nz,radar_3d,pprt,iprlvl(klev),tem1) CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1), & cdf3d(1,1,klev,idxrr) ) END DO END IF ! Good return from read invlen(1)=ncdflvls invlen(2)=1 vecistr(2)=idxtime volistart(4)=idxtime DO ivar=1,n3dvar WRITE(6,'(a,i4,2x,a)') ' Writing 3d ivar:',ivar, & v3dlngnam(ivar) IF(idxtime == 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)=idxtime DO ivar=1,n2dvar WRITE(6,'(a,i4,2x,a)') ' Writing 2d ivar:',ivar, & v2dlngnam(ivar) IF(adddefs) THEN IF(idxtime == 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 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 IF(adddefs) THEN ! ! 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 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) ! ! Write 1-d variables if idxtime is 1 (i.e. the first history file) ! IF (idxtime == 1) THEN ! ! set up the first element of the time arrays ! if more than 1 more valid time will be written to the file ! foreach time ! set valtime, vtmref, and reftime using ncdtinterval ! write out the inventory arrays with blanks ! AWIPS preprocessor must handle updating inventory arrays if appending ! valtime(idxtime)=FLOAT(itime) vtmref(idxtime)=nint(time) reftime(idxtime)=FLOAT(itimref) IF (ncdntime>idxtime) THEN DO ntime=idxtime+1,ncdntime valtime(ntime)=valtime(idxtime)+FLOAT((ntime-1)*ncdtinterval) vtmref(ntime)=vtmref(idxtime)+(ntime-1)*ncdtinterval reftime(ntime)=reftime(idxtime) ! ! write the 3d var inventory with the fill character ! start array: vecistr(1)=1 vecistr(2)=ntime ! count array: invlen(1)=1 (1 character) invlen(2)=1 ! vecistr(1)=1 vecistr(2)=ntime invlen(1)=1 invlen(2)=1 DO ivar=1,n3dvar istatus=nf_put_vara_text(ncid,v3dinvid(ivar), & vecistr,invlen,nf_fill_char) IF (istatus /= nf_noerr) & CALL handle_err('nf_put_vara_text 3d inv',istatus) END DO ! ! write the 2d var inventory with the fill character ! start array: vecistr(1)=1 vecistr(2)=ntime ! count array: invlen(1)=1 invlen(2)=1 DO ivar=1,n2dvar istatus=nf_put_vara_text(ncid,v2dinvid(ivar), & vecistr,invlen,nf_fill_char) IF (istatus /= nf_noerr) & CALL handle_err('nf_put_vara_text 2d inv',istatus) END DO ! END DO ENDIF istatus=nf_put_vara_int(ncid,vtmrefid, & 1,ncdntime,vtmref) ! istatus=nf_put_vara_double(ncid,reftimid, & ! 1,ncdntime,reftime) istatus=nf_put_vara_double(ncid,reftimid, & 1,1,reftime) istatus=nf_put_vara_double(ncid,valtimid, & 1,ncdntime,valtime) ! ELSE ! WRITE(6,'(a,i4)') ' Writing st ivar:',ivar, & END IF END IF ! end idxtime=1 ntwrt=ntwrt+nhisfile ! ! 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 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE extrph(nx,ny,nz,zps,t,pr,iprlvl,hgtcdf) 3 ! !----------------------------------------------------------------------- ! ! 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, INTENT(IN) :: nx,ny,nz REAL, INTENT(IN) :: zps(nx,ny,nz) REAL, INTENT(IN) :: t(nx,ny,nz) REAL, INTENT(IN) :: pr(nx,ny,nz) INTEGER, INTENT(IN) :: iprlvl REAL, INTENT(OUT) :: hgtcdf(nx,ny) INCLUDE 'phycst.inc' REAL, PARAMETER :: gamma = 0.0065 ! degrees/m lapse rate REAL, PARAMETER :: 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 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE extrpt(nx,ny,nz,t,pr,zps,iprlvl,tcdf) 5 ! !----------------------------------------------------------------------- ! ! 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, INTENT(IN) :: nx,ny,nz REAL, INTENT(IN) :: t(nx,ny,nz) REAL, INTENT(IN) :: pr(nx,ny,nz) REAL, INTENT(IN) :: zps(nx,ny,nz) INTEGER, INTENT(IN) :: iprlvl REAL, INTENT(OUT) :: tcdf(nx,ny) INCLUDE 'phycst.inc' REAL, PARAMETER :: gamma = 0.0065 ! degrees/m lapse rate REAL, PARAMETER :: 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 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE extrpq(nx,ny,nz,q,pr,iprlvl,qcdf) 17 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Assign a zero-gradient value to scalars below ground. ! !----------------------------------------------------------------------- ! ! ! AUTHOR: Keith Brewster, from arps2gem ! 2/19/99 ! ! MODIFICATION HISTORY: ! ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN) :: nx,ny,nz REAL, INTENT(IN) :: q(nx,ny,nz) REAL, INTENT(IN) :: pr(nx,ny,nz) INTEGER, INTENT(IN) :: iprlvl REAL, INTENT(OUT) :: 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 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE extrpuv(nx,ny,nz,us,vs,pr,zps,iprlvl,ucdf,vcdf) 3 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Assign a constant value to sub-terrainian winds ! !----------------------------------------------------------------------- ! ! ! AUTHOR: Keith Brewster, from arps2gem ! 2/19/99 ! ! MODIFICATION HISTORY: ! ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN) :: nx,ny,nz REAL, INTENT(IN) :: us(nx,ny,nz) REAL, INTENT(IN) :: vs(nx,ny,nz) REAL, INTENT(IN) :: pr(nx,ny,nz) REAL, INTENT(IN) :: zps(nx,ny,nz) INTEGER, INTENT(IN) :: iprlvl REAL, INTENT(OUT) :: ucdf(nx,ny) REAL, INTENT(OUT) :: 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 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE handle_err(string,istatus) 151 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Report netCDF error in plain English. Exit. ! ! !----------------------------------------------------------------------- ! ! ! AUTHOR: Keith Brewster ! 2/19/99 ! ! MODIFICATION HISTORY: ! ! !----------------------------------------------------------------------- ! IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: string INTEGER, INTENT(IN) :: istatus CHARACTER (LEN=80) :: nf_strerror WRITE(6,'(a,a,/a)') 'netCDF error: ',string, nf_strerror(istatus) STOP END SUBROUTINE handle_err ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE CDF_FILL ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE cdf_fill(nx,ny,nxcdf,nycdf,arpsvar,cdfvar) ! !----------------------------------------------------------------------- ! ! 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, INTENT(IN) :: nx,ny INTEGER, INTENT(IN) :: nxcdf,nycdf REAL, INTENT(IN) :: arpsvar(nx,ny) REAL, INTENT(OUT) :: 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