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