PROGRAM arpsplt,547
!
!
!##################################################################
!##################################################################
!###### ######
!###### 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=256) :: grdbasfn
CHARACTER (LEN=256) :: hisfile(nhisfile_max)
COMMON /init2_hisf/ hinfmt,nhisfile, grdbasfn, hisfile
INTEGER :: first_frame
COMMON /frstfrm/ first_frame
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'
INCLUDE 'mp.inc'
!
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
REAL, DIMENSION(:), POINTER :: x ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, DIMENSION(:), POINTER :: y ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, DIMENSION(:), POINTER :: z ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, DIMENSION(:,:,:), POINTER :: zp ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, DIMENSION(:,:,:), POINTER :: zpsoil ! The physical height coordinate defined at
! w-point of the staggered grid for soil model.
REAL, DIMENSION(:,:,:), POINTER :: uprt ! Perturbation u-velocity (m/s)
REAL, DIMENSION(:,:,:), POINTER :: vprt ! Perturbation v-velocity (m/s)
REAL, DIMENSION(:,:,:), POINTER :: wprt ! Perturbation w-velocity (m/s)
REAL, DIMENSION(:,:,:), POINTER :: ptprt ! Perturbation potential temperature
! from that of base state atmosphere (K)
REAL, DIMENSION(:,:,:), POINTER :: pprt ! Perturbation pressure from that
! of base state atmosphere (Pascal)
REAL, DIMENSION(:,:,:), POINTER :: qvprt
REAL, DIMENSION(:,:,:), POINTER :: qv ! Water vapor specific humidity (kg/kg)
REAL, DIMENSION(:,:,:), POINTER :: qc ! Cloud water mixing ratio (kg/kg)
REAL, DIMENSION(:,:,:), POINTER :: qr ! Rain water mixing ratio (kg/kg)
REAL, DIMENSION(:,:,:), POINTER :: qi ! Cloud ice mixing ratio (kg/kg)
REAL, DIMENSION(:,:,:), POINTER :: qs ! Snow mixing ratio (kg/kg)
REAL, DIMENSION(:,:,:), POINTER :: qh ! Hail mixing ratio (kg/kg)
REAL, DIMENSION(:,:,:), POINTER :: tke ! Turbulent Kinetic Energy ((m/s)**2)
REAL, DIMENSION(:,:,:), POINTER :: kmh ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, DIMENSION(:,:,:), POINTER :: kmv ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, DIMENSION(:,:,:), POINTER :: ubar ! Base state u-velocity (m/s)
REAL, DIMENSION(:,:,:), POINTER :: vbar ! Base state v-velocity (m/s)
REAL, DIMENSION(:,:,:), POINTER :: wbar ! Base state w-velocity (m/s)
REAL, DIMENSION(:,:,:), POINTER :: ptbar ! Base state potential temperature (K)
REAL, DIMENSION(:,:,:), POINTER :: pbar ! Base state pressure (Pascal)
REAL, DIMENSION(:,:,:), POINTER :: rhobar ! Base state density rhobar
REAL, DIMENSION(:,:,:), POINTER :: qvbar ! Base state water vapor specific humidity
! (kg/kg)
INTEGER, DIMENSION(:,:,:), POINTER :: soiltyp ! Soil type
INTEGER, DIMENSION(:,:), POINTER :: vegtyp ! Vegetation type
REAL, DIMENSION(:,:,:), POINTER :: stypfrct ! Soil type fraction
REAL, DIMENSION(:,:), POINTER :: lai ! Leaf Area Index
REAL, DIMENSION(:,:), POINTER :: roufns ! Surface roughness
REAL, DIMENSION(:,:), POINTER :: veg ! Vegetation fraction
REAL, DIMENSION(:,:,:,:), POINTER :: tsoil ! soil temperature (K)
REAL, DIMENSION(:,:,:,:), POINTER :: qsoil ! soil moisture
REAL, DIMENSION(:,:,:), POINTER :: wetcanp ! Canopy water amount
REAL, DIMENSION(:,:), POINTER :: snowdpth ! Snow depth (m)
REAL, DIMENSION(:,:), POINTER :: raing ! Grid supersaturation rain
REAL, DIMENSION(:,:), POINTER :: rainc ! Cumulus convective rain
REAL, DIMENSION(:,:), POINTER :: raint ! Total rain (rainc+raing)
REAL, DIMENSION(:,:,:), POINTER :: 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, DIMENSION(:,:,:), POINTER :: radfrc ! Radiation forcing (K/s)
REAL, DIMENSION(:,:), POINTER :: radsw ! Solar radiation reaching the surface
REAL, DIMENSION(:,:), POINTER :: rnflx ! Net radiation flux absorbed by surface
REAL, DIMENSION(:,:), POINTER :: radswnet ! Net shortwave radiation
REAL, DIMENSION(:,:), POINTER :: radlwin ! Incoming longwave radiation
REAL, DIMENSION(:,:), POINTER :: usflx ! Surface flux of u-momentum (kg/(m*s**2))
REAL, DIMENSION(:,:), POINTER :: vsflx ! Surface flux of v-momentum (kg/(m*s**2))
REAL, DIMENSION(:,:), POINTER :: ptsflx ! Surface heat flux (K*kg/(m*s**2))
REAL, DIMENSION(:,:), POINTER :: qvsflx ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
! Arrays derived from the read-in arrays
!
!-----------------------------------------------------------------------
!
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 :: pt (:,:,:) ! Total poten
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 :: acc_raing(:,:,:) ! Accumulated grid supersaturation rain
! between 2 successive time levels
REAL, ALLOCATABLE :: acc_rainc(:,:,:) ! Accumulated cumulus convective rain
! between 2 successive time levels
REAL, ALLOCATABLE :: acc_raint(:,:,:) ! Accumulated rain (rainc+raing)
! between 2 successive time levels
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 :: t00,p00
REAL, 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, vtpunit, 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 :: rainicminc, rainicmaxc, rainigminc, rainigmaxc, &
rainitminc,rainitmaxc
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 :: raicovr,raigovr,raitovr
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
INTEGER :: raichlf,raighlf,raithlf
!
!-----------------------------------------------------------------------
!
! 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
INTEGER :: raiczro,raigzro,raitzro
!
!-----------------------------------------------------------------------
!
! 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 :: raicsty,raigsty,raitsty
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
REAL :: ctinc,ctmin,ctmax,vtunt ! contour interval and vector unit
INTEGER :: icolor,icolor1,lbcolor,trcolor ! required color
CHARACTER (LEN=12) :: varname ! save the plot variable
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 /laypar/ layover
COMMON /incunt/ ctinc,ctmin,ctmax,vtunt
COMMON /recolor/ icolor,icolor1,lbcolor,trcolor
COMMON /varplt1/ varname
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=256) :: filename
CHARACTER (LEN=256) :: outfilename
!
!-----------------------------------------------------------------------
!
INTEGER :: length
CHARACTER (LEN=6) :: stem2
CHARACTER (LEN=1) :: stem1
INTEGER :: flagsin
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 :: xr1,yr1,x11,y11
INTEGER :: nxpic,nypic,islice,jslice,kslice,layout
INTEGER :: idxfmt,filetm
INTEGER :: isoilslice,jsoilslice,ksoilslice ! Zuwen 05/31/2002
REAL :: time,angl,aspect
REAL :: time1c, time2c, time1g, time2g, time1t, time2t, timediff
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
LOGICAL :: rdallhist
INTEGER :: onvf
LOGICAL, PARAMETER :: back=.true.
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 :: rainicplt,rainigplt,rainitplt
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
REAL :: rainicinc,rainiginc,rainitinc
!
! 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=256) :: 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=256) :: 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
CHARACTER (LEN=256) :: stalofl
COMMON /sta_par/ ovrstaopt,ovrstam,staset,ovrstan,ovrstav,stacol, &
markprio,nsta_typ,sta_typ,sta_marktyp, &
sta_markcol,sta_marksz,stalofl,wrtstax,wrtstad
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=256) :: 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, &
raiccol1,raigcol1,raitcol1
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, &
raiccol2,raigcol2,raitcol2
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, &
raicprio,raigprio,raitprio
INTEGER :: col_table,smooth
CHARACTER (LEN=80) :: color_map
INTEGER :: number_of_boxes, boxcol
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
! trajectory namelists and variables -Dan Dawson 12/03/2004
INTEGER :: trajfileflag
INTEGER :: trajopt
INTEGER :: ntimes
INTEGER, PARAMETER :: nmax_times=20
INTEGER :: nunit(nmax_times)
INTEGER, PARAMETER :: npoints_max=10000
INTEGER, PARAMETER :: ntrajcs_max=30
INTEGER :: traj_col(ntrajcs_max)
INTEGER :: labelopt
INTEGER :: labelfreq
REAL :: labelmag
REAL :: tzero, tend, tstart_calc, tend_calc, tinc_calc, reftime(nmax_times)
CHARACTER(LEN=256) :: trajc_fn_in(nmax_times), trajc_fn_out
COMMON /trajectopt/ trajopt,tstart_calc,tend_calc,tinc_calc,traj_col,reftime,trajc_fn_in,ntimes,labelopt,labelfreq,labelmag
INTEGER :: ntrajc0,ntrajcs, ntrajcs_in, npoints, npoints_in, npoints_cur,itrajc1,jtrajc1
INTEGER :: istat
REAL, allocatable :: xtrajc(:,:,:),ytrajc(:,:,:),ztrajc(:,:,:),ttrajc(:)
REAL, allocatable :: xtrajc1(:,:),ytrajc1(:,:),ztrajc1(:,:)
REAL :: xlow, xhigh, ylow, yhigh, zlow, zhigh, pttem
INTEGER :: nstart, nend, npnt, ntrj, ninc
CHARACTER(LEN=20) :: CH
INTEGER :: LCH
REAL :: chsize
!
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 :: raicunit,raigunit,raitunit
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
COMMON /init5_coltab/ color_map
INTEGER :: arbvaropt ! plot arbitrary variable
INTEGER :: istatus
INTEGER :: finfmt3d(maxarbvar), finfmt2d(maxarbvar)
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
COMMON /init8_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, &
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
COMMON /init9_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
COMMON /init10_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, &
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
COMMON /init810_char_units/ tunits, tdunits
COMMON /init11_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
COMMON /init12_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
COMMON /init13_cntl_vctr/ istride,jstride,kstride, &
vtrplt, vtrunit,vtrovr,vtrcol1,vtrcol2,vtrprio,vtrunits,vtrtype, &
vtpplt, vtpunit, 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
COMMON /init14_cntl_strm/ &
vtrstrm, vtrstmovr, vtrstmcol1, vtrstmcol2, vtrstmprio, &
vtpstrm, vtpstmovr, vtpstmcol1, vtpstmcol2, vtpstmprio
COMMON /init4_plotset/ iorig,zbgn,zend,zsoilbgn,zsoilend, &
yxstrch,zxstrch,zystrch, &
zhstrch,zsoilxstrch,zsoilystrch, &
margnx,margny
COMMON /init6_style/ lnmag, ctrlbopt, ctrstyle
COMMON /init7_slice/ nslice_xy, nslice_xz, nslice_yz, nslice_h, &
nslice_v, nslice_p, nslice_pt, &
nslice_xy_soil, nslice_xz_soil, nslice_yz_soil, &
slice_xy, slice_xz, slice_yz, slice_h, &
slice_p, slice_pt, xpnt1, ypnt1, xpnt2, ypnt2, &
slice_xy_soil, slice_xz_soil, slice_yz_soil, &
imove
COMMON /init15_sfc/ &
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, &
rainicplt,rainicinc,rainicminc,rainicmaxc,raicovr,raiccol1, &
raiccol2,raichlf,raicprio,raiczro,raicsty,raicunit, &
rainigplt,rainiginc,rainigminc,rainigmaxc,raigovr,raigcol1, &
raigcol2,raighlf,raigprio,raigzro,raigsty,raigunit, &
rainitplt,rainitinc,rainitminc,rainitmaxc,raitovr,raitcol1, &
raitcol2,raithlf,raitprio,raitzro,raitsty,raitunit
COMMON /init19_soil/ &
tsoilplt,tsoilinc,tsoilminc,tsoilmaxc,tsoilovr, &
tsoilcol1,tsoilcol2,tsoilhlf,tsoilprio,tsoilzro, &
qsoilplt,qsoilinc,qsoilminc,qsoilmaxc,qsoilovr, &
qsoilcol1,qsoilcol2,qsoilhlf,qsoilprio,qsoilzro
COMMON /init16_sfc/ &
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
COMMON /init17_sfc/ &
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
COMMON /init18_sfc/ &
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
COMMON /init20_sfccha/ &
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
COMMON /init23_wirfrm/ w3dplt, wisosf, q3dplt, qisosf
REAL :: paprlnth
COMMON /init3_page/ layout, inwfrm, paprlnth
COMMON /init25_prof/ 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
COMMON /init24_obs/ sfcobfl
COMMON /init22_ovrlay/ovrlaymulopt,ovrname,ovrmul_num,ovrmulname
COMMON /init21_cntl_arbvar/arbvaropt, &
var3dnum,dirname3d,finfmt3d, &
var3d,var3dplot, var3dinc, var3dminc,var3dmaxc, &
var3dovr, var3dhlf, var3dzro,var3dsty,var3dcol1, var3dcol2, &
var3dprio, var2dnum,dirname2d,finfmt2d, &
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
!
!-----------------------------------------------------------------------
!
! Variables for mpi jobs
!
!-----------------------------------------------------------------------
INTEGER, PARAMETER :: fzone = 3
INTEGER :: ncompressx, ncompressy ! compression in x and y direction:
! ncompressx=nprocx_in/nproc_x
! ncompressy=nprocy_in/nproc_y
INTEGER :: nproc_node
COMMON /init1_mpi/ ncompressx, ncompressy, nproc_node
INTEGER :: nxsm, nysm ! smaller domain
INTEGER :: nxlg, nylg ! global domain
INTEGER :: ii,jj,ia,ja
REAL, DIMENSION(:), POINTER :: xsm,ysm
REAL, DIMENSION(:,:,:), POINTER :: zpsm,zpsoilsm
REAL, DIMENSION(:,:,:), POINTER :: uprtsm, vprtsm, wprtsm, &
ptprtsm, pprtsm, qvprtsm
REAL, DIMENSION(:,:,:), POINTER :: qcsm, qrsm, qism, qssm, qhsm
REAL, DIMENSION(:,:,:), POINTER :: tkesm, kmhsm, kmvsm
REAL, DIMENSION(:,:,:), POINTER :: ubarsm, vbarsm, wbarsm, &
ptbarsm,pbarsm, qvbarsm,rhobarsm
REAL, DIMENSION(:,:,:), POINTER :: prcratesm
INTEGER, DIMENSION(:,:,:), POINTER :: soiltypsm
INTEGER, DIMENSION(:,:), POINTER :: vegtypsm
REAL, DIMENSION(:,:,:), POINTER :: stypfrctsm
REAL, DIMENSION(:,:,:,:),POINTER :: tsoilsm, qsoilsm
REAL, DIMENSION(:,:,:), POINTER :: wetcanpsm
REAL, DIMENSION(:,:), POINTER :: laism, roufnssm, vegsm, &
snowdpthsm, raingsm, raincsm
REAL, DIMENSION(:,:,:), POINTER :: radfrcsm(:,:,:)
REAL, DIMENSION(:,:), POINTER :: radswsm, rnflxsm, radswnetsm, radlwinsm
REAL, DIMENSION(:,:), POINTER :: usflxsm, vsflxsm, ptsflxsm, qvsflxsm
INTEGER :: procbgn ! processor rank holds (ibgn,jbgn)
INTEGER :: xpbgn, xpend, ypbgn, ypend
INTEGER :: xpbgn1,xpend1,ypbgn1,ypend1
INTEGER :: ibgnl,iendl,jbgnl,jendl
INTEGER :: ibgnl1, iendl1
INTEGER :: ibgnl2, jbgnl2 ! hold the value of ibgnl and jbgnl
! for vector plot
INTEGER :: ibgnl21
COMMON /processors/ xpbgn,xpend,ypbgn,ypend
INTEGER :: lenfilb ! the length of base file name, not including
! processor information
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directives for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
REAL :: amin, amax, tem, pref
REAL :: xr2, yr2
INTEGER :: nch
COMMON /XOUTCH/ nch
INTEGER :: pen_status
!fpp$ expand (f_qvsat)
!!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
cpu0 = f_cputime
()
xi1= 50.0
yi1= 50.0
xi2= 50.0
yi2= 50.0
trcolor = 183 ! xz slice terrain color
plotovr = .false.
sovrlay = 0
lbcolor = 1 ! laber color 1-black
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
first_frame = 1
!
!-----------------------------------------------------------------------
!
! 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.
time1c = 0.0
time2c = 0.0
time1g = 0.0
time2g = 0.0
time1t = 0.0
time2t = 0.0
timediff = 0.0
trajfileflag = 0
CALL initpltpara(nx,ny,nz,nzsoil,nstyps,outfilename)
!-----------------------------------------------------------------------
!
! Set some parameters
!
!-----------------------------------------------------------------------
trcolor = trncol1
nxsm = (nx-3)/ncompressx + 3
nysm = (ny-3)/ncompressy + 3
nxlg = (nx-3)*nproc_x + 3
nylg = (ny-3)*nproc_y + 3
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
ipriority(44)=qsoilprio
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
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
ipriority(111)=raicprio
ipriority(112)=raigprio
ipriority(113)=raitprio
DO i=1, var3dnum
ipriority(150+i) = var3dprio(i)
END DO
DO i=1, var2dnum
ipriority(170+i) = var2dprio(i)
END DO
nslice(1)=nslice_xy
nslice(2)=nslice_xz
nslice(3)=nslice_yz
nslice(4)=nslice_h
nslice(5)=nslice_v
nslice(6)=nslice_p
nslice(7)=nslice_pt
nslice(8)= 1 ! surface plots
iplot(8) = 1
iplot(1:7)=1
IF(nslice(1) == 0) iplot(1)=0
IF(nslice(2) == 0) iplot(2)=0
IF(nslice(3) == 0) iplot(3)=0
IF(nslice(4) == 0) iplot(4)=0
IF(nslice(5) == 0) iplot(5)=0
IF(nslice(6) == 0) iplot(6)=0
IF(nslice(7) == 0) iplot(7)=0
nslice(9) =nslice_xy_soil
nslice(10)=nslice_xz_soil
nslice(11)=nslice_yz_soil
iplot(9:11)=1
IF(nslice(9) == 0) iplot(9)=0
IF(nslice(10) == 0) iplot(10)=0
IF(nslice(11) == 0) iplot(11)=0
DO indxslic = 1,nslice(5)
xi1(indxslic)=xpnt1(indxslic)
yi1(indxslic)=ypnt1(indxslic)
xi2(indxslic)=xpnt2(indxslic)
yi2(indxslic)=ypnt2(indxslic)
END DO
IF(arbvaropt < 10) THEN
rdallhist=.true.
ELSE
arbvaropt=mod(arbvaropt,10)
rdallhist=.false.
END IF
umove_readin = umove
vmove_readin = vmove
!
!-----------------------------------------------------------------------
!
! 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
!
!-----------------------------------------------------------------------
!
! 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(acc_rainc(nx,ny,2),STAT=istatus)
ALLOCATE(acc_raing(nx,ny,2),STAT=istatus)
ALLOCATE(acc_raint(nx,ny,2),STAT=istatus)
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)
CALL check_alloc_status(istatus, "arpsplt:mapfct")
allocate(xtrajc(npoints_max,ntrajcs_max,nmax_times),stat=istatus)
CALL check_alloc_status(istatus, "arpsplt:xtrajc")
allocate(ytrajc(npoints_max,ntrajcs_max,nmax_times),stat=istatus)
CALL check_alloc_status(istatus, "arpsplt:ytrajc")
allocate(ztrajc(npoints_max,ntrajcs_max,nmax_times),stat=istatus)
CALL check_alloc_status(istatus, "arpsplt:ztrajc")
allocate(ttrajc(npoints_max),stat=istatus)
CALL check_alloc_status(istatus, "arpsplt:ttrajc")
allocate(xtrajc1(ntrajcs_max,nmax_times),stat=istatus)
allocate(ytrajc1(ntrajcs_max,nmax_times),stat=istatus)
allocate(ztrajc1(ntrajcs_max,nmax_times),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
acc_rainc = 0.0
acc_raing = 0.0
acc_raint = 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
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
!-----------------------------------------------------------------------
!
! Allocate arrays for smaller subdomain if needed
!
!-----------------------------------------------------------------------
IF(ncompressx > 1 .OR. ncompressy > 1) THEN ! allocate arrays
! otherwise, just link pointers
ALLOCATE(xsm (nxsm), STAT=istatus)
ALLOCATE(ysm (nysm), STAT=istatus)
ALLOCATE(zpsm (nxsm,nysm,nz), STAT=istatus)
ALLOCATE(zpsoilsm(nxsm,nysm,nzsoil),STAT=istatus)
ALLOCATE(uprtsm(nxsm,nysm,nz), STAT=istatus)
ALLOCATE(vprtsm(nxsm,nysm,nz), STAT=istatus)
ALLOCATE(wprtsm(nxsm,nysm,nz), STAT=istatus)
ALLOCATE(ptprtsm(nxsm,nysm,nz),STAT=istatus)
ALLOCATE(pprtsm(nxsm,nysm,nz), STAT=istatus)
ALLOCATE(qvprtsm(nxsm,nysm,nz),STAT=istatus)
ALLOCATE(qcsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(qrsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(qism (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(qssm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(qhsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(tkesm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(kmhsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(kmvsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(ubarsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(vbarsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(wbarsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(ptbarsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(pbarsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(rhobarsm(nxsm,nysm,nz),STAT=istatus)
ALLOCATE(qvbarsm (nxsm,nysm,nz),STAT=istatus)
ALLOCATE(soiltypsm (nxsm,nysm,nstyps),STAT=istatus)
ALLOCATE(stypfrctsm(nxsm,nysm,nstyps),STAT=istatus)
ALLOCATE(vegtypsm (nxsm,nysm),STAT=istatus)
ALLOCATE(laism (nxsm,nysm),STAT=istatus)
ALLOCATE(roufnssm (nxsm,nysm),STAT=istatus)
ALLOCATE(vegsm (nxsm,nysm),STAT=istatus)
ALLOCATE(tsoilsm (nxsm,nysm,nzsoil,0:nstyps),STAT=istatus)
ALLOCATE(qsoilsm (nxsm,nysm,nzsoil,0:nstyps),STAT=istatus)
ALLOCATE(wetcanpsm (nxsm,nysm,0:nstyps),STAT=istatus)
ALLOCATE(snowdpthsm(nxsm,nysm),STAT=istatus)
ALLOCATE(raingsm (nxsm,nysm),STAT=istatus)
ALLOCATE(raincsm (nxsm,nysm),STAT=istatus)
ALLOCATE(prcratesm (nxsm,nysm,4),STAT=istatus)
ALLOCATE(radfrcsm(nxsm,nysm,nz),STAT=istatus)
ALLOCATE(radswsm (nxsm,nysm),STAT=istatus)
ALLOCATE(rnflxsm (nxsm,nysm),STAT=istatus)
ALLOCATE(radswnetsm (nxsm,nysm),STAT=istatus)
ALLOCATE(radlwinsm (nxsm,nysm),STAT=istatus)
ALLOCATE(usflxsm (nxsm,nysm),STAT=istatus)
ALLOCATE(vsflxsm (nxsm,nysm),STAT=istatus)
ALLOCATE(ptsflxsm(nxsm,nysm),STAT=istatus)
ALLOCATE(qvsflxsm(nxsm,nysm),STAT=istatus)
CALL check_alloc_status(istatus, "arpsplt:qvsflxsm")
xsm = 0.0
ysm = 0.0
zpsm = 0.0
zpsoilsm= 0.0
uprtsm = 0.0
vprtsm = 0.0
wprtsm = 0.0
ptprtsm = 0.0
pprtsm = 0.0
qvprtsm = 0.0
qcsm = 0.0
qrsm = 0.0
qism = 0.0
qssm = 0.0
qhsm = 0.0
tkesm = 0.0
kmhsm = 0.0
kmvsm = 0.0
ubarsm = 0.0
vbarsm = 0.0
wbarsm = 0.0
ptbarsm = 0.0
pbarsm = 0.0
rhobarsm= 0.0
qvbarsm = 0.0
soiltypsm = 0.0
stypfrctsm= 0.0
vegtypsm = 0.0
laism = 0.0
roufnssm = 0.0
vegsm = 0.0
tsoilsm = 0.0
qsoilsm = 0.0
wetcanpsm = 0.0
snowdpthsm= 0.0
raingsm = 0.0
raincsm = 0.0
prcratesm = 0.0
radfrcsm = 0.0
radswsm = 0.0
rnflxsm = 0.0
radswnetsm = 0.0
radlwinsm = 0.0
usflxsm = 0.0
vsflxsm = 0.0
ptsflxsm = 0.0
qvsflxsm = 0.0
ELSE
xsm => x
ysm => y
zpsm => zp
zpsoilsm => zpsoil
uprtsm => uprt
vprtsm => vprt
wprtsm => wprt
ptprtsm => ptprt
pprtsm => pprt
qvprtsm => qvprt
qcsm => qc
qrsm => qr
qism => qi
qssm => qs
qhsm => qh
tkesm => tke
kmhsm => kmh
kmvsm => kmv
ubarsm => ubar
vbarsm => vbar
wbarsm => wbar
ptbarsm => ptbar
pbarsm => pbar
rhobarsm => rhobar
qvbarsm => qvbar
soiltypsm => soiltyp
stypfrctsm => stypfrct
vegtypsm => vegtyp
laism => lai
roufnssm => roufns
vegsm => veg
tsoilsm => tsoil
qsoilsm => qsoil
wetcanpsm => wetcanp
snowdpthsm => snowdpth
raingsm => raing
raincsm => rainc
prcratesm => prcrate
radfrcsm => radfrc
radswsm => radsw
rnflxsm => rnflx
radswnetsm => radswnet
radlwinsm => radlwin
usflxsm => usflx
vsflxsm => vsflx
ptsflxsm => ptsflx
qvsflxsm => qvsflx
END IF
IF(myproc == 0) THEN
!-----------------------------------------------------------------------
!
! Overlaying observations
!
!-----------------------------------------------------------------------
IF(ovrobs == 1) THEN
lsfcobfl=LEN(sfcobfl)
CALL strlnth(sfcobfl,lsfcobfl)
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
!-----------------------------------------------------------------------
!
! Overlaying airport location
!
!-----------------------------------------------------------------------
!
IF(ovrstaopt /= 0) THEN
lstalofl=LEN(stalofl)
CALL strlnth(stalofl,lstalofl)
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
END IF ! myproc == 0
CALL mpupdatei(n_obs_b,1)
CALL mpupdatei(staset,1)
!-----------------------------------------------------------------------
!
! Plot beginning below
!
!-----------------------------------------------------------------------
zxplot_called=0
DO nf=1, nhisfile
sigplt(:) = 0
filename=hisfile(nf)
lenfilb = len_trim(filename)
lengbf = len_trim(grdbasfn)
second1= f_walltime()
cpu1 = f_cputime()
IF(nf==1 .OR. rdallhist) THEN
DO jj = 1, ncompressy
DO ii = 1, ncompressx
IF (mp_opt > 0 .AND. readsplit <= 0) THEN
WRITE(grdbasfn(lengbf-5+1:lengbf) ,'(a,i2.2,i2.2)') '_', &
ncompressx*(loc_x-1)+ii,ncompressy*(loc_y-1)+jj
lengbf = lengbf
WRITE(filename(lenfilb+1:lenfilb+5),'(a,i2.2,i2.2)') '_', &
ncompressx*(loc_x-1)+ii,ncompressy*(loc_y-1)+jj
lenfil= lenfilb + 5
ELSE
lenfil = lenfilb
END IF
! 15 CONTINUE ! also continue to read another time recode
! ! from GrADS file
!
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
CALL setgbrd(0) ! read grid/base state file
IF (nproc_node <= 1) THEN ! the first readstride processes read then
! the next readstride processes and so on
! blocking inserted for ordering i/o for message passing
! DO i=nprocs-1,0,-1*readstride
! IF(myproc <= i .AND. myproc >= i-readstride+1) THEN
DO i=0,nprocs-1,readstride
IF(myproc >= i .AND. myproc <= i+readstride-1) THEN
IF (myproc == 0 .OR. readsplit < 1) WRITE(6,'(a,I5,a,a/)') &
'processor: ',myproc, ' reading file: ', filename(1:lenfil)
CALL dtaread(nxsm,nysm,nz,nzsoil, nstyps, &
hinfmt, nchin,grdbasfn(1:lengbf),lengbf, &
filename(1:lenfil),lenfil,time, &
xsm,ysm,z,zpsm,zpsoilsm, &
uprtsm,vprtsm,wprtsm,ptprtsm, pprtsm, qvprtsm, &
qcsm,qrsm,qism,qssm,qhsm,tkesm,kmhsm,kmvsm, &
ubarsm,vbarsm,wbarsm,ptbarsm,pbarsm,rhobarsm,qvbarsm, &
soiltypsm,stypfrctsm,vegtypsm,laism,roufnssm,vegsm, &
tsoilsm,qsoilsm,wetcanpsm,snowdpthsm, &
raingsm,raincsm,prcratesm, &
radfrcsm,radswsm,rnflxsm,radswnetsm,radlwinsm, &
usflxsm,vsflxsm,ptsflxsm,qvsflxsm, &
ireturn, tem1,tem2, tem3)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
ELSE ! Only one process read at one node
DO i=0,nproc_node-1
IF (MOD(myproc,nproc_node) == i) THEN
IF (myproc == 0 .OR. readsplit < 1) WRITE(6,'(a,I5,a,a/)') &
'processor: ',myproc, ' reading file: ', filename(1:lenfil)
CALL dtaread(nxsm,nysm,nz,nzsoil, nstyps, &
hinfmt, nchin,grdbasfn(1:lengbf),lengbf, &
filename(1:lenfil),lenfil,time, &
xsm,ysm,z,zpsm,zpsoilsm, &
uprtsm,vprtsm,wprtsm,ptprtsm, pprtsm, qvprtsm, &
qcsm,qrsm,qism,qssm,qhsm,tkesm,kmhsm,kmvsm, &
ubarsm,vbarsm,wbarsm,ptbarsm,pbarsm,rhobarsm,qvbarsm, &
soiltypsm,stypfrctsm,vegtypsm,laism,roufnssm,vegsm, &
tsoilsm,qsoilsm,wetcanpsm,snowdpthsm, &
raingsm,raincsm,prcratesm, &
radfrcsm,radswsm,rnflxsm,radswnetsm,radlwinsm, &
usflxsm,vsflxsm,ptsflxsm,qvsflxsm, &
ireturn, tem1,tem2, tem3)
END IF
IF (mp_opt > 0) CALL mpbarrier
END DO
END IF
!
!-----------------------------------------------------------------------
!
! 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
IF(ncompressx > 1 .OR. ncompressy > 1) THEN ! need join
DO j = 1, nysm
ja = (jj-1)*(nysm-3)+j
DO i = 1, nxsm
ia = (ii-1)*(nxsm-3)+i
x(ia) = xsm(i)
vegtyp(ia,ja) = vegtypsm(i,j)
lai(ia,ja) = laism(i,j)
roufns(ia,ja) = roufnssm(i,j)
veg(ia,ja) = vegsm(i,j)
snowdpth(ia,ja) = snowdpthsm(i,j)
raing(ia,ja) = raingsm(i,j)
rainc(ia,ja) = raincsm(i,j)
radsw(ia,ja) = radswsm(i,j)
radswnet(ia,ja) = radswnetsm(i,j)
radlwin(ia,ja) = radlwinsm(i,j)
rnflx(ia,ja) = rnflxsm(i,j)
usflx(ia,ja) = usflxsm(i,j)
vsflx(ia,ja) = vsflxsm(i,j)
ptsflx(ia,ja) = ptsflxsm(i,j)
qvsflx(ia,ja) = qvsflxsm(i,j)
prcrate(ia,ja,:) = prcratesm(i,j,:)
zp(ia,ja,:) = zpsm(i,j,:)
zpsoil(ia,ja,:) = zpsoilsm(i,j,:)
uprt(ia,ja,:) = uprtsm(i,j,:)
vprt(ia,ja,:) = vprtsm(i,j,:)
wprt(ia,ja,:) = wprtsm(i,j,:)
ptprt(ia,ja,:) = ptprtsm(i,j,:)
pprt(ia,ja,:) = pprtsm(i,j,:)
qvprt(ia,ja,:) = qvprtsm(i,j,:)
qc(ia,ja,:) = qcsm(i,j,:)
qr(ia,ja,:) = qrsm(i,j,:)
qi(ia,ja,:) = qism(i,j,:)
qs(ia,ja,:) = qssm(i,j,:)
qh(ia,ja,:) = qhsm(i,j,:)
tke(ia,ja,:) = tkesm(i,j,:)
kmh(ia,ja,:) = kmhsm(i,j,:)
kmv(ia,ja,:) = kmvsm(i,j,:)
ubar(ia,ja,:) = ubarsm(i,j,:)
vbar(ia,ja,:) = vbarsm(i,j,:)
wbar(ia,ja,:) = wbarsm(i,j,:)
ptbar(ia,ja,:) = ptbarsm(i,j,:)
pbar(ia,ja,:) = pbarsm(i,j,:)
rhobar(ia,ja,:) = rhobarsm(i,j,:)
qvbar(ia,ja,:) = qvbarsm(i,j,:)
radfrc(ia,ja,:) = radfrcsm(i,j,:)
soiltyp(ia,ja,:) = soiltypsm(i,j,:)
stypfrct(ia,ja,:) = stypfrctsm(i,j,:)
tsoil(ia,ja,:,:) = tsoilsm(i,j,:,:)
qsoil(ia,ja,:,:) = qsoilsm(i,j,:,:)
wetcanp(ia,ja,:) = wetcanpsm(i,j,:)
END DO !i
y(ja) = ysm(j)
END DO !j
END IF ! need join
END DO !ii
END DO ! jj
! IF(ncompressx > 1 .OR. ncompressy > 1) THEN ! DEALLOCATE arrays
!
! DEALLOCATE(xsm, ysm, zpsm, zpsoilsm)
!
! DEALLOCATE(uprtsm,vprtsm,wprtsm,ptprtsm,pprtsm,qvprtsm)
! DEALLOCATE(qcsm, qrsm, qism, qssm, qhsm)
! DEALLOCATE(tkesm, kmhsm, kmvsm)
! DEALLOCATE(ubarsm, vbarsm, wbarsm, ptbarsm, pbarsm)
! DEALLOCATE(rhobarsm, qvbarsm)
!
! DEALLOCATE(soiltypsm, stypfrctsm, vegtypsm)
! DEALLOCATE(laism, roufnssm, vegsm)
! DEALLOCATE(tsoilsm, qsoilsm, wetcanpsm, snowdpthsm)
!
! DEALLOCATE(raingsm, raincsm, prcratesm)
! DEALLOCATE(radfrcsm, radswsm, rnflxsm, radswnetsm, radlwinsm)
! DEALLOCATE(usflxsm, vsflxsm, ptsflxsm, qvsflxsm)
! END IF
cpu2 = f_cputime()
second2 = f_walltime()
! print*,'!!!! total clock time for read data:',second2-second1
! print*,'!!!! total cpu time for read data :', cpu2-cpu1
!
!-----------------------------------------------------------------------
!
! 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 .AND. myproc == 0) THEN
!
!-----------------------------------------------------------------------
!
! Initialize ZXPLOT plotting package
!
!-----------------------------------------------------------------------
!
! The following two ZXPLOT routines, if called, should be called before xdevic
!
zxout_unit=2
IF (LEN_TRIM(outfilename) < 1) THEN
WRITE(outfilename,'(a)') TRIM(dirname)//runname(1:lfnkey)//'.ps'
ELSE
WRITE(outfilename,'(a)') TRIM(dirname)//TRIM(outfilename)//'.ps'
END IF
CALL xpsfn(TRIM(outfilename), 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,'(a,a1,a)') '(f8.',stem1,')'
CALL xaxfmt( stem2 )
END IF
layover = 0
CALL xpmagn(margnx*xlimit, margny*ylimit)
END IF ! Initalization of plotting package and parameters
!
! Somewhere in one of the above calls, NCH gets set. Pass it around.
!
CALL mpupdatei(nch,1)
zxplot_called=1
!
!-----------------------------------------------------------------------
!
! 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
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)
! IF (mp_opt > 0) THEN ! If we know the read-in data is valid
! ! there is no need to call message passing
! CALL acc_interrupt(mp_acct)
! CALL mpsendrecv2dew(u,nx,ny,nz,1,1,1,tem4)
! CALL mpsendrecv2dns(u,nx,ny,nz,1,1,1,tem5)
!
! CALL mpsendrecv2dew(v,nx,ny,nz,1,1,2,tem4)
! CALL mpsendrecv2dns(v,nx,ny,nz,1,1,2,tem5)
!
! CALL mpsendrecv2dew(qv,nx,ny,nz,1,1,0,tem4)
! CALL mpsendrecv2dns(qv,nx,ny,nz,1,1,0,tem5)
!
! CALL mpsendrecv2dew(pprt,nx,ny,nz,1,1,0,tem4)
! CALL mpsendrecv2dns(pprt,nx,ny,nz,1,1,0,tem5)
!
! CALL mpsendrecv2dew(pbar,nx,ny,nz,1,1,0,tem4)
! CALL mpsendrecv2dns(pbar,nx,ny,nz,1,1,0,tem5)
! END IF
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 .OR. &
strmplt /= 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)
CALL mpupdatei(imax,1) ! only processor 0 contains the right indices
CALL mpupdatei(jmax,1)
CALL mpupdatei(kmax,1)
CALL mpupdatei(imin,1)
CALL mpupdatei(jmin,1)
CALL mpupdatei(kmin,1)
! 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)
!ibgn, iend, jbgn, jend, kbgn, kend etc. computed above are over global domain
!Introduced ibgnl,iendl, jbgnl,jendl below for message passing
!version, which are relative to the local domain. - WYH
xpbgn = (ibgn-2)/(nx-3) + 1
xpend = (iend-2)/(nx-3) + 1
ypbgn = (jbgn-2)/(ny-3) + 1
ypend = (jend-2)/(ny-3) + 1
xpbgn1 = xpbgn
xpend1 = xpend
ypbgn1 = ypbgn
ypend1 = ypend
IF (loc_x == xpbgn) THEN
ibgnl = MOD((ibgn-2),(nx-3)) + 2
ELSE IF (loc_x < xpbgn .OR. loc_x > xpend) THEN
ibgnl = 2 ! can be any value
ELSE
ibgnl = 1 ! should be 2, 1 for overlay
END IF
IF (loc_x == xpend) THEN
iendl = MOD((iend-2),(nx-3)) + 2
ELSE IF (loc_x > xpend .OR. loc_x < xpbgn) THEN
iendl = nx-2 ! can be any value
ELSE
iendl = nx - 1 ! should be nx-2, nx-1 for overlay
END IF
IF (loc_y == ypbgn) THEN
jbgnl = MOD((jbgn-2),(ny-3)) + 2
ELSE IF (loc_y < ypbgn .OR. loc_y > ypend) THEN
jbgnl = 2
ELSE
jbgnl = 1
END IF
IF (loc_y == ypend) THEN
jendl = MOD((jend-2),(ny-3)) + 2
ELSE IF (loc_y > ypend .OR. loc_y < ypbgn) THEN
jendl = ny-2
ELSE
jendl = ny - 1
END IF
ibgnl1 = ibgnl
iendl1 = iendl
! END -- WYH
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
ibgnl2 = (loc_x-1)*(nx-3)+ibgnl ! hold global index for ibgnl temporarily
jbgnl2 = (loc_y-1)*(ny-3)+jbgnl ! hold global index for jbgnl
IF(MOD(ibgnl2-ibgn,ist) /= 0) THEN
ibgnl2 = ibgnl + (ist- MOD(ibgnl2 - ibgn,ist)) ! for vector plot
ELSE
ibgnl2 = ibgnl
END IF
IF(MOD(jbgnl2-jbgn,jst) /= 0) THEN
jbgnl2 = jbgnl + (jst- MOD(jbgnl2 - jbgn,jst)) ! for vector plot
ELSE
jbgnl2 = jbgnl
END IF
ibgnl21 = ibgnl2
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 = (iend-ibgn)*dxkm + dxkm
yr = (jend-jbgn)*dykm + dykm
procbgn = (ypbgn-1)*nproc_x+xpbgn-1 ! processor rank which holds (ibgn,jbgn)
IF (myproc == procbgn) THEN
x1 = xc(ibgnl,2,2)-dxkm/2
y1 = yc(2,jbgnl,2)-dykm/2
END IF
CALL mpbcastr(x1,procbgn)
CALL mpbcastr(y1,procbgn)
x2 = x1 + xr
y2 = y1 + yr
IF( zbgn /= zend ) THEN
zmin=zbgn ! these values are already global
zmax=zend
ELSE
IF (loc_x >= xpbgn .AND. loc_x <= xpend .AND. &
loc_y >= ypbgn .AND. loc_y <= ypend) THEN
zmax=zp(ibgnl,jbgnl,kend+1)
zmin=zp(ibgnl,jbgnl,kbgn)
DO j=jbgnl,jendl
DO i=ibgnl,iendl
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.
ELSE ! we do not care their values
zmin = 1.0E6 ! should be large enough
zmax = -1.0E6 ! should be small enough
END IF
CALL mpmax0(zmax,zmin) ! get global values
END IF
zr = zmax-zmin
z1 = zmin
z2 = zmax
ibgn1 = ibgn
iend1 = iend
xr2 = xr ! save the window to be reseted later.
yr2 = yr
x11 = x1
y11 = y1
!
! 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
IF (loc_x >= xpbgn .AND. loc_x <= xpend .AND. &
loc_y >= ypbgn .AND. loc_y <= ypend) THEN
zsoilmax=zpsoilc(ibgnl,jbgnl,ksoilbgn)
zsoilmin=zpsoilc(ibgnl,jbgnl,ksoilend)
DO j=jbgnl,jendl
DO i=ibgnl,iendl
zsoilmax=AMAX1(zsoilmax,zpsoilc(i,j,ksoilbgn))
zsoilmin=AMIN1(zsoilmin,zpsoilc(i,j,ksoilend))
END DO
END DO
ELSE ! we do not care their values
zsoilmin = 1.0E6 ! should be large enough
zsoilmax = -1.0E6 ! should be small enough
END IF
CALL mpmax0(zsoilmax,zsoilmin)
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,'(a,a1,a)') '(f8.',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 = (nxlg-3)*dxkm
yl = (nylg-3)*dykm
IF (myproc == 0) THEN
xorig_0 = (xc(1,2,2)+xc(2,2,2))*0.5
yorig_0 = (yc(2,1,2)+yc(2,2,2))*0.5
END IF
CALL mpupdater(xorig_0,1)
CALL mpupdater(yorig_0,1)
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,ibgnl,iendl,nx,jbgnl,jendl,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
CALL mpmax0(ztmax,ztmin)
!
!-----------------------------------------------------------------------
!
! 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))
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
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)
ELSE
idxfmt=0
IF(hinfmt==1) THEN
idxfmt=INDEX(filename,'bin',back)
ELSE IF(hinfmt==2) THEN
idxfmt=INDEX(filename,'asc',back)
ELSE IF(hinfmt==3) THEN
idxfmt=INDEX(filename,'hdf',back)
ELSE IF(hinfmt==4) THEN
idxfmt=INDEX(filename,'pak',back)
ELSE IF(hinfmt==5) THEN
idxfmt=INDEX(filename,'svi',back)
ELSE IF(hinfmt==6) THEN
idxfmt=INDEX(filename,'bn2',back)
ELSE IF(hinfmt==7) THEN
idxfmt=INDEX(filename,'net',back)
ELSE IF(hinfmt==8) THEN
idxfmt=INDEX(filename,'net',back)
ELSE IF(hinfmt==9) THEN
idxfmt=INDEX(filename,'gad',back)
ELSE IF(hinfmt==10) THEN
idxfmt=INDEX(filename,'grb',back)
ELSE IF(hinfmt==11) THEN
idxfmt=INDEX(filename,'v5d',back)
END IF
IF(idxfmt > 0) THEN
READ(filename((idxfmt+3):(idxfmt+8)),'(i6)') filetm
WRITE(6,'(/a,i7)') ' Read file time as: ',filetm
time=FLOAT(filetm)
curtim=time
END IF
END IF !rdallhist
!
!-----------------------------------------------------------------------
!
! 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
first_frame = 1
!
! 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: surface xy slice
! 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 = xr2
yr = yr2
x1 = x11
x2 = x1+xr
y1 = y11
y2 = y1+yr
z1 = zmin
z2 = zmax
zsoil1 = zsoilmin
zsoil2 = zsoilmax
xpbgn = xpbgn1
xpend = xpend1
ypbgn = ypbgn1
ypend = ypend1
ibgnl = ibgnl1
iendl = iendl1
ibgnl2 = ibgnl21
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)
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 = (nylg-2)/2+1
IF (jslice == -2) jslice = jmax
IF (islice == -1) islice = (nxlg-2)/2+1
IF (islice == -2) islice = imax
IF (slicopt == 2) THEN
ypbgn = (jslice-2)/(ny-3) + 1
ypend = ypbgn
jslice = MOD((jslice-2),(ny-3)) + 2
END IF
IF (slicopt == 3) THEN
xpbgn = (islice-2)/(nx-3) + 1
xpend = xpbgn
islice = MOD((islice-2),(nx-3)) + 2
END IF
! 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
IF (n > nx+ny) THEN
WRITE(6,*) 'ERROR: Working arrays are too small for ', &
' arbitrary vertical slice. '
WRITE(6,*) ' See restriction for MPI mode in arpsplt.input'
CALL arpsstop('Local arrays too small',1)
END IF
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
xpbgn = 1
xpend = 1
xpbgn = 1
ypend = 1
ibgnl = 1
iendl = n
ibgnl2 = 1
xr = (n-1)*sqrtdxy
z1 = zmin
z2 = zmax
IF (zbgn /= zend) THEN
z1 = MAX(zbgn,zmin)
z2 = MIN(zend,zmax)
END IF
zr = z2-z1
x1 = 0.0
x2 = x1+xr
IF (zhstrch /= 0.0) THEN
aspratio = zhstrch
ELSE
aspratio = xr/zr
END IF
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 = (nylg-2)/2+1
IF (isoilslice == -1) isoilslice = (nxlg-2)/2+1
IF (slicopt == 10) THEN
ypbgn = (jsoilslice-2)/(ny-3) + 1
ypend = (jsoilslice-2)/(ny-3) + 1
jsoilslice = MOD((jsoilslice-2),(ny-3)) + 2
END IF
IF (slicopt == 11) THEN
xpbgn = (isoilslice-2)/(nx-3) + 1
xpend = (isoilslice-2)/(nx-3) + 1
isoilslice = MOD((isoilslice-2),(nx-3)) + 2
END IF
END IF
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,ibgnl,iendl,ny,jbgnl,jendl,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,ibgnl,iendl,ny, &
jbgnl,jendl, label(1:27),time, runname,1.0, &
tem1,tem2,tem3,tem4,tem5,hterain,slicopt,thkplt)
! size of tem5 must be >= 6*nx*ny
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
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,ibgnl2,iendl,ist, &
ny,jbgnl2,jendl,jst,label(1:15),time, runname,1.0, &
slicopt,tem1,tem2,tem3,tem4,tem5,tem6,hterain)
! size of tem6 must be >= 5*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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 == 2) 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,ibgnl,iendl, ny,jbgnl,jendl, 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 == 2) 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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:19),time, runname, 1.0,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,rfcplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,tem1)
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,ibgnl,iendl, ny,jbgnl,jendl, 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,tem1)
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,ibgnl,iendl, ny,jbgnl,jendl, 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,tem1)
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,ibgnl,iendl, ny,jbgnl,jendl, 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 ( vtpunit )
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,ibgnl2,iendl,ist,ny,jbgnl2,jendl,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,ibgnl,iendl,ist, &
ny,jbgnl,jendl,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,ibgnl,iendl,ist, &
ny,jbgnl,jendl,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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl,iendl, ny,jbgnl,jendl, 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,ibgnl2,iendl,ist, ny,jbgnl2,jendl,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) THEN
CALL vtr3d( tem7,tem8,tem9,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, nx,ibgnl2,iendl,ist,&
ny,jbgnl2,jendl,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
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,ibgnl,iendl, ny,jbgnl,jendl, 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
IF(mp_opt > 0) THEN
CALL mpsendrecv2dew(tem9,nx,ny,nz,1,1,0,tem4)
CALL mpsendrecv2dns(tem9,nx,ny,nz,1,1,0,tem5)
END IF
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:24),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,ipvplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:18),time0, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,trnplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,wetcanplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl,&
label(1:28),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,raincplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl,&
label(1:30),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,raingplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:22),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,raintplt)
! size of tem5 must be >= 6*nx*ny
sigplt(59) = 1
END IF
END IF
!
! Added Accumulated rainfall plots for priority 111, 112, 113
!
IF(ipriority(111) == ip) THEN
IF( (rainicplt == 1 .OR. rainicplt == 2 .OR. &
rainicplt == 4 .OR. rainicplt == 5 ) .AND. &
sigplt(111) == 0 ) THEN
CALL get_mulovrlay('rainicplt',8,ovrmul_num,ovrmulname, &
sovrlay)
DO j=1,ny
DO i=1,nx
acc_rainc(i,j,2)=rainc(i,j)
time2c=time
END DO
END DO
timediff=time2c-time1c ! time in sec
IF(raicunit == 0) THEN
WRITE(label,'(f8.1,a)') timediff, &
' s Accumulated Cumulus Rainfall (mm)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = acc_rainc(i,j,2)-acc_rainc(i,j,1)
IF (tem9(i,j,1) < 0.001) tem9(i,j,1)=0.0
END DO
END DO
ELSE IF( raicunit == 1) THEN ! unit is inch
WRITE(label,'(f8.1,a)') timediff, &
' s Accumulated Cumulus Rainfall (in)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1)=(acc_rainc(i,j,2)-acc_rainc(i,j,1)) &
* 0.039370079 ! (1/25.4)
IF (tem9(i,j,1) < 0.0001) tem9(i,j,1)=0.0
END DO
END DO
END IF
CALL ctrsetup(rainicinc,rainicminc,rainicmaxc, &
raicovr,raichlf,raiczro,raiccol1,raiccol2, &
'rainicplt ')
CALL ctrsfc(tem9(1,1,1),xc(1,1,2),yc(1,1,2),x1,x2,dxkm, &
y1,y2,dykm, nx,ibgn,iend, ny,jbgnl,jendl, &
TRIM(label),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,rainicplt)
! size of tem5 must be >= 6*nx*ny
sigplt(111) = 1
DO j=1,ny
DO i=1,nx
acc_rainc(i,j,1)=rainc(i,j)
END DO
END DO
time1c=time
END IF
END IF
IF(ipriority(112) == ip) THEN
IF( (rainigplt == 1 .OR. rainigplt == 2 .OR. &
rainigplt == 4 .OR. rainigplt == 5 ) .AND. &
sigplt(112) == 0 ) THEN
CALL get_mulovrlay('rainigplt',8,ovrmul_num,ovrmulname, &
sovrlay)
DO j=1,ny
DO i=1,nx
acc_raing(i,j,2)=raing(i,j)
time2g=time
END DO
END DO
timediff=time2g-time1g ! time in sec
IF(raigunit == 0) THEN
WRITE(label,'(f8.1,a)') timediff, &
' s Accumulated Gridscale Rainfall (mm)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = acc_raing(i,j,2)-acc_raing(i,j,1)
IF (tem9(i,j,1) < 0.001) tem9(i,j,1)=0.0
END DO
END DO
ELSE IF( raigunit == 1) THEN ! unit is inch
WRITE(label,'(f8.1,a)') timediff, &
' s Accumulated Gridscale Rainfall (in)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1)=(acc_raing(i,j,2)-acc_raing(i,j,1)) &
* 0.039370079 ! (1/25.4)
IF (tem9(i,j,1) < 0.0001) tem9(i,j,1)=0.0
END DO
END DO
END IF
CALL ctrsetup(rainiginc,rainigminc,rainigmaxc, &
raigovr,raighlf,raigzro,raigcol1,raigcol2, &
'rainigplt ')
CALL ctrsfc(tem9(1,1,1),xc(1,1,2),yc(1,1,2),x1,x2,dxkm, &
y1,y2,dykm, nx,ibgn,iend, ny,jbgnl,jendl, &
TRIM(label),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,rainigplt)
! size of tem5 must be >= 6*nx*ny
sigplt(112) = 1
DO j=1,ny
DO i=1,nx
acc_raing(i,j,1)=raing(i,j)
END DO
END DO
time1g=time
END IF
END IF
IF(ipriority(113) == ip) THEN
IF( (rainitplt == 1 .OR. rainitplt == 2 .OR. &
rainitplt == 4 .OR. rainitplt == 5 ) .AND. &
sigplt(113) == 0 ) THEN
CALL get_mulovrlay('rainitplt',8,ovrmul_num,ovrmulname, &
sovrlay)
DO j=1,ny
DO i=1,nx
acc_raint(i,j,2)=rainc(i,j) + raing(i,j)
time2t=time
END DO
END DO
timediff=time2t-time1t ! time in sec
IF(raitunit == 0) THEN
WRITE(label,'(f8.1,a)') timediff, &
' s Accumulated Rainfall (mm)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1) = acc_raint(i,j,2)-acc_raint(i,j,1)
IF (tem9(i,j,1) < 0.001) tem9(i,j,1)=0.0
END DO
END DO
ELSE IF( raitunit == 1) THEN ! unit is inch
WRITE(label,'(f8.1,a)') timediff, &
' s Accumulated Rainfall (in)'
DO j = 1,ny
DO i=1,nx
tem9(i,j,1)=(acc_raint(i,j,2)-acc_raint(i,j,1)) &
* 0.039370079 ! (1/25.4)
IF (tem9(i,j,1) < 0.0001) tem9(i,j,1)=0.0
END DO
END DO
END IF
CALL ctrsetup(rainitinc,rainitminc,rainitmaxc, &
raitovr,raithlf,raitzro,raitcol1,raitcol2, &
'rainitplt ')
CALL ctrsfc(tem9(1,1,1),xc(1,1,2),yc(1,1,2),x1,x2,dxkm, &
y1,y2,dykm, nx,ibgn,iend, ny,jbgnl,jendl, &
TRIM(label),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,rainitplt)
! size of tem5 must be >= 6*nx*ny
sigplt(113) = 1
DO j=1,ny
DO i=1,nx
acc_raint(i,j,1)=rainc(i,j) + raing(i,j)
END DO
END DO
time1t=time
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:23),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,pslplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,liplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,capeplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:18),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,cinplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,thetplt)
! size of tem5 must be >= 6*nx*ny
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)'
! IF(mp_opt >0) THEN
! CALL mpsend1dew(heli,nx,ny,1,1,0,mptag1,tem5)
! CALL mpsend1dns(heli,nx,ny,1,1,0,mptag2,tem5)
!
! CALL mprecv1dew(heli,nx,ny,1,1,0,mptag1,tem5)
! CALL mprecv1dns(heli,nx,ny,1,1,0,mptag2,tem5)
! END IF
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:24),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,heliplt)
! size of tem5 must be >= 6*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:19),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,ctcplt)
! size of tem5 must be >= 6*nx*ny
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)'
! IF(mp_opt >0) THEN
! CALL mpsend1dew(brn,nx,ny,1,1,0,mptag1,tem5)
! CALL mpsend1dns(brn,nx,ny,1,1,0,mptag2,tem5)
!
! CALL mprecv1dew(brn,nx,ny,1,1,0,mptag1,tem5)
! CALL mprecv1dns(brn,nx,ny,1,1,0,mptag2,tem5)
! END IF
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,ibgnl,iendl, ny,jbgnl,jendl, &
label(1:23),time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,slicopt,brnplt)
! size of tem5 must be >= 6*nx*ny
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)'
! IF(mp_opt >0) THEN
! CALL mpsend1dew(brnu,nx,ny,1,1,0,mptag1,tem5)
! CALL mpsend1dns(brnu,nx,ny,1,1,0,mptag2,tem5)
!
! CALL mprecv1dew(brnu,nx,ny,1,1,0,mptag1,tem5)
! CALL mprecv1dns(brnu,nx,ny,1,1,0,mptag2,tem5)
! END IF
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,ibgnl,iendl, ny,jbgnl,jendl, &
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)'
! IF(mp_opt >0) THEN
! CALL mpsend1dew(srlfl,nx,ny,1,1,0,mptag1,tem5)
! CALL mpsend1dns(srlfl,nx,ny,1,1,0,mptag2,tem5)
!
! CALL mprecv1dew(srlfl,nx,ny,1,1,0,mptag1,tem5)
! CALL mprecv1dns(srlfl,nx,ny,1,1,0,mptag2,tem5)
! END IF
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,ibgnl,iendl, ny,jbgnl,jendl, &
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)'
! IF(mp_opt >0) THEN
! CALL mpsend1dew(srmfl,nx,ny,1,1,0,mptag1,tem5)
! CALL mpsend1dns(srmfl,nx,ny,1,1,0,mptag2,tem5)
!
! CALL mprecv1dew(srmfl,nx,ny,1,1,0,mptag1,tem5)
! CALL mprecv1dns(srmfl,nx,ny,1,1,0,mptag2,tem5)
! END IF
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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)'
! IF(mp_opt >0) THEN
! CALL mpsend1dew(blcon,nx,ny,1,1,0,mptag1,tem5)
! CALL mpsend1dns(blcon,nx,ny,1,1,0,mptag2,tem5)
!
! CALL mprecv1dew(blcon,nx,ny,1,1,0,mptag1,tem5)
! CALL mprecv1dns(blcon,nx,ny,1,1,0,mptag2,tem5)
! END IF
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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 )
! IF(mp_opt >0) THEN
! CALL mpsend1dew(ustrm,nx,ny,1,1,0,mptag1,tem4)
! CALL mpsend1dns(ustrm,nx,ny,1,1,0,mptag2,tem5)
! CALL mprecv1dew(ustrm,nx,ny,1,1,0,mptag1,tem4)
! CALL mprecv1dns(ustrm,nx,ny,1,1,0,mptag2,tem5)
! CALL mpsend1dew(vstrm,nx,ny,1,1,0,mptag3,tem2)
! CALL mpsend1dns(vstrm,nx,ny,1,1,0,mptag4,tem3)
! CALL mprecv1dew(vstrm,nx,ny,1,1,0,mptag3,tem2)
! CALL mprecv1dns(vstrm,nx,ny,1,1,0,mptag4,tem3)
! END IF
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,ibgnl,iendl,ist, &
ny,jbgnl,jendl,jst, &
label(1:length),time, runname,1.0,slicopt, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
! size of tem6 must be >= 5*nx*ny
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,ibgnl,iendl, ny,jbgnl,jendl, &
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,'' (index)'')') soiltpn
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,ibgnl,iendl, ny,jbgnl,jendl, &
TRIM(label),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 (index)'
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,ibgnl,iendl, ny,jbgnl,jendl, &
TRIM(label),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 (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,ibgnl,iendl, ny,jbgnl,jendl, &
TRIM(label),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 (nounit)'
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,ibgnl,iendl, ny,jbgnl,jendl, &
TRIM(label),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 (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,ibgnl,iendl, ny,jbgnl,jendl, &
TRIM(label),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,ibgnl,iendl, ny,jbgnl,jendl, &
TRIM(label),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,ibgnl,iendl,ist, &
ny,jbgnl,jendl,jst, &
label(1:length),time, runname,1.0,slicopt, &
tem1,tem2,tem3,tem4,tem5,tem6,hterain)
! size of tem6 must be >= 5*nx*ny
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)
CALL readvar2(nx,ny,nz, var3dv, var3d(i),var_name, &
varunits, time, runname, dirname3d(i), &
finfmt3d(i), mp_opt*readsplit, istatus)
label=var3d(i)(1:6)//'('//varunits//')'
CALL ctr3d(var3dv,xc,yc,zps3d, &
x1,x2,dxkm,y1,y2,dykm,z1,z2,dzkm, &
nx,ibgnl,iendl, ny,jbgnl,jendl, nz,kbgn,kend, &
label,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)
CALL readvar2(nx,ny,1, var2dv, var2d(i),var_name, &
varunits, time, runname, dirname2d(i), &
finfmt2d(i), mp_opt*readsplit, istatus)
label=var2d(i)(1:6)//'('//varunits//')'
CALL ctrsfc(var2dv,xc(1,1,2),yc(1,1,2), &
x1,x2,dxkm,y1,y2,dykm, &
nx,ibgnl,iendl, ny,jbgnl,jendl, &
label,time, runname,1.0 ,tem1,tem2,tem3, &
tem4,tem5,hterain,1,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,ibgnl,iendl,ny,jbgnl,jendl,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,ibgnl,iendl, ny,jbgnl,jendl,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
! Plotting trajectories (Dan Dawson 12/03/2004) -should probably be
! plotted in priority loop like everything else
IF (trajopt == 1) THEN
DO k=1,ntimes
print*,'k = ',k
!-----------------------------------------------------------------------
! Read trajectory data
!-----------------------------------------------------------------------
!print*,trajc_fn_in(k)
IF(curtim >= tstart_calc .and. curtim < tend_calc) THEN
IF(trajfileflag == 0) THEN
trajfileflag = 1
CALL getunit (nunit(k))
write(6,'(a)') 'Opening ',trim(trajc_fn_in(k))
OPEN(UNIT=nunit(k),FILE=trim(trajc_fn_in(k)),STATUS='old', &
FORM='formatted',IOSTAT= istat )
IF(istat == 0) THEN
print*,'Trajectory file ',trim(trajc_fn_in(k)),' successfully read.'
READ(nunit(k),'(a)') runname
READ(nunit(k),'(6e17.6)') xlow, xhigh, ylow, yhigh, zlow, zhigh
write(6,'(6e17.6)') xlow, xhigh, ylow, yhigh, zlow, zhigh
READ(nunit(k),'(3e17.6)') dx, dy, dz
write(6,'(3e17.6)') dx, dy, dz
READ(nunit(k),'(3e17.6)') tstart, tzero, tend
READ(nunit(k),'(i10)') npoints
READ(nunit(k),'(i10)') ntrajcs
DO j=1,npoints
READ(nunit(k),'(4e17.6)',err=115,end=115) ttrajc(j)
READ(nunit(k),'(i10)',err=115,end=115) ntrajcs_in
IF( ntrajcs_in /= ntrajcs ) then
print*,'ntrajcs read in .ne. ntrajcs in program.'
print*,'Job stopped'
STOP
ENDIF
READ(nunit(k),'(6e17.6)',err=115,end=115) ((xtrajc(j,i,k),ytrajc(j,i,k),ztrajc(j,i,k)),i=1,ntrajcs)
ENDDO
CLOSE(UNIT=nunit(k))
CALL retunit(nunit(k))
ELSE
CALL retunit(nunit(k))
GOTO 135
END IF
END IF
npoints_in = npoints
GOTO 125
115 continue
npoints_in = max(1,j-1)
125 continue
! Plot the trajectories up to the current time
CALL xthick(2)
IF(labelopt == 1) THEN
CALL XQCHSZ(chsize)
CALL XCHSIZ(chsize*labelmag)
END IF
npoints_cur = INT((curtim-tstart)/tinc_calc) + 1
IF(npoints_cur >= npoints) THEN
npoints_cur = npoints
END IF
DO i=1,ntrajcs
CALL xcolor(traj_col(i))
print*,'Plotting trajectory ',i
IF(slicopt == 1) THEN
pen_status = 0
DO j=1,npoints_cur
IF( ( abs( xtrajc(j,i,k) + 99999.0 ) .lt. 0.001 .or. &
abs( ytrajc(j,i,k) + 99999.0 ) .lt. 0.001 ) ) then
pen_status = 0
ELSE
IF( pen_status == 0 ) then
CALL XPENUP(xtrajc(j,i,k)/1000.,ytrajc(j,i,k)/1000.)
ELSE
CALL XPENDN(xtrajc(j,i,k)/1000.,ytrajc(j,i,k)/1000.)
ENDIF
pen_status = 1
IF(labelopt == 1) THEN
IF(MOD(j,labelfreq) == 0) THEN
itrajc1 = MAX(1, MIN(nx-2, INT((xtrajc(j,i,k)-xs(1))/(xs(2)-xs(1)))+1 ))
jtrajc1 = MAX(1, MIN(ny-2, INT((ytrajc(j,i,k)-ys(1))/(ys(2)-ys(1)))+1 ))
CALL XRNUMB(xtrajc(j,i,k)/1000.,ytrajc(j,i,k)/1000., 0.001*(ztrajc(j,i,k)-zp(itrajc1,jtrajc1,2)),'*')
! pen_status = 0
CALL XPENUP(xtrajc(j,i,k)/1000.,ytrajc(j,i,k)/1000.)
END IF
END IF
Endif
END DO
ELSE IF(slicopt == 2) THEN
pen_status = 0
DO j=1,npoints_cur
IF( ( abs( xtrajc(j,i,k) + 99999.0 ) .lt. 0.001 .or. &
abs( ztrajc(j,i,k) + 99999.0 ) .lt. 0.001 ) ) then
pen_status = 0
ELSE
IF( pen_status == 0 ) then
CALL XPENUP(xtrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000.)
ELSE
CALL XPENDN(xtrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000.)
ENDIF
pen_status = 1
IF(labelopt == 1) THEN
IF(MOD(j,labelfreq) == 0) THEN
CALL XRNUMB(xtrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000., ytrajc(j,i,k)*0.001,'*')
! pen_status = 0
CALL XPENUP(xtrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000.)
END IF
END IF
ENDIF
END DO
ELSE IF(slicopt == 3) THEN
pen_status = 0
DO j=1,npoints_cur
IF( ( abs( ytrajc(j,i,k) + 99999.0 ) .lt. 0.001 .or. &
abs( ztrajc(j,i,k) + 99999.0 ) .lt. 0.001 ) ) then
pen_status = 0
ELSE
IF( pen_status == 0 ) then
CALL XPENUP(ytrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000.)
ELSE
CALL XPENDN(ytrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000.)
ENDIF
pen_status = 1
IF(labelopt == 1) THEN
IF(MOD(j,labelfreq) == 0) THEN
CALL XRNUMB(ytrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000., xtrajc(j,i,k)*0.001,'*')
! pen_status = 0
CALL XPENUP(ytrajc(j,i,k)/1000.,ztrajc(j,i,k)/1000.)
END IF
END IF
ENDIF
END DO
END IF
END DO ! i=1,ntrajcs
END IF ! curtim > tstart .and. curtim < tend
END DO ! k=1,ntimes
CALL XCHSIZ(chsize)
CALL xthick(1)
CALL xfull
END IF !trajopt == 1
GOTO 145
135 write(6,'(a)') 'Cannot read trajectory file, continuing'
145 continue
!
!-----------------------------------------------------------------------
!
! 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)
IF(myproc == 0) THEN
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
!
! Plot w=-wisosf isosurface.
!
CALL wirfrm(w,nx,ibgn,iend,ny,jbgn,jend,nz,kbgn,kend, &
-wisosf, tem2)
IF(myproc == 0) THEN
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
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)
IF(myproc == 0) THEN
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
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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) CALL runlab (runname)
END IF
IF (rfprof == 1) THEN
IF (rfopt == 2) THEN
CALL reflec_ferrier(nx,ny,nz, rhobar, qr, qs, qh, tz, tem5)
ELSE
CALL reflec(nx,ny,nz, rhobar, qr, qs, qh, tem5)
END IF
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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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)
IF(myproc == 0) 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
END DO ! data file loop
GO TO 800
700 CONTINUE
WRITE(6,'(a,/a,i3,a)') &
' Bad return status from data file reading', &
' ireturn = ', ireturn,' program stopped.'
STOP
800 CONTINUE
IF(myproc ==0 ) THEN
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
END IF
IF (mp_opt > 0) CALL mpexit(0)
STOP
END PROGRAM ARPSPLT