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