PROGRAM arpsplt,457
!
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSPLT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
! This is a graphic analysis/plotting program for display ARPS
! history format data set.
!
! It is based on Ming Xue's graphics package ZXPLOT, which is an
! independent package from the ARPS code. The object code library
! of ZXPLOT for most platforms are freely available from
! ftp://ftp.caps.ou.edu/pub/ZXPLOT3. Documentation and other info
! on ZXPLOT can be found at http://www.caps.ou.edu/ZXPLOT.
!
! ZXPLOT can be interfaced with NCAR graphics to produce metafile
! output consistent with NCAR graphics viewing facilities
! or linked to the Postscript driver (pure fortran program) to
! produce Postscript output directly (i.e., no NCAR graphics needed).
!
! Documentation on the control parameters for plotting can be found
! in input/arpsplt.input.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue, CAPS/OU.
! 2/19/1992.
!
! MODIFICATION HISTORY:
! See HISTORY file.
!
!-----------------------------------------------------------------------
!
! 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)
!
! 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
! wetcanp Canopy water amount
! raing Grid supersaturation rain (mm)
! rainc Cumulus convective rain(mm)
! raint Total rain (rainc+raing)(mm)
!
! radfrc Radiation forcing (K/s)
! radsw Solar radiation reaching the surface
! rnflx Net radiation flux absorbed by surface
! radswnet Net shortwave radiation, SWin - SWout
! radlwin Incoming longwave radiation
!
! usflx Surface flux of u-momentum (kg/(m*s**2))
! vsflx Surface flux of v-momentum (kg/(m*s**2))
! ptsflx Surface heat flux (K*kg/(m**2 * s ))
! qvsflx Surface moisture flux of (kg/(m**2 * s))
!
! psl Sea level pressure (mb)
!
! CALCULATED DATA ARRAYS:
!
! u x-component of velocity (m/s)
! v y-component of velocity (m/s)
! w z-component of velocity (m/s)
! pt potential temperature (K)
! qv water vapor mixing ratio (kg/kg)
! td dew-point temperature (C)
! cape CAPE (J/kg)
! cin CIN (J/kg)
! thet theta_E (K)
! heli helicity (m2/s2)
! srlfl storm-relative low-level flow (0-2km AGL)
! srmfl storm-relative mid-level flow (2-9km AGL)
! shr37 7km - 3km wind shear
! ustrm Estimated storm motion (Bob Johns)
! vstrm Estimated storm motion (Bob Johns)
! capst CAPE strength
! blcon boundary layer convergence
! ct convective temperature
!
! 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.
! tem7 Temporary work array.
! tem8 Temporary work array.
! tem9 Temporary work array.
!
! (These arrays are defined and used locally (i.e. inside this
! subroutine), they may also be passed into routines called by
! this one. Exiting the call to this subroutine, these temporary
! work arrays may be used for other purposes, and therefore their
! contents may be overwritten. Please examine the usage of work
! arrays before you make any change to the code.)
!
!-----------------------------------------------------------------------
!
! Arrays for plots on constant pressure levels
!
!-----------------------------------------------------------------------
!
! tz Temperature (K) on computational grids
! t700 Temperature (K) on 700mb pressure grids
! zps3d negative log pressure(Pascal) at ARPS grid points
! algpzc -log(pressure) at scalar grid points
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Grid dimensions.
INTEGER :: nzprofc ! the maximum vertical index in height (zpc/zpsoilc)
! and variables to be profiled when calling vprofil
! subroutine. 06/10/2002, Zuwen He
! In the atmosphere model, the vertical index is
! typically nz-1, while in the soil model, it's nzsoil.
INTEGER :: nzsoil ! levels of soil model
INTEGER :: nstyps ! Maximum number of soil types.
INTEGER :: hinfmt
INTEGER :: nhisfile_max,nhisfile
PARAMETER (nhisfile_max=200)
CHARACTER (LEN=132) :: grdbasfn
CHARACTER (LEN=132) :: hisfile(nhisfile_max)
INTEGER :: lengbf,nf,lenfil
INTEGER, PARAMETER :: max_dim=200
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'indtflg.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.inc'
INCLUDE 'arpsplt.inc'
INCLUDE 'alloc.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 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 staggered grid for soil model.
REAL, ALLOCATABLE :: u (:,:,:) ! Total u-velocity (m/s)
REAL, ALLOCATABLE :: v (:,:,:) ! Total v-velocity (m/s)
REAL, ALLOCATABLE :: w (:,:,:) ! Total w-velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature
! from that of base state atmosphere (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure from that
! of base state atmosphere (Pascal)
REAL, ALLOCATABLE :: qv (:,:,:) ! Water vapor specific humidity (kg/kg)
REAL, ALLOCATABLE :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qh (:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: 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
INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type fraction
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 :: raing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: raint(:,:) ! Total rain (rainc+raing)
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))
!
!-----------------------------------------------------------------------
!
! Arrays derived from the read-in arrays
!
!-----------------------------------------------------------------------
!
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 :: pt (:,:,:) ! Total poten
REAL, ALLOCATABLE :: qvprt (:,:,:)
REAL, ALLOCATABLE :: xc (:,:,:) ! x-coor of scalar point (km)
REAL, ALLOCATABLE :: yc (:,:,:) ! y-coor of scalar point (km)
REAL, ALLOCATABLE :: zc (:,:,:) ! z-coor of scalar point in computational
! space (km)
REAL, ALLOCATABLE :: zpc (:,:,:) ! z-coor of scalar point in physical
! space (km)
REAL, ALLOCATABLE :: zpsoilc (:,:,:) ! zsoil-coor of scalar point in physical
! space (m)
REAL, ALLOCATABLE :: hterain(:,:) ! The height of the terrain.
REAL, ALLOCATABLE :: psl (:,:) ! Sea level pressure (mb)
REAL, ALLOCATABLE :: td (:,:,:) ! dew-point temperature (C)
!
!-----------------------------------------------------------------------
!
! Arrays for plots on constant pressure levels
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tz (:,:,:) ! Temperature (K) on computational grids
REAL, ALLOCATABLE :: t700 (:,:) ! Temperature (K) on 700mb pressure grids
REAL, ALLOCATABLE :: algpzc(:,:,:) ! -log(pressure) at scalar grid points
REAL, ALLOCATABLE :: zps3d(:,:,:) !
REAL, ALLOCATABLE :: zpsoils3d(:,:,:) !
!-----------------------------------------------------------------------
!
! Array for CAPE , CIN
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: wrk1(:),wrk2(:),wrk3(:),wrk4(:),wrk5(:),wrk6(:)
REAL, ALLOCATABLE :: wrk7(:),wrk8(:),wrk9(:),wrk10(:),wrk11(:),wrk12(:)
REAL, ALLOCATABLE :: lcl(:,:) ! The lifting condensation level
REAL, ALLOCATABLE :: lfc(:,:) ! Level of free convection
REAL, ALLOCATABLE :: el(:,:) ! Equilibrium level
REAL, ALLOCATABLE :: twdf(:,:) ! Max wet-bulb potential temperature difference
REAL, ALLOCATABLE :: mbe(:,:) ! Moist Convective Potential Energy
! (MCAPE, includes water loading)
REAL, ALLOCATABLE :: li(:,:) ! Lifted Index (K)
REAL, ALLOCATABLE :: cape(:,:) ! CAPE (J/kg)
REAL, ALLOCATABLE :: cin(:,:) ! CIN (J/kg)
REAL, ALLOCATABLE :: thet(:,:) ! theta_E (K)
REAL, ALLOCATABLE :: heli(:,:) ! helicity
REAL, ALLOCATABLE :: brn(:,:) ! Bulk Richardson Number (Weisman and Klemp)
REAL, ALLOCATABLE :: brnu(:,:) ! Shear parameter of BRN, "U-squared"
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 (Bob Johns)
REAL, ALLOCATABLE :: vstrm(:,:) ! Estimated storm motion (Bob Johns)
REAL, ALLOCATABLE :: capst(:,:) ! cap strength
REAL, ALLOCATABLE :: blcon(:,:) ! boundary llayer convergence
REAL, ALLOCATABLE :: ct(:,:) ! convective temperature
REAL, ALLOCATABLE :: sinlat(:,:) ! Sin of latitude at each grid point
REAL, ALLOCATABLE :: xs(:),ys(:)
REAL :: dxkm, dykm, dzkm
REAL :: dzsoilcm ! Zuwen He, 05/31/2002 in cm.
REAL :: alttostpr
!
!-----------------------------------------------------------------------
!
! Temporary work arrays for general use
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
REAL, ALLOCATABLE :: tem4(:,:,:)
REAL, ALLOCATABLE :: tem5(:,:,:)
REAL, ALLOCATABLE :: tem6(:,:,:)
REAL, ALLOCATABLE :: tem7(:,:,:)
REAL, ALLOCATABLE :: tem8(:,:,:)
REAL, ALLOCATABLE :: tem9(:,:,:)
!
!-----------------------------------------------------------------------
!
! Work arrays used by profile plotting.
!
!-----------------------------------------------------------------------
!
REAL :: xprof(max_dim)
REAL :: yprof(max_dim)
!
!-----------------------------------------------------------------------
!
! Work arrays to be used in interpolation subroutines
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: b1(:,:),b2(:,:)
REAL, ALLOCATABLE :: u1(:,:),v1(:,:)
REAL, ALLOCATABLE :: u2(:,:),v2(:,:),w2(:,:),zs2(:,:)
REAL, ALLOCATABLE :: xp(:),yp(:)
!
!-----------------------------------------------------------------------
!
! Some universal constants
!
!-----------------------------------------------------------------------
!
REAL :: kappa,gamma,ex1,ex2,t00,p00,mbtopa
PARAMETER(kappa=287.053/cp, &
gamma=6.5, & ! 6.5 K/km
ex1=0.1903643, & ! R*gamma/g
ex2=5.2558774, & ! g/R/gamma
mbtopa=100.)
!-----------------------------------------------------------------------
!
! Plotting control parameters entered by user
!
!-----------------------------------------------------------------------
!
INTEGER :: nchin ! input file unit number
INTEGER :: n, iorig
INTEGER :: vtrplt, vtpplt, vagplt,vtrstrm,vtpstrm, xuvplt, strmplt
INTEGER :: uplot, vplot, wplot, ptplot, pplot, hplot, tplot, vhplot, &
vsplot
INTEGER :: qvplot,qcplot,qrplot,qiplot,qsplot,qhplot,kmhplt, &
kmvplt,tkeplt,rhplot,rfplot,pteplt,tdplot,qwplot, &
rfcplt,qtplot, rhiplot
INTEGER :: rfopt
INTEGER :: upplot,vpplot,wpplot,ptpplt,ppplot,qvpplt,vorpplt,divpplt
INTEGER :: divqplt,gricplt, avorplt
INTEGER :: viqcplt,viqrplt,viqiplt,viqhplt,viqsplt,vilplt,viiplt, &
vicplt, vitplt, pwplt,tprplt, gprplt, cprplt
!
!-----------------------------------------------------------------------
!
! Contour intervals
!
!-----------------------------------------------------------------------
!
REAL :: vtrunit, vtpunt, vagunit,xuvunit, strmunit
REAL :: uinc, vinc, winc, ptinc, pinc, hinc, tinc, vhinc,vsinc
REAL :: qvinc,qcinc,qrinc,qiinc,qsinc,qhinc,qwinc,qtinc
REAL :: kmhinc,kmvinc,tkeinc,rhinc,rfinc,pteinc,tdinc,rfcinc, rhiinc
REAL :: upinc,vpinc,wpinc,ptpinc,ppinc,qvpinc,vorpinc,divpinc,divqinc
REAL :: gricinc,avorinc
REAL :: viqcinc,viqrinc,viqiinc,viqhinc,viqsinc,vilinc,viiinc,vicinc
REAL :: vitinc, pwinc,tprinc, gprinc, cprinc
!
!-----------------------------------------------------------------------
!
! Limited variable minimum and maximum for color contour shade
!
!-----------------------------------------------------------------------
!
REAL :: wcpminc, wcpmaxc, trnminc, trnmaxc
REAL :: raincminc, raincmaxc, raingminc, raingmaxc, &
raintminc,raintmaxc
REAL :: capeminc, capemaxc, cinminc, cinmaxc, thetminc, thetmaxc, &
heliminc, helimaxc, capsminc, capsmaxc, blcominc, blcomaxc, &
ctcminc, ctcmaxc
REAL :: brnminc, brnmaxc, bruminc, brumaxc, srlminc, srlmaxc, &
srmminc, srmmaxc, liminc, limaxc
REAL :: uminc, umaxc, vminc, vmaxc, wminc, wmaxc,ptminc, ptmaxc
REAL :: vhminc, vhmaxc,vsminc, vsmaxc
REAL :: pminc, pmaxc, qvminc, qvmaxc, qcminc, qcmaxc
REAL :: qrminc, qrmaxc, qiminc, qimaxc, qsminc, qsmaxc
REAL :: qwminc, qwmaxc, qtminc, qtmaxc
REAL :: qhminc, qhmaxc, rhminc, rhmaxc, rhiminc, rhimaxc
REAL :: kmhminc, kmhmaxc,kmvminc, kmvmaxc,tkeminc, tkemaxc
REAL :: rfminc, rfmaxc, pteminc, ptemaxc, upminc, upmaxc
REAL :: rfcminc, rfcmaxc
REAL :: vpminc, vpmaxc, wpminc, wpmaxc, ptpminc, ptpmaxc
REAL :: ppminc, ppmaxc, qvpminc, qvpmaxc, vorpminc, vorpmaxc
REAL :: divpminc, divpmaxc, divqminc, divqmaxc
REAL :: gricminc, gricmaxc,avorminc, avormaxc
REAL :: hminc, hmaxc, tminc, tmaxc,pslmaxc,pslminc
REAL :: tdminc, tdmaxc
REAL :: soiltpminc,soiltpmaxc,vegtpminc,vegtpmaxc,laiminc,laimaxc, &
rouminc,roumaxc,vegminc,vegmaxc,snowdminc,snowdmaxc
REAL :: viqcminc,viqrminc,viqiminc,viqhminc,viqsminc,vilminc,viiminc, &
vicminc, vitminc, pwminc, tprminc, gprminc, cprminc
REAL :: viqcmaxc,viqrmaxc,viqimaxc,viqhmaxc,viqsmaxc,vilmaxc,viimaxc, &
vicmaxc, vitmaxc, pwmaxc, tprmaxc, gprmaxc, cprmaxc
!-----------------------------------------------------------------------
!
! Overlay control parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: vtrovr,vtpovr,vagovr,vtrstmovr,vtpstmovr,xuvovr,strmovr
INTEGER :: uovr , vovr , wovr , ptovr , povr, hovr, tovr,vhovr,vsovr
INTEGER :: qvovr ,qcovr ,qrovr ,qiovr ,qsovr ,qhovr ,kmhovr , &
kmvovr,tkeovr,rhovr ,rfovr ,pteovr,tdovr, qwovr, qtovr, &
rfcovr, rhiovr
INTEGER :: upovr ,vpovr ,wpovr ,ptpovr,ppovr ,qvpovr,vorpovr,divpovr
INTEGER :: trnovr,wcpovr,racovr,ragovr, &
ratovr,pslovr,capovr,cinovr,theovr,helovr,brnovr,brnuovr, &
srlfovr,srmfovr,liovr,capsovr,blcoovr,ctcovr
INTEGER :: divqovr,gricovr,avorovr
INTEGER :: styovr,vtyovr,laiovr,rouovr,vegovr,snowdovr
INTEGER :: viqcovr,viqrovr,viqiovr,viqhovr,viqsovr,vilovr,viiovr, &
vicovr, vitovr, pwovr, tprovr, gprovr, cprovr
!
!-----------------------------------------------------------------------
!
! highlighting frequency for contour parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: uhlf , vhlf , whlf , pthlf , phlf, hhlf, thlf,vhhlf,vshlf
INTEGER :: qvhlf ,qchlf ,qrhlf ,qihlf ,qshlf ,qhhlf ,kmhhlf , &
kmvhlf,tkehlf,rhhlf ,rfhlf ,ptehlf,tdhlf, qwhlf, qthlf, &
rfchlf, rhihlf
INTEGER :: uphlf ,vphlf ,wphlf ,ptphlf,pphlf ,qvphlf,vorphlf,divphlf
INTEGER :: trnhlf,wcphlf,rachlf,raghlf, &
rathlf,pslhlf,caphlf,cinhlf,thehlf,helhlf,brnhlf,brnuhlf, &
srlfhlf,srmfhlf,lihlf,capshlf,blcohlf,ctchlf
INTEGER :: divqhlf,grichlf,avorhlf
INTEGER :: styhlf,vtyhlf,laihlf,rouhlf,veghlf,snowdhlf
INTEGER :: viqchlf,viqrhlf,viqihlf,viqhhlf,viqshlf,vilhlf,viihlf, &
vichlf, vithlf, pwhlf,tprhlf, gprhlf, cprhlf
!
!-----------------------------------------------------------------------
!
! define the attributes of zero contour to be plotted parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: uzro , vzro , wzro , ptzro , pzro, hzro, tzro,vhzro,vszro
INTEGER :: qvzro ,qczro ,qrzro ,qizro ,qszro ,qhzro ,kmhzro , &
kmvzro,tkezro,rhzro ,rfzro ,ptezro,tdzro, qwzro, qtzro, &
rfczro, rhizro
INTEGER :: upzro ,vpzro ,wpzro ,ptpzro,ppzro ,qvpzro,vorpzro,divpzro
INTEGER :: trnzro,wcpzro,raczro,ragzro, &
ratzro,pslzro,capzro,cinzro,thezro,helzro,brnzro,brnuzro, &
srlfzro,srmfzro,lizro,capszro,blcozro,ctczro
INTEGER :: divqzro,griczro,avorzro
INTEGER :: styzro,vtyzro,laizro,rouzro,vegzro,snowdzro
INTEGER :: viqczro,viqrzro,viqizro,viqhzro,viqszro,vilzro,viizro, &
viczro, vitzro, pwzro, tprzro, gprzro, cprzro
!
!-----------------------------------------------------------------------
!
! Define the option for contour line stypes.
!
!-----------------------------------------------------------------------
!
INTEGER :: usty , vsty , wsty , ptsty , psty, hsty, tsty,vhsty,vssty
INTEGER :: qvsty ,qcsty ,qrsty ,qisty ,qssty ,qhsty ,kmhsty , &
kmvsty,tkesty,rhsty ,rfsty ,ptesty,tdsty, qwsty, qtsty, &
rfcsty, rhisty
INTEGER :: upsty ,vpsty ,wpsty ,ptpsty,ppsty ,qvpsty,vorpsty,divpsty
INTEGER :: trnsty,wcpsty,racsty,ragsty, &
ratsty,pslsty,capsty,cinsty,thesty,helsty,brnsty,brnusty, &
srlfsty,srmfsty,listy,capssty,blcosty,ctcsty
INTEGER :: divqsty,gricsty,avorsty
INTEGER :: stysty,vtysty,laisty,rousty,vegsty,snowdsty
INTEGER :: viqcsty,viqrsty,viqisty,viqhsty,viqssty,vilsty,viisty, &
vicsty, vitsty, pwsty, tprsty, gprsty, cprsty
INTEGER :: msfplt,msfovr,msfcol1,msfcol2,msfprio,msfhlf,msfzro,msfsty
REAL :: msfinc, msfminc, msfmaxc
INTEGER :: thkplt,thkovr,thkcol1,thkcol2,thkprio,thkhlf,thkzro,thksty
REAL :: thkinc, thkminc, thkmaxc
INTEGER :: ipvplt,ipvovr,ipvcol1,ipvcol2,ipvprio,ipvhlf,ipvzro,ipvsty
REAL :: ipvinc, ipvminc, ipvmaxc
!
!-----------------------------------------------------------------------
!
! Profile Plotting control parameters entered by user
!
!-----------------------------------------------------------------------
!
INTEGER :: uprof, vprof, wprof, ptprof, pprof
INTEGER :: qvprof,qcprof,qrprof,qiprof,qsprof,qhprof,kmhprof, &
kmvprof,tkeprof,rhprof,rfprof,pteprf
INTEGER :: upprof,vpprof,wpprof,ptpprf,ppprof,qvpprf,vorpprf,divpprf
INTEGER :: tsoilprof,qsoilprof ! Zuwen He 05/31/2002
!
!-----------------------------------------------------------------------
!
! Profile plot lower bound
!
!-----------------------------------------------------------------------
!
REAL :: uprmin, vprmin, wprmin, ptprmin, pprmin
REAL :: qvprmin,qcpmin,qrpmin,qipmin,qspmin,qhpmin
REAL :: kmhpmin,kmvpmin,tkepmin
REAL :: rhpmin,rfpmin,ptepmin,vorppmin,divppmin,avormin
REAL :: uppmin,vppmin,wppmin,ptppmin,pppmin,qvppmin
REAL :: tsoilprofmin,qsoilprofmin ! Zuwen He 06/03/2002
!
!-----------------------------------------------------------------------
!
! Profile plot upper bound
!
!-----------------------------------------------------------------------
!
REAL :: uprmax, vprmax, wprmax, ptprmax, pprmax
REAL :: qvprmax,qcpmax,qrpmax,qipmax,qspmax,qhpmax
REAL :: kmhpmax, kmvpmax, tkepmax
REAL :: rhpmax,rfpmax,ptepmax,vorppmax,divppmax, avormax
REAL :: uppmax,vppmax,wppmax,ptppmax,pppmax,qvppmax
REAL :: tsoilprofmax,qsoilprofmax ! Zuwen He 06/03/2002
!
!-----------------------------------------------------------------------
!
! 3-D wireframe plotting
!
!-----------------------------------------------------------------------
!
INTEGER :: w3dplt, q3dplt
REAL :: wisosf,qisosf
INTEGER :: idisplay
INTEGER :: imove, inwfrm
!
!-----------------------------------------------------------------------
!
! Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: layover ! set by OVERLAY, used in ctr2d,vtr2d
COMMON /laypar/ layover
!
REAL :: ctinc,ctmin,ctmax,vtunt ! contour interval and vector unit
COMMON /incunt/ ctinc,ctmin,ctmax,vtunt
!
INTEGER :: icolor,icolor1,lbcolor,trcolor ! required color
COMMON /recolor/icolor,icolor1,lbcolor,trcolor
!
CHARACTER (LEN=12) :: varname ! save the plot variable
COMMON /varplt1/ varname
!
REAL :: x01,y01 ! first interpolation point for a
! vertical slice specified by two pts
REAL :: x02,y02 ! second interpolation point for a
! vertical slice specified by two pts
REAL :: zlevel ! altitude (meters) of a horizontal
! slice so specified
REAL :: sinaf,cosaf,dist,sqrtdxy
COMMON /slicev/x01,y01,x02,y02,sinaf,cosaf,dist,sqrtdxy
COMMON /sliceh/zlevel
REAL :: yxstrch ! Stretching factor for x-y plots
REAL :: zxstrch ! Stretching factor for x-z plots
REAL :: zystrch ! Stretching factor for y-z plots
REAL :: zhstrch ! Stretching factor for arbitrary vertical slices
REAL :: zsoilxstrch ! Stretching factor for x-z plots for the soil model
REAL :: zsoilystrch ! Stretching factor for y-z plots for the soil model
REAL :: aspratio
COMMON /yratio/ aspratio ! scaling factor: the y/x ratio.
!
! Pass the plotting window parameters into subroutine encodwd.
!
COMMON /pltwdw/ xbgn,xend,ybgn,yend
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: label
CHARACTER (LEN=132) :: filename
!
!-----------------------------------------------------------------------
!
INTEGER :: length
! character*12 pltvar
CHARACTER (LEN=6) :: stem2
CHARACTER (LEN=1) :: stem1
INTEGER :: flagsin
!wdt update
REAL :: f_cputime,cpu0,cpu1,cpu2
DOUBLE PRECISION :: f_walltime,second1, second2
INTEGER :: ireturn
INTEGER :: slicopt
INTEGER :: i,j,k,ibgn,iend,jbgn,jend,kbgn,kend,ist,jst,kst,iob
INTEGER :: ksoilbgn,ksoilend ! Zuwen He, 05/31/2002
INTEGER :: ibgn1, iend1
INTEGER :: nxpic,nypic,islice,jslice,kslice,layout
INTEGER :: isoilslice,jsoilslice,ksoilslice ! Zuwen 05/31/2002
REAL :: time,angl,aspect
REAL :: xr,yr,zr,x1,x2,y1,y2,z1,z2,xlimit,ylimit
REAL :: zmax,zmin
REAL :: zsoilr,zsoil1,zsoil2
REAL :: zsoilmax,zsoilmin ! Zuwen He, 05/31/2002
REAL :: pmb,drot
INTEGER :: imax,jmax,kmax,imin,jmin,kmin
REAL :: wmin, wmax, xorold, yorold, tk, tdk, psta
REAL :: xbgn,xend,ybgn,yend,zbgn,zend
REAL :: zsoilbgn,zsoilend ! Zuwen He, 05/31/2002 <0
LOGICAL :: fexist
INTEGER :: onvf
REAL :: utmax,utmin,vtmax,vtmin,wtmax,wtmin,vsmax,vsmin
REAL :: qvmax,qvmin,qcmax,qcmin,qrmax,qrmin,qimax,qimin,qsmax,qsmin
REAL :: qwmax,qwmin
REAL :: qhmax,qhmin,rhmax,rhmin,rfmax,rfmin,ptemax,ptemin,upmax,upmin
REAL :: vpmax,vpmin,wpmax,wpmin,ptpmax,ptpmin,ppmax,ppmin
REAL :: qvpmax,qvpmin,vormax,vormin,divmax,divmin
REAL :: hmin, hmax, tdmin,tdmax, rhimin, rhimax
INTEGER :: nprof, profopt, nxprpic, nyprpic, npicprof
REAL :: zprofbgn, zprofend
REAL :: zsoilprofbgn, zsoilprofend ! Zuwen He, 05/31/2002
REAL :: time0
INTEGER :: trnplt,wetcanplt
INTEGER :: soiltpplt,vegtpplt,laiplt,rouplt,vegplt,snowdplt
INTEGER :: soiltpn ! number of soil type 1 to 4
INTEGER :: pslplt
INTEGER :: raincplt,raingplt,raintplt
INTEGER :: capeplt, cinplt, thetplt, heliplt
INTEGER :: brnplt, brnuplt, srlfplt, srmfplt, liplt
INTEGER :: capsplt, blcoplt, ctcplt
INTEGER :: ip,ipp,ipriority(nprio),iptemp(nprio),sigplt(nprio)
REAL :: trninc,wcpinc
REAL :: soiltpinc,vegtpinc,laiinc,rouinc,veginc,snowdinc
REAL :: pslinc
REAL :: raincinc, rainginc,raintinc
REAL :: capeinc, cininc, thetinc, heliinc
REAL :: brninc, brnuinc, srlfinc, srmfinc, liinc
REAL :: capsinc, blcoinc, ctcinc
!
! 05/30/2002 Zuwen He
!
! soil plot options and parameters
!
INTEGER :: tsoilplt
REAL :: tsoilinc,tsoilminc,tsoilmaxc
INTEGER :: tsoilovr,tsoilhlf,tsoilzro,tsoilcol1,tsoilcol2,tsoilprio
INTEGER :: qsoilplt
REAL :: qsoilinc,qsoilminc,qsoilmaxc
INTEGER :: qsoilovr,qsoilhlf,qsoilzro,qsoilcol1,qsoilcol2,qsoilprio
!
! END Zuwen He
!
INTEGER :: ovrmap,mapgrid,mapcol(maxmap),mapline_style(maxmap), &
mapgridcol,nmapfile
REAL :: latgrid,longrid
CHARACTER (LEN=132) :: mapfile(maxmap)
INTEGER :: lmapfile
COMMON /mappar / ovrmap
COMMON /mappar1/ nmapfile,mapcol,mapline_style,mapfile
COMMON /mappar2/ mapgrid,mapgridcol, latgrid,longrid
REAL :: ztmin,ztmax
INTEGER :: ovrtrn
COMMON /trnpar/trnplt,ovrtrn,trninc,trnminc,trnmaxc, &
ztmin,ztmax
INTEGER :: ovrobs,obsset,obscol,obs_marktyp
REAL :: obs_marksz
COMMON /obspar/ ovrobs,obsset,obscol,obs_marktyp, obs_marksz
CHARACTER (LEN=132) :: sfcobfl
INTEGER :: obunit,lsfcobfl
INTEGER :: ovrstaopt
INTEGER :: ovrstam,staset,ovrstan,ovrstav,stacol,markprio,wrtstax
INTEGER :: nsta_typ,sta_typ(10),sta_marktyp(10),sta_markcol(10)
REAL :: sta_marksz(10)
REAL :: wrtstad
COMMON /sta_par/ ovrstaopt,ovrstam,staset,ovrstan,ovrstav,stacol, &
markprio,nsta_typ,sta_typ,sta_marktyp, &
sta_markcol,sta_marksz,wrtstax,wrtstad
CHARACTER (LEN=132) :: stalofl
INTEGER :: stunit,lstalofl
!
!-----------------------------------------------------------------------
!
! Surface (single-level) read-in observation variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=24) :: atime
CHARACTER (LEN=5) :: stn(mxsfcob)
INTEGER :: n_meso_g, &
n_meso_pos,n_sao_g,n_sao_pos_g,n_sao_b,n_sao_pos_b, &
n_obs_g,n_obs_pos_g,n_obs_b,n_obs_pos_b,nobs
CHARACTER (LEN=8) :: obstype(mxsfcob)
CHARACTER (LEN=8) :: wx(mxsfcob)
CHARACTER (LEN=1) :: store_emv(mxsfcob,5)
CHARACTER (LEN=4) :: store_amt(mxsfcob,5)
REAL :: latob(mxsfcob),lonob(mxsfcob),elevob(mxsfcob)
REAL :: tob(mxsfcob),tdob(mxsfcob),dd(mxsfcob),ff(mxsfcob)
REAL :: ddg(mxsfcob),ffg(mxsfcob)
REAL :: pstn(mxsfcob),pmsl(mxsfcob),alt(mxsfcob)
REAL :: ceil(mxsfcob),lowcld(mxsfcob),cover(mxsfcob)
REAL :: vis(mxsfcob),rad(mxsfcob)
REAL :: store_hgt(mxsfcob,5)
INTEGER :: kloud(mxsfcob),idp3(mxsfcob)
INTEGER :: obstime(mxsfcob)
REAL :: obs1(mxsfcob),obs2(mxsfcob)
!
COMMON /sfc_obs1/ nobs
COMMON /sfc_obs2/ latob,lonob,obs1,obs2
!
INTEGER :: nsta,nstapro(mxstalo),nstatyp(mxstalo)
REAL :: latsta(mxstalo), lonsta(mxstalo)
CHARACTER (LEN=2 ) :: s_state(mxstalo)
CHARACTER (LEN=5 ) :: s_name(mxstalo)
CHARACTER (LEN=20) :: s_site(mxstalo)
INTEGER :: s_elev(mxstalo)
COMMON /sta_loc/latsta,lonsta,nstatyp,nstapro,nsta
COMMON /sta_loc1/s_name
!
REAL :: lblmag ! A global magnification factor for labels.
REAL :: ctrlbsiz, axlbsiz
COMMON /labmag/ lblmag, ctrlbsiz, axlbsiz
REAL :: winsiz ! A global factor for window size
COMMON /windows/ winsiz
REAL :: margnx, margny ! margin
INTEGER :: pcolbar ! position of color bar
INTEGER :: tickopt
INTEGER :: axlbfmt ! format for the axis's label
INTEGER :: fontopt ! the font of character
INTEGER :: ctrlbopt ! Contour labelling option
INTEGER :: ctrstyle
INTEGER :: ctrlbfrq
COMMON /clb_frq/ ctrlbfrq
INTEGER :: lbmaskopt
INTEGER :: haxisu, vaxisu
INTEGER :: lbaxis
INTEGER :: lnmag
REAL :: hmintick,vmajtick,vmintick,hmajtick
COMMON /var_par/ fontopt,haxisu, vaxisu,lbaxis,tickopt, &
hmintick,vmajtick,vmintick,hmajtick,axlbfmt
INTEGER :: presaxis_no
REAL :: pres_val(20), pres_z(20)
COMMON /pressbar_par/presaxis_no,pres_val,pres_z
!
INTEGER :: ntitle,titcol, wpltime
REAL :: titsiz
CHARACTER (LEN=132) :: title(3), footer_l, footer_c, footer_r
COMMON /titpar1/title, footer_l, footer_c, footer_r
COMMON /titpar2/ntitle,titcol,wpltime, nxpic, nypic
COMMON /titpar3/titsiz
COMMON /tmphc1/x1,x2,y1,y2,z1,z2
!
! Work arrays used by contouring and contour color fill routines
! m*n=<50000, 8*m=<m if the array to be contours is demensioned m x n.
!
INTEGER :: sovrlay
LOGICAL :: plotovr
!
!-----------------------------------------------------------------------
!
! Common block used by the Ncar Graphic streamline routine
! INITA Used to precondition grid boxes to be eligible to
! start a streamline. For example, a value of 4
! means that every fourth grid box is eligible ;
! a value of 2 means that every other grid box is eligible.
! (see INITB)
!
! INITB Used to precondition grid boxes to be eligible for
! direction arrows. If the user changes the default values
! of INITA and/or INITB, it should be done such that
! MOD(INITA,INITB) = 0. For a dense grid try INITA=4 and
! INITB=2 to reduce the CPU time.
!
!-----------------------------------------------------------------------
!
INTEGER :: inita , initb , iterp , iterc , igflg, imsg , icyc
REAL :: arowl , uvmsg , displ , dispc , cstop
COMMON /str03/ inita , initb , arowl , iterp , iterc , igflg &
, imsg , uvmsg , icyc , displ , dispc , cstop
!-----------------------------------------------------------------------
!
! Namelist for plot input
!
!-----------------------------------------------------------------------
!
INTEGER :: maxslicopt
PARAMETER (maxslicopt=11)
!
! slicopt = 1: xy slice
! slicopt = 2: xz slice
! slicopt = 3: yz slice
! slicopt = 4: constant height
! slicopt = 5: vertical slice through two points
! slicopt = 6: constant pressure
! slicopt = 7: isentropic level
! slicopt = 8: ??
! slicopt = 9: xy-soil-slice for the soil model
! slicopt = 10: xz-soil-slice for the soil model
! slicopt = 11: yz-soil-slice for the soil model
INTEGER :: nslice(maxslicopt), indxslic
INTEGER :: iplot(maxslicopt)
INTEGER :: nslice_xy,nslice_xz, nslice_yz, nslice_h, nslice_v, &
nslice_p,nslice_pt
INTEGER :: slice_xy(max_dim),slice_xz(max_dim), slice_yz(max_dim)
!
! 05/30/2002 Zuwen He
! Slice for soil variables
!
INTEGER :: nslice_xy_soil
INTEGER :: slice_xy_soil(max_dim)
INTEGER :: nslice_xz_soil
INTEGER :: slice_xz_soil(max_dim)
INTEGER :: nslice_yz_soil
INTEGER :: slice_yz_soil(max_dim)
!
! END Zuwen
!
REAL :: xi1(max_dim),yi1(max_dim),xi2(max_dim),yi2(max_dim)
REAL :: slice_h(max_dim), slice_p(max_dim), slice_pt(max_dim)
REAL :: xpnt1(max_dim),ypnt1(max_dim), xpnt2(max_dim),ypnt2(max_dim)
INTEGER :: ucol1,vcol1,wcol1,ptcol1,pcol1,qvcol1,qccol1, &
qrcol1,qicol1,qscol1,qhcol1,kmhcol1,kmvcol1,tkecol1, &
rhcol1,rfcol1,rfccol1,ptecol1,upcol1,vpcol1,wpcol1, &
ptpcol1,ppcol1,qvpcol1,vorpcol1,divpcol1,vtrcol1, &
vtpcol1,vagcol1,vtrstmcol1,vtpstmcol1,hcol1,tcol1, &
divqcol1, &
tdcol1,qwcol1,qtcol1,trncol1, &
wcpcol1,raccol1,ragcol1,ratcol1,pslcol1,stycol1, &
vtycol1,laicol1,roucol1,vegcol1,snowdcol1, &
capcol1,cincol1,thecol1, &
helcol1, vhcol1, vscol1,brncol1,brnucol1,srlfcol1, &
srmfcol1,licol1, capscol1, blcocol1,griccol1,avorcol1, &
viqccol1,viqrcol1,viqicol1,viqhcol1,viqscol1,vilcol1, &
viicol1,viccol1, xuvcol1, strmcol1, ctccol1, rhicol1, &
vitcol1, pwcol1, tprcol1, gprcol1, cprcol1
INTEGER :: ucol2,vcol2,wcol2,ptcol2,pcol2,qvcol2,qccol2, &
qrcol2,qicol2,qscol2,qhcol2,kmhcol2,kmvcol2,tkecol2, &
rhcol2,rfcol2,rfccol2,ptecol2,upcol2,vpcol2,wpcol2, &
ptpcol2,ppcol2,qvpcol2,vorpcol2,divpcol2,vtrcol2, &
vtpcol2,vagcol2,vtrstmcol2,vtpstmcol2,hcol2,tcol2, &
divqcol2, &
tdcol2,qwcol2,qtcol2,trncol2, &
wcpcol2,raccol2,ragcol2,ratcol2,pslcol2,stycol2, &
vtycol2,laicol2,roucol2,vegcol2,snowdcol2, &
capcol2,cincol2,thecol2, &
helcol2, vhcol2, vscol2,brncol2,brnucol2,srlfcol2, &
srmfcol2,licol2, capscol2, blcocol2,griccol2,avorcol2, &
viqccol2,viqrcol2,viqicol2,viqhcol2,viqscol2,vilcol2, &
viicol2,viccol2, xuvcol2, strmcol2, ctccol2, rhicol2, &
vitcol2, pwcol2, tprcol2, gprcol2, cprcol2
INTEGER :: uprio,vprio,wprio,ptprio,pprio,qvprio,qcprio, &
qrprio,qiprio,qsprio,qhprio,kmhprio,kmvprio,tkeprio, &
rhprio,rfprio,rfcprio,rhiprio, &
pteprio,upprio,vpprio,wpprio,ptpprio,ppprio,qvpprio, &
vorpprio,divpprio,vtrprio,vtpprio,vagprio,vtrstmprio, &
vtpstmprio,hprio,tprio,divqprio,tdprio,qwprio,qtprio, &
trnprio,wcpprio,racprio, &
ragprio,ratprio,pslprio, &
styprio,vtyprio,laiprio,rouprio,vegprio,snowdprio, &
capprio,cinprio,theprio,helprio, vhprio,vsprio, &
brnprio,bruprio,srlprio,srmprio,liprio, &
capsprio,blcoprio,gricprio,avorprio, &
viqcprio,viqrprio,viqiprio,viqhprio,viqsprio,vilprio, &
viiprio,vicprio, xuvprio, strmprio, ctcprio, vitprio, &
pwprio, tprprio, gprprio, cprprio
INTEGER :: col_table,smooth
CHARACTER (LEN=80) :: color_map
INTEGER :: number_of_boxes, boxcol
REAL :: bctrx(10),bctry(10),blengx(10),blengy(10)
REAL :: bx1(10), bx2(10),by1(10),by2(10)
INTEGER :: number_of_polys, polycol
REAL :: vertx(max_verts,max_polys), verty(max_verts,max_polys)
INTEGER :: istride,jstride,kstride
COMMON /coltable/col_table,pcolbar
COMMON /smoothopt/smooth
COMMON /boxesopt/number_of_boxes,boxcol,bx1,bx2,by1,by2
COMMON /polysopt/number_of_polys,polycol,vertx,verty
CHARACTER (LEN=1) :: tunits ! units for temperature F or C
CHARACTER (LEN=1) :: tdunits ! units for dew-point temp F or C
INTEGER :: vhunits,vtrunits,vtpunits,vagunits,xuvunits,strmunits
INTEGER :: vtrtype,vtptype,vagtype,xuvtype, strmtype
INTEGER :: tprunits, gprunits, cprunits
INTEGER :: iunits, itype
COMMON /windvtr/iunits, itype
INTEGER :: racunit, ragunit, ratunit
INTEGER :: ovrlaymulopt, ovrmul_num
CHARACTER (LEN=12) :: ovrmulname(50)
CHARACTER (LEN=12) :: ovrname
INTEGER :: setcontopt ,setcontnum
CHARACTER (LEN=12) :: setcontvar(maxuneva)
REAL :: setconts(maxunevm,maxuneva)
COMMON /setcont_var/setcontvar
COMMON /setcon_par/setcontopt,setcontnum,setconts
NAMELIST /col_table_cntl/ col_table,color_map
NAMELIST /smooth_cntl/ smooth
INTEGER :: arbvaropt ! plot arbitrary variable
INTEGER :: istatus
CHARACTER (LEN=40) :: dirname3d(maxarbvar),dirname2d(maxarbvar)
CHARACTER (LEN=40) :: varunits, var_name
CHARACTER (LEN=6) :: var3d(maxarbvar)
INTEGER :: var3dnum
INTEGER :: var3dplot(maxarbvar)
REAL :: var3dinc(maxarbvar), var3dminc(maxarbvar), &
var3dmaxc(maxarbvar)
INTEGER :: var3dovr(maxarbvar),var3dcol1(maxarbvar), &
var3dcol2(maxarbvar),var3dprio(maxarbvar), &
var3dhlf(maxarbvar),var3dzro(maxarbvar), &
var3dsty(maxarbvar)
REAL, ALLOCATABLE :: var3dv(:,:,:)
CHARACTER (LEN=6) :: var2d(maxarbvar)
INTEGER :: var2dnum
INTEGER :: var2dplot(maxarbvar)
REAL :: var2dinc(maxarbvar), var2dminc(maxarbvar), &
var2dmaxc(maxarbvar)
INTEGER :: var2dovr(maxarbvar),var2dcol1(maxarbvar), &
var2dcol2(maxarbvar), var2dprio(maxarbvar), &
var2dhlf(maxarbvar),var2dzro(maxarbvar), &
var2dsty(maxarbvar)
REAL, ALLOCATABLE :: var2dv(:,:)
INTEGER :: missfill_opt,missval_colind ! miss value color index
COMMON /multi_value/ missfill_opt, missval_colind
INTEGER :: xnwpic_called
COMMON /callnwpic/xnwpic_called
NAMELIST /sclrplt_cntl1/ &
hplot, hinc, hminc, hmaxc, hovr, hcol1,hcol2, hprio, &
hhlf, hzro, hsty, &
msfplt,msfinc,msfminc, msfmaxc,msfovr,msfcol1,msfcol2,msfprio, &
msfhlf, msfzro, msfsty, &
thkplt,thkinc,thkminc, thkmaxc,thkovr,thkcol1,thkcol2,thkprio, &
thkhlf, thkzro, thksty, &
tplot, tinc, tminc, tmaxc, tovr, tcol1,tcol2, tprio, &
tunits, thlf, tzro, tsty, &
uplot, uinc, uminc, umaxc, uovr, ucol1,ucol2, uprio, &
uhlf, uzro, usty, &
vplot, vinc, vminc, vmaxc, vovr, vcol1,vcol2, vprio, &
vhlf, vzro, vsty, &
vhplot, vhinc, vhminc, vhmaxc, vhovr, vhcol1,vhcol2,vhprio, &
vhunits, vhhlf, vhzro, vhsty, &
vsplot, vsinc, vsminc, vsmaxc, vsovr, vscol1,vscol2,vsprio, &
vshlf, vszro, vssty, &
wplot, winc, wminc, wmaxc, wovr, wcol1,wcol2, wprio, &
whlf, wzro, wsty, &
ptplot, ptinc, ptminc, ptmaxc, ptovr, ptcol1,ptcol2,ptprio, &
pthlf, ptzro, ptsty, &
pplot , pinc, pminc, pmaxc, povr, pcol1,pcol2, pprio, &
phlf, pzro, psty, &
ipvplt,ipvinc,ipvminc, ipvmaxc,ipvovr,ipvcol1,ipvcol2,ipvprio, &
ipvhlf, ipvzro, ipvsty
NAMELIST /sclrplt_cntl2/ &
qvplot, qvinc, qvminc, qvmaxc, qvovr, qvcol1,qvcol2,qvprio, &
qvhlf, qvzro, qvsty, &
qcplot, qcinc, qcminc, qcmaxc, qcovr, qccol1,qccol2,qcprio, &
qchlf, qczro, qcsty, &
qrplot, qrinc, qrminc, qrmaxc, qrovr, qrcol1,qrcol2,qrprio, &
qrhlf, qrzro, qrsty, &
qiplot, qiinc, qiminc, qimaxc, qiovr, qicol1,qicol2,qiprio, &
qihlf, qizro, qisty, &
qsplot, qsinc, qsminc, qsmaxc, qsovr, qscol1,qscol2,qsprio, &
qshlf, qszro, qssty, &
qhplot, qhinc, qhminc, qhmaxc, qhovr, qhcol1,qhcol2,qhprio, &
qhhlf, qhzro, qhsty, &
qwplot, qwinc, qwminc, qwmaxc, qwovr, qwcol1,qwcol2,qwprio, &
qwhlf, qwzro, qwsty, &
qtplot, qtinc, qtminc, qtmaxc, qtovr, qtcol1,qtcol2,qtprio, &
qthlf, qtzro, qtsty
NAMELIST /sclrplt_cntl3/ &
kmhplt, kmhinc, kmhminc,kmhmaxc,kmhovr,kmhcol1,kmhcol2,kmhprio, &
kmhhlf, kmhzro, kmhsty, &
kmvplt, kmvinc, kmvminc,kmvmaxc,kmvovr,kmvcol1,kmvcol2,kmvprio, &
kmvhlf, kmvzro, kmvsty, &
tkeplt, tkeinc, tkeminc, tkemaxc,tkeovr,tkecol1,tkecol2,tkeprio, &
tkehlf, tkezro, tkesty, &
rhplot, rhinc, rhminc, rhmaxc, rhovr, rhcol1,rhcol2,rhprio, &
rhhlf, rhzro, rhsty, &
tdplot, tdinc, tdminc, tdmaxc, tdovr, tdcol1,tdcol2,tdprio, &
tdunits, tdhlf, tdzro, tdsty, &
rfopt, &
rfplot, rfinc, rfminc, rfmaxc, rfovr, rfcol1,rfcol2,rfprio, &
rfhlf, rfzro, rfsty, &
rfcplt, rfcinc, rfcminc, rfcmaxc,rfcovr,rfccol1,rfccol2,rfcprio, &
rfchlf, rfczro, rfcsty, &
pteplt, pteinc, pteminc, ptemaxc,pteovr,ptecol1,ptecol2,pteprio, &
ptehlf, ptezro, ptesty
NAMELIST /sclrplt_cntl_prt1/ &
upplot, upinc, upminc, upmaxc, upovr,upcol1,upcol2,upprio, &
uphlf, upzro, upsty, &
vpplot, vpinc, vpminc, vpmaxc, vpovr,vpcol1,vpcol2,vpprio, &
vphlf, vpzro, vpsty, &
wpplot, wpinc, wpminc, wpmaxc, wpovr,wpcol1,wpcol2,wpprio, &
wphlf, wpzro, wpsty, &
ptpplt, ptpinc, ptpminc,ptpmaxc,ptpovr,ptpcol1,ptpcol2,ptpprio, &
ptphlf, ptpzro, ptpsty, &
ppplot, ppinc, ppminc, ppmaxc, ppovr, ppcol1,ppcol2,ppprio, &
pphlf, ppzro, ppsty, &
qvpplt, qvpinc, qvpminc,qvpmaxc,qvpovr,qvpcol1,qvpcol2,qvpprio, &
qvphlf, qvpzro, qvpsty, &
vorpplt,vorpinc,vorpminc, vorpmaxc, vorpovr, vorpcol1,vorpcol2, &
vorphlf, vorpprio, vorpzro, vorpsty, &
divpplt,divpinc,divpminc, divpmaxc, divpovr, divpcol1,divpcol2, &
divphlf, divpprio, divpzro, divpsty, &
divqplt,divqinc,divqminc, divqmaxc, divqovr, divqcol1,divqcol2, &
divqhlf,divqprio, divqzro,divqsty
NAMELIST /sclrplt_cntl_prt2/ &
gricplt,gricinc,gricminc, gricmaxc, gricovr, griccol1,griccol2, &
grichlf,gricprio, griczro, gricsty, &
avorplt,avorinc,avorminc, avormaxc, avorovr, avorcol1,avorcol2, &
avorhlf,avorprio, avorzro, avorsty, &
rhiplot, rhiinc, rhiminc, rhimaxc, rhiovr, rhicol1,rhicol2, &
rhiprio, rhihlf, rhizro , rhisty
NAMELIST /vctrplt_cntl/istride,jstride,kstride, &
vtrplt, vtrunit,vtrovr,vtrcol1,vtrcol2,vtrprio,vtrunits,vtrtype, &
vtpplt, vtpunt, vtpovr,vtpcol1,vtpcol2,vtpprio,vtpunits,vtptype, &
xuvplt, xuvunit,xuvovr,xuvcol1,xuvcol2,xuvprio,xuvunits,xuvtype, &
strmplt,strmunit,strmovr,strmcol1,strmcol2,strmprio,strmunits, &
strmtype, &
vagplt, vagunit,vagovr,vagcol1,vagcol2,vagprio,vagunits,vagtype
NAMELIST /strmplt_cntl/ &
vtrstrm, vtrstmovr, vtrstmcol1, vtrstmcol2, vtrstmprio, &
vtpstrm, vtpstmovr, vtpstmcol1, vtpstmcol2, vtpstmprio
NAMELIST /plotting_setup/ iorig, xorig, yorig, &
xbgn, xend, ybgn, yend, zbgn, zend, zsoilbgn, zsoilend, &
yxstrch, zxstrch, zystrch, zhstrch, zsoilxstrch, zsoilystrch, &
winsiz, &
margnx, margny,pcolbar
NAMELIST /style_tuning/ lblmag,lnmag, fontopt, &
lbaxis,axlbfmt,axlbsiz, haxisu, vaxisu, &
tickopt,hmintick,vmajtick,vmintick,hmajtick, &
presaxis_no,pres_val, &
ctrlbopt,ctrstyle,ctrlbfrq,ctrlbsiz,lbmaskopt
NAMELIST /title_setup/ntitle,titcol,titsiz,title
NAMELIST /footer_setup/wpltime,footer_l,footer_c, footer_r
NAMELIST /domain_move/ imove, umove, vmove
NAMELIST /xy_slice_cntl/ nslice_xy, slice_xy
NAMELIST /xz_slice_cntl/ nslice_xz, slice_xz
NAMELIST /yz_slice_cntl/ nslice_yz, slice_yz
NAMELIST /h_slice_cntl/ nslice_h, slice_h
!
! 05/30/2002 Zuwen He
! Slice controls for soil model
!
NAMELIST /xy_soil_slice_cntl/ nslice_xy_soil, slice_xy_soil
NAMELIST /xz_soil_slice_cntl/ nslice_xz_soil, slice_xz_soil
NAMELIST /yz_soil_slice_cntl/ nslice_yz_soil, slice_yz_soil
!
! END Zuwen
!
NAMELIST /v_slice_cntl/ nslice_v, xpnt1,ypnt1,xpnt2,ypnt2
NAMELIST /p_slice_cntl/ nslice_p, slice_p
NAMELIST /pt_slice_cntl/ nslice_pt, slice_pt
NAMELIST /map_plot/ ovrmap,mapgrid,latgrid,longrid,mapgridcol, &
nmapfile,mapfile,mapcol,mapline_style
NAMELIST /multi_setup/ missfill_opt,missval_colind
NAMELIST /ovr_terrain/ ovrtrn
NAMELIST /sfc_plot1/ &
trnplt,trninc,trnminc, trnmaxc,trnovr,trncol1,trncol2,trnprio, &
trnhlf, trnzro, trnsty, &
wetcanplt,wcpinc,wcpminc,wcpmaxc,wcpovr,wcpcol1,wcpcol2,wcpprio, &
wcphlf, wcpzro, wcpsty, &
raincplt,raincinc,raincminc,raincmaxc,racovr,raccol1,raccol2, &
rachlf, racprio, raczro, racsty, racunit, &
raingplt,rainginc,raingminc,raingmaxc,ragovr,ragcol1,ragcol2, &
raghlf, ragprio, ragzro, ragsty, ragunit, &
raintplt,raintinc,raintminc,raintmaxc,ratovr,ratcol1,ratcol2, &
rathlf, ratprio, ratzro, ratsty, ratunit
!
! 05/30/2002 Zuwen He
!
NAMELIST /soil_plot/ &
tsoilplt,tsoilinc,tsoilminc,tsoilmaxc,tsoilovr, &
tsoilcol1,tsoilcol2,tsoilhlf,tsoilprio,tsoilzro, &
qsoilplt,qsoilinc,qsoilminc,qsoilmaxc,qsoilovr, &
qsoilcol1,qsoilcol2,qsoilhlf,qsoilprio,qsoilzro
!
! END Zuwen
!
NAMELIST /sfc_plot2/ &
pslplt ,pslinc, pslminc, pslmaxc,pslovr,pslcol1,pslcol2,pslprio, &
pslhlf, pslzro, pslsty, &
capeplt,capeinc,capeminc,capemaxc,capovr,capcol1,capcol2,capprio, &
caphlf, capzro, capsty, &
cinplt, cininc, cinminc, cinmaxc, cinovr,cincol1,cincol2,cinprio, &
cinhlf, cinzro, cinsty, &
thetplt,thetinc,thetminc,thetmaxc,theovr,thecol1,thecol2,theprio, &
thehlf, thezro, thesty, &
heliplt,heliinc,heliminc,helimaxc,helovr,helcol1,helcol2,helprio, &
helhlf, helzro, helsty, &
brnplt, brninc, brnminc, brnmaxc, brnovr,brncol1,brncol2,brnprio, &
brnhlf, brnzro, brnsty, &
brnuplt, brnuinc, bruminc, brumaxc, brnuovr, brnucol1,brnucol2, &
brnuhlf, brnuzro, brnusty, bruprio, &
srlfplt, srlfinc, srlminc, srlmaxc, srlfovr, srlfcol1,srlfcol2, &
srlfhlf, srlfzro, srlfsty, srlprio, &
srmfplt, srmfinc, srmminc, srmmaxc, srmfovr, srmfcol1,srmfcol2, &
srmfhlf, srmfzro, srmfsty, srmprio
NAMELIST /sfc_plot3/ &
liplt, liinc, liminc, limaxc, liovr, licol1,licol2,liprio, &
lihlf, lizro, listy, &
capsplt, capsinc, capsminc, capsmaxc, capsovr, capscol1,capscol2, &
capshlf, capszro, capssty, capsprio, &
blcoplt, blcoinc, blcominc, blcomaxc, blcoovr, blcocol1,blcocol2, &
blcohlf, blcozro, blcosty, blcoprio, &
viqcplt, viqcinc, viqcminc, viqcmaxc, viqcovr, viqccol1,viqccol2, &
viqchlf, viqczro, viqcsty, viqcprio, &
viqiplt, viqiinc, viqiminc, viqimaxc, viqiovr, viqicol1,viqicol2, &
viqihlf, viqizro, viqisty, viqiprio, &
viqrplt, viqrinc, viqrminc, viqrmaxc, viqrovr, viqrcol1,viqrcol2, &
viqrhlf, viqrzro, viqrsty, viqrprio, &
viqsplt, viqsinc, viqsminc, viqsmaxc, viqsovr, viqscol1,viqscol2, &
viqshlf, viqszro, viqssty,viqsprio, &
viqhplt, viqhinc, viqhminc, viqhmaxc, viqhovr, viqhcol1,viqhcol2, &
viqhhlf, viqhzro, viqhsty,viqhprio, &
vilplt, vilinc, vilminc, vilmaxc, vilovr, vilcol1,vilcol2, &
vilhlf, vilzro, vilsty, vilprio
NAMELIST /sfc_plot4/ &
viiplt, viiinc, viiminc, viimaxc, viiovr, viicol1,viicol2, &
viihlf, viizro, viisty, viiprio, &
vicplt, vicinc, vicminc, vicmaxc, vicovr, viccol1,viccol2, &
vichlf, viczro, vicsty, vicprio, &
ctcplt, ctcinc, ctcminc, ctcmaxc, ctcovr, ctccol1,ctccol2, &
ctchlf, ctczro, ctcsty, ctcprio, &
vitplt, vitinc, vitminc, vitmaxc, vitovr, vitcol1,vitcol2, &
vithlf, vitzro, vitsty, vitprio, &
pwplt, pwinc, pwminc, pwmaxc, pwovr, pwcol1,pwcol2, &
pwhlf, pwzro, pwsty, pwprio, &
tprplt, tprinc, tprminc, tprmaxc, tprovr, tprcol1,tprcol2, &
tprhlf, tprzro, tprsty, tprprio, tprunits, &
gprplt, gprinc, gprminc, gprmaxc, gprovr, gprcol1,gprcol2, &
gprhlf, gprzro, gprsty, gprprio, gprunits, &
cprplt, cprinc, cprminc, cprmaxc, cprovr, cprcol1,cprcol2, &
cprhlf, cprzro, cprsty, cprprio, cprunits
NAMELIST /sfc_cha_plot/ &
soiltpplt,soiltpinc,soiltpminc,soiltpmaxc,styovr,stycol1,stycol2, &
styhlf, styzro, stysty,styprio,soiltpn, &
vegtpplt,vegtpinc,vegtpminc,vegtpmaxc,vtyovr,vtycol1,vtycol2, &
vtyhlf, vtyzro, vtysty,vtyprio, &
laiplt,laiinc,laiminc,laimaxc,laiovr,laicol1,laicol2,laiprio, &
laihlf, laizro, laisty, &
rouplt,rouinc,rouminc,roumaxc,rouovr,roucol1,roucol2,rouprio, &
rouhlf, rouzro, rousty, &
vegplt,veginc,vegminc,vegmaxc,vegovr,vegcol1,vegcol2,vegprio, &
veghlf, vegzro, vegsty, &
snowdplt,snowdinc,snowdminc,snowdmaxc,snowdovr,snowdcol1, &
snowdcol2, snowdprio,snowdhlf, snowdzro, snowdsty
NAMELIST /wirfrm_plot/ w3dplt, wisosf, q3dplt, qisosf
REAL :: paprlnth
NAMELIST /page_setup/ layout, nxpic, nypic, inwfrm,paprlnth
NAMELIST /profile_cntl/ profopt, nprof, xprof, yprof, &
npicprof, uprof, uprmin, uprmax, vprof, vprmin, vprmax, &
wprof,wprmin,wprmax, ptprof,ptprmin,ptprmax, &
pprof,pprmin,pprmax, qvprof,qvprmin,qvprmax, &
qcprof,qcpmin,qcpmax, qrprof,qrpmin,qrpmax, &
qiprof,qipmin,qipmax, qsprof,qspmin,qspmax, &
qhprof,qhpmin,qhpmax, rhprof,rhpmin,rhpmax, &
kmhprof,kmhpmin,kmhpmax, kmvprof,kmvpmin,kmvpmax, &
tkeprof,tkepmin,tkepmax, &
rfprof,rfpmin,rfpmax, pteprf,ptepmin,ptepmax, &
upprof,uppmin,uppmax, vpprof,vppmin,vppmax, &
wpprof,wppmin,wppmax, ptpprf,ptppmin,ptppmax, &
ppprof,pppmin,pppmax, qvpprf,qvppmin,qvppmax, &
vorpprf, vorppmin, vorppmax, divpprf, divppmin, divppmax, &
zprofbgn,zprofend, &
tsoilprof,tsoilprofmin,tsoilprofmax, &
qsoilprof,qsoilprofmin,qsoilprofmax, &
zsoilprofbgn,zsoilprofend, &
nxprpic, nyprpic
NAMELIST /plot_boxes/ number_of_boxes,boxcol, &
bctrx,bctry,blengx,blengy
NAMELIST /plot_polylines/ number_of_polys,polycol,vertx,verty
NAMELIST /plot_obs/ ovrobs,sfcobfl,obscol,obs_marktyp,obs_marksz
NAMELIST /plot_sta/ ovrstaopt,ovrstam,ovrstan,ovrstav,wrtstax, &
wrtstad, stacol, markprio, nsta_typ, sta_typ, sta_marktyp, &
sta_markcol,sta_marksz,stalofl
NAMELIST /ovrlay_mul/ovrlaymulopt,ovrname,ovrmul_num,ovrmulname
NAMELIST /setcont_cntl/setcontopt,setcontnum,setcontvar,setconts
NAMELIST /arbvar_cntl/arbvaropt, &
var3dnum,dirname3d, &
var3d,var3dplot, var3dinc, var3dminc,var3dmaxc, &
var3dovr, var3dhlf, var3dzro,var3dsty,var3dcol1, var3dcol2, &
var3dprio, var2dnum,dirname2d, &
var2d,var2dplot, var2dinc, var2dminc,var2dmaxc, &
var2dovr, var2dhlf, var2dzro, var2dsty, var2dcol1, var2dcol2, &
var2dprio
INTEGER :: zxplot_called,zxout_unit
REAL :: xorig_0, yorig_0
REAL :: omega2, r2d
REAL :: relvort
REAL, ALLOCATABLE :: fcorio(:,:), mapfct(:,:)
REAL :: umove_readin, vmove_readin,uadd,vadd
!
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directives for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
REAL :: amin, amax, tem, pref
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
WRITE(6,'(/ 16(/5x,a)//)')
WRITE(6,'(/ 16(/5x,a)//)') &
'###############################################################', &
'###############################################################', &
'#### ####', &
'#### Welcome to ARPSPLT ####', &
'#### ####', &
'#### A graphic analysis program for model ARPS 5.0 ####', &
'#### ####', &
'#### The graphic plotting is based ####', &
'#### on graphic package ZXPLOT ####', &
'#### by Ming Xue CAPS/SOM/OU ####', &
'#### (http://www.caps.ou.edu/ZXPLOT) ####', &
'#### ####', &
'###############################################################', &
'###############################################################'
!
!-----------------------------------------------------------------------
!
! Get the names of the input data files.
!
!-----------------------------------------------------------------------
!
CALL get_input_file_names
(hinfmt,grdbasfn,hisfile,nhisfile)
lengbf = len_trim(grdbasfn)
CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf), &
nx,ny,nz,nzsoil,nstyps, ireturn)
nstyp = nstyps ! Copy to global variable
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,', nzsoil=',nzsoil
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,', nzsoil=',nzsoil
print*,'nstyps =', nstyps
cpu0 = f_cputime()
!
!-----------------------------------------------------------------------
!
! Allocate arrays.
!
!-----------------------------------------------------------------------
!
ALLOCATE(x (nx),STAT=istatus)
ALLOCATE(y (ny),STAT=istatus)
ALLOCATE(z (nz),STAT=istatus)
ALLOCATE(zp (nx,ny,nz),STAT=istatus)
ALLOCATE(zpsoil(nx,ny,nzsoil),STAT=istatus)
ALLOCATE(u (nx,ny,nz),STAT=istatus)
ALLOCATE(v (nx,ny,nz),STAT=istatus)
ALLOCATE(w (nx,ny,nz),STAT=istatus)
ALLOCATE(ptprt (nx,ny,nz),STAT=istatus)
ALLOCATE(pprt (nx,ny,nz),STAT=istatus)
ALLOCATE(qv (nx,ny,nz),STAT=istatus)
ALLOCATE(qc (nx,ny,nz),STAT=istatus)
ALLOCATE(qr (nx,ny,nz),STAT=istatus)
ALLOCATE(qi (nx,ny,nz),STAT=istatus)
ALLOCATE(qs (nx,ny,nz),STAT=istatus)
ALLOCATE(qh (nx,ny,nz),STAT=istatus)
ALLOCATE(tke (nx,ny,nz),STAT=istatus)
ALLOCATE(kmh (nx,ny,nz),STAT=istatus)
ALLOCATE(kmv (nx,ny,nz),STAT=istatus)
ALLOCATE(ubar (nx,ny,nz),STAT=istatus)
ALLOCATE(vbar (nx,ny,nz),STAT=istatus)
ALLOCATE(wbar (nx,ny,nz),STAT=istatus)
ALLOCATE(ptbar (nx,ny,nz),STAT=istatus)
ALLOCATE(pbar (nx,ny,nz),STAT=istatus)
ALLOCATE(rhobar(nx,ny,nz),STAT=istatus)
ALLOCATE(qvbar (nx,ny,nz),STAT=istatus)
ALLOCATE(soiltyp (nx,ny,nstyps),STAT=istatus)
ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
ALLOCATE(vegtyp (nx,ny),STAT=istatus)
ALLOCATE(lai (nx,ny),STAT=istatus)
ALLOCATE(roufns (nx,ny),STAT=istatus)
ALLOCATE(veg (nx,ny),STAT=istatus)
ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
ALLOCATE(wetcanp (nx,ny,0:nstyps),STAT=istatus)
ALLOCATE(snowdpth(nx,ny),STAT=istatus)
ALLOCATE(raing (nx,ny),STAT=istatus)
ALLOCATE(rainc (nx,ny),STAT=istatus)
ALLOCATE(raint (nx,ny),STAT=istatus)
ALLOCATE(prcrate (nx,ny,4),STAT=istatus)
ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
ALLOCATE(radsw (nx,ny),STAT=istatus)
ALLOCATE(rnflx (nx,ny),STAT=istatus)
ALLOCATE(radswnet (nx,ny),STAT=istatus)
ALLOCATE(radlwin (nx,ny),STAT=istatus)
ALLOCATE(usflx (nx,ny),STAT=istatus)
ALLOCATE(vsflx (nx,ny),STAT=istatus)
ALLOCATE(ptsflx(nx,ny),STAT=istatus)
ALLOCATE(qvsflx(nx,ny),STAT=istatus)
ALLOCATE(uprt (nx,ny,nz),STAT=istatus)
ALLOCATE(vprt (nx,ny,nz),STAT=istatus)
ALLOCATE(wprt (nx,ny,nz),STAT=istatus)
ALLOCATE(pt (nx,ny,nz),STAT=istatus)
ALLOCATE(qvprt (nx,ny,nz),STAT=istatus)
ALLOCATE(xc (nx,ny,max(nz,nzsoil)),STAT=istatus) ! shared with soil
ALLOCATE(yc (nx,ny,max(nz,nzsoil)),STAT=istatus) ! shared with soil
ALLOCATE(zc (nx,ny,nz),STAT=istatus)
ALLOCATE(zpc (nx,ny,nz),STAT=istatus)
ALLOCATE(zpsoilc(nx,ny,nzsoil),STAT=istatus) !Zuwen
ALLOCATE(hterain(nx,ny),STAT=istatus)
ALLOCATE(psl (nx,ny),STAT=istatus)
ALLOCATE(td (nx,ny,nz),STAT=istatus)
ALLOCATE(tz (nx,ny,nz),STAT=istatus)
ALLOCATE(t700 (nx,ny),STAT=istatus)
ALLOCATE(algpzc(nx,ny,nz),STAT=istatus)
ALLOCATE(zps3d (nx,ny,nz),STAT=istatus)
ALLOCATE(zpsoils3d (nx,ny,nzsoil),STAT=istatus)
ALLOCATE(wrk1 (nz),wrk2 (nz),wrk3 (nz),STAT=istatus)
ALLOCATE(wrk4 (nz),wrk5 (nz),wrk6 (nz),STAT=istatus)
ALLOCATE(wrk7 (nz),wrk8 (nz),wrk9 (nz),STAT=istatus)
ALLOCATE(wrk10(nz),wrk11(nz),wrk12(nz),STAT=istatus)
ALLOCATE(lcl (nx,ny),STAT=istatus)
ALLOCATE(lfc (nx,ny),STAT=istatus)
ALLOCATE(el (nx,ny),STAT=istatus)
ALLOCATE(twdf (nx,ny),STAT=istatus)
ALLOCATE(mbe (nx,ny),STAT=istatus)
ALLOCATE(li (nx,ny),STAT=istatus)
ALLOCATE(cape (nx,ny),STAT=istatus)
ALLOCATE(cin (nx,ny),STAT=istatus)
ALLOCATE(thet (nx,ny),STAT=istatus)
ALLOCATE(heli (nx,ny),STAT=istatus)
ALLOCATE(brn (nx,ny),STAT=istatus)
ALLOCATE(brnu (nx,ny),STAT=istatus)
ALLOCATE(srlfl(nx,ny),STAT=istatus)
ALLOCATE(srmfl(nx,ny),STAT=istatus)
ALLOCATE(shr37(nx,ny),STAT=istatus)
ALLOCATE(ustrm(nx,ny),STAT=istatus)
ALLOCATE(vstrm(nx,ny),STAT=istatus)
ALLOCATE(capst(nx,ny),STAT=istatus)
ALLOCATE(blcon(nx,ny),STAT=istatus)
ALLOCATE(ct (nx,ny),STAT=istatus)
ALLOCATE(sinlat(nx,ny),STAT=istatus)
ALLOCATE(xs(nx),STAT=istatus)
ALLOCATE(ys(ny),STAT=istatus)
ALLOCATE(tem1(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem2(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem3(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem4(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem5(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem6(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem7(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem8(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(tem9(nx,ny,max(nz,nzsoil)),STAT=istatus)
ALLOCATE(b1(nx,ny),STAT=istatus)
ALLOCATE(b2(nx,ny),STAT=istatus)
ALLOCATE(u1(nx,ny),STAT=istatus)
ALLOCATE(v1(nx,ny),STAT=istatus)
ALLOCATE(u2 (nx+ny,nz),STAT=istatus)
ALLOCATE(v2 (nx+ny,nz),STAT=istatus)
ALLOCATE(w2 (nx+ny,nz),STAT=istatus)
ALLOCATE(zs2(nx+ny,nz),STAT=istatus)
ALLOCATE(xp (nx+ny),STAT=istatus)
ALLOCATE(yp (nx+ny),STAT=istatus)
ALLOCATE(var3dv(nx,ny,nz),STAT=istatus)
ALLOCATE(var2dv(nx,ny),STAT=istatus)
ALLOCATE(fcorio(nx,ny),STAT=istatus)
ALLOCATE(mapfct(nx,ny),STAT=istatus)
x =0.0
y =0.0
z =0.0
zp =0.0
zpsoil=0.0
u =0.0
v =0.0
w =0.0
ptprt =0.0
pprt =0.0
qv =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
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
raing =0.0
rainc =0.0
raint =0.0
prcrate =0.0
radfrc=0.0
radsw =0.0
rnflx =0.0
usflx =0.0
vsflx =0.0
ptsflx=0.0
qvsflx=0.0
uprt =0.0
vprt =0.0
wprt =0.0
pt =0.0
qvprt =0.0
xc =0.0
yc =0.0
zc =0.0
zpc =0.0
zpsoilc =0.0
hterain=0.0
psl =0.0
td =0.0
tz =0.0
t700 =0.0
algpzc=0.0
zps3d =0.0
zpsoils3d =0.0
wrk1 =0.0
wrk4 =0.0
wrk7 =0.0
wrk10=0.0
lcl =0.0
lfc =0.0
el =0.0
twdf =0.0
mbe =0.0
li =0.0
cape =0.0
cin =0.0
thet =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
capst=0.0
blcon=0.0
ct =0.0
sinlat=0.0
xs=0.0
ys=0.0
tem1=0.0
tem2=0.0
tem3=0.0
tem4=0.0
tem5=0.0
tem6=0.0
tem7=0.0
tem8=0.0
tem9=0.0
xprof=0.0
yprof=0.0
b1=0.0
b2=0.0
u1=0.0
v1=0.0
u2 =0.0
v2 =0.0
w2 =0.0
zs2=0.0
xp =0.0
yp =0.0
b1=0.0
b2=0.0
u1=0.0
v1=0.0
u2 =0.0
v2 =0.0
w2 =0.0
zs2=0.0
xp =0.0
yp =0.0
var3dv=0.0
var2dv=0.0
fcorio=0.0
!-----------------------------------------------------------------------
! Set certain defaul options / values
!-----------------------------------------------------------------------
msfplt = 0
ipvplt = 0
vagplt = 0
thkplt = 0
lbcolor = 1 ! laber color 1-black
trcolor = 183 ! xz slice terrain color
imsg = 1 ! Activate missing value checking for streamline routine
uvmsg = -9999.0 ! Set the missing data value for streamline routine
inita = 8
initb = 8
layover = 0
ctinc = 0.0
vtunt = 0.0
aspratio = 1.0
xnwpic_called =0
!
!-----------------------------------------------------------------------
!
! Set the max. and min. values for producing HDF images.
!
!-----------------------------------------------------------------------
!
utmax = 40.0
utmin = -40.0
vtmax = 40.0
vtmin = -40.0
wtmax = 40.0
wtmin = -40.0
vsmax = 40.0
vsmin = -40.0
qvmax = 15.0
qvmin = 0.0
qcmax = 5.0
qcmin = 0.0
qrmax = 10.0
qrmin = 0.0
qimax = 10.0
qimin = 0.0
qsmax = 10.0
qsmin = 0.0
qhmax = 10.0
qhmin = 0.0
qwmax = 10.0
qwmin = 0.0
rhmax = 1.0
rhmin = 0.0
rhimax = 1.0
rhimin = 0.0
rfmax = 75.0
rfmin = 0.0
ptemax= 400.0
ptemin= 250.0
upmax = 40.0
upmin = -40.0
vpmax = 40.0
vpmin = -40.0
wpmax = 5.0
wpmin = -5.0
ptpmax= 0.0
ptpmin= -8.0
ppmax = 200.0
ppmin =-200.0
qvpmax= 10.0
qvpmin= 0.0
vormax= 0.02
vormin= -0.02
avormax= 0.02
avormin= -0.02
divmax= 0.02
divmin= -0.02
hmax = 2000.0
hmin = 0.0
tdmin = 200.
tdmax = 400.
xi1= 50.0
yi1= 50.0
xi2= 50.0
yi2= 50.0
paprlnth = 1.5 ! default value
lnmag = 1 ! default value
lblmag = 1.0
winsiz = 1.0
margnx = 0.1
margny = 0.1
pcolbar = 1
axlbfmt = -1
axlbsiz = 0.025
tickopt=0
ctrlbopt = 1
ctrstyle = 1
ctrlbfrq = 2
ctrlbsiz = 0.02
lbmaskopt = 0
istride = 0
jstride = 0
kstride = 0
!
!-----------------------------------------------------------------------
!
! Read in plotting control parameters
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Page set-up parameters
!
!-----------------------------------------------------------------------
!
READ(5,page_setup,ERR=100)
WRITE(6,'(a)')'Namelist page_setup was successfully read.'
READ(5,plotting_setup,ERR=100)
WRITE(6,'(a)')'Namelist plotting_setup was successfully read.'
READ(5,col_table_cntl, ERR=100)
WRITE(6,'(a)')'namelist col_table_cntl was successfully read.'
WRITE(6,'(a,i3)') 'Color table is : ',col_table
READ(5,style_tuning,ERR=100)
WRITE(6,'(a)')'Namelist style_tuning was successfully read.'
READ(5,smooth_cntl, ERR=100)
WRITE(6,'(a)')'Namelist ,smooth_cntl was successfully read.'
WRITE(6,'(a,i3)') 'Smoothing option is : ',smooth
READ(5,title_setup, ERR=100)
WRITE(6,'(a)')'Namelist title_setup was successfully read.'
READ(5,footer_setup, ERR=100)
WRITE(6,'(a)')'Namelist footer_setup was successfully read.'
!
!-----------------------------------------------------------------------
!
! Input control parameters map plotting
!
!-----------------------------------------------------------------------
!
mapgrid = 0 ! no longer used.
READ(5,map_plot,ERR=100)
WRITE(6,'(a)')'Namelist map_plot was successfully read.'
READ(5,multi_setup,ERR=100)
WRITE(6,'(a)')'Namelist multi_setup was successfully read.'
IF(nmapfile > maxmap) &
WRITE(6,'(a)')'Warning: the maximum map files should be ',maxmap
DO i=1,nmapfile
lmapfile=132
CALL strlnth(mapfile(i), lmapfile)
WRITE(6,'(1x,a,a)') 'Input was ',mapfile(i)(1:lmapfile)
IF(ovrmap == 1) THEN
INQUIRE(FILE=mapfile(i)(1:lmapfile), EXIST = fexist )
IF( .NOT.fexist) THEN
WRITE(6,'(a)') 'Warning: Map file '//mapfile(i)(1:lmapfile) &
//' not found. Program will be continue'
END IF
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Input control parameters plotting type
!
!-----------------------------------------------------------------------
!
READ(5,xy_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist xy_slice_cntl was successfully read.'
nslice(1)=nslice_xy
iplot(1)=1
IF(nslice(1) == 0) iplot(1)=0
READ(5,xz_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist xz_slice_cntl was successfully read.'
nslice(2)=nslice_xz
iplot(2)=1
IF(nslice(2) == 0) iplot(2)=0
READ(5,yz_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist yz_slice_cntl was successfully read.'
nslice(3)=nslice_yz
iplot(3)=1
IF(nslice(3) == 0) iplot(3)=0
READ(5,h_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist h_slice_cntl was successfully read.'
nslice(4)=nslice_h
iplot(4)=1
IF(nslice(4) == 0) iplot(4)=0
IF(nslice(4) > max_dim ) THEN
WRITE(6,'(a,i3,a)') 'Please give value smaller than ', &
nz,'. Program stopped !'
STOP
END IF
!
! 05/30/2002 Zuwen He
!
! xy_soil_slice_cntl is different from xy_slice_cntl.
! nslice(8) has been used for other purpose.
!
READ(5,xy_soil_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist xy_soil_slice_cntl was successfully read.'
nslice(9)=nslice_xy_soil
iplot(9)=1
IF(nslice(9) == 0) iplot(9)=0
READ(5,xz_soil_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist xz_soil_slice_cntl was successfully read.'
nslice(10)=nslice_xz_soil
iplot(10)=1
IF(nslice(10) == 0) iplot(10)=0
READ(5,yz_soil_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist yz_soil_slice_cntl was successfully read.'
nslice(11)=nslice_yz_soil
iplot(11)=1
IF(nslice(11) == 0) iplot(11)=0
!
! END Zuwen
!
READ(5,v_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist v_slice_cntl was successfully read.'
nslice(5)=nslice_v
iplot(5)=1
IF(nslice(5) == 0) iplot(5)=0
IF(nslice(5) > nx+ny) THEN
WRITE(6,'(1x,a,i3,a)') ' Please give a value smaller than ', &
nx+ny,'. Program stopped !'
STOP
END IF
DO indxslic = 1,nslice(5)
xi1(indxslic)=xpnt1(indxslic)
yi1(indxslic)=ypnt1(indxslic)
xi2(indxslic)=xpnt2(indxslic)
yi2(indxslic)=ypnt2(indxslic)
END DO
READ(5,p_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist p_slice_cntl was successfully read.'
WRITE(6,'(a)') 'Pressure slices to be plotted are at p=: '
nslice(6)=nslice_p
iplot(6)=1
IF(nslice(6) == 0) iplot(6)=0
IF(nslice(6) > nz) THEN
WRITE(6,'(1x,a,i3,a)') 'Please give value smaller than ', &
nz,'. Program stopped !'
STOP
END IF
DO indxslic = 1,nslice(6)
WRITE(6,'(1x,f10.3)') slice_p(indxslic)
END DO
READ(5,pt_slice_cntl,ERR=100)
WRITE(6,'(a)')'Namelist pt_slice_cntl was successfully read.'
IF(nslice_pt > max_dim) THEN
WRITE(6,'(1x,2(a,i3))') &
'Warning: Maximum number of PT slices allowed is max_dim= ', &
max_dim,'nslice_pt reset t= ',max_dim
nslice_pt = max_dim
END IF
nslice(7)=nslice_pt
iplot(7)=1
IF(nslice(7) == 0) iplot(7)=0
WRITE(6,'(a)')'Isentropic slices to be plotted are at theta=:'
DO indxslic = 1,nslice(7)
WRITE(6,'(1x,f10.3)') slice_pt(indxslic)
END DO
READ(5,domain_move,ERR=100)
WRITE(6,'(a)')'Namelist domain_move was successfully read.'
umove_readin = umove
vmove_readin = vmove
!
READ(5,sclrplt_cntl1, ERR=100)
WRITE(6,'(a)')'Namelist sclrplt_cntl1 was successfully read.'
READ(5,sclrplt_cntl2, ERR=100)
WRITE(6,'(a)')'Namelist sclrplt_cntl2 was successfully read.'
READ(5,sclrplt_cntl3, ERR=100)
WRITE(6,'(a)')'Namelist sclrplt_cntl3 was successfully read.'
READ(5,sclrplt_cntl_prt1, ERR=100)
WRITE(6,'(a)')'Namelist sclrplt_cntl_prt1 was successfully read.'
READ(5,sclrplt_cntl_prt2, ERR=100)
WRITE(6,'(a)')'Namelist sclrplt_cntl_prt2 was successfully read.'
READ(5,vctrplt_cntl, ERR=100)
WRITE(6,'(a)')'Namelist vctrplt_cntl was successfully read.'
READ(5,strmplt_cntl,ERR=100)
WRITE(6,'(a)')'Namelist strmplt_cntl was successfully read.'
ipriority(1)=hprio
ipriority(2)=tprio
ipriority(3)=uprio
ipriority(4)=vprio
ipriority(5)=vhprio
ipriority(6)=wprio
ipriority(7)=ptprio
ipriority(8)=pprio
ipriority(9)=qvprio
ipriority(10)=qcprio
ipriority(11)=qrprio
ipriority(12)=qiprio
ipriority(13)=qsprio
ipriority(14)=qhprio
ipriority(15)=qwprio
ipriority(16)=kmhprio
ipriority(17)=kmvprio
ipriority(18)=tkeprio
ipriority(19)=rhprio
ipriority(20)=tdprio
ipriority(21)=rfprio
ipriority(22)=pteprio
ipriority(23)=upprio
ipriority(24)=vpprio
ipriority(25)=wpprio
ipriority(26)=ptpprio
ipriority(27)=ppprio
ipriority(28)=qvpprio
ipriority(29)=vorpprio
ipriority(30)=divpprio
ipriority(31)=divqprio
ipriority(32)=vtrprio
ipriority(33)=vtpprio
ipriority(34)=vtrstmprio
ipriority(35)=vtpstmprio
ipriority(36)=rfcprio
ipriority(37)=vsprio
ipriority(38)=gricprio
ipriority(39)=avorprio
ipriority(40)=xuvprio
ipriority(41)=qtprio
ipriority(42)=rhiprio
ipriority(43)=tsoilprio ! Zuwen He 05/31/2002
ipriority(44)=qsoilprio ! Zuwen He 05/31/2002
!
!-----------------------------------------------------------------------
!
! Input control parameters for 2-d surface feild plotting
!
!-----------------------------------------------------------------------
!
READ(5,sfc_plot1,ERR=100)
WRITE(6,'(a)')'Namelist sfc_plot1 was successfully read.'
READ(5,soil_plot,ERR=100)
WRITE(6,'(a)')'Namelist sfc_plot1 was successfully read.'
READ(5,sfc_plot2,ERR=100)
WRITE(6,'(a)')'Namelist sfc_plot2 was successfully read.'
READ(5,sfc_plot3,ERR=100)
WRITE(6,'(a)')'Namelist sfc_plot3 was successfully read.'
READ(5,sfc_plot4,ERR=100)
WRITE(6,'(a)')'Namelist sfc_plot4 was successfully read.'
ipriority(51)=trnprio
ipriority(56)=wcpprio
ipriority(57)=racprio
ipriority(58)=ragprio
ipriority(59)=ratprio
ipriority(60)=pslprio
ipriority(61)=liprio
ipriority(62)=capprio
ipriority(63)=cinprio
ipriority(64)=theprio
ipriority(65)=helprio
ipriority(66)=ctcprio
ipriority(67)=brnprio
ipriority(68)=bruprio
ipriority(69)=srlprio
ipriority(70)=srmprio
ipriority(71)=capsprio
ipriority(72)=blcoprio
ipriority(73)=viqcprio
ipriority(74)=viqrprio
ipriority(75)=viqiprio
ipriority(76)=viqsprio
ipriority(77)=viqhprio
ipriority(78)=vilprio
ipriority(79)=viiprio
ipriority(80)=vicprio
ipriority(81)=strmprio
ipriority(82)=vitprio
ipriority(83)=pwprio
ipriority(84)=tprprio
ipriority(85)=gprprio
ipriority(86)=cprprio
!-----------------------------------------------------------------------
!
! Input control parameters for 2-d surface characteristics plotting
!
!-----------------------------------------------------------------------
!
READ(5,sfc_cha_plot,ERR=100)
WRITE(6,'(a)') &
'Namelist sfc_cha_plot was successfully read.'
ipriority(101)=styprio
ipriority(102)=vtyprio
ipriority(103)=laiprio
ipriority(104)=rouprio
ipriority(105)=vegprio
ipriority(106)=snowdprio
ipriority(107)=msfprio
ipriority(108)=ipvprio
ipriority(109)=vagprio
ipriority(110)=thkprio
!-----------------------------------------------------------------------
!
! Set some parameters
!
!-----------------------------------------------------------------------
trcolor = trncol1
nslice(8)=1
!-----------------------------------------------------------------------
!
! Input control parameter for uneven contour interval
!
!-----------------------------------------------------------------------
DO i=1,maxuneva
setcontvar(i)(1:12) = ' '
DO k=1,maxunevm
setconts(k,i) = -9999.
END DO
END DO
READ(5,setcont_cntl,ERR=100)
WRITE(6,'(a)')'Namelist setcont_cntl was successfully read.'
READ(5,arbvar_cntl,ERR=100)
WRITE(6,'(a)')'Namelist arbvar_cntl was successfully read.'
DO i=1, var3dnum
ipriority(150+i) = var3dprio(i)
END DO
DO i=1, var2dnum
ipriority(170+i) = var2dprio(i)
END DO
!
!-----------------------------------------------------------------------
!
! Input control parameters plotting boxes
!
!-----------------------------------------------------------------------
!
READ(5,plot_boxes,ERR=100)
WRITE(6,'(a)')'Namelist plot_box was successfully read.'
IF(number_of_boxes /= 0) THEN
DO k=1,number_of_boxes
WRITE(6,'(1x,a,i3,a,2f10.5)') &
'Center of box No.',k,' is at ',bctrx(k),bctry(k)
WRITE(6,'(1x,a,i3,a,2f10.5)') &
'The size of box No.',k,' is ',blengx(k),blengy(k)
END DO
DO k=1,number_of_boxes
bx1(k)=bctrx(k) - blengx(k)*0.5
bx2(k)=bctrx(k) + blengx(k)*0.5
by1(k)=bctry(k) - blengy(k)*0.5
by2(k)=bctry(k) + blengy(k)*0.5
END DO
END IF
!
!
!-----------------------------------------------------------------------
!
! Input control parameters plotting polylines
!
!-----------------------------------------------------------------------
!
DO j=1,max_polys
DO i=1,max_verts
vertx(i,j) = -9999.
verty(i,j) = -9999.
END DO
END DO
READ(5,plot_polylines,ERR=100)
WRITE(6,'(a)')'Namelist plot_polylines was successfully read.'
IF(number_of_polys /= 0) THEN
DO k=1,number_of_polys
WRITE(6,'(1x,a,i2)')'The number of polyline is : ',k
DO j = 1, max_verts
IF(vertx(j,k) /= -9999. .AND. verty(j,k) /= -9999.) &
WRITE(6,'(1x,a,2f10.5)') &
'The position of vertices are: ',vertx(j,k),verty(j,k)
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Input control parameters for overlay one filed to many fields
!
!-----------------------------------------------------------------------
!
READ(5,ovrlay_mul,ERR=100)
WRITE(6,'(a)')'Namelist ovrlay_mul was successfully read.'
plotovr = .false.
sovrlay = 0
IF(ovrlaymulopt == 0) THEN
ovrname(1:12)=' '
ovrmul_num = 0
DO i = 1,50
ovrmulname(i)(1:12) =' '
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Input control parameters for terrain overlay
!
!-----------------------------------------------------------------------
!
READ(5,ovr_terrain,ERR=100)
WRITE(6,'(a)')'Namelist ovr_terrain was successfully read.'
!
!-----------------------------------------------------------------------
!
! Input control parameters for 3-D wireframe plotting
!
!-----------------------------------------------------------------------
!
READ(5,wirfrm_plot,ERR=100)
WRITE(6,'(a)')'Namelist wirfrm_plot was successfully read.'
!
!-----------------------------------------------------------------------
!
! Parameters for overlaying observations
!
!-----------------------------------------------------------------------
!
ovrobs=0
READ(5,plot_obs,ERR=71)
WRITE(6,'(a)')'Namelist plot_obs was successfully read.'
lsfcobfl=LEN(sfcobfl)
CALL strlnth(sfcobfl,lsfcobfl)
WRITE(6,'(1x,a,a)') 'The surface observation file name was ', &
sfcobfl(1:lsfcobfl)
71 CONTINUE
nobs=0
IF(ovrobs == 1) THEN
INQUIRE(FILE=sfcobfl(1:lsfcobfl), EXIST = fexist )
IF( .NOT.fexist) THEN
WRITE(6,'(a)') 'Surface obs file '//sfcobfl(1:lsfcobfl)// &
' not found. Program will continue.'
END IF
CALL getunit(obunit)
CALL read_surface_obs(sfcobfl,mxsfcob,atime,n_meso_g, &
n_meso_pos,n_sao_g,n_sao_pos_g,n_sao_b,n_sao_pos_b,n_obs_g, &
n_obs_pos_g,n_obs_b,n_obs_pos_b, &
stn,obstype,latob,lonob,elevob,wx, &
tob,tdob,dd,ff,ddg,ffg,pstn,pmsl,alt, &
kloud,ceil,lowcld,cover,rad, &
idp3,store_emv,store_amt,store_hgt,vis,obstime,istatus)
nobs=n_obs_b
END IF
!-----------------------------------------------------------------------
!
! Parameters for overlaying airport location
!
!-----------------------------------------------------------------------
!
ovrstam=0
ovrstan=0
ovrstav=0
wrtstax=0
READ(5,plot_sta, ERR=72)
WRITE(6,'(a)') &
'Namelist plot_sta was successfully read.'
lstalofl=LEN(stalofl)
CALL strlnth(stalofl,lstalofl)
WRITE(6,'(1x,a,a)') 'Station file name was ',stalofl(1:lstalofl)
72 CONTINUE
nsta=0
IF(ovrstaopt /= 0) THEN
INQUIRE(FILE=stalofl(1:lstalofl), EXIST = fexist )
IF( .NOT.fexist) THEN
WRITE(6,'(a)') 'WARNING: Surface obs file ' &
//stalofl(1:lstalofl)// &
' not found. Program continue in ARPSPLT.'
ELSE
CALL getunit(stunit)
CALL read_station(stalofl,mxstalo,latsta,lonsta,nstatyp, &
nstapro,nsta,s_name,s_state,s_site,s_elev)
IF(nsta > 0) staset=1
END IF
END IF
!
!-----------------------------------------------------------------------
!
! When priority equal to zero , the order of the plotting is default.
!
!-----------------------------------------------------------------------
!
DO i=1,nprio
iptemp(i)=i
sigplt(i)=0
END DO
iptemp(nprio)=0
!-----------------------------------------------------------------------
!
! Input control parameter for profile
!
!-----------------------------------------------------------------------
READ(5,profile_cntl,ERR=100)
WRITE(6,'(a)')'Namelist profile_cntl was successfully read.'
IF (nprof > max_dim) THEN
WRITE (6,'(1x,a,i4)') 'Too many profiles. Limited to ',max_dim
nprof = max_dim
END IF
zxplot_called=0
lengbf = len_trim(grdbasfn)
DO nf=1, nhisfile
DO i=1,nprio
sigplt(i)=0
END DO
filename=hisfile(nf)
lenfil = len_trim(filename)
WRITE(6,'(a,a/)') &
' The file name of data set is ', filename(1:lenfil)
15 CONTINUE ! also continue to read another time recode
! from GrADS file
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
second1= f_walltime()
cpu1 = f_cputime()
CALL dtaread(nx,ny,nz,nzsoil, nstyps, &
hinfmt, nchin,grdbasfn(1:lengbf),lengbf, &
filename(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)
cpu2 = f_cputime()
second2 = f_walltime()
! print*,'!!!! total clock time for read data:',second2-second1
! print*,'!!!! total cpu time for read data :', cpu2-cpu1
!
!-----------------------------------------------------------------------
!
! 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 700 ! unsuccessful read
!
!-----------------------------------------------------------------------
!
! Set ice to be icein and then ice is used to control if snow and
! graupel/hail affect reflectivity calculation in REFLEC.
!
!-----------------------------------------------------------------------
!
ice = icein
CALL gtlfnkey( runname, lfnkey )
mgrid = 1 ! Grid number 1
nestgrd = 0 ! Not nested grid data
IF( zxplot_called == 0) THEN
!
!-----------------------------------------------------------------------
!
! Initialize ZXPLOT plotting package
!
!-----------------------------------------------------------------------
!
! The following two ZXPLOT routines, if called, should be called before xdevic
!
zxout_unit=2
CALL xpsfn(runname(1:lfnkey)//'.ps', zxout_unit)
CALL xpaprlnth( paprlnth )
CALL xdevic
CALL xstctfn(color_map)
CALL xsetclrs(col_table)
! CALL xafstyl(1)
CALL xcolor(lbcolor)
CALL xartyp(2)
CALL xlnmag(lnmag)
CALL xcfont(fontopt)
CALL xlabmask(lbmaskopt)
IF(ctrlbopt == 0) CALL xcltyp(0)
IF(ctrlbopt == 1) CALL xclfmt('(f10.1)')
IF(ctrlbopt == 2) CALL xclfmt('(I10)')
CALL xcmixl
IF( ctrstyle == 2) THEN
CALL xcfull
ELSE IF( ctrstyle == 3) THEN
CALL xcdash
END IF
CALL xclfrq( ctrlbfrq )
CALL xctrbadv( 1 ) ! Turn on missing value checking for contouring
CALL xvtrbadv( 1 ) ! Turn on missing value checking for vector plotting
!
! Default value of -9999.0 for the flag is used.
!
CALL xbadval ( -9999.0 ) ! Set the missing value flag to -9999.0
!
! above three routines were not available with older version of ZXPLOT.
!
CALL xdspac(0.9*winsiz)
IF( nxpic == 1 .AND. nypic == 1) CALL xdspac( 0.85*winsiz)
IF( layout == 1 ) THEN
angl = 00.0
ELSE
angl = 90.0
END IF
CALL xspace(nxpic, nypic, angl , xlimit,ylimit)
IF(axlbfmt == (-1)) THEN
CALL xaxfmt( '*' )
ELSE IF(axlbfmt == 0) THEN
CALL xaxfmt( '(i5)' )
ELSE
WRITE(stem1,'(i1)') axlbfmt
WRITE(stem2,43)stem1
CALL xaxfmt( stem2 )
END IF
43 FORMAT('(f8.',a1,')')
layover = 0
CALL xpmagn(margnx*xlimit, margny*ylimit)
zxplot_called=1
END IF ! Initalization of plotting package and parameters
!
!-----------------------------------------------------------------------
!
! Add the domain translation/grid movement speed back to the wind fields,
! so that the wind is the Earth relative.
!
!-----------------------------------------------------------------------
!
IF( imove == 1 ) THEN
uadd = umove_readin
vadd = vmove_readin
ELSEIF( imove == 2 ) THEN
uadd = umove
vadd = vmove
ELSE
uadd = 0.0
vadd = 0.0
END IF
DO k=1,nz
DO j=1,ny
DO i=1,nx
u(i,j,k)=uprt(i,j,k)+ubar(i,j,k)+uadd
v(i,j,k)=vprt(i,j,k)+vbar(i,j,k)+vadd
w(i,j,k)=wprt(i,j,k)+wbar(i,j,k)
qv(i,j,k)=MAX(0.0,qvprt(i,j,k)+qvbar(i,j,k))
pt(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
tem1(i,j,k)=0.0
tem2(i,j,k)=0.0
tem3(i,j,k)=0.0
tem4(i,j,k)=0.0
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Set coordinates at the grid volume center
!
!-----------------------------------------------------------------------
!
DO k=1,max(nz,nzsoil)
DO j=1,ny
DO i=1,nx-1
xc(i,j,k) = (x(i)+x(i+1))*0.5 *0.001
END DO
END DO
END DO
DO k=1,max(nz,nzsoil)
DO j=1,ny-1
DO i=1,nx
yc(i,j,k) = (y(j)+y(j+1))*0.5 *0.001
END DO
END DO
END DO
DO k=1,nz-1
DO j=1,ny
DO i=1,nx
zc(i,j,k) = (z(k)+z(k+1))*0.5 *0.001
zpc(i,j,k)= (zp(i,j,k)+zp(i,j,k+1))*0.5 *0.001 ! km
END DO
END DO
END DO
!
! 05/30/2002 Zuwen He
!
! zpsoilc is the vertical distance of the soil model
! from the terrain surface in meters.
!
DO k=1,nzsoil
DO j=1,ny
DO i=1,nx
zpsoilc(i,j,k)=(zpsoil(i,j,k)-zpsoil(i,j,1))*100.
! cm for soil model, Zuwen
END DO
END DO
END DO
!
! END Zuwen
!
!-----------------------------------------------------------------------
!
! Set negative logrithm pressure (Pascal) coordinates
! for scalar and w grid points
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
algpzc(i,j,k) = -ALOG(pbar(i,j,k)+pprt(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate temperature (K) at ARPS grid points and 700mb
! pressure level
!
!-----------------------------------------------------------------------
!
CALL temper (nx,ny,nz,pt, pprt ,pbar,tz)
zlevel=-ALOG(70000.0)
CALL hintrp(nx,ny,nz,tz,algpzc,zlevel,t700)
!
!----------------------------------------------------------------------
!
! Calculate sea level pressure (mb)
! Reduction method: Benjamin and Miller: 1990, MWR, vol.118, No.10,
! Page: 2100-2101
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
DO j=1,ny-1
p00 = 0.01*(pbar(i,j,2)+pprt(i,j,2))
IF(p00 <= 700.0) THEN
t00=tz(i,j,2)
ELSE
t00 = t700(i,j)*(p00/700.0)**ex1
END IF
psl(i,j) = p00*((t00+gamma*zpc(i,j,2))/t00)**ex2
IF(mod(i,5) == 0 .AND. mod(j,5) == 0 .AND. zpc(i,j,2) > 0.5 ) &
print *, &
' zpc,p0,t0,t,psl:',zpc(i,j,2),p00,t00,tz(i,j,2),psl(i,j)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate dew-point temperature td using enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = pbar(i,j,k )+pprt (i,j,k)
tem3(i,j,k) = ptbar(i,j,k)+ptprt(i,j,k)
END DO
END DO
END DO
CALL temper(nx,ny,nz,tem3, pprt ,pbar,tem2)
CALL getdew(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, tem1, tem2, qv, td)
!
!-----------------------------------------------------------------------
!
! Calculate CAPE and CIN
! change qv to mixing ratio
!
!-----------------------------------------------------------------------
!
CALL edgfill(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem6(i,j,k)=0.5*(u(i,j,k)+u(i+1,j,k))
tem7(i,j,k)=0.5*(v(i,j,k)+v(i,j+1,k))
tem8(i,j,k)=qv(i,j,k)/(1.-qv(i,j,k))
tem9(i,j,k)=pprt(i,j,k)+pbar(i,j,k)
END DO
END DO
END DO
CALL edgfill(tem6,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(tem7,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(tem8,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(tem9,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(tz,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(td,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
IF (capeplt /= 0 .OR. cinplt /= 0 .OR. capsplt /= 0 .OR. &
liplt /= 0 .OR. brnplt /= 0) THEN
second1= f_walltime()
cpu1 = f_cputime()
CALL arps_be(nx,ny,nz, &
tem9,zpc,tz,tem8, &
lcl,lfc,el,twdf,li,cape,mbe,cin,capst, &
wrk1,wrk2,wrk3,wrk4,wrk5,wrk6, &
wrk7,wrk8,wrk9,wrk10,wrk11,wrk12,tem5)
cpu2 = f_cputime()
second2 = f_walltime()
! print*,'!!!! total clock time for arps_be:',second2-second1
! print*,'!!!! total cpu time for arps_be :', cpu2-cpu1
! print *, ' back from arps_be'
END IF
!-----------------------------------------------------------------------
!
! Calculate the convective temperature (celsius)
!
!-----------------------------------------------------------------------
IF(ctcplt /= 0) THEN
CALL arps_ct(ct,nx,ny,nz,tem9,tz,td,tem1,wrk1,wrk2,wrk3 )
END IF
!
!-----------------------------------------------------------------------
!
! Store k=2 theta-e and dewpt
! level 1 and 2 respectively.
!
!-----------------------------------------------------------------------
!
CALL pt2pte(nx,ny,1, 1,nx,1,ny,1,1, &
tem9(1,1,2),tem3(1,1,2),qv(1,1,2),thet)
!
!-----------------------------------------------------------------------
!
! Calculate helicity and other shear-related paramaters
!
!-----------------------------------------------------------------------
!
IF (heliplt /= 0 .OR. brnplt /= 0 .OR. brnuplt /= 0 .OR. &
srlfplt /= 0 .OR. srmfplt /= 0 .OR. blcoplt /= 0 ) THEN
CALL xytomf(nx,ny,x,y,mapfct)
CALL calcshr(nx,ny,nz,x,y,zp,mapfct, &
tem9,tz,tem6,tem7,cape, &
shr37,ustrm,vstrm,srlfl,srmfl,heli,brn,brnu,blcon, &
tem1,u1,v1,tem2,tem5)
END IF
!
!-----------------------------------------------------------------------
!
! Modify the coordinate arrays to make certain that the origin is
! the same as specified above.
!
!-----------------------------------------------------------------------
!
dx = x(2)-x(1) ! in meters
dy = y(2)-y(1) ! in meters
dz = z(2)-z(1) ! in meters
dxinv = 1.0/(x(2)-x(1))
dyinv = 1.0/(y(2)-y(1))
dzinv = 1.0/(z(2)-z(1))
IF( iorig == 1) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
xc(i,j,k) = xc(i,j,k) + xgrdorg * 0.001
yc(i,j,k) = yc(i,j,k) + ygrdorg * 0.001
END DO
END DO
END DO
ELSE IF( iorig == 2) THEN
xorold = (xc(1,2,2)+xc(2,2,2))/2 + xgrdorg * 0.001
yorold = (yc(2,1,2)+yc(2,2,2))/2 + ygrdorg * 0.001
dxkm = xc(3,2,2)-xc(2,2,2) ! dx in km
dykm = yc(2,3,2)-yc(2,2,2) ! dy in km
WRITE(6,'(/1x,4(a,f8.3),a/)') &
'Model grid origin (',xorold,',',yorold,') was reset to (' &
,xorig,',',yorig,').'
IF( ABS(xorig-xorold) > 0.01*dxkm .OR. &
ABS(yorig-yorold) > 0.01*dykm ) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
xc(i,j,k) = xc(i,j,k) - xorold + xorig
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny
DO i=1,nx
yc(i,j,k) = yc(i,j,k) - yorold + yorig
END DO
END DO
END DO
END IF
END IF
curtim = time
!
!-----------------------------------------------------------------------
!
! If islice=-2, plot the y-z cross-sections through w maximum.
! If jslice=-2, plot the x-z cross-sections through w maximum.
! If kslice=-2, plot the x-y cross-sections through w maximum.
!
!-----------------------------------------------------------------------
!
CALL a3dmax(wprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,2,nz-1, &
wmax,wmin, imax,jmax,kmax, imin,jmin,kmin)
! write(6,'(/2(1x,a,f9.3,3(a,i3))/)')
! : 'wmin =',wmin,' at i=',imin,', j=',jmin,', k=',kmin,
! : 'wmax =',wmax,' at i=',imax,', j=',jmax,', k=',kmax
!-----------------------------------------------------------------------
!
! Get index bounds of the domain to be plotted
!
!-----------------------------------------------------------------------
!
CALL indxbnds(xc,yc,zpc,zpsoilc,nx,ny,nz,nzsoil, &
xbgn,xend,ybgn,yend,zbgn,zend,zsoilbgn,zsoilend, &
ibgn,iend,jbgn,jend,kbgn,kend,ksoilbgn,ksoilend)
IF( istride /= 0 ) THEN
ist = istride
ELSE IF ( (iend-ibgn+1) > 110) THEN
ist = 4
ELSE IF( (iend-ibgn+1) > 30) THEN
ist = 2
ELSE
ist = 1
END IF
IF( jstride /= 0 ) THEN
jst = jstride
ELSE IF ( (jend-jbgn+1) > 110) THEN
jst = 4
ELSE IF( (jend-jbgn+1) > 30) THEN
jst = 2
ELSE
jst = 1
END IF
IF( kstride /= 0 ) THEN
kst = kstride
ELSE IF ( (kend-kbgn+1) > 90 ) THEN
kst = 3
ELSE IF( (kend-kbgn+1) > 45) THEN
kst = 2
ELSE
kst = 1
END IF
dxkm = xc(3,2,2)-xc(2,2,2)
dykm = yc(2,3,2)-yc(2,2,2)
dzkm = zc(2,2,3)-zc(2,2,2)
dzsoilcm = abs(zpsoilc(1,1,2)-zpsoilc(1,1,1))/2.
xr = xc(iend,2,2)-xc(ibgn,2,2)+dxkm
yr = yc(2,jend,2)-yc(2,jbgn,2)+dykm
zmax=zp(ibgn,jbgn,kend+1)
zmin=zp(ibgn,jbgn,kbgn)
DO j=jbgn,jend
DO i=ibgn,iend
zmax=AMAX1(zmax,zp(i,j,kend+1))
zmin=AMIN1(zmin,zp(i,j,kbgn) )
END DO
END DO
IF( kbgn == 2 ) zmin=z(kbgn)
zmin=zmin / 1000.
zmax=zmax / 1000.
IF( zbgn /= zend ) THEN
zmin=zbgn
zmax=zend
ELSE
zmax=zp(ibgn,jbgn,kend+1)
zmin=zp(ibgn,jbgn,kbgn)
DO j=jbgn,jend
DO i=ibgn,iend
zmax=AMAX1(zmax,zp(i,j,kend+1))
zmin=AMIN1(zmin,zp(i,j,kbgn) )
END DO
END DO
IF( kbgn == 2 ) zmin=z(kbgn)
zmin=zmin / 1000.
zmax=zmax / 1000.
END IF
zr = zmax-zmin
x1 = xc(ibgn,2,2)-dxkm/2
x2 = x1+xr
y1 = yc(2,jbgn,2)-dykm/2
y2 = y1+yr
z1 = zmin
z2 = zmax
ibgn1 = ibgn
iend1 = iend
!
! 05/31/2002 Zuwen He
!
! For soil model, the zsoilmax and zsoilmin are the
! maximum and minimum of zpsoilc.
!
IF( zsoilbgn /= zsoilend ) THEN
zsoilmax=zsoilbgn
zsoilmin=zsoilend
ELSE
zsoilmax=zpsoilc(ibgn,jbgn,ksoilbgn)
zsoilmin=zpsoilc(ibgn,jbgn,ksoilend)
DO j=jbgn,jend
DO i=ibgn,iend
zsoilmax=AMAX1(zsoilmax,zpsoilc(i,j,ksoilbgn))
zsoilmin=AMIN1(zsoilmin,zpsoilc(i,j,ksoilend))
END DO
END DO
END IF
zsoilr = zsoilmax-zsoilmin
zsoil1 = zsoilmin
zsoil2 = zsoilmax
!
! END Zuwen
!
!-----------------------------------------------------------------------
!
! Re-Initialize the plotting space parameters, which might have
! been set differently for the vertical profile plotting.
!
!-----------------------------------------------------------------------
!
IF (profopt == 1) THEN
CALL xdspac(0.9*winsiz)
IF( nxpic == 1 .AND. nypic == 1) CALL xdspac(0.85*winsiz)
CALL xspace(nxpic, nypic, angl , xlimit,ylimit)
IF(axlbfmt == -1) THEN
CALL xaxfmt( '*' )
ELSE IF(axlbfmt == 0) THEN
CALL xaxfmt( '(i5)' )
ELSE
WRITE(stem1,'(i1)')axlbfmt
WRITE(stem2,43)stem1
CALL xaxfmt( stem2 )
END IF
layover = 0
CALL xpmagn(margnx*xlimit, margny*ylimit)
END IF
!
!-----------------------------------------------------------------------
!
! Set map
!
!-----------------------------------------------------------------------
!
! IF( ovrmap .eq. 1 .or. ovrobs .eq.1 .or. blcoplt.ne.0 .or.
! : ovrstaopt.ne.0 ) THEN
xl = (nx-3)*dxkm
yl = (ny-3)*dykm
xorig_0 = (xc(1,2,2)+xc(2,2,2))*0.5
yorig_0 = (yc(2,1,2)+yc(2,2,2))*0.5
CALL xstpjgrd(mapproj,trulat1,trulat2,trulon, &
ctrlat,ctrlon,xl,yl,xorig_0,yorig_0)
! ENDIF
!
!----------------------------------------------------------------------
!
! Calculate the coordinate value of pres_z to pres_val
! if presaxis_no >0
!
!----------------------------------------------------------------------
!
IF( presaxis_no > 0 ) &
CALL interp_p (pbar,zpc,ibgn,iend,nx,jbgn,jend,ny,kbgn,kend,nz)
!
!
!-----------------------------------------------------------------------
!
! Set terrain data
!
!-----------------------------------------------------------------------
!
DO j=1,ny
DO i=1,nx
hterain(i,j)=zp(i,j,2) ! Give zp(i,j,2) to hterain(i,j)
END DO
END DO
ztmin=hterain(1,1)
ztmax=ztmin
DO j=1,ny
DO i=1,nx
ztmax=MAX(ztmax,hterain(i,j))
ztmin=MIN(ztmin,hterain(i,j))
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate Coriolis parameter and map factor
!
!-----------------------------------------------------------------------
!
DO i=1,nx-1
xs(i) = (x(i)+x(i+1))*0.5
END DO
xs(nx) = -xs(nx-2)+2.0*xs(nx-1)
DO j=1,ny-1
ys(j) = (y(j)+y(j+1))*0.5
END DO
ys(ny) = -ys(ny-2)+2.0*ys(ny-1)
CALL xxytoll(nx,ny,xs,ys,tem8(1,1,1),tem8(1,1,2))
! print*,'sample lat, long = ', tem8(1,1,1), tem8(1,1,2)
r2d = ATAN(1.0)*4.0/180.0
omega2 = 2.0* omega
DO j=1,ny-1
DO i=1,nx-1
fcorio(i,j) = omega2*SIN( r2d * tem8(i,j,1) )
END DO
END DO
! print*,'sample fcorio =', fcorio(1,1)
DO i=1,nx-1
xp(i) = 0.5*(x(i)+x(i+1))
END DO
DO j=1,ny-1
yp(j) = 0.5*(y(j)+y(j+1))
END DO
CALL xytomf(nx-1,ny-1,xp,yp,mapfct)
!mx
! print*,'sample map factor =', mapfct(1,1)
!
!-----------------------------------------------------------------------
!
! Loop through slicopts, for different types of
! 2-d slices each time. If iplot(slicopt)=0, plotting for the
! given slicopt is switched off.
!
!-----------------------------------------------------------------------
!
DO slicopt = 1,maxslicopt ! Loop over different types of slices
!
! slicopt = 1: xy slice
! slicopt = 2: xz slice
! slicopt = 3: yz slice
! slicopt = 4: constant height
! slicopt = 5: vertical slice through two points
! slicopt = 6: constant pressure
! slicopt = 7: isentropic level
! slicopt = 8: ??
! slicopt = 9: xy-soil-slice for the soil model
! slicopt = 10: xz-soil-slice for the soil model
! slicopt = 11: yz-soil-slice for the soil model
! reset the windows
iend = iend1
ibgn = ibgn1
xr = xc(iend,2,2)-xc(ibgn,2,2)+dxkm
yr = yc(2,jend,2)-yc(2,jbgn,2)+dykm
x1 = xc(ibgn,2,2)-dxkm/2
x2 = x1+xr
y1 = yc(2,jbgn,2)-dykm/2
y2 = y1+yr
z1 = zmin
z2 = zmax
zsoil1 = zsoilmin
zsoil2 = zsoilmax
IF(slicopt == 1.OR.slicopt == 4.OR.slicopt == 6.OR. &
slicopt == 7.OR.slicopt == 8.OR.slicopt == 9) THEN
aspect=xr/yr
ELSE IF (slicopt == 2) THEN
aspect=xr/zr
ELSE IF( slicopt == 5) THEN
! aspect=
ELSE IF (slicopt == 10) THEN
aspect=xr/zsoilr
ELSE IF (slicopt == 11) THEN
aspect=yr/zsoilr
ELSE
aspect=yr/zr
END IF
aspratio = aspect
IF(slicopt == 1.AND.yxstrch /= 0.0) aspratio=yxstrch
IF(slicopt == 2.AND.zxstrch /= 0.0) aspratio=zxstrch
IF(slicopt == 3.AND.zystrch /= 0.0) aspratio=zystrch
IF(slicopt == 4.AND.yxstrch /= 0.0) aspratio=yxstrch
IF(slicopt == 5.AND.zhstrch /= 0.0) aspratio=zhstrch
IF(slicopt == 6.AND.yxstrch /= 0.0) aspratio=yxstrch
IF(slicopt == 7.AND.yxstrch /= 0.0) aspratio=yxstrch
IF(slicopt == 8.AND.yxstrch /= 0.0) aspratio=yxstrch
IF(slicopt == 9.AND.yxstrch /= 0.0) aspratio=yxstrch
!
!-----------------------------------------------------------------------
!
! Loop over all slices (nslice(slicopt)) for each slice type
!
!-----------------------------------------------------------------------
!
DO indxslic = 1, nslice(slicopt)
! print*,'slicopt=',slicopt
! print*,'indxslic =',indxslic ,'nslice(slicopt)',nslice(slicopt)
!
! 06/03/2002 Zuwen He Reorganized slice indexes.
!
IF (iplot(slicopt) /= 0) THEN
!
! Regular vertical cross-sections, i.e., i or j is fixed
!
IF (slicopt == 2) jslice=slice_xz(indxslic)
IF (slicopt == 3) islice=slice_yz(indxslic)
IF (jslice == -1) jslice = (ny-2)/2+1
IF (jslice == -2) jslice = jmax
IF (islice == -1) islice = (nx-2)/2+1
IF (islice == -2) islice = imax
!
! Horizontal plots, i.e., k is fixed
!
IF( slicopt == 1) kslice=slice_xy(indxslic)
IF( kslice == -1) kslice = (nz-2)/2+1
IF( kslice == -2) kslice = kmax
!
! Other plots
!
IF( slicopt == 4) zlevel=slice_h (indxslic)
IF( slicopt <= 5) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
zps3d(i,j,k) = zpc(i,j,k)
END DO
END DO
END DO
END IF
IF( slicopt == 6) THEN
zlevel=-ALOG(100.0*slice_p(indxslic))
DO k=1,nz
DO j=1,ny
DO i=1,nx
zps3d(i,j,k) = algpzc(i,j,k)
END DO
END DO
END DO
END IF
IF( slicopt == 7) THEN
zlevel=slice_pt(indxslic)
DO k=1,nz
DO j=1,ny
DO i=1,nx
zps3d(i,j,k) = ptbar(i,j,k)+ptprt(i,j,k)
END DO
END DO
END DO
END IF
IF( slicopt == 5) THEN
x01=xi1(indxslic)
y01=yi1(indxslic)
x02=xi2(indxslic)
y02=yi2(indxslic)
CALL clipwd(x01,y01,x02,y02, idisplay )
IF(idisplay == 0) THEN
WRITE(6,'(2(/1x,a))') &
'The specified vertical slice was outside the plotting ', &
'window, no cross-section is plotted.'
CYCLE
END IF
sqrtdxy=SQRT(dxkm*dykm)
dist=SQRT((x01-x02)**2+(y01-y02)**2)
n=INT(dist/sqrtdxy+0.5)+1
sinaf=(y02-y01)/dist
cosaf=(x02-x01)/dist
DO i=1,n
xp(i)=x01+(FLOAT(i-1))*sqrtdxy*cosaf
yp(i)=y01+(FLOAT(i-1))*sqrtdxy*sinaf
END DO
ibgn = 1
iend = n
xr = (n-1)*sqrtdxy
z1=zp(ibgn,jbgn,kbgn)
z2=zp(ibgn,jbgn,kend+1)
DO j=jbgn,jend
DO i=ibgn,iend
z1=AMIN1(z1,zp(i,j,kbgn))
z2=AMAX1(z2,zp(i,j,kend+1))
END DO
END DO
z1=z1 / 1000.
z2=z2 / 1000.
IF( zbgn /= zend ) THEN
z1 = MAX(zbgn,z1)
z2 = MIN(zend,z2)
END IF
zr = z2-z1
x1 = 0.0
x2 = x1+xr
aspratio = xr/zr
END IF
!
! slices for soil variables, tsoil and qsoil
!
IF (slicopt == 9) THEN
ksoilslice=slice_xy_soil(indxslic)
IF (ksoilslice == -1) ksoilslice = (nzsoil-2)/2+1
DO k=1,nzsoil
DO j=1,ny
DO i=1,nx
zpsoils3d(i,j,k) = zpsoilc(i,j,k)
END DO
END DO
END DO
END IF
IF (slicopt == 10) jsoilslice=slice_xz_soil(indxslic)
IF (slicopt == 11) isoilslice=slice_yz_soil(indxslic)
IF (jsoilslice == -1) jsoilslice = (ny-2)/2+1
IF (isoilslice == -1) isoilslice = (nx-2)/2+1
END IF
!
! END Zuwen
!
CALL xlbint(3)
CALL xczero(3)
!
!-----------------------------------------------------------------------
!
! Loop for every priority. small number is high priority, doing high
! priority first , then low priority, finally doing the zero priority.
! ipriority(): 1-50 3D plot, 51-100 2D plot, 101-150 surface character.
! 151-170 arbitrary 3D plot, 171-189 arbitrary 2D plot
!
!-----------------------------------------------------------------------
!
DO ipp = 1,nprio ! do loop priority
sovrlay = 0
plotovr = .false.
ip=iptemp(ipp)
!
!-----------------------------------------------------------------------
!
! Start plotting for all slicopt
!
!-----------------------------------------------------------------------
!
IF( iplot(slicopt) /= 0 .AND. slicopt <= 7 ) THEN
!
!-----------------------------------------------------------------------
!
! Plot height(10m) on constant pressure level
! z-coor of sacalar point in physical space (10m)
!
!-----------------------------------------------------------------------
!
IF(ipriority(1) == ip ) THEN
IF(hplot == 1 .OR. hplot == 2 .OR. hplot == 4 .OR. &
hplot == 5) THEN
CALL get_mulovrlay('hplot',5,ovrmul_num,ovrmulname,sovrlay)
IF(slicopt == 6 .OR. slicopt == 7 ) THEN
CALL filzero(tem9,nx*ny*nz)
DO k=1,nz-1
DO j=1,ny
DO i=1,nx
tem9(i,j,k)=100.0*zpc(i,j,k)
END DO
END DO
END DO
CALL ctrsetup(hinc,hminc,hmaxc, &
hovr,hhlf,hzro,hcol1,hcol2,'hplot ')
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'h (10m)',time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,hplot)
END IF
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot thickness in 10meters between a given pressure level and
! 850mb level
!
!-----------------------------------------------------------------------
!
IF(ipriority(110) == ip ) THEN
IF(thkplt == 1 .OR. thkplt == 2 .OR. thkplt == 4 .OR. &
thkplt == 5) THEN
CALL get_mulovrlay('thkplt',5,ovrmul_num,ovrmulname,sovrlay)
IF(slicopt == 6 .OR. slicopt == 7 ) THEN
CALL hintrp1 &
(nx,ny,nz,2,nz-2,zpc ,zps3d,zlevel,tem9(1,1,1))
pref = 800.0 ! mb
IF( slicopt == 6) THEN
tem=-ALOG(100.0*pref)
END IF
IF( slicopt == 7) THEN
tem = 300.0 ! K
END IF
CALL hintrp1 &
(nx,ny,nz,2,nz-2,zpc ,zps3d,tem,tem9(1,1,2))
DO j=1,ny-1
DO i=1,nx-1
IF( ABS(tem9(i,j,1)+9999.0) < 0.1 .OR. &
ABS(tem9(i,j,2)+9999.0) < 0.1 ) THEN
tem9(i,j,1)=-9999.0
ELSE
tem9(i,j,1)=100.0*(tem9(i,j,1)-tem9(i,j,2))
END IF
END DO
END DO
CALL ctrsetup(thkinc,thkminc,thkmaxc,thkovr, &
thkhlf,thkzro,thkcol1,thkcol2,'thkplt ')
WRITE(label,'(i4,''-'',I4,''MB THICKNESS (10M)'')') &
nint(slice_p(indxslic)),nint(pref)
CALL ctrsfc(tem9(1,1,1),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm,nx,ibgn,iend,ny,jbgn,jend, &
label(1:27),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,thkplt)
END IF
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot ageostrophic winds on pressure surfaces
!
!-----------------------------------------------------------------------
!
IF(ipriority(109) == ip ) THEN
IF(vagplt == 1 .OR. vagplt == 2 .OR. vagplt == 4 .OR. &
vagplt == 5) THEN
CALL get_mulovrlay('vagplt',5,ovrmul_num,ovrmulname,sovrlay)
IF(slicopt == 6 ) THEN ! Pressue surface plot
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=(u(i+1,j,k)+u(i,j,k))*0.5
tem2(i,j,k)=(v(i,j+1,k)+v(i,j,k))*0.5
END DO
END DO
END DO
CALL hintrp1 &
(nx,ny,nz,2,nz-2,tem1,zps3d,zlevel,tem9(1,1,1))
CALL hintrp1 &
(nx,ny,nz,2,nz-2,tem2,zps3d,zlevel,tem9(1,1,2))
CALL hintrp1 &
(nx,ny,nz,2,nz-2,zpc ,zps3d,zlevel,tem9(1,1,3))
! do i=1,20
! CALL smooth9pmv(tem9(1,1,3),nx,ny,2,nx-2,2,ny-2,tem5)
! enddo
! print*,'sample p-level height =', tem9(nx/2,ny/2,3)
! print*,'dxkm, dykm =', dxkm, dykm
DO j=2,ny-2
DO i=1,nx-1
IF(ABS(tem9(i,j ,1)+9999.0) < 0.1.OR. &
ABS(tem9(i,j+1,3)+9999.0) < 0.1.OR. &
ABS(tem9(i,j-1,3)+9999.0) < 0.1 )THEN
tem8(i,j,1)=-9999.0
ELSE
tem8(i,j,1)= &
tem9(i,j,1)+ mapfct(i,j)* &
g*(tem9(i,j+1,3)-tem9(i,j-1,3))/(fcorio(i,j)*2*dykm)
END IF
END DO
END DO
DO j=1,ny-1
DO i=2,nx-2
IF(ABS(tem9(i,j ,1)+9999.0) < 0.1.OR. &
ABS(tem9(i+1,j,3)+9999.0) < 0.1.OR. &
ABS(tem9(i-1,j,3)+9999.0) < 0.1)THEN
tem8(i,j,2)=-9999.0
ELSE
tem8(i,j,2)= &
tem9(i,j,2)- mapfct(i,j)* &
g*(tem9(i+1,j,3)-tem9(i-1,j,3))/(fcorio(i,j)*2*dxkm)
END IF
END DO
END DO
CALL vtrunt ( vagunit )
IF(plotovr) CALL overlay ( 1 )
IF(.NOT.plotovr) CALL overlay( vagovr )
CALL ctrcol(vagcol1,vagcol2)
CALL varplt( 'vagplt ' )
CALL ctrvtr( vagunits, vagtype )
WRITE(label,'(''AG AT '',I4,''(HPA)'')') &
nint(slice_p(indxslic))
CALL vtrsfc( tem8(1,1,1),tem8(1,1,2),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, &
label(1:15),time, runname,1.0,slicopt, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
sigplt(109) = 1
END IF
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot temperature (C)
!
!-----------------------------------------------------------------------
!
IF(ipriority(2) == ip ) THEN
IF( tplot > 0) THEN
CALL cal_t ( tem9,tz,nx, ny, nz,tob, label,length,tunits )
END IF
IF(tplot == 1 .OR. tplot == 2 .OR. tplot == 4 .OR. &
tplot == 5 ) THEN
CALL get_mulovrlay('tplot',5,ovrmul_num,ovrmulname,sovrlay)
! pltvar = 'tplot'
! CALL spltpara(tinc, tminc, tmaxc, tovr, thlf, tzro,
! : tcol1,tcol2,pltvar)
CALL ctrsetup(tinc,tminc,tmaxc, &
tovr,thlf,tzro,tcol1,tcol2,'tplot ')
CALL xcfull
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:length),time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,tplot)
CALL xcmixl
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot wind components
!
!-----------------------------------------------------------------------
!
IF(ipriority(3) == ip ) THEN
IF(uplot == 1 .OR. uplot == 2 .OR. uplot == 4 .OR. &
uplot == 5) THEN
CALL get_mulovrlay('uplot',5,ovrmul_num,ovrmulname,sovrlay)
IF(nobs > 0 .AND. slicopt == 1 .AND. kslice == 2) THEN
DO iob=1,nobs
IF(dd(iob) >= 0. .AND. dd(iob) < 360. .AND. &
ff(iob) >= 0. .AND. ff(iob) < 60.) THEN
CALL ddrotuv(1,lonob(iob),dd(iob),ff(iob), &
drot,obs1(iob),obs2(iob))
obs1(iob)=0.51444*obs1(iob)
ELSE
obs1(iob)=-999.
END IF
END DO
obsset=1
END IF
onvf = 0
CALL avgx(u , onvf, &
nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)
CALL ctrsetup(uinc,uminc,umaxc, &
uovr,uhlf,uzro,ucol1,ucol2,'uplot ')
CALL ctr3d( tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'u (m/s)',time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,uplot)
END IF
END IF
IF(ipriority(4) == ip ) THEN
IF(vplot == 1 .OR. vplot == 2 .OR. vplot == 4 .OR. &
vplot == 5 ) THEN
CALL get_mulovrlay('vplot',5,ovrmul_num,ovrmulname,sovrlay)
IF(nobs > 0 .AND. slicopt == 1 .AND. kslice == 2) THEN
DO iob=1,nobs
IF(dd(iob) >= 0. .AND. dd(iob) < 360. .AND. &
ff(iob) >= 0. .AND. ff(iob) < 60.) THEN
CALL ddrotuv(1,lonob(iob),dd(iob),ff(iob), &
drot,obs1(iob),obs2(iob))
obs1(iob)=0.51444*obs2(iob)
ELSE
obs1(iob)=-999.
END IF
END DO
obsset=1
END IF
onvf = 0
CALL avgy(v,onvf,nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem9)
CALL ctrsetup(vinc,vminc,vmaxc, &
vovr,vhlf,vzro,vcol1,vcol2,'vplot ')
CALL ctr3d(tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'v (m/s)',time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,vplot)
END IF
END IF
IF(ipriority(5) == ip ) THEN
IF(vhplot > 0) THEN
CALL cal_vh(tem9,u,v,nx,ny,nz,vhunits,label,length,tem8)
END IF
IF(vhplot == 1 .OR. vhplot == 2 .OR. vhplot == 4 .OR. &
vhplot == 5 .OR.vhplot == 6 ) THEN
CALL get_mulovrlay('vhplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(vhinc,vhminc,vhmaxc, &
vhovr,vhhlf,vhzro,vhcol1,vhcol2,'vhplot ')
CALL ctr3d(tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:length),time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,vhplot)
END IF
END IF
IF(ipriority(6) == ip ) THEN
IF(wplot == 1 .OR. wplot == 2 .OR. wplot == 4 .OR. &
wplot == 5 ) THEN
CALL get_mulovrlay('wplot',5,ovrmul_num,ovrmulname,sovrlay)
onvf = 0
CALL avgz(w,onvf,nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem9)
CALL ctrsetup(winc,wminc,wmaxc, &
wovr,whlf,wzro,wcol1,wcol2,'wplot ')
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'w (m/s)',time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,wplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot scalars
!
!-----------------------------------------------------------------------
!
IF(ipriority(7) == ip ) THEN
IF(ptplot == 1 .OR. ptplot == 2 .OR. ptplot == 4 .OR. &
ptplot == 5 ) THEN
CALL get_mulovrlay('ptplot',6,ovrmul_num,ovrmulname,sovrlay)
IF(ovrobs == 1 .AND. nobs > 0) THEN
DO iob=1,nobs
!
! Station pressure (mb)
!
IF(pstn(iob) > 0.) THEN
psta=mbtopa*pstn(iob)
ELSE IF(alt(iob) > 0.) THEN
psta=mbtopa*alttostpr(alt(iob),elevob(iob))
ELSE
psta=mbtopa*alttostpr(1013.,elevob(iob))
END IF
!
! Potential temperature (K)
!
IF(tob(iob) > -98.) THEN
tk=(5.*(tob(iob)-32.)/9.) + 273.15
obs1(iob)=tk*((p0/psta)**rddcp)
ELSE
obs1(iob)=-999.
END IF
END DO
obsset=1
END IF
CALL ctrsetup(ptinc,ptminc,ptmaxc, &
ptovr,pthlf,ptzro,ptcol1,ptcol2,'ptplot ')
CALL ctr3d( pt ,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'pt (K)',time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,ptplot)
END IF
END IF
IF(ipriority(8) == ip ) THEN
IF(pplot == 1 .OR. pplot == 2 .OR. pplot == 4 .OR. &
pplot == 5 ) THEN
CALL get_mulovrlay('pplot',5,ovrmul_num,ovrmulname,sovrlay)
IF( ovrobs == 1 .AND. nobs > 0) THEN
!
DO iob=1,nobs
!
! Station pressure (mb)
!
IF(pstn(iob) > 0.) THEN
obs1(iob)=pstn(iob)-100.*(MOD(INT(pstn(iob)),100))
ELSE IF(alt(iob) > 0.) THEN
pmb=alttostpr(alt(iob),elevob(iob))
obs1(iob)=pmb -100.*(MOD(INT(pmb),100))
ELSE
obs1(iob)=-999.
END IF
END DO
obsset=1
END IF
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem9(i,j,k)=pbar(i,j,k)+pprt (i,j,k)
END DO
END DO
END DO
CALL ctrsetup(pinc,pminc,pmaxc, &
povr,phlf,pzro,pcol1,pcol2,'pplot ')
CALL ctr3d( tem9 ,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'p (Pa)',time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,pplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot scalars for qv, qc,qr,qi,qs,qh,qw
!
!-----------------------------------------------------------------------
!
IF(ipriority(9) == ip) THEN
IF(qvplot == 1 .OR. qvplot == 2 .OR. qvplot == 4 .OR. &
qvplot == 5 ) THEN
CALL get_mulovrlay('qvplot',6,ovrmul_num,ovrmulname,sovrlay)
IF( ovrobs == 1 .AND. nobs > 0 ) THEN
DO iob=1,nobs
!
! Station pressure
!
IF(pstn(iob) > 0.) THEN
psta=mbtopa*pstn(iob)
ELSE IF(alt(iob) > 0.) THEN
psta=mbtopa*alttostpr(alt(iob),elevob(iob))
ELSE
psta=mbtopa*alttostpr(1013.,elevob(iob))
END IF
!
! Specific humidity
!
IF(tdob(iob) > -90.) THEN
tdk=(5.*(tdob(iob)-32.)/9.) + 273.15
obs1(iob) = f_qvsat( psta,tdk ) *1000.
ELSE
obs1(iob)=-999.
END IF
END DO
obsset=1
END IF
CALL ctrsetup(qvinc,qvminc,qvmaxc, &
qvovr,qvhlf,qvzro,qvcol1,qvcol2,'qvplot ')
CALL ctr3d( qv,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'qv (g/kg)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qvplot)
END IF
END IF
IF(ipriority(10) == ip) THEN
IF(qcplot == 1 .OR. qcplot == 2 .OR. qcplot == 4 .OR. &
qcplot == 5 ) THEN
CALL get_mulovrlay('qcplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(qcinc,qcminc,qcmaxc, &
qcovr,qchlf,qczro,qccol1,qccol2,'qcplot ')
CALL ctr3d( qc,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'qc (g/kg)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qcplot)
END IF
END IF
!
IF(ipriority(11) == ip) THEN
IF(qrplot == 1 .OR. qrplot == 2 .OR. qrplot == 4 .OR. &
qrplot == 5 ) THEN
CALL get_mulovrlay('qrplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(qrinc,qrminc,qrmaxc, &
qrovr,qrhlf,qrzro,qrcol1,qrcol2,'qrplot ')
CALL ctr3d( qr,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'qr (g/kg)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qrplot)
END IF
END IF
IF(ipriority(12) == ip ) THEN
IF(qiplot == 1 .OR. qiplot == 2 .OR. qiplot == 4 .OR. &
qiplot == 5 ) THEN
CALL get_mulovrlay('qiplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(qiinc,qiminc,qimaxc, &
qiovr,qihlf,qizro,qicol1,qicol2,'qiplot ')
CALL ctr3d( qi,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'qi (g/kg)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qiplot)
END IF
END IF
IF(ipriority(13) == ip ) THEN
IF(qsplot == 1 .OR. qsplot == 2 .OR. qsplot == 4 .OR. &
qsplot == 5 ) THEN
CALL get_mulovrlay('qsplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(qsinc,qsminc,qsmaxc, &
qsovr,qshlf,qszro,qscol1,qscol2,'qsplot ')
CALL ctr3d( qs,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'qs (g/kg)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qsplot)
END IF
END IF
IF(ipriority(14) == ip ) THEN
IF(qhplot == 1 .OR. qhplot == 2 .OR. qhplot == 4 .OR. &
qhplot == 5 ) THEN
CALL get_mulovrlay('qhplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(qhinc,qhminc,qhmaxc, &
qhovr,qhhlf,qhzro,qhcol1,qhcol2,'qhplot ')
CALL ctr3d( qh,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'qh (g/kg)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qhplot)
END IF
END IF
IF(ipriority(15) == ip ) THEN
IF(qwplot > 0) THEN
CALL cal_qw(tem9,qc,qr,qi,qs,qh, nx,ny,nz)
END IF
IF(qwplot == 1 .OR. qwplot == 2 .OR. qwplot == 4 .OR. &
qwplot == 5 ) THEN
CALL get_mulovrlay('qwplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(qwinc,qwminc,qwmaxc, &
qwovr,qwhlf,qwzro,qwcol1,qwcol2,'qwplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'Total water (g/kg)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qwplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate relative humidity, where tem1 = temperature,
! tem2 = saturation qv.
!
!-----------------------------------------------------------------------
!
IF(ipriority(16) == ip ) THEN
IF(kmhplt == 1 .OR. kmhplt == 2 .OR. kmhplt == 4 .OR. &
kmhplt == 5 ) THEN
CALL get_mulovrlay('kmhplt',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(kmhinc,kmhminc,kmhmaxc, &
kmhovr,kmhhlf,kmhzro,kmhcol1,kmhcol2,'kmhplot ')
CALL ctr3d( kmh,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'kmh (m**2/s)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,kmhplt)
END IF
END IF
IF(ipriority(17) == ip ) THEN
IF(kmvplt == 1 .OR. kmvplt == 2 .OR. kmvplt == 4 .OR. &
kmvplt == 5 ) THEN
CALL get_mulovrlay('kmvplt',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(kmvinc,kmvminc,kmvmaxc, &
kmvovr,kmvhlf,kmvzro,kmvcol1,kmvcol2,'kmvplot ')
CALL ctr3d( kmv,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'kmv (m**2/s)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,kmvplt)
END IF
END IF
IF(ipriority(18) == ip ) THEN
IF(tkeplt == 1 .OR. tkeplt == 2 .OR. tkeplt == 4 .OR. &
tkeplt == 5 ) THEN
CALL get_mulovrlay('tkeplt',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(tkeinc,tkeminc,tkemaxc, &
tkeovr,tkehlf,tkezro,tkecol1,tkecol2,'tkeplot ')
CALL ctr3d( tke,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'tke ((m/s)**2)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,tkeplt)
END IF
END IF
IF(ipriority(19) == ip) THEN
IF(rhplot > 0) THEN
CALL cal_rh(tem9,pt, pprt ,pbar,qv,tem1,tem2,nx,ny,nz)
END IF
IF(rhplot == 1 .OR. rhplot == 2 .OR. rhplot == 4 .OR. &
rhplot == 5 ) THEN
CALL get_mulovrlay('rhplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL ctrsetup(rhinc,rhminc,rhmaxc, &
rhovr,rhhlf,rhzro,rhcol1,rhcol2,'rhplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'RH',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,rhplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot dew-point temperature td tdplot
!
!-----------------------------------------------------------------------
IF(ipriority(20) == ip ) THEN
IF(tdplot > 0) THEN
CALL cal_td(tem9,td,nx,ny,nz,tdunits,label, length)
END IF
IF(tdplot == 1 .OR. tdplot == 2 .OR. tdplot == 4 .OR. &
tdplot == 5 ) THEN
CALL get_mulovrlay('tdplot',6,ovrmul_num,ovrmulname,sovrlay)
CALL cal_tdobs(tdob,tdunits)
CALL ctrsetup(tdinc,tdminc,tdmaxc, &
tdovr,tdhlf,tdzro,tdcol1,tdcol2,'tdplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:length),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,tdplot)
END IF
END IF
IF(ipriority(21) == ip) THEN
IF(rfplot == 1 .OR. rfplot == 2 .OR. rfplot == 4 .OR. &
rfplot == 5 ) THEN
CALL get_mulovrlay('rfplot',6,ovrmul_num,ovrmulname,sovrlay)
IF (rfopt == 1) THEN
CALL reflec_ferrier(nx,ny,nz, rhobar, qr, qs, qh, tz, tem9)
ELSE
CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, tem9)
ENDIF
CALL ctrsetup(rfinc,rfminc,rfmaxc, &
rfovr,rfhlf,rfzro,rfcol1,rfcol2,'rfplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
'Ref (dBZ)',time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,rfplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate composite reflectivity
!
!-----------------------------------------------------------------------
!
IF(ipriority(36) == ip) THEN
IF(rfcplt > 0) THEN
IF (rfopt == 1) THEN
CALL reflec_ferrier(nx,ny,nz, rhobar, qr, qs, qh, tz, tem8)
ELSE
CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, tem8)
ENDIF
CALL cal_rfc(nx, ny, nz, tem8, tem9)
END IF
IF(rfcplt == 1 .OR. rfcplt == 2 .OR. rfcplt == 4 .OR. &
rfcplt == 5 ) THEN
CALL get_mulovrlay('rfcplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Composite Ref (dBZ)'
CALL ctrsetup(rfcinc,rfcminc,rfcmaxc, &
rfcovr,rfchlf,rfczro,rfccol1,rfccol2,'rfcplot ')
CALL ctrsfc(tem9,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time, runname, 1.0,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,rfcplt)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate equivalent potential temperature.
!
!-----------------------------------------------------------------------
!
IF(ipriority(22) == ip ) THEN
IF(pteplt == 1 .OR. pteplt == 2 .OR. pteplt == 4 .OR. &
pteplt == 5 ) THEN
CALL get_mulovrlay('pteplt',6,ovrmul_num,ovrmulname,sovrlay)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=pprt(i,j,k)+pbar(i,j,k)
END DO
END DO
END DO
CALL pt2pte(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem1,pt,qv,tem9)
label = 'pte (K)'
CALL ctrsetup(pteinc,pteminc,ptemaxc, &
pteovr,ptehlf,ptezro,ptecol1,ptecol2,'pteplot ')
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:7),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,pteplt)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot perturbation of wind components
!
!-----------------------------------------------------------------------
!
IF(ipriority(23) == ip ) THEN
IF(upplot > 0) THEN
onvf = 0
CALL avgx(uprt , onvf, &
nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)
END IF
IF(upplot == 1 .OR. upplot == 2 .OR. upplot == 4 .OR. &
upplot == 5 ) THEN
CALL get_mulovrlay('upplot',6,ovrmul_num,ovrmulname,sovrlay)
label = 'uprt (m/s)'
CALL ctrsetup(upinc,upminc,upmaxc, &
upovr,uphlf,upzro,upcol1,upcol2,'upplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:10),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,upplot)
END IF
END IF
IF(ipriority(24) == ip ) THEN
IF(vpplot > 0) THEN
onvf = 0
CALL avgy(vprt , onvf, &
nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)
END IF
IF(vpplot == 1 .OR. vpplot == 2 .OR. vpplot == 4 .OR. &
vpplot == 5 ) THEN
CALL get_mulovrlay('vpplot',6,ovrmul_num,ovrmulname,sovrlay)
label = 'vprt (m/s)'
CALL ctrsetup(vpinc,vpminc,vpmaxc, &
vpovr,vphlf,vpzro,vpcol1,vpcol2,'vpplot ')
CALL ctr3d( tem9 ,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:10),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,vpplot)
END IF
END IF
IF(ipriority(25) == ip ) THEN
IF(wpplot > 0) THEN
onvf = 0
CALL avgz(wprt , onvf, &
nx,ny,nz, 1,nx-1, 1,ny-1, 1,nz-1, tem9)
END IF
IF(wpplot == 1 .OR. wpplot == 2.OR. wpplot == 4 .OR. &
wpplot == 5) THEN
CALL get_mulovrlay('wpplot',6,ovrmul_num,ovrmulname,sovrlay)
label = 'wprt (m/s)'
CALL ctrsetup(wpinc,wpminc,wpmaxc, &
wpovr,wphlf,wpzro,wpcol1,wpcol2,'wpplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:10),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,wpplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot perturbation scalars, calculating and storing
! perturbations in tem4, where necessary
!
!-----------------------------------------------------------------------
!
IF(ipriority(26) == ip ) THEN
IF(ptpplt == 1 .OR. ptpplt == 2.OR. ptpplt == 4 .OR. &
ptpplt == 5 ) THEN
if(.true.) then ! plot ptprt or buoyancy including water loading
CALL get_mulovrlay('ptpplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'ptprt (K)'
CALL ctrsetup(ptpinc,ptpminc,ptpmaxc, &
ptpovr,ptphlf,ptpzro,ptpcol1,ptpcol2,'ptpplot ')
CALL ctr3d(ptprt,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:10),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,ptpplt)
else
call BUOYCY_plt(nx,ny,nz,ptprt,pprt,qv,qc,qr,qi,qs,qh, &
ptbar,pbar,rhobar,qvbar, tem6, tem1)
do k=1,nz-1
do j=1,ny-1
do i=1,nx-1
tem6(i,j,k) = tem6(i,j,k)/(rhobar(i,j,k)*g) *ptbar(i,j,k)
enddo
enddo
enddo
CALL get_mulovrlay('Buoyancy',6,ovrmul_num,ovrmulname,sovrlay)
label = 'T-Equivalent Buoyancy (K)'
CALL ctrsetup(ptpinc,ptpminc,ptpmaxc, &
ptpovr,ptphlf,ptpzro,ptpcol1,ptpcol2,'ptpplot ')
CALL ctr3d(tem6,xc,yc,zps3d,x1,x2,dx,y1,y2,dy,z1,z2,dz, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:17),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,ptpplt)
endif
END IF
END IF
IF(ipriority(27) == ip ) THEN
IF(ppplot == 1 .OR. ppplot == 2 .OR. ppplot == 4 .OR. &
ppplot == 5 ) THEN
CALL get_mulovrlay('ppplot',6,ovrmul_num,ovrmulname,sovrlay)
label = 'pprt (Pa)'
CALL ctrsetup(ppinc,ppminc,ppmaxc, &
ppovr,pphlf,ppzro,ppcol1,ppcol2,'ppplot ')
CALL ctr3d(pprt ,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:10),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,ppplot)
END IF
END IF
IF(ipriority(28) == ip ) THEN
IF(qvpplt == 1 .OR. qvpplt == 2 .OR. qvpplt == 4 .OR. &
qvpplt == 5 ) THEN
CALL get_mulovrlay('qvpplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'qvprt (g/kg)'
CALL ctrsetup(qvpinc,qvpminc,qvpmaxc, &
qvpovr,qvphlf,qvpzro,qvpcol1,qvpcol2,'qvpplt ')
CALL ctr3d(qvprt,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:12),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qvpplt)
END IF
END IF
IF(ipriority(29) == ip ) THEN
IF(vorpplt > 0) THEN
CALL cal_vorp(tem9,u,v,x,y,nx,ny,nz)
END IF
IF(vorpplt == 1 .OR. vorpplt == 2 .OR. vorpplt == 4 .OR. &
vorpplt == 5 ) THEN
CALL get_mulovrlay('vorpplt',7,ovrmul_num,ovrmulname, &
sovrlay)
label = 'Vort*10^5 (1/s)'
CALL ctrsetup(vorpinc,vorpminc,vorpmaxc, &
vorpovr,vorphlf,vorpzro,vorpcol1,vorpcol2,'vorpplt ')
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:15),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3,tem4, &
tem5,tem6,hterain,vorpplt)
END IF
END IF
IF(ipriority(30) == ip ) THEN
IF(divpplt > 0) THEN
CALL cal_div(tem9,u,v,x,y,nx,ny,nz)
END IF
IF( divpplt == 1 .OR. divpplt == 2 .OR. divpplt == 4 .OR. &
divpplt == 5 ) THEN
label = 'Div.*1000 (1/s)'
CALL ctrsetup(divpinc,divpminc,divpmaxc, &
divpovr,divphlf,divpzro,divpcol1,divpcol2,'divpplt ')
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:15),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3,tem4, &
tem5,tem6,hterain,divpplt)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot divergence of moist (qv*u )
!
!-----------------------------------------------------------------------
!
IF(ipriority(31) == ip ) THEN
IF(divqplt > 0) THEN
CALL cal_divq(tem9,u,v,qv,x,y,nx,ny,nz)
END IF
IF (divqplt == 1 .OR. divqplt == 2 .OR. divqplt == 4 .OR. &
divqplt == 5 ) THEN
CALL get_mulovrlay('divqplt',7,ovrmul_num,ovrmulname, &
sovrlay)
label = 'Moist Conv.*1000. (g/kg/s)'
CALL ctrsetup(divqinc,divqminc,divqmaxc, &
divqovr,divqhlf,divqzro,divqcol1,divqcol2,'divqplt ')
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:26),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3,tem4, &
tem5,tem6,hterain,divqplt)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate and plot perturbation wind components
!
!-----------------------------------------------------------------------
!
IF(ipriority(33) == ip ) THEN
IF(vtpplt == 1) THEN
CALL get_mulovrlay('vtpplt',6,ovrmul_num,ovrmulname,sovrlay)
CALL vtrunt ( vtpunt )
CALL overlay(vtpovr )
CALL ctrcol(vtpcol1,vtpcol2)
CALL varplt( 'vtpplt ' )
CALL ctrvtr( vtpunits, vtptype )
CALL cal_vtp(tem7,tem8,tem9,uprt,vprt,wprt,nx,ny,nz, &
vtpunits,label,length)
CALL vtr3d(tem7, tem8, tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, nz,kbgn,kend,kst, &
kslice,jslice,islice, label(1:length),time, runname,1.0, &
slicopt,n,xp,yp,zs2,u1,v1,u2,v2,w2, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot streamline field
!
!-----------------------------------------------------------------------
!
IF(ipriority(34) == ip ) THEN
IF(vtrstrm == 1) THEN
CALL get_mulovrlay('vtrstrm',7,ovrmul_num,ovrmulname, &
sovrlay)
CALL cal_vtrstrm(tem7,tem8,tem9,u,v,w,nx,ny,nz,aspratio)
CALL overlay( vtrstmovr )
CALL ctrcol(vtrstmcol1,vtrstmcol2)
CALL varplt( 'vtrstrm ' )
CALL strm3d( tem7 , tem8 , tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, nz,kbgn,kend,kst, &
kslice,jslice,islice, time, runname,1.0, slicopt, &
n,xp,yp,zs2,u1,v1,u2,v2,w2, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Calculate and plot the perturbation of wind streamlines
!
!-----------------------------------------------------------------------
!
IF(ipriority(35) == ip ) THEN
IF(vtpstrm == 1) THEN
CALL get_mulovrlay('vtpstrm',7,ovrmul_num,ovrmulname, &
sovrlay)
CALL cal_vtpstrm(tem7,tem8,tem9,uprt,vprt,wprt,nx,ny,nz, &
aspratio)
CALL overlay( vtpstmovr )
CALL varplt( 'vtpstrm ' )
CALL strm3d( tem7 , tem8 , tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, nz,kbgn,kend,kst, &
kslice,jslice,islice, time, runname,1.0, slicopt, &
n,xp,yp,zs2,u1,v1,u2,v2,w2, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot vertical wind shear, SQRT[(du/dz)^2 + (dv/dz)^2].
!
!-----------------------------------------------------------------------
!
!
IF(ipriority(37) == ip ) THEN
IF(vsplot > 0) THEN
CALL cal_vs(tem9,u,v,zp,tem7,tem8,nx,ny,nz)
END IF
IF(vsplot == 1 .OR. vsplot == 2 .OR. vsplot == 4 .OR. &
vsplot == 5 ) THEN
CALL get_mulovrlay('vsplot',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Vertical wind shear*1000(1/s)'
CALL ctrsetup(vsinc,vsminc,vsmaxc, &
vsovr,vshlf,vszro,vscol1,vscol2,'vsplot ')
CALL ctr3d( tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:30),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,vsplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot gradient Richardson Number g / theta * d(theta)/dz
! -----------------------
! (dV/dz)^2g
!
!-----------------------------------------------------------------------
!
!
IF(ipriority(38) == ip ) THEN
IF(gricplt > 0) THEN
CALL cal_gric(tem9,u,v,zp,pt,tem7,tem8,nx,ny,nz)
END IF
IF(gricplt == 1 .OR. gricplt == 2 .OR. gricplt == 4 .OR. &
gricplt == 5 ) THEN
CALL get_mulovrlay('gricplt',7,ovrmul_num,ovrmulname, &
sovrlay)
label = 'Richardson Number'
CALL ctrsetup(gricinc,gricminc,gricmaxc, &
gricovr,grichlf,griczro,griccol1,griccol2,'gricplt ')
CALL ctr3d( tem9, xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:17),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,gricplt)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot absolute vorticity
!
!-----------------------------------------------------------------------
!
IF(ipriority(39) == ip ) THEN
IF(avorplt > 0) THEN
CALL cal_avor(tem9,u,v,x,y,nx,ny,nz,slicopt,flagsin,omega, &
sinlat,tem1,tem2,tem3)
END IF
IF( avorplt == 1 .OR. avorplt == 2 .OR. avorplt == 3 &
.OR. avorplt == 4 .OR. avorplt == 5 ) THEN
CALL get_mulovrlay('avorplt',7,ovrmul_num,ovrmulname, &
sovrlay)
label = 'Absolute Vort*10^5 (1/s)'
CALL ctrsetup(avorinc,avorminc,avormaxc, &
avorovr,avorhlf,avorzro,avorcol1,avorcol2,'avorplt ')
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:24),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3,tem4, &
tem5,tem6,hterain,avorplt)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot scalars for qtplot ,similar to qwplot, except that the
! vapor component is included.
!
!-----------------------------------------------------------------------
!
IF(ipriority(41) == ip ) THEN
IF(qtplot > 0) THEN
CALL cal_qt(tem9,qc,qr,qi,qs,qh,qv,nx,ny,nz)
END IF
IF( qtplot == 1 .OR. qtplot == 2 .OR. qtplot == 4 &
.OR. qtplot == 5 ) THEN
CALL get_mulovrlay('qtplot',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Total water & vapor (g/kg)'
CALL ctrsetup(qtinc,qtminc,qtmaxc, &
qtovr,qthlf,qtzro,qtcol1,qtcol2,'qtplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:26),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1000.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qtplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot RHI relative humidity with ice phase
!
!-----------------------------------------------------------------------
!
IF(ipriority(42) == ip) THEN
IF(rhiplot > 0) THEN
CALL cal_rhi(tem9,pt, pprt ,pbar,qv,tem1,tem2,nx,ny,nz)
END IF
IF(rhiplot == 1 .OR. rhiplot == 2 .OR. rhiplot == 4 .OR. &
rhiplot == 5 ) THEN
CALL get_mulovrlay('rhiplot',7,ovrmul_num,ovrmulname, &
sovrlay)
label = 'RHI'
CALL ctrsetup(rhiinc,rhiminc,rhimaxc, &
rhiovr,rhihlf,rhizro,rhicol1,rhicol2,'rhiplot ')
CALL ctr3d( tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:3),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,rhiplot)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot vector field
!
!-----------------------------------------------------------------------
!
CALL varplt( 'vtrplt ' )
IF(varname(1:6) == ovrname(1:6) .AND. sovrlay == 1) plotovr=.true.
IF( (ipriority(32) == ip .AND. vtrplt == 0) &
.AND. (.NOT.plotovr)) GOTO 1100
IF( ipriority(32) == ip .OR. plotovr) THEN
IF( vtrplt == 1 .OR. (ovrlaymulopt == 0 .OR. plotovr .OR. &
(ovrlaymulopt == 1 .AND. .NOT.plotovr)) )THEN
CALL vtrunt ( vtrunit )
IF(plotovr) CALL overlay ( 1 )
IF(.NOT.plotovr) CALL overlay( vtrovr )
CALL ctrcol(vtrcol1,vtrcol2)
CALL varplt( 'vtrplt ' )
CALL ctrvtr( vtrunits, vtrtype )
!
IF(nobs > 0 .AND. slicopt == 1 .AND. kslice == 2) THEN
CALL cal_vtrobs(dd,ff,drot, vtrunits)
END IF
CALL cal_vtr(tem7,tem8,tem9,u,v,w,nx,ny,nz,vtrunits, &
label,length)
CALL vtr3d( tem7,tem8,tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, nz,kbgn,kend,kst, &
kslice,jslice,islice, label(1:length),time, runname,1.0, &
slicopt,n,xp,yp,zs2,u1,v1,u2,v2,w2, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
END IF
IF(plotovr) THEN
plotovr=.false.
sovrlay = 0
END IF
END IF
1100 CONTINUE
!
!-----------------------------------------------------------------------
!
! Plot (u,v) in cross-section .
!
!-----------------------------------------------------------------------
!
CALL varplt( 'xuvplt ' )
IF(varname(1:6) == ovrname(1:6) .AND. sovrlay == 1) plotovr=.true.
IF( ipriority(40) == ip .OR. plotovr) THEN
IF( xuvplt == 1 .AND. (ovrlaymulopt == 0 .OR. plotovr .OR. &
(ovrlaymulopt == 1 .AND. .NOT.plotovr) ) )THEN
CALL vtrunt ( xuvunit )
IF(plotovr) CALL overlay ( 1 )
IF(.NOT.plotovr) CALL overlay( xuvovr )
CALL ctrcol(xuvcol1,xuvcol2)
CALL varplt( 'xuvplt ' )
CALL ctrvtr( xuvunits, xuvtype )
CALL cal_xuv(tem7,tem8,tem9,u,v,w,nx,ny,nz,xuvunits,label, &
length)
IF(slicopt == 2 .OR. slicopt == 3 .OR. slicopt == 5) &
CALL vtr3d( tem7,tem8,tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, nz,kbgn,kend,kst, &
kslice,jslice,islice, label(1:length),time, runname,1.0, &
slicopt,n,xp,yp,zs2,u1,v1,u2,v2,w2, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
END IF
IF(plotovr) THEN
plotovr=.false.
sovrlay = 0
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot Mongtgomery Streamfunction (10m) on constant theta surfaces
!
!-----------------------------------------------------------------------
!
IF(slicopt == 7) THEN
IF(ipriority(107) == ip ) THEN
IF(msfplt == 1.OR.msfplt == 2.OR.msfplt == 4.OR.msfplt == 5) &
THEN
CALL get_mulovrlay('msfplt',6,ovrmul_num,ovrmulname,sovrlay)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem9(i,j,k)=(cp*tz(i,j,k)+g*zpc(i,j,k)*1000.0)/g*0.1
END DO
END DO
END DO
CALL ctrsetup(msfinc,msfminc,msfmaxc,msfovr, &
msfhlf,msfzro,msfcol1,msfcol2,'msfplt ')
label = 'MSF (10m)'
CALL ctr3d(tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
label(1:9), &
time,slicopt, kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,msfplt)
END IF
END IF
IF(ipriority(108) == ip ) THEN
IF( ipvplt == 1 .OR. ipvplt == 2 .OR. ipvplt == 4 .OR. &
ipvplt == 5 ) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)=(u(i+1,j,k)+u(i,j,k))*0.5
tem2(i,j,k)=(v(i,j+1,k)+v(i,j,k))*0.5
END DO
END DO
END DO
DO k=2,nz-2
DO j=1,ny-1
DO i=1,nx-1
tem3(i,j,k)=-g*(zps3d(i,j,k+1)-zps3d(i,j,k-1))/ &
(pbar(i,j,k+1)+pprt(i,j,k+1)-pbar(i,j,k-1)-pprt(i,j,k-1))
END DO
END DO
END DO
CALL hintrp1(nx,ny,nz,2,nz-2,tem1,zps3d,zlevel,tem9(1,1,1))
CALL hintrp1(nx,ny,nz,2,nz-2,tem2,zps3d,zlevel,tem9(1,1,2))
CALL hintrp1(nx,ny,nz,2,nz-2,tem3,zps3d,zlevel,tem9(1,1,3))
DO j=2,ny-2
DO i=2,nx-2
IF(ABS(tem9(i+1,j,2)+9999.0) < 0.1.OR. &
ABS(tem9(i-1,j,2)+9999.0) < 0.1.OR. &
ABS(tem9(i,j+1,1)+9999.0) < 0.1.OR. &
ABS(tem9(i,j-1,1)+9999.0) < 0.1.OR. &
ABS(tem9(i,j,3 )+9999.0) < 0.1) THEN
tem9(i,j,4)= -9999.0
relvort = -9999.0
ELSE
relvort = (tem9(i+1,j,2)-tem9(i-1,j,2))*dxinv &
-(tem9(i,j+1,1)-tem9(i,j-1,1))*dyinv
tem9(i,j,4)=(relvort &
+fcorio(i,j))*tem9(i,j,3)*1.0E6
END IF
END DO
END DO
DO j=2,ny-2
tem9( 1,j,4)=tem9( 2,j,4)
tem9(nx-1,j,4)=tem9(nx-2,j,4)
END DO
DO j=1,nx-1
tem9(i, 1,4)=tem9(i, 2,4)
tem9(i,ny-1,4)=tem9(i,ny-2,4)
END DO
CALL get_mulovrlay('ipvplt',6,ovrmul_num,ovrmulname,sovrlay)
WRITE(label,'(''IPV (PVU) AT THETA='',F5.1)') zlevel
CALL ctrsetup(ipvinc,ipvminc,ipvmaxc, &
ipvovr,ipvhlf,ipvzro,ipvcol1,ipvcol2,'ipvplt ')
CALL ctrsfc(tem9(1,1,4),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:24),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,ipvplt)
END IF
END IF
END IF ! For slicopt=7 only
END IF ! For slicopt=1 - 7
!
!-----------------------------------------------------------------------
!
! Start plotting of 2-D surface fields...
!
!-----------------------------------------------------------------------
!
IF (slicopt == 8 .OR. slicopt == 1) THEN
IF (slicopt == 8 ) THEN ! 2-D surface field
aspect=xr/yr
IF(yxstrch /= 0.0) THEN
aspratio=yxstrch
ELSE
aspratio= aspect
END IF
END IF
IF(ipriority(51) == ip) THEN
IF( (trnplt == 1 .OR. trnplt == 2 .OR. trnplt == 4 .OR. &
trnplt == 5 ) .AND. sigplt(51) == 0 ) THEN
CALL get_mulovrlay('trnplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Terrain height (m)'
time0 = 0.0
CALL ctrsetup(trninc,trnminc,trnmaxc, &
trnovr,trnhlf,trnzro,trncol1,trncol2,'trnplt ')
CALL ctrsfc(zp(1,1,2),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:18),time0, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,trnplt)
sigplt(51) = 1
END IF
END IF
IF(ipriority(56) == ip ) THEN
IF( (wetcanplt == 1 .OR. wetcanplt == 2 .OR. wetcanplt == 4 &
.OR. wetcanplt == 5 ) .AND. sigplt(56) == 0 ) THEN
CALL get_mulovrlay('wetcanplt',9,ovrmul_num,ovrmulname, &
sovrlay)
label = 'Canopy water amount'
CALL ctrsetup(wcpinc,wcpminc,wcpmaxc, &
wcpovr,wcphlf,wcpzro,wcpcol1,wcpcol2,'wetcanplt ')
CALL ctrsfc(wetcanp,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,wetcanplt)
sigplt(56) = 1
END IF
END IF
IF(ipriority(57) == ip) THEN
IF( (raincplt == 1 .OR.raincplt == 2 .OR.raincplt == 4 .OR. &
raincplt == 5 ) .AND. sigplt(57) == 0 ) THEN
CALL get_mulovrlay('raincplt',8,ovrmul_num,ovrmulname,sovrlay)
IF(racunit == 0) THEN
label = 'Cumulus Rainfall (mm)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = rainc(i,j)
END DO
END DO
ELSE IF( racunit == 1) THEN
label = 'Cumulus Rainfall (in)' ! unit is inch
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = rainc(i,j)*0.039370079 ! (1/25.4)
END DO
END DO
END IF
CALL ctrsetup(raincinc,raincminc,raincmaxc, &
racovr,rachlf,raczro,raccol1,raccol2,'raincplt ')
CALL ctrsfc(tem9(1,1,1),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, nx,ibgn,iend, ny,jbgn,jend, &
label(1:28),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,raincplt)
sigplt(57) = 1
END IF
END IF
IF(ipriority(58) == ip) THEN
IF( (raingplt == 1 .OR. raingplt == 2 .OR. raingplt == 4 .OR. &
raingplt == 5 ) .AND. sigplt(58) == 0 ) THEN
CALL get_mulovrlay('raingplt',8,ovrmul_num,ovrmulname,sovrlay)
IF(ragunit == 0 ) THEN
label = 'Gridscale Rainfall (mm) '
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = raing(i,j)
END DO
END DO
ELSE IF( ragunit == 1) THEN
label = 'Gridscale Rainfall (in) '
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = raing(i,j)*0.039370079 ! (1/25.4)
END DO
END DO
END IF
CALL ctrsetup(rainginc,raingminc,raingmaxc, &
ragovr,raghlf,ragzro,ragcol1,ragcol2,'raingplt ')
CALL ctrsfc(tem9(1,1,1),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, nx,ibgn,iend, ny,jbgn,jend, &
label(1:30),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,raingplt)
sigplt(58) = 1
END IF
END IF
IF(ipriority(59) == ip) THEN
IF( (raintplt == 1 .OR.raintplt == 2 .OR.raintplt == 4 &
.OR. raintplt == 5 ) &
.AND. sigplt(59) == 0 ) THEN
CALL get_mulovrlay('raintplt',8,ovrmul_num,ovrmulname,sovrlay)
IF(ratunit == 0 ) THEN
label = 'Total Rainfall (mm)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = rainc(i,j) + raing(i,j)
END DO
END DO
ELSE IF (ratunit == 1) THEN
label = 'Total Rainfall (in)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = (rainc(i,j)+raing(i,j))*0.039370079 ! (1/25.4)
END DO
END DO
END IF
CALL ctrsetup(raintinc,raintminc,raintmaxc, &
ratovr,rathlf,ratzro,ratcol1,ratcol2,'raintplt ')
CALL ctrsfc(tem9,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:22),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,raintplt)
sigplt(59) = 1
END IF
END IF
IF(ipriority(60) == ip) THEN
IF( (pslplt == 1 .OR. pslplt == 2 .OR. pslplt == 4 .OR. &
pslplt == 5 ) &
.AND. sigplt(60) == 0 ) THEN
CALL get_mulovrlay('pslplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Sea Level Pressure (mb)'
IF(ovrobs == 1 .AND. nobs > 0) THEN
DO iob=1,nobs
!
! Sea-level pressure
!
IF(pmsl(iob) > 0. ) THEN
obs1(iob)=mod(nint(pmsl(iob)),100)
print *, ' pmsl, obs1: ',pmsl(iob),obs1(iob)
! ELSE IF (pstn(iob) > 0. .AND. tob(iob) > -98. ) THEN
! tk=(5.*(tob(iob)-32.)/9.) + 273.15
! obs1(iob)=pstn(iob)*((tk+gamma*elevob(iob))/tk)**ex2
! obs1(iob)=mod(nint(obs1(iob)),100)
! print *, ' pstn, obs1: ',pstn(iob),obs1(iob)
! ELSE IF( alt(iob) > 0. ) THEN
! IF ( tob(iob) > -98.) THEN
! tk=(5.*(tob(iob)-32.)/9.) + 273.15
! obs1(iob)=alttostpr(alt(iob),elevob(iob))* &
! ((tk+gamma*elevob(iob))/tk)**ex2
! obs1(iob)=mod(nint(obs1(iob)),100)
! print *, ' alt1, obs1: ',alt(iob),obs1(iob)
! ELSE
! obs1(iob)=0.01*alt(iob)
! obs1(iob)=mod(nint(obs1(iob)),100)
! print *, ' alt2, obs1: ',alt(iob),obs1(iob)
! END IF
ELSE
obs1(iob)=-99.
END IF
END DO
obsset=1
END IF
CALL ctrsetup(pslinc,pslminc,pslmaxc, &
pslovr,pslhlf,pslzro,pslcol1,pslcol2,'pslplt ')
CALL ctrsfc(psl,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:23),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,pslplt)
sigplt(60) = 1
END IF
END IF
IF(ipriority(61) == ip) THEN
IF( (liplt == 1 .OR. liplt == 2 .OR. liplt == 4 .OR. liplt == 5 ) &
.AND. sigplt(61) == 0 ) THEN
CALL get_mulovrlay('liplt',5,ovrmul_num,ovrmulname,sovrlay)
label = 'Sfc-based LI (K)'
CALL ctrsetup(liinc,liminc,limaxc, &
liovr,lihlf,lizro,licol1,licol2,'liplt ')
CALL ctrsfc(li,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,liplt)
sigplt(61) = 1
END IF
END IF
IF(ipriority(62) == ip) THEN
IF( (capeplt == 1 .OR. capeplt == 2 .OR. capeplt == 4 .OR. &
capeplt == 5) &
.AND. sigplt(62) == 0 ) THEN
CALL get_mulovrlay('capeplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Surface CAPE (J/kg)'
CALL ctrsetup(capeinc,capeminc,capemaxc, &
capovr,caphlf,capzro,capcol1,capcol2,'capeplt ')
CALL ctrsfc(cape,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,capeplt)
sigplt(62) = 1
END IF
END IF
IF(ipriority(63) == ip) THEN
IF( (cinplt == 1 .OR. cinplt == 2 .OR. cinplt == 4 .OR. &
cinplt == 5 ) &
.AND. sigplt(63) == 0 )THEN
CALL get_mulovrlay('cinplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Surface CIN (J/kg)'
CALL ctrsetup(cininc,cinminc,cinmaxc, &
cinovr,cinhlf,cinzro,cincol1,cincol2,'cinplt ')
CALL ctrsfc(cin,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm,nx,ibgn,iend, ny,jbgn,jend, &
label(1:18),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,cinplt)
sigplt(63) = 1
END IF
END IF
IF(ipriority(64) == ip ) THEN
IF( (thetplt == 1 .OR. thetplt == 2 .OR. thetplt == 4 .OR. &
thetplt == 5 ) &
.AND. sigplt(64) == 0 ) THEN
CALL get_mulovrlay('thetplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Surface Theta-e (K)'
CALL ctrsetup(thetinc,thetminc,thetmaxc, &
theovr,thehlf,thezro,thecol1,thecol2,'thetplt ')
CALL ctrsfc(thet,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm,nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,thetplt)
sigplt(64) = 1
END IF
END IF
IF(ipriority(65) == ip ) THEN
IF( (heliplt == 1 .OR. heliplt == 2 .OR. heliplt == 4 .OR. &
heliplt == 5 ) &
.AND. sigplt(65) == 0 ) THEN
CALL get_mulovrlay('heliplt',7,ovrmul_num,ovrmulname,sovrlay)
label = '0-3km Helicity (m^2/s^2)'
CALL ctrsetup(heliinc,heliminc,helimaxc, &
helovr,helhlf,helzro,helcol1,helcol2,'heliplt ')
CALL ctrsfc(heli,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm,nx,ibgn,iend, ny,jbgn,jend, &
label(1:24),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,heliplt)
sigplt(65) = 1
END IF
END IF
IF(ipriority(66) == ip ) THEN
IF( (ctcplt == 1 .OR. ctcplt == 2 .OR. ctcplt == 4 .OR. &
ctcplt == 5 ) .AND. sigplt(66) == 0 ) THEN
CALL get_mulovrlay('ctcplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Convective temp (C)'
CALL ctrsetup(ctcinc,ctcminc,ctcmaxc, &
ctcovr,ctchlf,ctczro,ctccol1,ctccol2,'ctcplt ')
CALL ctrsfc(ct,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,ctcplt)
sigplt(66) = 1
END IF
END IF
IF(ipriority(67) == ip ) THEN
IF( (brnplt == 1 .OR. brnplt == 2 .OR. brnplt == 4 .OR. &
brnplt == 5 ) &
.AND. sigplt(67) == 0 ) THEN
CALL get_mulovrlay('brnplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'BRN Shear Denom (m2/s2)'
CALL ctrsetup(brninc,brnminc,brnmaxc, &
brnovr,brnhlf,brnzro,brncol1,brncol2,'brnplt ')
CALL ctrsfc(brn,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:22),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,brnplt)
sigplt(67) = 1
END IF
END IF
IF(ipriority(68) == ip ) THEN
IF( (brnuplt == 1 .OR. brnuplt == 2 .OR. brnuplt == 4 .OR. &
brnuplt == 5) &
.AND. sigplt(68) == 0 ) THEN
CALL get_mulovrlay('brnuplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Bulk Richardson Shear (1/s)'
CALL ctrsetup(brnuinc,bruminc,brumaxc, &
brnuovr,brnuhlf,brnuzro,brnucol1,brnucol2,'brnuplt ')
CALL ctrsfc(brnu,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:27),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,brnuplt)
sigplt(68) = 1
END IF
END IF
IF(ipriority(69) == ip ) THEN
IF( (srlfplt == 1 .OR. srlfplt == 2 .OR. srlfplt == 4 .OR. &
srlfplt == 5 ) .AND. sigplt(69) == 0 ) THEN
CALL get_mulovrlay('srlfplt',7,ovrmul_num,ovrmulname,sovrlay)
label = '0-2km StRel Flow (m/s)'
CALL ctrsetup(srlfinc,srlminc,srlmaxc, &
srlfovr,srlfhlf,srlfzro,srlfcol1,srlfcol2,'srlfplt ')
CALL ctrsfc(srlfl,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:22),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,srlfplt)
sigplt(69) = 1
END IF
END IF
IF(ipriority(70) == ip ) THEN
IF( (srmfplt == 1 .OR. srmfplt == 2 .OR. srmfplt == 4 .OR. &
srmfplt == 5) .AND. sigplt(70) == 0 ) THEN
CALL get_mulovrlay('srmfplt',7,ovrmul_num,ovrmulname,sovrlay)
label = '2-9km StRel Flow (m/s)'
CALL ctrsetup(srmfinc,srmminc,srmmaxc, &
srmfovr,srmfhlf,srmfzro,srmfcol1,srmfcol2,'srmfplt ')
CALL ctrsfc(srmfl,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:22),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,srmfplt)
sigplt(70) = 1
END IF
END IF
IF(ipriority(71) == ip ) THEN
IF( (capsplt == 1 .OR. capsplt == 2 .OR. capsplt == 4 .OR. &
capsplt == 5 ) .AND. sigplt(71) == 0 ) THEN
CALL get_mulovrlay('capsplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Cap Strength (K)'
CALL ctrsetup(capsinc,capsminc,capsmaxc, &
capsovr,capshlf,capszro,capscol1,capscol2,'capsplt ')
CALL ctrsfc(capst,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:16),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,capsplt)
sigplt(71) = 1
END IF
END IF
IF(ipriority(72) == ip ) THEN
IF( (blcoplt == 1 .OR. blcoplt == 2 .OR. blcoplt == 4 .OR. &
blcoplt == 5 ) .AND. sigplt(72) == 0 ) THEN
CALL get_mulovrlay('blcoplt',7,ovrmul_num,ovrmulname,sovrlay)
! label = '0-2km Layer Conv.*1000 (1/s)'
label = '0-2km Conv.*1000 (1/s)'
CALL ctrsetup(blcoinc,blcominc,blcomaxc, &
blcoovr,blcohlf,blcozro,blcocol1,blcocol2,'blcoplt ')
CALL ctrsfc(blcon,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:22),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,blcoplt)
sigplt(72) = 1
END IF
END IF
IF(ipriority(73) == ip ) THEN
IF( (viqcplt == 1 .OR. viqcplt == 2 .OR. viqcplt == 4 &
.OR. viqcplt == 5 ) &
.AND. sigplt(73) == 0 ) THEN
CALL get_mulovrlay('viqcplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integrated qc (kg/m2)'
CALL ctrsetup(viqcinc,viqcminc,viqcmaxc, &
viqcovr,viqchlf,viqczro,viqccol1,viqccol2,'viqcplt ')
CALL cal_viqc(tem7, qc, rhobar, zp, nx,ny,nz)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:27),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,viqcplt)
sigplt(73) = 1
END IF
END IF
IF(ipriority(74) == ip ) THEN
IF( (viqrplt == 1 .OR. viqrplt == 2 .OR. viqrplt == 4 &
.OR. viqrplt == 5 ) &
.AND. sigplt(74) == 0 ) THEN
CALL get_mulovrlay('viqrplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integrated qr (kg/m2)'
CALL ctrsetup(viqrinc,viqrminc,viqrmaxc, &
viqrovr,viqrhlf,viqrzro,viqrcol1,viqrcol2,'viqrplt ')
CALL cal_viqr(tem7, qr, rhobar, zp, nx,ny,nz)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:27),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,viqrplt)
sigplt(74) = 1
END IF
END IF
IF(ipriority(75) == ip ) THEN
IF( (viqiplt == 1 .OR. viqiplt == 2 .OR. viqiplt == 4 &
.OR. viqiplt == 5 ) &
.AND. sigplt(75) == 0 ) THEN
CALL get_mulovrlay('viqiplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integrated qi (kg/m2)'
CALL ctrsetup(viqiinc,viqiminc,viqimaxc, &
viqiovr,viqihlf,viqizro,viqicol1,viqicol2,'viqiplt ')
CALL cal_viqi(tem7, qi, rhobar, zp, nx,ny,nz)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:27),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,viqiplt)
sigplt(75) = 1
END IF
END IF
IF(ipriority(76) == ip ) THEN
IF( (viqsplt == 1 .OR. viqsplt == 2 .OR. viqsplt == 4 &
.OR. viqsplt == 5 ) &
.AND. sigplt(76) == 0 ) THEN
CALL get_mulovrlay('viqsplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integrated qs (kg/m2)'
CALL ctrsetup(viqsinc,viqsminc,viqsmaxc, &
viqsovr,viqshlf,viqszro,viqscol1,viqscol2,'viqsplt ')
CALL cal_viqs(tem7, qs, rhobar, zp, nx,ny,nz)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:27),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,viqsplt)
sigplt(76) = 1
END IF
END IF
IF(ipriority(77) == ip ) THEN
IF( (viqhplt == 1 .OR. viqhplt == 2 .OR. viqhplt == 4 &
.OR. viqhplt == 5 ) &
.AND. sigplt(77) == 0 ) THEN
CALL get_mulovrlay('viqhplt',7,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integrated qh (kg/m2)'
CALL ctrsetup(viqhinc,viqhminc,viqhmaxc, &
viqhovr,viqhhlf,viqhzro,viqhcol1,viqhcol2,'viqhplt ')
CALL cal_viqh(tem7, qh, rhobar, zp, nx,ny,nz)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:27),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,viqhplt)
sigplt(77) = 1
END IF
END IF
IF(ipriority(78) == ip ) THEN
IF( (vilplt == 1 .OR. vilplt == 2 .OR. vilplt == 4 &
.OR. vilplt == 5 ) &
.AND. sigplt(78) == 0 ) THEN
CALL get_mulovrlay('vilplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integ Liquid (kg/m2)'
CALL ctrsetup(vilinc,vilminc,vilmaxc, &
vilovr,vilhlf,vilzro,vilcol1,vilcol2,'vilplt ')
CALL cal_vil(tem7,qc,qr,rhobar,zp, nx,ny,nz,tem6)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:26),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,vilplt)
sigplt(78) = 1
END IF
END IF
IF(ipriority(79) == ip ) THEN
IF( (viiplt == 1 .OR. viiplt == 2 .OR. viiplt == 4 &
.OR. viiplt == 5 ) &
.AND. sigplt(79) == 0 ) THEN
CALL get_mulovrlay('viiplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integrated ice (kg/m2)'
CALL ctrsetup(viiinc,viiminc,viimaxc, &
viiovr,viihlf,viizro,viicol1,viicol2,'viiplt ')
CALL cal_vii(tem7,qi,qs,qh,rhobar,zp, nx,ny,nz,tem6)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:28),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,viiplt)
sigplt(79) = 1
END IF
END IF
IF(ipriority(80) == ip ) THEN
IF( (vicplt == 1 .OR. vicplt == 2 .OR. vicplt == 4 &
.OR. vicplt == 5 ) &
.AND. sigplt(80) == 0 ) THEN
CALL get_mulovrlay('vicplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integ Condensate (kg/m2)'
CALL ctrsetup(vicinc,vicminc,vicmaxc, &
vicovr,vichlf,viczro,viccol1,viccol2,'vicplt ')
CALL cal_vic(tem7,qc,qr,qi,qs,qh,rhobar,zp,nx,ny,nz,tem6)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:30),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,vicplt)
sigplt(80) = 1
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot storm motion vector
!
!-----------------------------------------------------------------------
!
IF(ipriority(81) == ip ) THEN
IF( strmplt == 1 .AND. sigplt(81) == 0 ) THEN
CALL get_mulovrlay('strmplt',7,ovrmul_num,ovrmulname,sovrlay)
CALL vtrunt ( strmunit )
CALL overlay( strmovr )
CALL ctrcol(strmcol1,strmcol2)
CALL varplt( 'strmplt ' )
CALL ctrvtr( strmunits, strmtype )
!
CALL cal_strm(tem7, tem8,ustrm,vstrm,strmunits,nx,ny,nz, &
label,length)
CALL vtrsfc( tem7(1,1,1),tem8(1,1,1),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, &
label(1:length),time, runname,1.0,slicopt, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
sigplt(81) = 1
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot vertically integrated total water (kg/m**2)
!
!-----------------------------------------------------------------------
!
IF(ipriority(82) == ip ) THEN
IF( (vitplt == 1 .OR. vitplt == 2 .OR. vitplt == 4 &
.OR. vitplt == 5 ) &
.AND. sigplt(82) == 0 ) THEN
CALL get_mulovrlay('vitplt',6,ovrmul_num,ovrmulname,sovrlay)
label = 'Vert. Integ Total Water (kg/m2)'
CALL ctrsetup(vitinc,vitminc,vitmaxc, &
vitovr,vithlf,vitzro,vitcol1,vitcol2,'vitplt ')
CALL cal_vit(tem7,qv,qc,qr,qi,qs,qh,rhobar,zp,nx,ny,nz,tem6)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:31),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,vitplt)
sigplt(82) = 1
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot precipitable water vapor, cm
!
!-----------------------------------------------------------------------
!
IF(ipriority(83) == ip ) THEN
IF( (pwplt == 1 .OR. pwplt == 2 .OR. pwplt == 4 .OR. pwplt == 5 ) &
.AND. sigplt(83) == 0 ) THEN
CALL get_mulovrlay('pwplt',5,ovrmul_num,ovrmulname,sovrlay)
label = 'Precipitable Water Vapor(cm)'
CALL ctrsetup(pwinc,pwminc,pwmaxc, &
pwovr,pwhlf,pwzro,pwcol1,pwcol2,'pwplt ')
CALL cal_pw(tem7,qv,rhobar,zp,nx,ny,nz,tem6)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:28),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,pwplt)
sigplt(83) = 1
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot total precipitation rate
!
!-----------------------------------------------------------------------
!
IF(ipriority(84) == ip ) THEN
IF( (tprplt == 1 .OR. tprplt == 2 .OR. tprplt == 4 &
.OR. tprplt == 5 ) &
.AND. sigplt(84) == 0 ) THEN
CALL get_mulovrlay('tprplt',6,ovrmul_num,ovrmulname,sovrlay)
! label = 'total Precipitation rate( mm/h)'
CALL ctrsetup(tprinc,tprminc,tprmaxc, &
tprovr,tprhlf,tprzro,tprcol1,tprcol2,'tprplt ')
CALL cal_tpr(tem7,prcrate(1,1,1),nx,ny,nz,tprunits, &
label,length)
! call xcltyp(0)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:length),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,tprplt)
sigplt(84) = 1
! call xcltyp(2)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot grid scale precip. rate
!
!-----------------------------------------------------------------------
!
IF(ipriority(85) == ip ) THEN
IF( (gprplt == 1 .OR. gprplt == 2 .OR. gprplt == 4 &
.OR. gprplt == 5 ) &
.AND. sigplt(85) == 0 ) THEN
CALL get_mulovrlay('gprplt',6,ovrmul_num,ovrmulname,sovrlay)
! label = 'grid scale precip. rate( mm/h)'
CALL ctrsetup(gprinc,gprminc,gprmaxc, &
gprovr,gprhlf,gprzro,gprcol1,gprcol2,'gprplt ')
CALL cal_gpr(tem7,prcrate(1,1,2),prcrate(1,1,4), nx,ny,nz, &
gprunits,label,length)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:length),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,gprplt)
sigplt(85) = 1
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot Convective precip. rate = prcrate(1,1,3)
!
!-----------------------------------------------------------------------
!
IF(ipriority(86) == ip ) THEN
IF( (cprplt == 1 .OR. cprplt == 2 .OR. cprplt == 4 &
.OR. cprplt == 5 ) &
.AND. sigplt(86) == 0 ) THEN
CALL get_mulovrlay('cprplt',6,ovrmul_num,ovrmulname,sovrlay)
! label = 'Convective precip. rate( mm/h)'
CALL ctrsetup(cprinc,cprminc,cprmaxc, &
cprovr,cprhlf,cprzro,cprcol1,cprcol2,'cprplt ')
CALL cal_cpr(tem7,prcrate(1,1,3),nx,ny,nz,cprunits, &
label,length)
CALL ctrsfc(tem7,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:length),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,cprplt)
sigplt(86) = 1
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Plot surface characteristics.
!
!-----------------------------------------------------------------------
!
IF(ipriority(101) == ip ) THEN
IF( (soiltpplt == 1 .OR. soiltpplt == 2 .OR. soiltpplt == 4 &
.OR. soiltpplt == 5 ) .AND. sigplt(101) == 0 ) THEN
WRITE(label,'(''SOIL TYPE '',I1 )') soiltpn
! label = 'Soil type'
time0=0.0
DO j=1,ny
DO i=1,nx
tem9(i,j,1)=soiltyp(i,j,soiltpn)
END DO
END DO
CALL ctrsetup(soiltpinc,soiltpminc,soiltpmaxc, &
styovr,styhlf,styzro,stycol1,stycol2,'soiltpplt ')
CALL ctrsfc(tem9,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:11),time0, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,soiltpplt)
sigplt(101) = 1
END IF
END IF
IF(ipriority(102) == ip ) THEN
IF( (vegtpplt == 1 .OR. vegtpplt == 2 .OR. vegtpplt == 4 &
.OR. vegtpplt == 5 ) .AND. sigplt(102) == 0) THEN
label = 'Vegetation type'
time0=0.0
DO j=1,ny
DO i=1,nx
tem9(i,j,1)=vegtyp(i,j)
END DO
END DO
CALL ctrsetup(vegtpinc,vegtpminc,vegtpmaxc, &
vtyovr,vtyhlf,vtyzro,vtycol1,vtycol2,'vegtpplt ')
CALL ctrsfc(tem9,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:15),time0, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,vegtpplt)
sigplt(102)=1
END IF
END IF
IF(ipriority(103) == ip ) THEN
IF( (laiplt == 1 .OR. laiplt == 2 .OR. laiplt == 4 .OR. &
laiplt == 5 ) .AND. sigplt(103) == 0) THEN
label = 'Leaf Area Index'
time0=0.0
CALL ctrsetup(laiinc,laiminc,laimaxc, &
laiovr,laihlf,laizro,laicol1,laicol2,'laiplt ')
CALL ctrsfc(lai,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:15),time0, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,laiplt)
sigplt(103) =1
END IF
END IF
IF(ipriority(104) == ip) THEN
IF( (rouplt == 1 .OR. rouplt == 2 .OR. rouplt == 4 .OR. &
rouplt == 5) .AND. sigplt(104) == 0) THEN
label = 'Surface roughness'
time0=0.0
CALL ctrsetup(rouinc,rouminc,roumaxc, &
rouovr,rouhlf,rouzro,roucol1,roucol2,'rouplt ')
CALL ctrsfc(roufns,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:17),time0, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,rouplt)
sigplt(104) = 1
END IF
END IF
IF(ipriority(105) == ip) THEN
IF( (vegplt == 1 .OR. vegplt == 2 .OR. vegplt == 4 .OR. &
vegplt == 5) .AND. sigplt(105) == 0) THEN
label = 'Vegetation fraction'
time0=0.0
CALL ctrsetup(veginc,vegminc,vegmaxc, &
vegovr,veghlf,vegzro,vegcol1,vegcol2,'vegplt ')
CALL ctrsfc(veg,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time0, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,vegplt)
sigplt(105) = 1
END IF
END IF
IF(ipriority(106) == ip) THEN
IF( (snowdplt == 1 .OR. snowdplt == 2 .OR. snowdplt == 4 .OR. &
snowdplt == 5) .AND. sigplt(106) == 0) THEN
label = 'Snow depth (m)'
CALL ctrsetup(snowdinc,snowdminc,snowdmaxc, &
snowdovr,snowdhlf,snowdzro,snowdcol1, &
snowdcol2,'snowdplt ')
DO j=1,ny
DO i=1,nx
tem9(i,j,1)=snowdpth(i,j)
END DO
END DO
CALL ctrsfc(tem9,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
label(1:19),time, runname, 1.0,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,snowdplt)
sigplt(106) = 1
END IF
END IF
!-----------------------------------------------------------------------
!
! Overlay wind vector/barb if necessary on 2-D surface fields.
!
!-----------------------------------------------------------------------
!
CALL varplt( 'vtrplt ' )
IF(varname(1:6) == ovrname(1:6) .AND. sovrlay == 1) plotovr=.true.
IF( vtrplt == 1 .AND. plotovr ) THEN
CALL vtrunt ( vtrunit )
IF(plotovr) CALL overlay ( 1 )
IF(.NOT.plotovr) CALL overlay( vtrovr )
CALL ctrcol(vtrcol1,vtrcol2)
CALL varplt( 'vtrplt ' )
CALL ctrvtr( vtrunits, vtrtype )
CALL cal_vtr(tem7,tem8,tem9,u,v,w,nx,ny,nz,vtrunits, &
label,length)
CALL vtrsfc( tem7(1,1,1),tem8(1,1,1),xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, nx,ibgn,iend,ist, &
ny,jbgn,jend,jst, &
label(1:length),time, runname,1.0,slicopt, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
IF(plotovr) THEN
plotovr=.false.
sovrlay = 0
END IF
END IF
END IF ! For slicopt=8 - two plots
!
!-----------------------------------------------------------------------
!
! plot arbitrary variables.
!
!-----------------------------------------------------------------------
!
!NEW
IF(arbvaropt == 1) THEN
IF(slicopt <= 7 .AND. var3dnum > 0) THEN
DO i = 1,var3dnum
IF(ipriority(150+i) == ip ) THEN
IF(var3dplot(i) /= 0) THEN
CALL ctrsetup(var3dinc(i),var3dminc(i),var3dmaxc(i), &
var3dovr(i),var3dhlf(i),var3dzro(i), &
var3dcol1(i),var3dcol2(i),var3d(i))
CALL readvar1(nx,ny,nz, var3dv, var3d(i),var_name, &
varunits, time, runname, dirname3d(i), istatus)
! print*, var_name, varunits
CALL ctr3d(var3dv,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgn,iend, ny,jbgn,jend, nz,kbgn,kend, &
var3d(i)(1:6),time,slicopt,kslice,jslice,islice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,var3dplot(i))
END IF
END IF
END DO
END IF
IF(slicopt == 8 .AND. var2dnum > 0) THEN
DO i = 1,var2dnum
IF(ipriority(170+i) == ip ) THEN
IF(var2dplot(i) /= 0 .AND. sigplt(170+i) == 0) THEN
CALL ctrsetup(var2dinc(i),var2dminc(i),var2dmaxc(i), &
var2dovr(i),var2dhlf(i),var2dzro(i), &
var2dcol1(i),var2dcol2(i),var2d(i))
CALL readvar1(nx,ny,1, var2dv, var2d(i),var_name, &
varunits, time, runname, dirname2d(i), istatus)
! print*, var_name, varunits
CALL ctrsfc(var2dv,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgn,iend, ny,jbgn,jend, &
var2d(i)(1:6),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,var2dplot(i))
sigplt(170+i) = 1
END IF
END IF
END DO
END IF
END IF
!
! 06/06/2002 Zuwen He
!
! plotting soil variables, tsoil and qsoil
!
IF (iplot(slicopt) /= 0 .AND. &
(slicopt == 9 .OR. slicopt == 10 .OR. slicopt == 11)) THEN
IF(ipriority(43) == ip ) THEN
IF(tsoilplt == 1 .OR. tsoilplt == 2 .OR. tsoilplt == 4 .OR. &
tsoilplt == 5) THEN
CALL get_mulovrlay('tsoilplt',8,ovrmul_num,ovrmulname,sovrlay)
label = 'Soil temperature (K)'
CALL ctrsetup(tsoilinc,tsoilminc,tsoilmaxc, &
tsoilovr,tsoilhlf,tsoilzro,tsoilcol1,tsoilcol2, &
'tsoilplt ')
CALL ctr3d(tsoil(1,1,1,0),xc,yc,zpsoils3d, &
x1,x2,dxkm,y1,y2,dykm,zsoil1,zsoil2,dzsoilcm, &
nx,ibgn,iend,ny,jbgn,jend,nzsoil,ksoilbgn,ksoilend, &
label(1:25), &
time,slicopt,ksoilslice,jsoilslice,isoilslice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,tsoilplt)
sigplt(43) = 1
END IF
END IF
IF(ipriority(44) == ip ) THEN
IF(qsoilplt == 1 .OR. qsoilplt == 2 .OR. qsoilplt == 4 .OR. &
qsoilplt == 5) THEN
CALL get_mulovrlay('qsoilplt',8,ovrmul_num,ovrmulname,sovrlay)
label = 'Soil moisture (m**3/m**3)'
CALL ctrsetup(qsoilinc,qsoilminc,qsoilmaxc, &
qsoilovr,qsoilhlf,qsoilzro,qsoilcol1,qsoilcol2, &
'qsoilplt ')
CALL ctr3d(qsoil(1,1,1,0),xc,yc,zpsoils3d, &
x1,x2,dxkm,y1,y2,dykm,zsoil1,zsoil2,dzsoilcm, &
nx,ibgn,iend, ny,jbgn,jend,nzsoil,ksoilbgn,ksoilend, &
label(1:25), &
time,slicopt,ksoilslice,jsoilslice,isoilslice, &
n,xp,yp,b1,b2,zs2, &
runname,1.0,tem1,tem2,tem3, &
tem4,tem5,tem6,hterain,qsoilplt)
sigplt(44) = 1
END IF
END IF
END IF
!
!-----------------------------------------------------------------------
!
! End of priority loop
!
!-----------------------------------------------------------------------
!
END DO
!
!-----------------------------------------------------------------------
!
! End of loop over slices
!
!-----------------------------------------------------------------------
!
END DO
!
!-----------------------------------------------------------------------
!
! End of loop over slicopt (i.e. x-y, x-z or y-z slice plotting).
!
!-----------------------------------------------------------------------
!
END DO
!
!-----------------------------------------------------------------------
!
! 3-D wireframe isosurface plots.
!
!-----------------------------------------------------------------------
!
IF( w3dplt == 1) THEN
!
! Plot w=-wisosf isosurface.
!
CALL wirfrm(w,nx,ibgn,iend,ny,jbgn,jend,nz,kbgn,kend, &
wisosf, tem2)
CALL xqpspc(x1,x2,y1,y2)
CALL xchmag(0.02*(y2-y1)*lblmag)
CALL xqmap (x1,x2,y1,y2)
CALL xcharl(x1,y1-0.02*(y2-y1),runname)
WRITE(label,'(a,f4.1,a,f6.1,a)') &
'w=',wisosf,' m/s isosurface at t=',time/60.0,' min'
CALL xcharl(x1,y1-0.04*(y2-y1),label )
!
! Plot w=-wisosf isosurface.
!
CALL wirfrm(w,nx,ibgn,iend,ny,jbgn,jend,nz,kbgn,kend, &
-wisosf, tem2)
CALL xqpspc(x1,x2,y1,y2)
CALL xchmag(0.02*(y2-y1)*lblmag)
CALL xqmap (x1,x2,y1,y2)
CALL xcharl(x1,y1-0.02*(y2-y1),runname)
WRITE(label,'(a,f4.1,a,f6.1,a)') &
'w=',-wisosf,' m/s isosurface at t=',time/60.0,' min'
CALL xcharl(x1,y1-0.04*(y2-y1),label )
END IF
IF( q3dplt == 1) THEN
DO k=1,nz
DO j=1,ny
DO i=1,nx
tem1(i,j,k)=(qr(i,j,k)+qc(i,j,k))*1000.0
END DO
END DO
END DO
CALL wirfrm(tem1, &
nx,ibgn,iend,ny,jbgn,jend,nz,kbgn,kend,qisosf,tem2)
CALL xqpspc(x1,x2,y1,y2)
CALL xchmag(0.02*(y2-y1)*lblmag)
CALL xqmap (x1,x2,y1,y2)
CALL xcharl(x1,y1-0.02*(y2-y1),runname)
WRITE(label,'(a,f4.1,a,f6.1,a)') &
'qr+qc=',qisosf,' kg/kg isosurface at t=',time/60.0,' min'
CALL xcharl(x1,y1-0.04*(y2-y1),label )
END IF
IF (profopt == 0) GO TO 900
!
!-----------------------------------------------------------------------
!
! Reinitialize plotting parameters for profile plotting
!
!-----------------------------------------------------------------------
!
CALL xdspac(0.9*winsiz)
IF( nxprpic == 1 .AND. nyprpic == 1) CALL xdspac(0.85*winsiz)
CALL xspace(nxprpic, nyprpic, angl , xlimit,ylimit)
CALL xpmagn(margnx*xlimit, margny*ylimit)
CALL xnwfrm
!
!-----------------------------------------------------------------------
!
! Draw all the profiles based on the user inputs
!
!-----------------------------------------------------------------------
!
IF (uprof == 1) THEN
onvf = 0
CALL avgx(u, onvf, nx,ny,nz, 1, nx-1, 1, ny-1, 1, nz-1, tem1)
onvf = 1
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,uprmin,uprmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'u (m/s)', 'height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (vprof == 1) THEN
onvf = 0
CALL avgy(v, onvf, nx,ny,nz, 1, nx-1, 1, ny-1, 1, nz-1, tem1)
onvf = 1
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,vprmin,vprmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'v (m/s)', 'height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (wprof == 1) THEN
onvf = 0
CALL avgz(w, onvf, nx,ny,nz, 1, nx-1, 1, ny-1, 1, nz-1, tem1)
onvf = 1
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,wprmin,wprmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'w (m/s)', 'height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (ptprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = ptprt(i,j,k)+ptbar(i,j,k)
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,ptprmin,ptprmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'pt (K)', 'height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (pprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = pprt(i,j,k)+pbar(i,j,k)
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,pprmin,pprmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'p (Pa)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qvprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = qv(i,j,k) * 1000.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,qvprmin,qvprmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'qv (g/kg)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qcprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = qc(i,j,k) * 1000.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,qcpmin,qcpmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'qc (g/kg)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qrprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = qr(i,j,k) * 1000.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,qrpmin,qrpmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'qr (g/kg)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qiprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = qi(i,j,k) * 1000.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,qipmin,qipmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'qi (g/kg)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qsprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = qs(i,j,k) * 1000.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,qspmin,qspmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'qs (g/kg)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qhprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = qh(i,j,k) * 1000.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,qhpmin,qhpmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'qh (g/kg)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (kmhprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = kmh(i,j,k)
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,kmhpmin,kmhpmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'kmh (m**2/s)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (kmvprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = kmv(i,j,k) * 1.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,kmvpmin,kmvpmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'kmv (m**2/s)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (tkeprof == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = tke(i,j,k) * 1.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,tkepmin,tkepmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'tke ((m/s)**2)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (rhprof == 1) THEN
CALL temper(nx,ny,nz,pt, pprt,pbar, tem1)
CALL getqvs(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, pbar,tem1,tem2)
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k) = qv(i,j,k)/tem2(i,j,k)
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,rhpmin,rhpmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'RH', 'height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (rfprof == 1) THEN
IF (rfopt == 1) THEN
CALL reflec_ferrier(nx,ny,nz, rhobar, qr, qs, qh, tz, tem5)
ELSE
CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, tem5)
ENDIF
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem5,xc,yc,zpc,rfpmin,rfpmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'ref (dBZ)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (pteprf == 1) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=pprt(i,j,k)+pbar(i,j,k)
END DO
END DO
END DO
CALL pt2pte(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem2,pt,qv,tem1)
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,ptepmin,ptepmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'pte (K)', 'height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (upprof == 1) THEN
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,uprt,xc,yc,zpc,uppmin,uppmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'uprt (m/s)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (vpprof == 1) THEN
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,vprt,xc,yc,zpc,vppmin,vppmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'vprt (m/s)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (wpprof == 1) THEN
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,wprt,xc,yc,zpc,wppmin,wppmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'wprt (m/s)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (ptpprf == 1) THEN
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,ptprt,xc,yc,zpc,ptppmin,ptppmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'ptprt (K)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (ppprof == 1) THEN
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,pprt,xc,yc,zpc,pppmin,pppmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'pprt (Pa)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qvpprf == 1) THEN
DO k = 1, nz
DO j = 1, ny
DO i = 1, nx
tem1(i,j,k) = qvprt(i,j,k) * 1000.0
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,qvppmin,qvppmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'qvprt (g/kg)','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (vorpprf == 1) THEN
DO k=2,nz-2
DO j=2,ny-2
DO i=2,nx-2
tem1(i,j,k)= 1.0E5*( &
(v(i+1,j,k)-v(i-1,j,k)+v(i+1,j+1,k)-v(i-1,j+1,k))/ &
(4*(x(i+1)-x(i)))- &
(u(i,j+1,k)-u(i,j-1,k)+u(i+1,j+1,k)-u(i+1,j-1,k))/ &
(4*(y(j+1)-y(j))) )
END DO
END DO
END DO
DO j=2,ny-2
DO i=2,nx-2
tem1(i,j, 1)=tem1(i,j, 2)
tem1(i,j,nz-1)=tem1(i,j,nz-2)
END DO
END DO
DO k=1,nz-1
DO j=2,ny-2
tem1( 1,j,k)=tem1( 2,j,k)
tem1(nx-1,j,k)=tem1(nx-2,j,k)
END DO
END DO
DO k=1,nz-1
DO i=1,nx-1
tem1(i, 1,k)=tem1(i, 2,k)
tem1(i,ny-1,k)=tem1(i,ny-2,k)
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,vorppmin,vorppmax, &
xprof, yprof, nprof, zprofbgn, zprofend, &
'Vort *10^5','height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (divpprf == 1) THEN
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem1(i,j,k)= 1000.0 * ( &
(u(i+1,j,k)-u(i,j,k))/(x(i+1)-x(i))+ &
(v(i,j+1,k)-v(i,j,k))/(y(j+1)-y(j)) )
END DO
END DO
END DO
nzprofc = nz-1
CALL vprofil (nx,ny,nz,nzprofc,tem1,xc,yc,zpc,divppmin,divppmax, &
xprof,yprof,nprof,zprofbgn,zprofend, &
'Div*1000', 'height (km)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (tsoilprof == 1) THEN
DO k=1,nzsoil
DO j=1,ny
DO i=1,nx
tem1(i,j,k)= tsoil(i,j,k,0)
END DO
END DO
END DO
nzprofc = nzsoil
CALL vprofil (nx,ny,nzsoil,nzprofc,tem1,xc,yc,zpsoilc, &
tsoilprofmin,tsoilprofmax, &
xprof,yprof,nprof,zsoilprofbgn,zsoilprofend, &
'tsoil', 'height (cm)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
IF (qsoilprof == 1) THEN
DO k=1,nzsoil
DO j=1,ny
DO i=1,nx
tem1(i,j,k)= qsoil(i,j,k,0)
END DO
END DO
END DO
nzprofc = nzsoil
CALL vprofil (nx,ny,nzsoil,nzprofc,tem1,xc,yc,zpsoilc, &
qsoilprofmin,qsoilprofmax, &
xprof,yprof,nprof,zsoilprofbgn,zsoilprofend, &
'qsoil', 'height (cm)',npicprof,tem2,tem3)
CALL runlab (runname)
END IF
!
!-----------------------------------------------------------------------
!
! Set up a new page if inwfrm = 1, i.e., flush the plot buffer and
! move to the next plotting page.
!
!-----------------------------------------------------------------------
!
900 CONTINUE
IF( inwfrm == 1 ) CALL xnwfrm
IF ( ireturn == 0 .AND. hinfmt == 9 ) THEN
GO TO 15 ! read one more time in the same GrADS file
ELSE
CLOSE ( nchin ) ! read another data file
END IF ! read another data file
END DO
GO TO 800
700 CONTINUE
WRITE(6,'(a,/a,i3,a)') &
' Bad return status from data file reading', &
' ireturn = ', ireturn,' program stopped.'
800 CALL xgrend
WRITE(6,'(a)') 'Program completed.'
WRITE(6,'(/1x,a,e14.6,a)') &
'Total CPU time used =',f_cputime()-cpu0,'s.'
WRITE (6,*) "Maxumum memory allocation (in words):", &
max_memory_use
GO TO 10
100 WRITE(6,'(a)') &
'Error reading NAMELIST file. Job stopped in ARPSPLT.'
STOP
10 CONTINUE
STOP
END PROGRAM ARPSPLT