PROGRAM arpspost,31 ! !################################################################## !################################################################## !###### ###### !###### PROGRAM ARPSPOST ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! ! !----------------------------------------------------------------------- ! PURPOSE: ! Read in a series of ARPS history files, generate 2-D pruducts data ! in ASCII, binary, and/or GEMPAK format (based on input option) ! for display and/or post analysis. The output is in the same grid ! as the input data. ! ! Author : Fanyou Kong ! History: January 30, 2007: First code implementation ! Februry 28, 2007: MPI version implemented ! !----------------------------------------------------------------------- ! ! MODIFIED: ! !----------------------------------------------------------------------- ! 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 (km) ! zp z-coordinate of grid points in computational space (m) ! zpsoil z-coordinate of grid points in computational space (m) ! ! uprt x-component of perturbation velocity (m/s) ! vprt y-component of perturbation velocity (m/s) ! wprt vertical component of perturbation velocity in Cartesian ! coordinates (m/s). ! ! ptprt perturbation potential temperature (K) ! pprt perturbation pressure (Pascal) ! ! qvprt perturbation 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) ! tke Turbulent Kinetic Energy ((m/s)**2) ! kmh Horizontal turb. mixing coef. for momentum ( m**2/s ) ! kmv Vertical turb. mixing coef. for momentum ( m**2/s ) ! ! 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 air density (kg/m**3) ! qvbar Base state water vapor mixing ratio (kg/kg) ! ! soiltyp Soil type ! stypfrct Soil type fraction ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! snowdpth ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain(mm) ! raing Total rain (rainc+raing)(mm) ! ! WORK ARRAYS: ! ! tem1 Temporary work array. ! tem2 Temporary work array. ! tem3 Temporary work array. ! tem4 Temporary work array. ! tem5 Temporary work array. ! tem6 Temporary work array. !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! ! IMPLICIT NONE !----------------------------------------------------------------------- ! Include files: !----------------------------------------------------------------------- INCLUDE 'indtflg.inc' INCLUDE 'globcst.inc' INCLUDE 'grid.inc' ! INCLUDE 'phycst.inc' INCLUDE 'enscv.inc' ! INCLUDE 'GEMINC:GEMPRM.PRM' ! INCLUDE 'GEMPRM.PRM' INCLUDE 'mp.inc' ! !----------------------------------------------------------------------- ! ! Define the dimensions nx, ny, nz, nzsoil ! !----------------------------------------------------------------------- INTEGER :: nx, ny, nz, nzsoil INTEGER :: nstyps ! Maximum number of soil types in each ! grid box !----------------------------------------------------------------------- ! ! 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 on the staggered grid. REAL, ALLOCATABLE :: zp (:,:,:) ! The physical height coordinate defined at ! w-point of the staggered grid. REAL, ALLOCATABLE :: zpsoil(:,:,:) ! The physical height coordinate defined at ! w-point of the soil grid. 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 REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal) REAL, ALLOCATABLE :: qvprt (:,:,:) REAL, ALLOCATABLE :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg) REAL, ALLOCATABLE :: qr (:,:,:) ! Rain water mixing ratio (kg/kg) REAL, ALLOCATABLE :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg) REAL, ALLOCATABLE :: qs (:,:,:) ! Snow mixing ratio (kg/kg) REAL, ALLOCATABLE :: qh (:,:,:) ! Hail mixing ratio (kg/kg) REAL, ALLOCATABLE :: tke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2) REAL, ALLOCATABLE :: kmh (:,:,:) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) 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 density rhobar REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific humidity ! (kg/kg) INTEGER, ALLOCATABLE :: soiltyp(:,:,:) ! Soil type REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type fraction INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! Soil temperature (K) REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! Soil moisture (m**3/m**3) REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m) REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s) REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface REAL, ALLOCATABLE :: radswnet(:,:) ! Net shortwave radiation, SWin - SWout 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)) ! !----------------------------------------------------------------------- ! ! Temporary work arrays for general use ! !----------------------------------------------------------------------- ! REAL, ALLOCATABLE :: tem1(:,:,:) REAL, ALLOCATABLE :: tem2(:,:,:) REAL, ALLOCATABLE :: tem3(:,:,:) REAL, ALLOCATABLE :: tem4(:,:,:) REAL, ALLOCATABLE :: tem5(:,:,:) REAL, ALLOCATABLE :: tem6(:,:,:) ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! CHARACTER (LEN=256) :: filename CHARACTER (LEN=256) :: grdbasfn CHARACTER (LEN=256) :: hisfile(200) INTEGER :: hinfmt,nhisfile,nchin INTEGER :: lengbf,lenfil,lenbin,lengem INTEGER :: ireturn, iret, ierr INTEGER :: i,j,k, itags,itagr REAL :: ctime REAL :: truelat(2) INTEGER :: nf INTEGER :: lenstag, linstit, lmodel ! !----------------------------------------------------------------------- INTEGER :: icape,iaccu,iascii,ibinary,igempak CHARACTER (LEN=256) :: outheader,gemoutheader INTEGER :: nprocx_in, nprocy_in INTEGER :: dmp_out_joined NAMELIST /message_passing/nproc_x, nproc_y, readsplit, nprocx_in, nprocy_in NAMELIST /output_data/dmp_out_joined,outheader,gemoutheader, & icape,iaccu,iascii,ibinary,igempak NAMELIST /output_grid/enstag,instit,model INTEGER :: istatus ! !----------------------------------------------------------------------- ! ! Variables for mpi jobs ! !----------------------------------------------------------------------- INTEGER, PARAMETER :: fzone = 3 INTEGER :: nxlg, nylg ! global domain INTEGER :: ii,jj,ia,ja ! !----------------------------------------------------------------------- ! ! Function f_qvsat and inline directives for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_qvsat !fpp$ expand (f_qvsat) !dir$ inline always f_qvsat ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! CALL mpinit_proc IF(myproc == 0) WRITE(6,'(/ 16(/5x,a)//)') & '###############################################################', & '###############################################################', & '#### ####', & '#### Welcome to ARPSPOST ####', & '#### A post-processing program for ARPS 5.2.6 that reads ####', & '#### in history files generated by ARPS and calculate ####', & '#### and write out required 2D fields in binary and ####', & '#### GEMPAK formats ####', & '#### ####', & '###############################################################', & '###############################################################' !----------------------------------------------------------------------- ! ! First, initialize MPI jobs ! !----------------------------------------------------------------------- IF(myproc == 0) THEN READ(5,message_passing,ERR=100) WRITE(6,'(a)')'Namelist message_passing was successfully read.' GO TO 10 100 WRITE(6,'(a)') & 'Error reading NAMELIST block message_passing. ', & 'Default values used.' 10 CONTINUE END IF CALL mpupdatei(nproc_x,1) CALL mpupdatei(nproc_y,1) CALL mpupdatei(readsplit,1) IF (mp_opt == 0) THEN nproc_x = 1 nproc_y = 1 nprocx_in = 1 nprocy_in = 1 readsplit = 0 nprocs = 1 max_fopen = 1 ELSE max_fopen = nproc_x * nproc_y END IF CALL mpinit_var !----------------------------------------------------------------------- ! ! Read grid/base name, then get the dimensions ! !----------------------------------------------------------------------- IF(myproc == 0) THEN CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile) lengbf = len_trim(grdbasfn) IF(mp_opt > 0 .AND. readsplit <= 0) THEN WRITE(grdbasfn,'(a,a,2i2.2)') grdbasfn(1:lengbf),'_',loc_x,loc_y lengbf = lengbf + 5 END IF WRITE(6,'(a,a/)')' The grid/base name is ', grdbasfn(1:lengbf) CALL get_dims_from_data(hinfmt,trim(grdbasfn), & nx,ny,nz,nzsoil,nstyps, ireturn) IF (mp_opt > 0 .AND. readsplit > 0) THEN ! ! Fiddle with nx/ny, which apparently are wrong. ! nx = (nx - 3) / nproc_x + 3 ny = (ny - 3) / nproc_y + 3 END IF IF (nstyps <= 0) nstyps = 1 nstyp = nstyps ! Copy to global variabl print*,'nstyps =', nstyps IF( ireturn /= 0 ) THEN PRINT*,'Problem occured when trying to get dimensions from data.' PRINT*,'Program stopped.' CALL arpsstop('Problem with data.',1) END IF END IF CALL mpupdatei(hinfmt,1) CALL mpupdatec(grdbasfn,256) CALL mpupdatei(nhisfile,1) CALL mpupdatec(hisfile,256*nhisfile) CALL mpupdatei(nx,1) CALL mpupdatei(ny,1) CALL mpupdatei(nz,1) CALL mpupdatei(nzsoil,1) CALL mpupdatei(nstyps,1) CALL mpupdatei(nstyp,1) IF(myproc == 0) & WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,' nzsoil =',nzsoil !----------------------------------------------------------------------- ! ! Allocate nx, ny, nz dependent arrays and initiliaze to zero ! !----------------------------------------------------------------------- ALLOCATE(x(nx),STAT=istatus) x=0 ALLOCATE(y(ny),STAT=istatus) y=0 ALLOCATE(z(nz),STAT=istatus) z=0 ALLOCATE(zp(nx,ny,nz),STAT=istatus) zp=0 ALLOCATE(zpsoil(nx,ny,nzsoil),STAT=istatus) zpsoil=0 ALLOCATE(ptprt(nx,ny,nz),STAT=istatus) ptprt=0 ALLOCATE(pprt(nx,ny,nz),STAT=istatus) pprt=0 ALLOCATE(qc(nx,ny,nz),STAT=istatus) qc=0 ALLOCATE(qr(nx,ny,nz),STAT=istatus) qr=0 ALLOCATE(qi(nx,ny,nz),STAT=istatus) qi=0 ALLOCATE(qs(nx,ny,nz),STAT=istatus) qs=0 ALLOCATE(qh(nx,ny,nz),STAT=istatus) qh=0 ALLOCATE(tke(nx,ny,nz),STAT=istatus) tke=0 ALLOCATE(kmh(nx,ny,nz),STAT=istatus) kmh=0 ALLOCATE(kmv(nx,ny,nz),STAT=istatus) kmv=0 ALLOCATE(ubar(nx,ny,nz),STAT=istatus) ubar=0 ALLOCATE(vbar(nx,ny,nz),STAT=istatus) vbar=0 ALLOCATE(wbar(nx,ny,nz),STAT=istatus) wbar=0 ALLOCATE(ptbar(nx,ny,nz),STAT=istatus) ptbar=0 ALLOCATE(pbar(nx,ny,nz),STAT=istatus) pbar=0 ALLOCATE(rhobar(nx,ny,nz),STAT=istatus) rhobar=0 ALLOCATE(qvbar(nx,ny,nz),STAT=istatus) qvbar=0 ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus) soiltyp=0 ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus) stypfrct=0 ALLOCATE(vegtyp(nx,ny),STAT=istatus) vegtyp=0 ALLOCATE(lai(nx,ny),STAT=istatus) lai=0 ALLOCATE(roufns(nx,ny),STAT=istatus) roufns=0 ALLOCATE(veg(nx,ny),STAT=istatus) veg=0 ALLOCATE(tsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus) tsoil=0 ALLOCATE(qsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus) qsoil=0 ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus) wetcanp=0 ALLOCATE(snowdpth(nx,ny),STAT=istatus) snowdpth=0 ALLOCATE(raing(nx,ny),STAT=istatus) raing=0 ALLOCATE(rainc(nx,ny),STAT=istatus) rainc=0 ALLOCATE(prcrate(nx,ny,4),STAT=istatus) prcrate=0 ALLOCATE(radfrc(nx,ny,nz),STAT=istatus) radfrc=0 ALLOCATE(radsw(nx,ny),STAT=istatus) radsw=0 ALLOCATE(rnflx(nx,ny),STAT=istatus) rnflx=0 ALLOCATE(radswnet(nx,ny),STAT=istatus) radswnet=0 ALLOCATE(radlwin(nx,ny),STAT=istatus) radlwin=0 ALLOCATE(usflx(nx,ny),STAT=istatus) usflx=0 ALLOCATE(vsflx(nx,ny),STAT=istatus) vsflx=0 ALLOCATE(ptsflx(nx,ny),STAT=istatus) ptsflx=0 ALLOCATE(qvsflx(nx,ny),STAT=istatus) qvsflx=0 ALLOCATE(uprt(nx,ny,nz),STAT=istatus) uprt=0 ALLOCATE(vprt(nx,ny,nz),STAT=istatus) vprt=0 ALLOCATE(wprt(nx,ny,nz),STAT=istatus) wprt=0 ALLOCATE(qvprt(nx,ny,nz),STAT=istatus) qvprt=0 ALLOCATE(tem1(nx,ny,nz),STAT=istatus) tem1=0 ALLOCATE(tem2(nx,ny,nz),STAT=istatus) tem2=0 ALLOCATE(tem3(nx,ny,nz),STAT=istatus) tem3=0 ALLOCATE(tem4(nx,ny,nz),STAT=istatus) tem4=0 ALLOCATE(tem5(nx,ny,nz),STAT=istatus) tem5=0 ALLOCATE(tem6(nx,ny,nz),STAT=istatus) tem6=0 nxlg = (nx-fzone)*nproc_x + fzone nylg = (ny-fzone)*nproc_y + fzone !----------------------------------------------------------------------- ! ! default output grid values (SAMEX settings): ! !----------------------------------------------------------------------- ! IF(myproc == 0) THEN READ(5,output_data,END=200) write(6,output_data) GO TO 20 200 WRITE(6,'(a)') & 'Error reading NAMELIST block output_data. ', & 'Default values used.' 20 CONTINUE END IF CALL mpupdatei(dmp_out_joined,1) CALL mpupdatec(outheader,256) CALL mpupdatec(gemoutheader,256) CALL mpupdatei(icape,1) CALL mpupdatei(iaccu,1) CALL mpupdatei(iascii,1) CALL mpupdatei(ibinary,1) CALL mpupdatei(igempak,1) joindmp = dmp_out_joined IF (mp_opt > 0) THEN ! should moved into mpinit_var later dumpstride = max_fopen IF (joindmp > 0) dumpstride = nprocs ! join and dump END IF enstag = 'ARPSENS' instit = 'CAPS_OU' model = 'ARPS5.0.0' ! Read remaining namelist IF(myproc == 0) THEN READ(5,output_grid,END=250) write(6,output_grid) GO TO 30 250 WRITE(6,'(a)') & 'Error reading NAMELIST block output_grid. ', & 'Default values used.' 30 CONTINUE END IF CALL mpupdatec(enstag,256) CALL mpupdatec(instit,256) CALL mpupdatec(model,256) ! Begin time loop notopn = .true. DO nf=1, nhisfile filename=hisfile(nf) lenfil = len_trim(filename) IF (mp_opt > 0 .AND. readsplit <= 0) THEN WRITE(filename,'(a,i2.2,i2.2)') filename(1:lenfil),'_',loc_x,loc_y lenfil = lenfil +5 END IF IF (myproc == 0) THEN WRITE(6,'(/2a/)') ' Reading file: ',filename(1:lenfil) WRITE(6,'(/2a/)') ' Reading file: ',grdbasfn(1:lengbf) END IF 15 CONTINUE ! also continue to read another time recode ! from GrADS file !----------------------------------------------------------------------- ! ! Read all input data arrays ! !----------------------------------------------------------------------- ! CALL dtaread(nx,ny,nz,nzsoil, nstyps, & hinfmt, nchin,grdbasfn(1:lengbf),lengbf, & filename(1:lenfil),lenfil,ctime, & 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 ! For hinfmt=9, i.e. the GraDs format data, ireturn is used as a ! flag indicating if there is any data at more time level to be read. ! !----------------------------------------------------------------------- ! IF( ireturn /= 0 .AND. hinfmt /= 9 ) GO TO 7000 !unsuccessful read ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! Subroutine call: postcore for calculating and writing out 2D fields ! !----------------------------------------------------------------------- lenbin = len_trim(outheader) lengem = len_trim(gemoutheader) CALL postcore(nx,ny,nz,nzsoil, nstyps,ctime, & filename(1:lenfil),lenfil,outheader(1:lenbin), & lenbin,gemoutheader(1:lengem),lengem, & icape,iaccu,iascii,ibinary,igempak, & x,y,z,zp,uprt ,vprt ,wprt ,ptprt, pprt , & qvprt, qc, qr, qi, qs, qh, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp, & raing,rainc, & tem1,tem2,tem3,tem4,tem5) 7000 CONTINUE ENDDO ! file-loop complete IF(myproc == 0) WRITE(6,'(a)') 'Program completed.' STOP 950 CONTINUE IF(myproc == 0) WRITE(6,'(a)') ' Error setting up GEMPAK file' STOP END PROGRAM arpspost