PROGRAM arpsenscv,110
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM ARPSENSCV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
! PURPOSE:
! Read in a series of ARPS history files, generate 2-D pruducts data
! for display. The output may be in the same grid as the input
! data, or can be different from it.
! The output data is on a qni x qnj grid with its southwest corner
! at qswcorn(lat degree),qswcorw(lon degree) specified in enscv.inc,
! where qnecorn and qnecorw are also set. However, the mapproj,
! trulon, trulat(1),trulat(2) parameters should be specified in this
! program. Please check the NEWMAP subroutine.
!
! Author : Dingchen Hou
! History: Apr. 30 1998, developmed from the framework of ARPSPLT.
! Sep. 15, 1999, modified to include perturbations in soil
! variables.
! Feb 2002 (F. KONG), modified to read nx,ny,nz directly from
! input file; add 'iaccu' into output_data namelist; & fix bugs
!-----------------------------------------------------------------------
!
! MODIFIED:
!
! 05/28/2002 (J. Brotzge)
! Added tsoil/qsoil to accomodate new soil scheme.
!
! 1 June 2002 Eric Kemp
! Soil variable updates.
!
!-----------------------------------------------------------------------
! DATA ARRAYS READ IN:
!
! x x-coordinate of grid points in physical/comp. space (m)
! y y-coordinate of grid points in physical/comp. space (m)
! z z-coordinate of grid points in computational space (km)
! zp z-coordinate of grid points in computational space (m)
! zpsoil z-coordinate of grid points in computational space (m)
!
! uprt x-component of perturbation velocity (m/s)
! vprt y-component of perturbation velocity (m/s)
! wprt vertical component of perturbation velocity in Cartesian
! coordinates (m/s).
!
! ptprt perturbation potential temperature (K)
! pprt perturbation pressure (Pascal)
!
! qvprt perturbation water vapor mixing ratio (kg/kg)
! qc Cloud water mixing ratio (kg/kg)
! qr Rainwater mixing ratio (kg/kg)
! qi Cloud ice mixing ratio (kg/kg)
! qs Snow mixing ratio (kg/kg)
! qh Hail mixing ratio (kg/kg)
! tke Turbulent Kinetic Energy ((m/s)**2)
! kmh Horizontal turb. mixing coef. for momentum ( m**2/s )
! kmv Vertical turb. mixing coef. for momentum ( m**2/s )
!
! ubar Base state x-velocity component (m/s)
! vbar Base state y-velocity component (m/s)
! wbar Base state z-velocity component (m/s)
! ptbar Base state potential temperature (K)
! pbar Base state pressure (Pascal)
! rhobar Base state air density (kg/m**3)
! qvbar Base state water vapor mixing ratio (kg/kg)
!
! soiltyp Soil type
! stypfrct Soil type fraction
! vegtyp Vegetation type
! lai Leaf Area Index
! roufns Surface roughness
! veg Vegetation fraction
!
! tsoil Soil temperature (K)
! qsoil Soil moisture (m**3/m**3)
! wetcanp Canopy water amount
! raing Grid supersaturation rain
! rainc Cumulus convective rain(mm)
! raing Total rain (rainc+raing)(mm)
!
! 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)
! cape1 CAPE (J/kg)
! cin1 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
!
!-----------------------------------------------------------------------
!
! hgt z-coor of scalar point in physical space (10m)
! 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
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
INCLUDE 'indtflg.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.inc'
INCLUDE 'enscv.inc'
!
!-----------------------------------------------------------------------
!
! Some universal constants
!
!-----------------------------------------------------------------------
!
REAL*4 kappa,gamma,ex1,ex2,t00,p00,mbtopa
PARAMETER (kappa=287.053/cp, &
gamma=6.5, & ! 6.5 K/km
ex1=0.1903643, & ! R*gamma/g
ex2=5.2558774, & ! g/R/gamma
mbtopa=100.)
!-----------------------------------------------------------------------
!
! Define the dimensions nx, ny, nz, nzsoil
!
!-----------------------------------------------------------------------
INTEGER :: nx, ny, nz, nzsoil
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x (:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, ALLOCATABLE :: y (:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, ALLOCATABLE :: z (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, ALLOCATABLE :: zp (:,:,:) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, ALLOCATABLE :: zpsoil(:,:,:) ! The physical height coordinate defined at
! w-point of the soil grid.
REAL, ALLOCATABLE :: u (:,:,:) ! Total u-velocity (m/s)
REAL, ALLOCATABLE :: v (:,:,:) ! Total v-velocity (m/s)
REAL, ALLOCATABLE :: w (:,:,:) ! Total w-velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature
! from that of base state atmosphere (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure from that
! of base state atmosphere (Pascal)
REAL, ALLOCATABLE :: qv (:,:,:) ! Water vapor specific humidity (kg/kg)
REAL, ALLOCATABLE :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qh (:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: tke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, ALLOCATABLE :: kmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: wbar (:,:,:) ! Base state w-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: rhobar(:,:,:) ! Base state density rhobar
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific humidity
! (kg/kg)
INTEGER :: nstyps ! Number of soil type
! PARAMETER ( nstyps = 4 )
INTEGER, ALLOCATABLE :: soiltyp(:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type fraction
INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! Soil temperature (K)
REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! Soil moisture (m**3/m**3)
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: raint(:,:) ! Total rain (rainc+raing)
REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface
REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface
REAL, ALLOCATABLE :: radswnet(:,:) ! Net shortwave radiation, SWin - SWout
REAL, ALLOCATABLE :: radlwin(:,:) ! Incoming longwave radiation
REAL, ALLOCATABLE :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: ptsflx(:,:) ! Surface heat flux (K*kg/(m*s**2))
REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
! Arrays derived from the read-in arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: uprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: wprt (:,:,:) ! Perturbation w-velocity (m/s)
REAL, ALLOCATABLE :: pt (:,:,:) ! Total poten
REAL, ALLOCATABLE :: qvprt (:,:,:)
REAL, ALLOCATABLE :: xc (:,:,:) ! x-coor of sacalar point (km)
REAL, ALLOCATABLE :: yc (:,:,:) ! y-coor of sacalar point (km)
REAL, ALLOCATABLE :: zc (:,:,:) ! z-coor of sacalar point in computational
! space (km)
REAL, ALLOCATABLE :: zpc (:,:,:) ! z-coor of sacalar point in physical
! space (km)
REAL, ALLOCATABLE :: hterain(:,:) ! The height of the terrain.
!-----------------------------------------------------------------------
!
! Arrays for plots on constant pressure levels
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: hgt (:,:,:) ! z-coor of sacalar point in physical
! space (10m)
REAL, ALLOCATABLE :: t700 (:,:) ! Temperature (K) on 700mb pressure grids
REAL, ALLOCATABLE :: qs700 (:,:) ! Temperature (K) on 700mb pressure grids
REAL, ALLOCATABLE :: algpzc(:,:,:) ! -log(pressure) at scalar grid points
!-----------------------------------------------------------------------
!
! 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 :: li1(:,:) ! Lifted Index (K)
REAL, ALLOCATABLE :: cape1(:,:) ! CAPE (J/kg)
REAL, ALLOCATABLE :: cin1(:,:) ! CIN (J/kg)
REAL, ALLOCATABLE :: thet(:,:) ! theta_E (K)
REAL, ALLOCATABLE :: heli(:,:) ! helicity
REAL, ALLOCATABLE :: brn1(:,:) ! Bulk Richardson Number (Weisman and Klemp)
REAL, ALLOCATABLE :: brnu(:,:) ! Shear parameter of BRN, "U"
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 :: oe ! function for caculate theta_e
REAL :: alttostpr
REAL, ALLOCATABLE :: sinlat(:,:) ! Sin of latitude at each grid point
SAVE sinlat
REAL :: lat1,lon1,lat2,lon2
REAL :: pw1d,totw
REAL :: pi
!
!-----------------------------------------------------------------------
! OUTPUT arrays
!-----------------------------------------------------------------------
INTEGER :: fileun
REAL :: resltn
INTEGER :: houroday, minohour, dayoweek, dateomon, month1, year1
REAL, ALLOCATABLE :: mslp(:,:)
REAL, ALLOCATABLE :: hgt250(:,:), vort250(:,:), uwind250(:,:), vwind250(:,:)
REAL, ALLOCATABLE :: hgt500(:,:), vort500(:,:), uwind500(:,:), vwind500(:,:)
REAL, ALLOCATABLE :: hgt850(:,:), vort850(:,:), uwind850(:,:), vwind850(:,:)
REAL, ALLOCATABLE :: temp850(:,:)
REAL, ALLOCATABLE :: thck(:,:), rh700(:,:), omega700(:,:),temp2m(:,:), dewp2m(:,:)
REAL, ALLOCATABLE :: accppt(:,:), convppt(:,:)
REAL, ALLOCATABLE :: tpptold(:,:), cpptold(:,:)
REAL, ALLOCATABLE :: sreh(:,:),li(:,:),cape(:,:),brn(:,:),cin(:,:),llmc(:,:),pw(:,:)
REAL, ALLOCATABLE :: cmpref(:,:)
INTEGER :: ni,nj
REAL, ALLOCATABLE :: f2(:,:),x2(:,:),y2(:,:)
REAL, ALLOCATABLE :: x1(:,:),y1(:,:)
!
!-----------------------------------------------------------------------
! Output Grid
!-----------------------------------------------------------------------
!
INTEGER :: max2d_out2
! PARAMETER (MAX2D_OUT2=max((nx-3)*(ny-3),MAX2D_OUT))
PARAMETER (max2d_out2=max2d_out)
REAL :: mslpn(max2d_out2)
REAL :: hgt250n(max2d_out2), vort250n(max2d_out2), &
uwind250n(max2d_out2), vwind250n(max2d_out2)
REAL :: hgt500n(max2d_out2), vort500n(max2d_out2), &
uwind500n(max2d_out2), vwind500n(max2d_out2)
REAL :: hgt850n(max2d_out2), vort850n(max2d_out2), &
uwind850n(max2d_out2), vwind850n(max2d_out2)
REAL :: temp850n(max2d_out2)
REAL :: thckn(max2d_out2), rh700n(max2d_out2), omega700n(max2d_out2), &
temp2mn(max2d_out2), dewp2mn(max2d_out2)
REAL :: accpptn(max2d_out2), convpptn(max2d_out2)
REAL :: srehn(max2d_out2), lin(max2d_out2), capen(max2d_out2), &
brnn(max2d_out2), &
cinn(max2d_out2), llmcn(max2d_out2), pwn(max2d_out2)
REAL :: cmprefn(max2d_out2)
REAL :: xn(max2d_out2),yn(max2d_out2)
REAL :: x1n(max2d_out2),y1n(max2d_out2)
REAL :: latn(max2d_out2),lonn(max2d_out2)
!
!-----------------------------------------------------------------------
!
! Temporary work arrays for general use
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
REAL, ALLOCATABLE :: tem4(:,:,:)
REAL, ALLOCATABLE :: tem5(:,:,:)
REAL, ALLOCATABLE :: tem6(:,:,:)
INTEGER, ALLOCATABLE :: item(:,:)
REAL, ALLOCATABLE :: var2(:,:)
REAL, ALLOCATABLE :: pt2m(:,:),qv2m(:,:)
REAL, ALLOCATABLE :: vtwo1(:,:),vtwo2(:,:),vtwo3(:,:)
REAL, ALLOCATABLE :: vtwo4(:,:),vtwo5(:,:),vtwo6(:,:)
REAL, ALLOCATABLE :: qvsfc (:,:,:) ! Soil Skin effective qv
REAL, ALLOCATABLE :: psf(:,:),psl(:,:)
INTEGER :: ijq
INTEGER, ALLOCATABLE :: item1(:,:),item2(:,:),item3(:,:),item4(:,:), &
item5(:,:),item6(:,:),item7(:,:)
REAL, ALLOCATABLE :: rtem1(:,:),rtem2(:,:),rtem3(:,:),rtem4(:,:), &
rtem5(:,:),rtem6(:,:),rtem7(:,:)
!
!-----------------------------------------------------------------------
!
! Work arrays to be used in interpolation subroutines
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: b1(:,:)
REAL, ALLOCATABLE :: t1(:,:),t2(:,:),t3(:,:),t4(:,:)
!
!-----------------------------------------------------------------------
!
! Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
REAL :: z01 ! altitude (meters) of a horizontal
! slice so specified
COMMON /sliceh/z01
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
REAL :: swx,ctrx,swy,ctry,tem
REAL :: ctrx1, ctry1
REAL :: iorig
REAL :: rho
!
INTEGER,ALLOCATABLE :: slice(:,:)
REAL,ALLOCATABLE :: xi1(:), yi1(:),xi2(:),yi2(:),zi1(:),ppil(:)
CHARACTER (LEN=80) :: label
CHARACTER (LEN=80) :: filename
CHARACTER (LEN=80) :: grdbasfn
CHARACTER (LEN=80) :: hisfile(25)
CHARACTER (LEN=80) :: outfile(25)
INTEGER :: hinfmt,nchin
INTEGER :: nslice(7), indxslic
INTEGER :: ireturn
INTEGER :: mode,lengbf,lenfil
INTEGER :: i,j,k,ibgn,iend,jbgn,jend,kbgn,kend,ist,jst,kst,iob
INTEGER :: ibgn1, iend1
INTEGER :: nxpic,nypic,islice,jslice,kslice,layout
REAL :: ctime,angl,aspect
! real xr,yr,zr,x1,x2,y1,y2,z1,z2,xlimit,ylimit
REAL :: zmax,zmin
REAL :: factor,pmb,drot
INTEGER :: imax,jmax,kmax,imin,jmin,kmin
REAL :: wmin, wmax, xorold, yorold, tk, tdk, psta
REAL :: xbgn,xend,ybgn,yend,zbgn,zend
LOGICAL :: fexist
INTEGER :: onvf
INTEGER :: nf
INTEGER :: lenstag, linstit, lmodel
!
!-----------------------------------------------------------------------
INTEGER :: intrpl ! flag indicating whether to interpolate output
REAL :: truelat(2)
REAL :: trltnew(2),trlnnew,dxnew
INTEGER :: mptjnew
INTEGER :: nhisfile
INTEGER :: iout(25),icape,iplot,iterr,iaccu
NAMELIST /history_data/ hinfmt, nhisfile, grdbasfn,hisfile
NAMELIST /output_data/ outfile,iout,icape,iplot,iterr,iaccu
NAMELIST /output_grid/ mptjnew,trltnew,trlnnew,dxnew, &
qni,qnj,qswcorn,qswcorw,qnecorn,qnecorw,intrpl, &
enstag,instit,model
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directives for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
WRITE(6,'(/ 16(/5x,a)//)') &
'###############################################################', &
'###############################################################', &
'#### ####', &
'#### Welcome to ARPSENSCV ####', &
'#### ####', &
'#### A graphic analysis program for model ARPS 4.5 ####', &
'#### ####', &
'#### Program version 4.5 ####', &
'#### ####', &
'#### The graphic plotting is based ####', &
'#### on graphic package ZXPLOT ####', &
'#### by Ming Xue CAPS/OU ####', &
'#### ####', &
'###############################################################', &
'###############################################################'
!-----------------------------------------------------------------------
!
! Read grid/base name, then get the dimensions
!
!-----------------------------------------------------------------------
READ(5,history_data,END=100)
WRITE(6,'(/a/)')'Namelist block history_data sucessfully read in.'
IF ( hinfmt == 9 ) GO TO 10
lengbf = len_trim(grdbasfn)
WRITE(6,'(a,a/)')' The grid/base name is ', grdbasfn(1:lengbf)
DO nf=1,nhisfile
lenfil = len_trim(hisfile(nf))
WRITE(6,'(1x,a,a)') 'Input values for hisfile: ',trim(hisfile(nf))
END DO
GO TO 10
100 WRITE(6,'(a)') &
'Error reading NAMELIST block history_data. STOP'
STOP
10 CONTINUE
!
CALL get_dims_from_data
(hinfmt,trim(grdbasfn), &
nx,ny,nz,nzsoil,nstyps, ireturn)
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
STOP
END IF
ni=nx-1 ! F. KONG add to fix a bug for zero ni,nj
nj=ny-1
WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,' nzsoil =',nzsoil
print*,'nstyps =', nstyps
!-----------------------------------------------------------------------
!
! Allocate nx, ny, nz dependent arrays and initiliaze to zero
!
!-----------------------------------------------------------------------
ALLOCATE( x=0
ALLOCATE( y=0
ALLOCATE( z=0
ALLOCATE(zp(nx,ny,nz),STAT=istatus)
zp=0
ALLOCATE(zpsoil(nx,ny,nzsoil),STAT=istatus)
zpsoil=0
ALLOCATE( u=0
ALLOCATE( v=0
ALLOCATE( w=0
ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
ptprt=0
ALLOCATE(pprt(nx,ny,nz),STAT=istatus)
pprt=0
ALLOCATE(qv(nx,ny,nz),STAT=istatus)
qv=0
ALLOCATE(qc(nx,ny,nz),STAT=istatus)
qc=0
ALLOCATE(qr(nx,ny,nz),STAT=istatus)
qr=0
ALLOCATE(qi(nx,ny,nz),STAT=istatus)
qi=0
ALLOCATE(qs(nx,ny,nz),STAT=istatus)
qs=0
ALLOCATE(qh(nx,ny,nz),STAT=istatus)
qh=0
ALLOCATE(tke(nx,ny,nz),STAT=istatus)
tke=0
ALLOCATE(kmh(nx,ny,nz),STAT=istatus)
kmh=0
ALLOCATE(kmv(nx,ny,nz),STAT=istatus)
kmv=0
ALLOCATE(ubar(nx,ny,nz),STAT=istatus)
ubar=0
ALLOCATE(vbar(nx,ny,nz),STAT=istatus)
vbar=0
ALLOCATE(wbar(nx,ny,nz),STAT=istatus)
wbar=0
ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
ptbar=0
ALLOCATE(pbar(nx,ny,nz),STAT=istatus)
pbar=0
ALLOCATE(rhobar(nx,ny,nz),STAT=istatus)
rhobar=0
ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
qvbar=0
ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus)
soiltyp=0
ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
stypfrct=0
ALLOCATE(vegtyp(nx,ny),STAT=istatus)
vegtyp=0
ALLOCATE(lai(nx,ny),STAT=istatus)
lai=0
ALLOCATE(roufns(nx,ny),STAT=istatus)
roufns=0
ALLOCATE(veg(nx,ny),STAT=istatus)
veg=0
ALLOCATE(tsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus)
tsoil=0
ALLOCATE(qsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus)
qsoil=0
ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus)
wetcanp=0
ALLOCATE(snowdpth(nx,ny),STAT=istatus)
snowdpth=0
ALLOCATE(raing(nx,ny),STAT=istatus)
raing=0
ALLOCATE(rainc(nx,ny),STAT=istatus)
rainc=0
ALLOCATE(raint(nx,ny),STAT=istatus)
raint=0
ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
prcrate=0
ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
radfrc=0
ALLOCATE(radsw(nx,ny),STAT=istatus)
radsw=0
ALLOCATE(rnflx(nx,ny),STAT=istatus)
rnflx=0
ALLOCATE(radswnet(nx,ny),STAT=istatus)
radswnet=0
ALLOCATE(radlwin(nx,ny),STAT=istatus)
radlwin=0
ALLOCATE(usflx(nx,ny),STAT=istatus)
usflx=0
ALLOCATE(vsflx(nx,ny),STAT=istatus)
vsflx=0
ALLOCATE(ptsflx(nx,ny),STAT=istatus)
ptsflx=0
ALLOCATE(qvsflx(nx,ny),STAT=istatus)
qvsflx=0
ALLOCATE(uprt(nx,ny,nz),STAT=istatus)
uprt=0
ALLOCATE(vprt(nx,ny,nz),STAT=istatus)
vprt=0
ALLOCATE(wprt(nx,ny,nz),STAT=istatus)
wprt=0
ALLOCATE(pt(nx,ny,nz),STAT=istatus)
pt=0
ALLOCATE(qvprt(nx,ny,nz),STAT=istatus)
qvprt=0
ALLOCATE(xc(nx,ny,nz),STAT=istatus)
xc=0
ALLOCATE(yc(nx,ny,nz),STAT=istatus)
yc=0
ALLOCATE(zc(nx,ny,nz),STAT=istatus)
zc=0
ALLOCATE(zpc(nx,ny,nz),STAT=istatus)
zpc=0
ALLOCATE(hterain(nx,ny),STAT=istatus)
hterain=0
ALLOCATE(hgt(nx,ny,nz),STAT=istatus)
hgt=0
ALLOCATE(t700(nx,ny),STAT=istatus)
t700=0
ALLOCATE(qs700(nx,ny),STAT=istatus)
qs700=0
ALLOCATE(algpzc(nx,ny,nz),STAT=istatus)
algpzc=0
ALLOCATE(wrk1(nz),STAT=istatus)
wrk1=0
ALLOCATE(wrk2(nz),STAT=istatus)
wrk2=0
ALLOCATE(wrk3(nz),STAT=istatus)
wrk3=0
ALLOCATE(wrk4(nz),STAT=istatus)
wrk4=0
ALLOCATE(wrk5(nz),STAT=istatus)
wrk5=0
ALLOCATE(wrk6(nz),STAT=istatus)
wrk6=0
ALLOCATE(wrk7(nz),STAT=istatus)
wrk7=0
ALLOCATE(wrk8(nz),STAT=istatus)
wrk8=0
ALLOCATE(wrk9(nz),STAT=istatus)
wrk9=0
ALLOCATE(wrk10(nz),STAT=istatus)
wrk10=0
ALLOCATE(wrk11(nz),STAT=istatus)
wrk11=0
ALLOCATE(wrk12(nz),STAT=istatus)
wrk12=0
ALLOCATE(lcl(nx,ny),STAT=istatus)
lcl=0
ALLOCATE(lfc(nx,ny),STAT=istatus)
lfc=0
ALLOCATE(el(nx,ny),STAT=istatus)
el=0
ALLOCATE(twdf(nx,ny),STAT=istatus)
twdf=0
ALLOCATE(mbe(nx,ny),STAT=istatus)
mbe=0
ALLOCATE(li1(nx,ny),STAT=istatus)
li1=0
ALLOCATE(cape1(nx,ny),STAT=istatus)
cape1=0
ALLOCATE(cin1(nx,ny),STAT=istatus)
cin1=0
ALLOCATE(thet(nx,ny),STAT=istatus)
thet=0
ALLOCATE(heli(nx,ny),STAT=istatus)
heli=0
ALLOCATE(brn1(nx,ny),STAT=istatus)
brn1=0
ALLOCATE(brnu(nx,ny),STAT=istatus)
brnu=0
ALLOCATE(srlfl(nx,ny),STAT=istatus)
srlfl=0
ALLOCATE(srmfl(nx,ny),STAT=istatus)
srmfl=0
ALLOCATE(shr37(nx,ny),STAT=istatus)
shr37=0
ALLOCATE(ustrm(nx,ny),STAT=istatus)
ustrm=0
ALLOCATE(vstrm(nx,ny),STAT=istatus)
vstrm=0
ALLOCATE(capst(nx,ny),STAT=istatus)
capst=0
ALLOCATE(blcon(nx,ny),STAT=istatus)
blcon=0
ALLOCATE(ct(nx,ny),STAT=istatus)
ct=0
ALLOCATE(sinlat(nx,ny),STAT=istatus)
sinlat=0
ALLOCATE(mslp(nx-1,ny-1),STAT=istatus)
mslp=0
ALLOCATE(hgt250(nx-1,ny-1),STAT=istatus)
hgt250=0
ALLOCATE(vort250(nx-1,ny-1),STAT=istatus)
vort250=0
ALLOCATE(uwind250(nx-1,ny-1),STAT=istatus)
uwind250=0
ALLOCATE(vwind250(nx-1,ny-1),STAT=istatus)
vwind250=0
ALLOCATE(hgt500(nx-1,ny-1),STAT=istatus)
hgt500=0
ALLOCATE(vort500(nx-1,ny-1),STAT=istatus)
vort500=0
ALLOCATE(uwind500(nx-1,ny-1),STAT=istatus)
uwind500=0
ALLOCATE(vwind500(nx-1,ny-1),STAT=istatus)
vwind500=0
ALLOCATE(hgt850(nx-1,ny-1),STAT=istatus)
hgt850=0
ALLOCATE(vort850(nx-1,ny-1),STAT=istatus)
vort850=0
ALLOCATE(uwind850(nx-1,ny-1),STAT=istatus)
uwind850=0
ALLOCATE(vwind850(nx-1,ny-1),STAT=istatus)
vwind850=0
ALLOCATE(temp850(nx-1,ny-1),STAT=istatus)
temp850=0
ALLOCATE(thck(nx-1,ny-1),STAT=istatus)
thck=0
ALLOCATE(rh700(nx-1,ny-1),STAT=istatus)
rh700=0
ALLOCATE(omega700(nx-1,ny-1),STAT=istatus)
omega700=0
ALLOCATE(temp2m(nx-1,ny-1),STAT=istatus)
temp2m=0
ALLOCATE(dewp2m(nx-1,ny-1),STAT=istatus)
dewp2m=0
ALLOCATE(accppt(nx-1,ny-1),STAT=istatus)
accppt=0
ALLOCATE(convppt(nx-1,ny-1),STAT=istatus)
convppt=0
ALLOCATE(tpptold(nx-1,ny-1),STAT=istatus)
tpptold=0
ALLOCATE(cpptold(nx-1,ny-1),STAT=istatus)
cpptold=0
ALLOCATE(sreh(nx-1,ny-1),STAT=istatus)
sreh=0
ALLOCATE(li(nx-1,ny-1),STAT=istatus)
li=0
ALLOCATE(cape(nx-1,ny-1),STAT=istatus)
cape=0
ALLOCATE(brn(nx-1,ny-1),STAT=istatus)
brn=0
ALLOCATE(cin(nx-1,ny-1),STAT=istatus)
cin=0
ALLOCATE(llmc(nx-1,ny-1),STAT=istatus)
llmc=0
ALLOCATE(pw(nx-1,ny-1),STAT=istatus)
pw=0
ALLOCATE(cmpref(nx-1,ny-1),STAT=istatus)
cmpref=0
ALLOCATE(f2(nx-1,ny-1),STAT=istatus)
f2=0
ALLOCATE(x2(nx-1,ny-1),STAT=istatus)
x2=0
ALLOCATE(y2(nx-1,ny-1),STAT=istatus)
y2=0
ALLOCATE(x1(nx-1,ny-1),STAT=istatus)
x1=0
ALLOCATE(y1(nx-1,ny-1),STAT=istatus)
y1=0
ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
tem1=0
ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
tem2=0
ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
tem3=0
ALLOCATE(tem4(nx,ny,nz),STAT=istatus)
tem4=0
ALLOCATE(tem5(nx,ny,nz),STAT=istatus)
tem5=0
ALLOCATE(tem6(nx,ny,nz),STAT=istatus)
tem6=0
ALLOCATE(item(nx,ny),STAT=istatus)
ALLOCATE(var2(nx,ny),STAT=istatus)
var2=0
ALLOCATE(pt2m(nx,ny),STAT=istatus)
pt2m=0
ALLOCATE(qv2m(nx,ny),STAT=istatus)
qv2m=0
ALLOCATE(vtwo1(nx,ny),STAT=istatus)
vtwo1=0
ALLOCATE(vtwo2(nx,ny),STAT=istatus)
vtwo2=0
ALLOCATE(vtwo3(nx,ny),STAT=istatus)
vtwo3=0
ALLOCATE(vtwo4(nx,ny),STAT=istatus)
vtwo4=0
ALLOCATE(vtwo5(nx,ny),STAT=istatus)
vtwo5=0
ALLOCATE(vtwo6(nx,ny),STAT=istatus)
vtwo6=0
ALLOCATE(qvsfc(nx,ny,0:nstyps),STAT=istatus)
qvsfc=0
ALLOCATE(psf(nx,ny),STAT=istatus)
psf=0
ALLOCATE(psl(nx,ny),STAT=istatus)
psl=0
ALLOCATE(item1(nx,ny),STAT=istatus)
item1=0
ALLOCATE(item2(nx,ny),STAT=istatus)
item2=0
ALLOCATE(item3(nx,ny),STAT=istatus)
item3=0
ALLOCATE(item4(nx,ny),STAT=istatus)
item4=0
ALLOCATE(item5(nx,ny),STAT=istatus)
item5=0
ALLOCATE(item6(nx,ny),STAT=istatus)
item6=0
ALLOCATE(item7(nx,ny),STAT=istatus)
item7=0
ALLOCATE(rtem1(nx,ny),STAT=istatus)
rtem1=0
ALLOCATE(rtem2(nx,ny),STAT=istatus)
rtem2=0
ALLOCATE(rtem3(nx,ny),STAT=istatus)
rtem3=0
ALLOCATE(rtem4(nx,ny),STAT=istatus)
rtem4=0
ALLOCATE(rtem5(nx,ny),STAT=istatus)
rtem5=0
ALLOCATE(rtem6(nx,ny),STAT=istatus)
rtem6=0
ALLOCATE(rtem7(nx,ny),STAT=istatus)
rtem7=0
ALLOCATE(b1(nx,ny),STAT=istatus)
b1=0
ALLOCATE(t1(nx,ny),STAT=istatus)
t1=0
ALLOCATE(t2(nx,ny),STAT=istatus)
t2=0
ALLOCATE(t3(nx,ny),STAT=istatus)
t3=0
ALLOCATE(t4(nx,ny),STAT=istatus)
t4=0
ALLOCATE(slice(nx+ny+nz,3),STAT=istatus)
ALLOCATE(xi1(nx+ny),STAT=istatus)
ALLOCATE(yi1(nx+ny),STAT=istatus)
ALLOCATE(xi2(nx+ny),STAT=istatus)
ALLOCATE(yi2(nx+ny),STAT=istatus)
ALLOCATE(zi1(nz),STAT=istatus)
ALLOCATE(ppil(nz),STAT=istatus)
!-----------------------------------------------------------------------
!
! default output grid values (SAMEX settings):
!
!-----------------------------------------------------------------------
!
intrpl = 1
enstag = 'ARPSENS'
instit = 'CAPS_OU'
model = 'ARPS5.0.0'
qswcorn = 27.0844
qswcorw = -110.068
qnecorn = 45.1552
qnecorw = -68.10911
qni = 117
qnj = 81
mptjnew=2
trltnew(1)=39.0
trltnew(2)=39.0
trlnnew=-100.0
dxnew=30.0
! Read remaining namelist
READ(5,output_data,END=200)
WRITE(6,'(/a/)')'Namelist block output_data sucessfully read in.'
GO TO 20
200 WRITE(6,'(a)') &
'Error reading NAMELIST block output_data. ', &
'Default values used.'
20 CONTINUE
DO nf=1,nhisfile
WRITE(6,'(1x,a,a)') 'Input values for outfile: ',trim(outfile(nf))
END DO
READ(5,output_grid,END=250)
WRITE(6,'(/a/)')'Namelist block output_grid sucessfully read in.'
GO TO 30
250 WRITE(6,'(a)') &
'Error reading NAMELIST block output_grid. ', &
'Default values used.'
30 CONTINUE
lenstag = len_trim(enstag)
linstit = len_trim(instit)
lmodel = len_trim(model)
WRITE (6,*) 'Descriptor enstag: ', enstag(1:lenstag)
WRITE (6,*) 'Descriptor instit: ', instit(1:linstit)
WRITE (6,*) 'Descriptor model: ', model(1:lmodel)
WRITE(6,'(1x,a,3I5)') &
'Grid variables intrpl,qni,qnj:', intrpl,qni,qnj
WRITE(6,'(1x,a,4F11.5)') &
'Grid variables qswcorn, qswcorw, qnecorn, qnecorw:', &
qswcorn, qswcorw, qnecorn, qnecorw
WRITE(6,'(1x,a,I5,4F10.3)') &
'Grid variables mptjnew,trltnew,trlnnew,dxnew:', &
mptjnew,trltnew,trlnnew,dxnew
IF (intrpl == 1 .AND. qni*qnj > max2d_out2) THEN
WRITE (6,'(a)') &
'Error: qni*nqj > MAX2D_OUT (in enscv.inc). Stopping.'
STOP
END IF
IF (intrpl == 0 .AND. (nx-3)*(ny-3) > max2d_out2) THEN
WRITE (6,'(a)') &
'Error: (nx-3)*(ny-3) > MAX2D_OUT (in enscv.inc). Stopping.'
STOP
END IF
IF (intrpl == 0 .AND. (qni /= (nx-3).OR.qnj /= (ny-3) ) ) THEN
WRITE (6,'(a,a)') &
' WARNING: qni =/= nx-3 and/or qnj=/=ny-3. ', &
'correction is made TO make qni=nx-3 AND qnj=ny-3 '
qni=nx-3
qnj=ny-3
END IF
IF (iplot >= 1.OR.iterr > 0) CALL xdevic
DO nf=1, nhisfile
filename=hisfile(nf)
lenfil = len_trim(filename)
WRITE(6,'(a,a/)') &
' The file name of data set is ', filename(1:lenfil)
15 CONTINUE ! also continue to read another time recode
! from GrADS file
!-----------------------------------------------------------------------
!
! Read all input data arrays
!
!-----------------------------------------------------------------------
!
CALL dtaread(nx,ny,nz,nzsoil, nstyps, &
hinfmt, nchin,grdbasfn(1:lengbf),lengbf, &
filename(1:lenfil),lenfil,ctime, &
x,y,z,zp,zpsoil,uprt ,vprt ,wprt ,ptprt, pprt , &
qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate, &
radfrc,radsw,rnflx,radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
ireturn, tem1,tem2, tem3)
year1=year
month1=month
dateomon=day
houroday=hour
minohour=ctime
PRINT *,year1,month1,dateomon,houroday,ctime
!
!-----------------------------------------------------------------------
!
! ireturn = 0 for a successful read
! For hinfmt=9, i.e. the GraDs format data, ireturn is used as a
! flag indicating if there is any data at more time level to be read.
!
!-----------------------------------------------------------------------
!
IF( ireturn /= 0 .AND. hinfmt /= 9 ) GO TO 7000 !unsuccessful read
!
!-----------------------------------------------------------------------
!
! write(6,'(/5x,a,i4)')
! : 'mapproj in the data was ',mapproj
! write(6,'(/5x,a,f10.3,a)')
! : 'trulat1 in the data was ',trulat1,' degree North'
! write(6,'(/5x,a,f10.3,a)')
! : 'trulat2 in the data was ',trulat2,' degree North'
! write(6,'(/5x,a,f10.3,a)')
! : 'trulon in the data was ',trulon,' degree East'
! write(6,'(/5x,a,e15.8)')
! : 'sclfct in the data was ',sclfct
!
!-----------------------------------------------------------------------
! Convert the rainfal arries and store the old ones.
DO j=1,ny-1
DO i=1,nx-1
raint(i,j)=raing(i,j)+rainc(i,j)
! IF (nf == 1) THEN ! don't know what's the purpose (F.KONG)
! raint(i,j)=0.0E0
! rainc(i,j)=0.0E0
! END IF
END DO
END DO
! print *,raint(35,35),rainc(35,35),'rain'
CALL ary2dcv(nx,ny,raint,ni,nj,accppt)
CALL ary2dcv(nx,ny,rainc,ni,nj,convppt)
! print *,accppt(35,35),convppt(35,35),'rain'
IF (iaccu == 1) THEN
! convert accppt from 0h-now accul. rain to (Tnf-1->Tnf) accul. rain
! (by subtracting tpptold) and store the original value in tpptold
IF (nf == 1) THEN ! to clear the 1st time level values
CALL ary2dcv(nx,ny,raint,ni,nj,tpptold)
CALL ary2dcv(nx,ny,rainc,ni,nj,cpptold)
END IF
CALL pptsto(ni,nj,accppt,tpptold)
CALL pptsto(ni,nj,convppt,cpptold)
! print *,accppt(35,35),convppt(35,35),'rain'
END IF
IF (iout(nf) /= 1) CYCLE
!
!
!-----------------------------------------------------------------------
!
! Set coordinates at the grid volume center
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx-1
xc(i,j,k) = (x(i)+x(i+1))*0.5 /1000.0
END DO
END DO
END DO
DO k=1,nz
DO j=1,ny-1
DO i=1,nx
yc(i,j,k) = (y(j)+y(j+1))*0.5 /1000.0
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 /1000.0
zpc(i,j,k)= (zp(i,j,k)+zp(i,j,k+1))*0.5 /1000.0
hgt(i,j,k)=1000.0*zpc(i,j,k)
END DO
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
x2(i,j) = (x(i)+x(i+1))*0.5
y2(i,j) = (y(j)+y(j+1))*0.5
x1(i,j)=x2(i,j)/1000.0
y1(i,j)=y2(i,j)/1000.0
END DO
END DO
!-----------------------------------------------------------------------
! Find the lat/log of the points in the output grid, LatN, LonN
! if interpolation is necessary. Set the flag for INTERPOLATION
! Note that in NEWMAP, setmaps and setorig are used.
IF (intrpl == 1) THEN
CALL newmap(qni,qnj,latn,lonn,dxnew, &
qswcorn,qswcorw,trltnew,trlnnew,mptjnew)
intrpl=1
END IF
!-----------------------------------------------------------------------
! Set the map parameters and put the origin the same as used in defining
! xc and yc, x and y etc. The origin is actually point (2,2)
dx = xc(3,2,2)-xc(2,2,2)
dy = yc(2,3,2)-yc(2,2,2)
truelat(1) = trulat1
truelat(2) = trulat2
CALL setmapr( mapproj, 1.0 , truelat , trulon)
! set up parameters for map projection
CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
! print *,ctrlat,ctrlon,ctrx,ctry
! print *,dx,dy,nx,ny,xc(1,2,2),xc(2,2,2),yc(2,1,2),yc(2,2,2)
swx = ctrx- (nx-3)*dx*0.5*1000.0
swy = ctry- (ny-3)*dy*0.5*1000.0
! print *, swx,swy,'swx,swy of the original data's physical domain'
CALL setorig( 1, swx, swy)
!-----------------------------------------------------------------------
! Calculate the Coriolis para. f=2wsin(lat) at the data's scalar points
pi=ATAN(1.0)*4.0
DO j=1,ny-1
DO i=1,nx-1
CALL xytoll( 1,1, x2(i,j),y2(i,j),lat1,lon1)
f2(i,j)=SIN(lat1*pi/180.0)*2*omega
END DO
END DO
PRINT *,f2(1,1),f2(nx-1,ny-1), 'f values at data map corners'
!-----------------------------------------------------------------------
! Find the SW and NE corners of the outputgrid if they are not set
! (and intrpl=0)
! The second most SW or NE scalar point are chosen as these points
! and so qni=nx-3 and qnj=nx-3 should be satisfied
! *** for direct assignment instead of interpolation,set the new - grid
! corners and the flag for NO INTERPOLATION.also calculate latN,lonN
IF (intrpl == 0) THEN
CALL xytoll(1,1,x2(2,2),y2(2,2),qswcorn,qswcorw)
CALL xytoll(1,1,x2(nx-2,ny-2),y2(nx-2,ny-2),qnecorn,qnecorw)
intrpl=0
dxnew=dx
trltnew(1)=truelat(1)
trltnew(2)=truelat(2)
trlnnew=trulon
DO j=1,qnj
DO i=1,qni
CALL xytoll(1,1,x2(i+1,j+1),y2(i+1,j+1),lat1,lon1)
ijq = i + qni*(j-1)
latn(ijq)=lat1
lonn(ijq)=lon1
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
! Find the lat/lon of the SW and NE corners of the data grid
! and print out for comparison with the output grid's corners
CALL xytoll( 1,1, x2(1,1),y2(1,1),lat1,lon1)
CALL xytoll( 1,1, x2(nx-1,ny-1),y2(nx-1,ny-1),lat2,lon2)
PRINT *,lat1,lon1,lat2,lon2,' Data Map corners'
PRINT *,qswcorn,qswcorw,qnecorn,qnecorw,'Output Map corners'
!
!-----------------------------------------------------------------------
! Calculate the following
! 1. the x,y coordinates of the points in the output grid relatvie to
! its SW corner,unit being km, used for plotting and CLLMC,x1N + y1N.
! 2. the x,y coordinates of the points in the output grid in the old
! grid,unit being m, used for interpolation,xN + yN. comparable to
! x2,y2.
DO j=1,qnj
DO i=1,qni
ijq = i + qni*(j-1)
x1n(ijq)=(i-1)*dxnew
y1n(ijq)=(j-1)*dxnew
CALL lltoxy( 1,1, latn(ijq),lonn(ijq), xn(ijq), yn(ijq) )
IF (i == 1.AND.j == 1.OR.i == qni.AND.j == qnj) THEN
PRINT *,'SW/NE of output grid, lat,lon,x and y in input grid'
PRINT *,i,j,latn(ijq),lonn(ijq),xn(ijq),yn(ijq)
END IF
END DO
END DO
!----------------------------------------------------------------------
!
DO j=1,nj
DO i=1,ni
sreh(i,j)=0.0
li(i,j)=0.0
cape(i,j)=0.0
brn(i,j)=0.0
cin(i,j)=0.0
llmc(i,j)=0.0
pw(i,j)=0.0
END DO
END DO
! combine prt and bar
DO k=1,nz
DO j=1,ny
DO i=1,nx
u(i,j,k)=uprt(i,j,k)+ubar(i,j,k)
v(i,j,k)=vprt(i,j,k)+vbar(i,j,k)
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)=pprt(i,j,k)+pbar(i,j,k)
tem2(i,j,k)=0.0
tem3(i,j,k)=0.0
tem4(i,j,k)=0.0
END DO
END DO
END DO
! convert u,v,w to scalar points.
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem2(i,j,k)=0.5*(u(i,j,k)+u(i+1,j,k))
tem3(i,j,k)=0.5*(v(i,j,k)+v(i,j+1,k))
tem4(i,j,k)=0.5*(w(i,j,k)+w(i,j,k+1))
END DO
END DO
END DO
u=0
v=0
w=0
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
u(i,j,k)=tem2(i,j,k)
v(i,j,k)=tem3(i,j,k)
w(i,j,k)=tem4(i,j,k)
END DO
END DO
END DO
!
!----------------------------------------------------------------------
!
! Calculate temperature (K) at ARPS grid points
! Calculate dew-point temperature td using enhanced Teten's formula.
!
!-----------------------------------------------------------------------
!
CALL temper(nx,ny,nz,pt, pprt ,pbar,tem2)
CALL getdew(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, tem1, tem2, qv,tem3)
CALL getqvs(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, tem1,tem2,tem4)
ice=1
CALL reflec(nx,ny,nz,rhobar,qr,qs,qh,tem5)
!
!-----------------------------------------------------------------------
!
! 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(tem1(i,j,k))
END DO
END DO
END DO
!-----------------------------------------------------------------------
! calculate surface pressure (-log(ps)=psf),pmsl=psl,and moist. con. llmc
CALL psfsl(nx,ny,tem2(1,1,2),qv(1,1,2),tem1(1,1,2), &
hgt(1,1,2),zp(1,1,2),psf,psl)
CALL cllmc(nx,ny,nz,u,v,qv,psf,zp,rhobar,algpzc,llmc, &
x2,y2, vtwo1,vtwo2)
CALL edgfill(llmc,1,nx-1,2,nx-2,1,ny-1,2,ny-2, &
1,1,1,1)
! calculate precipitable water and composite reflectivity
DO j=1,ny-1
DO i=1,nx-1
cmpref(i,j)=0.0
pw1d=0.0
DO k=2,nz-1
! totw=qv(i,j,k)+qc(i,j,k)+qr(i,j,k)
! : +qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
totw=qv(i,j,k)
rho=tem1(i,j,k)/(rd*tem2(i,j,k)*(1.0+0.61*qv(i,j,k)))
pw1d=pw1d+totw*(zp(i,j,k+1)-zp(i,j,k))*rho
cmpref(i,j)=MAX(cmpref(i,j),tem5(i,j,k))
END DO
pw(i,j)=pw1d*1000.0
END DO
END DO
! calculate pt2m and qv2m by calling PTQV2M
CALL ptqv2m(nx,ny,nz,nstyps,zp, &
u ,v ,w ,ptprt, pprt, qvprt, &
ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil(:,:,1,:),qsoil(:,:,1,:),wetcanp, &
pt2m,qv2m, vtwo1,vtwo2, &
qvsfc,vtwo3, &
vtwo4,vtwo5,vtwo6)
!-----------------------------------------------------------------------
! Generate other 2-D data sets.
z01=-ALOG(25000.0)
! CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,hgt250)
! CALL secthrza(nx,ny,nz,u,algpzc,var2)
CALL pintp(nx,ny,nz,u,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,uwind250)
! CALL secthrza(nx,ny,nz,v,algpzc,var2)
CALL pintp(nx,ny,nz,v,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,vwind250)
CALL vrtcalc(uwind250,vwind250,ni,nj,vort250,f2,x2,y2)
z01=-ALOG(50000.0)
! CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,hgt500)
! CALL secthrz (nx,ny,nz,u,algpzc,var2)
CALL pintp(nx,ny,nz,u,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,uwind500)
! CALL secthrza(nx,ny,nz,v,algpzc,var2)
CALL pintp(nx,ny,nz,v,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,vwind500)
CALL vrtcalc(uwind500,vwind500,ni,nj,vort500,f2,x2,y2)
z01=-ALOG(85000.0)
! CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,hgt850)
! CALL secthrza(nx,ny,nz,u,algpzc,var2)
CALL pintp(nx,ny,nz,u,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,uwind850)
! CALL secthrza(nx,ny,nz,v,algpzc,var2)
CALL pintp(nx,ny,nz,v,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,vwind850)
CALL vrtcalc(uwind850,vwind850,ni,nj,vort850,f2,x2,y2)
! CALL secthrza(nx,ny,nz,tem2,algpzc,var2)
CALL pintp(nx,ny,nz,tem2,algpzc,var2,2,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,temp850)
z01=-ALOG(100000.0)
! CALL secthrza(nx,ny,nz,hgt,algpzc,var2)
CALL pintp(nx,ny,nz,hgt,algpzc,var2,3,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,thck)
z01=-ALOG(70000.0)
! CALL secthrza(nx,ny,nz,tem2,algpzc,t700)
CALL pintp(nx,ny,nz,tem2,algpzc,t700,2,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
! CALL secthrza(nx,ny,nz,tem4,algpzc,qs700)
CALL pintp(nx,ny,nz,tem4,algpzc,qs700,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
! CALL secthrza(nx,ny,nz,qv,algpzc,var2)
CALL pintp(nx,ny,nz,qv,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,rh700)
! CALL secthrza(nx,ny,nz,w,algpzc,var2)
CALL pintp(nx,ny,nz,w,algpzc,var2,1,z01, &
tem2,qv,hgt(1,1,2),zp(1,1,2),psf)
CALL ary2dcv(nx,ny,var2,ni,nj,omega700)
DO j=1,ny-1
DO i=1,nx-1
mslp(i,j) = psl(i,j)
thck(i,j) = hgt500(i,j)-thck(i,j)
rh700(i,j)= rh700(i,j)/qs700(i,j)*100.0
temp850(i,j)=temp850(i,j)-273.15
omega700(i,j)=omega700(i,j)*(-700000.0*g/287.053)/t700(i,j)
END DO
END DO
IF (minohour < 1800) THEN
! temp2m, dewp2m scheme 1, direct interpolation from grid values
z01=2.0
CALL secthrza(nx,ny,nz,tem3,zc,var2)
CALL ary2dcv(nx,ny,var2,ni,nj,dewp2m)
z01=2.0
CALL secthrza(nx,ny,nz,tem2,zc,var2)
CALL ary2dcv(nx,ny,var2,ni,nj,temp2m)
ELSE
!temp2m dewp2m scheme 2,boundary layer physics (PTQV2M's results are used)
z01=2.0
CALL secthrza(nx,ny,nz,tem1,zc,vtwo1) !vtwo1, pressure
DO j=1,ny-1
DO i=1,nx-1
vtwo2(i,j)=pt2m(i,j)*(vtwo1(i,j)/100000.0)**rddcp
END DO
END DO
CALL ary2dcv(nx,ny,vtwo2,ni,nj,temp2m)
CALL getdew(nx,ny,1,1,nx-1,1,ny-1,1,1, vtwo1, vtwo2, qv2m,vtwo3)
!vtwo3,dew point temperature
CALL ary2dcv(nx,ny,vtwo3,ni,nj,dewp2m)
!*** for the initial time,dewp at the level k=2 is used as dewp2m
! CALL ARY2DCV(nx,ny,tem3(1,1,2),ni,nj,dewp2m)
END IF
! unit conversion
DO j=1,ny-1
DO i=1,nx-1
temp2m(i,j)=temp2m(i,j)-273.15
dewp2m(i,j)=dewp2m(i,j)-273.15
IF (dewp2m(i,j) > temp2m(i,j)) dewp2m(i,j)=temp2m(i,j)
END DO
END DO
!----------------------------------------------------------------------
!
! Calculate sea level pressure (mb)
! Reduction method: Benjamin and Miller: 1990, MWR, vol.118, No.10,
! Page: 2100-2101
!
!-----------------------------------------------------------------------
!
! do 700 j=1,ny-1
! do 700 i=1,nx-1
! p00 = 0.01*(pbar(i,j,2)+pprt(i,j,2))
! if(p00.le.700.0) then
! t00=tem2(i,j,2)
! else
! t00 = t700(i,j)*(p00/700.0)**ex1
! end if
! mslp(i,j) = p00*(1.0+gamma*zpc(i,j,2)/t00)**ex2
! mslp(i,j) = psl(i,j)*0.01
! thck(i,j) = hgt500(i,j)-8.0*(mslp(i,j)-1000.0)
! thck(i,j) = hgt500(i,j)-thck(i,j)
700 CONTINUE
!-----------------------------------------------------------------------
!
! Calculate CAPE and CIN
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
tem4(i,j,k)=qv(i,j,k)/(1.-qv(i,j,k))
END DO
END DO
END DO
CALL edgfill(u,1,nx,1,nx-1,1,ny,1,ny-1, &
1,nz,1,nz-1)
CALL edgfill(v,1,nx,1,nx-1,1,ny,1,ny-1, &
1,nz,1,nz-1)
CALL edgfill(tem4,1,nx,1,nx-1,1,ny,1,ny-1, &
1,nz,1,nz-1)
CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1, &
1,nz,1,nz-1)
CALL edgfill(tem2,1,nx,1,nx-1,1,ny,1,ny-1, &
1,nz,1,nz-1)
CALL edgfill(tem3,1,nx,1,nx-1,1,ny,1,ny-1, &
1,nz,1,nz-1)
IF (icape == 2) THEN
PRINT *, ' calling cape_ncep'
CALL cape_ncep(nx,ny,nz, &
tem2,qv,tem1,cape1,cin1,li1, &
tem3,tem4,tem5,tem6, &
vtwo1,vtwo2,vtwo3,vtwo4,vtwo5,item, &
item1,item2,item3,item4,item5,item6,item7, & ! work arrays
rtem1,rtem2,rtem3,rtem4,rtem5,rtem6,rtem7) ! work arrays
PRINT *, ' back from cape_ncep'
ELSE IF (icape == 1) THEN
PRINT *, ' calling arps_be'
CALL arps_be(nx,ny,nz, &
tem1,zpc,tem2,tem4, &
lcl,lfc,el,twdf,li1,cape1,mbe,cin1,capst, &
wrk1,wrk2,wrk3,wrk4,wrk5,wrk6, &
wrk7,wrk8,wrk9,wrk10,wrk11,wrk12,t1)
PRINT *, ' back from arps_be'
ELSE
PRINT *,'icape option is incorrect, chose 1 or 2'
STOP
END IF
CALL ary2dcv(nx,ny,cape1,ni,nj,cape)
CALL ary2dcv(nx,ny,cin1,ni,nj,cin)
CALL ary2dcv(nx,ny,li1,ni,nj,li)
!
!-----------------------------------------------------------------------
!
! Calculate helicity and other shear-related paramaters
!
!-----------------------------------------------------------------------
!
PRINT *, ' calling calcshr'
CALL xytomf(nx,ny,x,y,b1)
CALL calcshr(nx,ny,nz,x,y,zp,b1, &
tem1,tem2,u,v,cape1, &
shr37,ustrm,vstrm,srlfl,srmfl,heli,brn1,brnu,blcon, &
t1,t2,t3,t4,tem4)
PRINT *, ' back from calcshr'
CALL ary2dcv(nx,ny,heli,ni,nj,sreh)
CALL ary2dcv(nx,ny,brn1,ni,nj,brn)
!-----------------------------------------------------------------------
! Interpolate the 2_d variables
!-----------------------------------------------------------------------
IF (iterr == 1) THEN
DO i=1,nx-1
DO j=1,ny-1
cmpref(i,j)=zp(i,j,2)
END DO
END DO
END IF
CALL d2intrpa(x2,y2,nx-1,ny-1, &
xn,yn,qni,qnj,intrpl, &
mslp,mslpn, &
hgt250,hgt250n, &
vort250,vort250n, &
uwind250,uwind250n, &
vwind250,vwind250n, &
hgt500,hgt500n, &
vort500,vort500n, &
uwind500,uwind500n, &
vwind500,vwind500n, &
hgt850,hgt850n, &
vort850,vort850n, &
uwind850,uwind850n, &
vwind850,vwind850n, &
temp850,temp850n, &
thck,thckn, &
rh700,rh700n, &
omega700,omega700n, &
temp2m,temp2mn, &
dewp2m,dewp2mn, &
accppt,accpptn, &
convppt,convpptn, &
sreh,srehn, &
li,lin, &
cape,capen, &
brn,brnn, &
cin,cinn, &
llmc,llmcn, &
pw,pwn, &
cmpref,cmprefn)
!-----------------------------------------------------------------------
! Find the information about the time and date of the datafile
! and add the lead time (minohour) to it, find the day of week
! read (filename(1:26),'(2x,I4,I2,I2,I2,8x,I6)')
! : year1, month1,dateomon,houroday,minohour
! print *, filename(1:26),'DATE VARIABLES'
PRINT *, year1,month1,dateomon,houroday,minohour
CALL calender(year1,month1,dateomon,houroday,minohour,dayoweek)
! print *, year1, month1,dateomon,houroday,minohour,dayoweek
! print *, 'year1, month1,dateomon,houroday,minohour,dayoweek'
! Write out the output 2-D fields
WRITE (0,*) 'About to call EnsWrtAll'
filename=outfile(nf)
lenfil = len_trim(filename)
WRITE(6,'(a,a/)') &
' The file name of output set is ', filename(1:lenfil)
CALL enswrtall (filename(1:lenfil), 2, dxnew, &
houroday,minohour,dayoweek,dateomon,month1,year1, &
mslpn, hgt250n, vort250n, uwind250n, vwind250n, &
hgt500n, vort500n, uwind500n, vwind500n, &
hgt850n, vort850n, uwind850n, vwind850n, &
temp850n, thckn, rh700n, omega700n, &
temp2mn, dewp2mn, accpptn, convpptn, &
srehn, lin, capen, brnn, cinn, llmcn,pwn,cmprefn)
!-----------------------------------------------------------------------
! Plotting of the relavent parameters if iplot=1
! Plotting of the terrain contour if iterr=1
IF (iterr == 1) THEN
CALL a2plt(cmprefn,x1n,y1n,qni,qnj,dxnew,100.0,'HTER')
ELSE IF (iplot == 1) THEN
CALL a2plt(accpptn,x1n,y1n,qni,qnj,dxnew,1.0,'TPPT')
CALL a2plt(convpptn,x1n,y1n,qni,qnj,dxnew,1.0,'CPPT')
ELSE IF (iplot == 2) THEN
CALL a2plt(mslpn,x1n,y1n,qni,qnj,dxnew,200.0,'MSLP')
CALL a2plt(hgt250n,x1n,y1n,qni,qnj,dxnew,50.0,'H250')
CALL a2plt(hgt500n,x1n,y1n,qni,qnj,dxnew,20.0,'H500')
CALL a2plt(hgt850n,x1n,y1n,qni,qnj,dxnew,20.0,'H850')
CALL a2plt(vort250n,x1n,y1n,qni,qnj,dxnew,0.00002,'V250')
CALL a2plt(vort500n,x1n,y1n,qni,qnj,dxnew,0.00002,'V500')
CALL a2plt(vort850n,x1n,y1n,qni,qnj,dxnew,0.00002,'V850')
CALL a2plt(temp850n,x1n,y1n,qni,qnj,dxnew,2.0,'T850')
CALL a2plt(thckn,x1n,y1n,qni,qnj,dxnew,20.0,'THCK')
CALL a2plt(rh700n,x1n,y1n,qni,qnj,dxnew,10.0,'RH70')
CALL a2plt(omega700n,x1n,y1n,qni,qnj,dxnew,2.0,'OM70')
CALL a2plt(temp2mn,x1n,y1n,qni,qnj,dxnew,2.0,'TC2M')
CALL a2plt(dewp2mn,x1n,y1n,qni,qnj,dxnew,2.0,'TD2M')
CALL a2plt(accpptn,x1n,y1n,qni,qnj,dxnew,2.0,'TPPT')
CALL a2plt(convpptn,x1n,y1n,qni,qnj,dxnew,2.0,'CPPT')
CALL a2plt(srehn,x1n,y1n,qni,qnj,dxnew,50.0,'SREH')
CALL a2plt(lin,x1n,y1n,qni,qnj,dxnew,2.0,'LIDX')
CALL a2plt(capen,x1n,y1n,qni,qnj,dxnew,500.0,'CAPE')
CALL a2plt(brnn,x1n,y1n,qni,qnj,dxnew,10.0,'BRN ')
CALL a2plt(cinn,x1n,y1n,qni,qnj,dxnew,50.0,'CIN ')
CALL a2plt(llmcn,x1n,y1n,qni,qnj,dxnew,0.0001,'LLMC')
CALL a2plt(pwn,x1n,y1n,qni,qnj,dxnew,2000.0,' PW ')
CALL a2plt(cmprefn,x1n,y1n,qni,qnj,dxnew,2000.0,'CPRF')
END IF
7000 CONTINUE
END DO
WRITE(6,'(a)') 'Program completed.'
IF (iplot >= 1) CALL xgrend
STOP
END PROGRAM arpsenscv
SUBROUTINE a2plt(a,x,y,nx,ny,dx,zinc,title) 26
IMPLICIT NONE
INTEGER :: nx,ny,ncl,mode
INTEGER :: i,j
REAL :: a(nx,ny),x(nx,ny),y(nx,ny),dx
INTEGER :: iw(212,152)
REAL :: cl(100),zmax,zmin,zinc
CHARACTER (LEN=4) :: title
CALL xpspac(0.1,0.9,0.05,0.15+80.0/116.0*0.8)
! call XMAP(0.0,2400.0,0.0,2080.0)
! call XMAP(0.0,3480.0,0.0,2400.0)
CALL xmap(0.0,(nx-1)*dx,0.0,(ny-1)*dx)
CALL xaxsca(0.0,(nx-1)*dx,dx,0.0,(ny-1)*dx,dx)
! call XAXSCA(0.0,3480.0,30.0,0.0,2400.0,30.0)
! call XAXSCA(0.0,2400.0,32.0,0.0,2080.0,32.0)
CALL xbordr
ncl=100
zmax=-1.0E-10
zmin=+1.0E+10
DO j=2,ny-1
DO i=2,nx-1
IF (a(i,j) > zmax) zmax=a(i,j)
IF (a(i,j) < zmin) zmin=a(i,j)
END DO
END DO
IF (ABS(zmax-zmin) < 1.0E-8) GO TO 200
cl(1)=INT(zmin/zinc)*zinc
cl(2)=cl(1)+zinc
mode=1
CALL xconta( zmax=cl(ncl)
zmin=cl(1)
zinc=cl(MIN(2,ncl))-cl(1)
200 CONTINUE
PRINT *,zmax,zmin,zinc
CALL xclimt(zmax,zmin,zinc,0.0)
CALL xchsiz(80.00)
CALL xcharc(1500.,-100.,title)
CALL xframe
RETURN
END SUBROUTINE a2plt
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C
SUBROUTINE pptsto(ni,nj,a,b) 3
INTEGER :: ni,nj,i,j
REAL :: a(ni,nj),b(ni,nj),c
DO j=1,nj
DO i=1,ni
c=a(i,j)-b(i,j)
b(i,j)=a(i,j)
a(i,j)=c
END DO
END DO
RETURN
END SUBROUTINE pptsto
!
SUBROUTINE ary2dcv(nx,ny,a,ni,nj,b) 27
IMPLICIT NONE
INTEGER :: nx,ny,ni,nj
REAL :: a(nx,ny),b(ni,nj)
INTEGER :: i,j
! print *,a(35,35),a(78,68)
DO j=1,nj
DO i=1,ni
b(i,j)=a(i,j)
END DO
END DO
RETURN
END SUBROUTINE ary2dcv
SUBROUTINE vrtcalc(u,v,ni,nj,b,f,x,y) 3,1
IMPLICIT NONE
INTEGER :: ni,nj
REAL :: u(ni,nj),v(ni,nj),b(ni,nj),f(ni,nj),x(ni,nj),y(ni,nj)
INTEGER :: i,j
DO j=1,nj
DO i=1,ni
b(i,j)=0.0
END DO
END DO
DO j=2,nj-1
DO i=2,ni-1
b(i,j)=(v(i+1,j)-v(i-1,j))/(x(i+1,j)-x(i-1,j)) &
-(u(i,j+1)-u(i,j-1))/(y(i,j+1)-y(i,j-1)) &
+f(i,j)
END DO
END DO
CALL edgfill(b,1,ni,2,ni-1,1,nj,2,nj-1,1,1,1,1)
RETURN
END SUBROUTINE vrtcalc
SUBROUTINE enswrthd (filenm, fileun, & 1,3
enstag, instit, model, resltn, resun, &
houroday, minohour, dayoweek, dateomon, &
month, year)
!
!======================================================================-
! &&&& G E N E R A L I N F O R M A T I O N &&&&
!======================================================================-
!
! PURPOSE: This output subroutine produces the header for SAMEX ASCII
! fields for conversion to images.
!
! AUTHOR: Henry Neeman, Univ. of Oklahoma, Project SAMEX
! (hneeman@tornado.caps.ou.edu)
!
! HISTORY:
! 1998/03/18 First written [HJN]
!
! INPUT ARGUMENTS:
! filenm: string for filename to output data to
! fileun: integer for Fortran file unit
! enstag: tag for ensemble (e.g., 'SAMEX')
! instit: string for name of institution
! (e.g., 'CAPS', 'NCEP', 'NCAR', 'AFWA', 'NSSL')
! model: string for name of model
! (e.g., 'ARPS4.3.3', 'MM5V2')
! Note: the model name should not contain any blanks or tabs.
! resltn: real for resolution (e.g., 30 for 30 km resolution)
! resun: string for name of resolution units (e.g., 'km')
! houroday: integer for hour of the day (e.g., 20)
! minohour: integer for minute of the hour (e.g., 30)
! dayoweek: integer for day the week
! 1=Sun, 2=Mon, 3=Tue, 4=Wed, 5=Thu, 6=Fri, 7=Sat
! dateomon: integer for date of the month
! month: integer for month (1=Jan, 2=Feb, etc.)
! year: integer for year (e.g., 1998)
!
! OUTPUT:
! ASCII file header for 2D data field
!
! I/O:
! Output data to a Fortran unformatted ASCII file.
!
! SPECIAL REQUIREMENTS:
! <none>
!
! OTHER INFORMATION:
! <none>
!
!======================================================================-
! %%%% D E C L A R A T I O N S %%%%
!======================================================================-
!
IMPLICIT NONE
! - - - - - - - - - Include files - - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - - Constant declarations - - - - - - - - - - - - - - -
!
! - - - - - - - - - Argument declarations - - - - - - - - - - - - - - -
CHARACTER (LEN=*) :: filenm
INTEGER :: fileun
CHARACTER (LEN=*) :: enstag, instit, model
REAL :: resltn
CHARACTER (LEN=*) :: resun
INTEGER :: houroday, minohour, dayoweek, dateomon, month, year
INTEGER :: lenstag, linstit, lmodel
! - - - - - - - - - Global/External declarations - - - - - - - - - - -
!
! - - - - - - - - - Local declarations - - - - - - - - - - - - - - - -
INTEGER :: iost
! - - - - - - - - - Data statements - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - - Statement functions - - - - - - - - - - - - - - - -
!
!======================================================================-
! @@@@ E X E C U T A B L E C O D E @@@@
!======================================================================-
!
100 FORMAT (a)
110 FORMAT (a, ' ', a, ' ', a, ' ', g10.3, ' ', a)
120 FORMAT (i2, i3, i2, i3, i3, i5)
lenstag = 0
CALL strlnth( enstag, lenstag)
linstit = 0
CALL strlnth( instit, linstit)
lmodel = 0
CALL strlnth( model, lmodel)
OPEN (UNIT=fileun,IOSTAT=iost,ERR=150,FILE=filenm)
WRITE (fileun,*) enstag(1:lenstag)
WRITE (fileun,*) instit(1:linstit), ' ', model(1:lmodel), &
' ', resltn, ' ', resun
WRITE (fileun,*) houroday, ' ', minohour, ' ', dayoweek, ' ', &
dateomon, ' ', month, ' ', year
RETURN
150 WRITE (0,160) filenm
RETURN
160 FORMAT ('Error: cannot open file ', a, &
' for writing ensemble header!')
END SUBROUTINE enswrthd
!
!
! ########################################
! ########################################
! ########################################
! ######## ########
! ######## ENSEMBLE ########
! ######## ########
! ########################################
! ########################################
! ########################################
!
!
SUBROUTINE enswrtq (filenm, fileun, & 29
qtag, qunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, qfield, &
islast)
!
!======================================================================-
! &&&& G E N E R A L I N F O R M A T I O N &&&&
!======================================================================-
!
! PURPOSE: This output subroutine produces each SAMEX ASCII quantity
! field for conversion to images.
!
! AUTHOR: Henry Neeman, Univ. of Oklahoma, Project SAMEX
! (hneeman@tornado.caps.ou.edu)
!
! HISTORY:
! 1998/03/18 First written [HJN]
!
! INPUT ARGUMENTS:
! filenm: string for filename to output data to
! fileun: integer for Fortran file unit
! qtag: string for quantity tag (e.g., 'MSLP')
! qunits: string for quantity units (e.g., 'km')
! qswcorn: real for quantity southwest corner north coordinate
! qswcore: real for quantity southwest corner west coordinate
! qnecorn: real for quantity northeast corner north coordinate
! qnecore: real for quantity northeast corner west coordinate
! qstag: string for quantity staggering (for now, CENTER)
! qni: integer for quantity number of cells along the major axis
! qnj: integer for quantity number of cells along the minor axis
! qfield: real array for the quantity field data
! islast: logical for the quantity being the last in the file
!
! OUTPUT:
! ASCII file for 2D data field
!
! I/O:
! Output data to a Fortran unformatted ASCII file.
!
! SPECIAL REQUIREMENTS:
! <none>
!
! OTHER INFORMATION:
! <none>
!
!======================================================================-
! %%%% D E C L A R A T I O N S %%%%
!======================================================================-
!
IMPLICIT NONE
! - - - - - - - - - Include files - - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - - Constant declarations - - - - - - - - - - - - - - -
!
! - - - - - - - - - Argument declarations - - - - - - - - - - - - - - -
CHARACTER (LEN=*) :: filenm
INTEGER :: fileun
CHARACTER (LEN=*) :: qtag, qunits
REAL :: qswcorn, qswcorw, qnecorn, qnecorw
CHARACTER (LEN=*) :: qstag
INTEGER :: qni, qnj
REAL :: qfield(qni,qnj)
LOGICAL :: islast
! - - - - - - - - - Global/External declarations - - - - - - - - - - -
!
! - - - - - - - - - Local declarations - - - - - - - - - - - - - - - -
INTEGER :: iost
INTEGER :: i, j
! - - - - - - - - - Data statements - - - - - - - - - - - - - - - - - -
! - - - - - - - - - Statement functions - - - - - - - - - - - - - - - -
!======================================================================-
! @@@@ E X E C U T A B L E C O D E @@@@
!======================================================================-
!
200 FORMAT (a, ' ', a)
210 FORMAT (g10.5, ' N ', g10.5, ' W ', g10.5, ' N ', g10.5, ' W ')
220 FORMAT (a)
230 FORMAT (i5, ' ', i5)
240 FORMAT (e16.8, e16.8, e16.8, e16.8, e16.8)
WRITE (fileun,*) qtag, ' ', qunits
WRITE (fileun,*) qswcorn, ' N ', qswcorw, ' W ', &
qnecorn, ' N ', qnecorw, ' W'
WRITE (fileun,*) qstag
WRITE (fileun,*) qni, ' ', qnj
DO j = 1, qnj
WRITE (fileun,240) (qfield(i,j),i=1,qni)
END DO
IF (islast) THEN
CLOSE (UNIT=fileun,IOSTAT=iost,ERR=250,STATUS='KEEP')
ELSE
WRITE (fileun,*) ' '
END IF
RETURN
250 WRITE (0,260) filenm
RETURN
260 FORMAT ('Error: cannot close file ', a, &
' after writing ensemble data!')
END SUBROUTINE enswrtq
!
!
! ########################################
! ########################################
! ########################################
! ######## ########
! ######## ENSEMBLE ########
! ######## ########
! ########################################
! ########################################
! ########################################
!
!
SUBROUTINE enswrtall (filenm, fileun, resltn, & 1,30
houroday, minohour, dayoweek, dateomon, &
month, year, &
mslp, hgt250, vort250, uwind250, vwind250, &
hgt500, vort500, uwind500, vwind500, &
hgt850, vort850, uwind850, vwind850, &
temp850, thck, rh700, omega700, &
temp2m, dewp2m, accppt, convppt, &
sreh, li, cape, brn, cin, llmc,pw,cmpref)
!
!======================================================================-
! &&&& G E N E R A L I N F O R M A T I O N &&&&
!======================================================================-
!
! PURPOSE: This output subroutine produces the header for SAMEX ASCII
! fields for conversion to images.
!
! AUTHOR: Henry Neeman, Univ. of Oklahoma, Project SAMEX
! (hneeman@tornado.caps.ou.edu)
!
! HISTORY:
! 1998/03/18 First written [HJN]
!
! INPUT ARGUMENTS:
! filenm: string for filename to output data to
! fileun: integer for Fortran file unit
! Note: the model name should not contain any blanks or tabs.
! resltn: real for resolution (e.g., 30 for 30 km resolution)
! houroday: integer for hour of the day (e.g., 20)
! minohour: integer for minute of the hour (e.g., 30)
! dayoweek: integer for day the week
! 1=Sun, 2=Mon, 3=Tue, 4=Wed, 5=Thu, 6=Fri, 7=Sat
! dateomon: integer for date of the month
! month: integer for month (1=Jan, 2=Feb, etc.)
! year: integer for year (e.g., 1998)
! mslp: real array for mean sea level pressure in pascals
! hgt250: real array for 250 mb height in m
! vort250: real array for 250 mb absolute vorticity in 1/s
! uwing250: real array for 250 mb u wind velocity in m/s
! vwing250: real array for 250 mb v wind velocity in m/s
! hgt500: real array for 500 mb height in m
! vort500: real array for 500 mb absolute vorticity in 1/s
! uwing500: real array for 500 mb u wind velocity in m/s
! vwing500: real array for 500 mb v wind velocity in m/s
! hgt850: real array for 850 mb height in m
! vort850: real array for 850 mb absolute vorticity in 1/s
! uwing850: real array for 850 mb u wind velocity in m/s
! vwing850: real array for 850 mb v wind velocity in m/s
! temp850: real array for 850 mb temperature in degrees C
! thck: real array for 1000-500 mb thickness in m
! rh700: real array for 700 mb relative humidity in %
! omega700: real array for 700 mb omega vertical velocity in -microb/s
! temp2m: real array for 2 m temperature in degrees C
! dewp2m: real array for 2 m dewpoint in degrees C
! accppt: real array for accumulated precipitation in mm
! convppt: real array for convective precipitation in mm
! sreh: real array for storm-relative environmental helicity
! in m/s^2
! li: real array for lifted index in degrees C
! cape: real array for convective available potential energy in J/kg
! brn: real array for bulk Richardson number (unitless)
! cin: real array for convective inhibition in J/kg
! llmc: real array for low-level moisture convergence in g/kg/s
! pw: real array for precipitable water in g/m^2
!
! OUTPUT:
! ASCII file header for 2D data field
!
! I/O:
! Output data to a Fortran unformatted ASCII file.
!
! SPECIAL REQUIREMENTS:
! <none>
!
! OTHER INFORMATION:
! <none>
!
!======================================================================-
! %%%% D E C L A R A T I O N S %%%%
!======================================================================-
!
IMPLICIT NONE
! - - - - - - - - - Include files - - - - - - - - - - - - - - - - - - -
INCLUDE 'enscv.inc'
! - - - - - - - - - Constant declarations - - - - - - - - - - - - - - -
!
! - - - - - - - - - Argument declarations - - - - - - - - - - - - - - -
CHARACTER (LEN=*) :: filenm
INTEGER :: fileun
REAL :: resltn
INTEGER :: houroday, minohour, dayoweek, dateomon, month, year
REAL :: mslp(qni,qnj)
REAL :: hgt250(qni,qnj), vort250(qni,qnj), &
uwind250(qni,qnj), vwind250(qni,qnj)
REAL :: hgt500(qni,qnj), vort500(qni,qnj), &
uwind500(qni,qnj), vwind500(qni,qnj)
REAL :: hgt850(qni,qnj), vort850(qni,qnj), &
uwind850(qni,qnj), vwind850(qni,qnj)
REAL :: temp850(qni,qnj)
REAL :: thck(qni,qnj), rh700(qni,qnj), omega700(qni,qnj), &
temp2m(qni,qnj), dewp2m(qni,qnj)
REAL :: accppt(qni,qnj), convppt(qni,qnj)
REAL :: sreh(qni,qnj), li(qni,qnj), cape(qni,qnj), brn(qni,qnj), &
cin(qni,qnj), llmc(qni,qnj), pw(qni,qnj)
REAL :: cmpref(qni,qnj)
! - - - - - - - - - Global/External declarations - - - - - - - - - - -
!
! - - - - - - - - - Local declarations - - - - - - - - - - - - - - - -
!
! - - - - - - - - - Data statements - - - - - - - - - - - - - - - - - -
!
! - - - - - - - - - Statement functions - - - - - - - - - - - - - - - -
!
!======================================================================-
! @@@@ E X E C U T A B L E C O D E @@@@
!======================================================================-
!
PRINT *, 'starting (inside) EnsWrt'
CALL enswrthd (filenm, fileun, &
enstag, instit, model, resltn, resun, &
houroday, minohour, dayoweek, dateomon, &
month, year)
CALL enswrtq (filenm, fileun, &
mslptag, mslpunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, mslp, &
.false.)
CALL enswrtq (filenm, fileun, &
hgt250tag, hgt250units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, hgt250, &
.false.)
CALL enswrtq (filenm, fileun, &
vort250tag, vort250units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, vort250, &
.false.)
CALL enswrtq (filenm, fileun, &
uwind250tag, uwind250units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, uwind250, &
.false.)
CALL enswrtq (filenm, fileun, &
vwind250tag, vwind250units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, vwind250, &
.false.)
CALL enswrtq (filenm, fileun, &
hgt500tag, hgt500units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, hgt500, &
.false.)
CALL enswrtq (filenm, fileun, &
vort500tag, vort500units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, vort500, &
.false.)
CALL enswrtq (filenm, fileun, &
uwind500tag, uwind500units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, uwind500, &
.false.)
CALL enswrtq (filenm, fileun, &
vwind500tag, vwind500units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, vwind500, &
.false.)
CALL enswrtq (filenm, fileun, &
hgt850tag, hgt850units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, hgt850, &
.false.)
CALL enswrtq (filenm, fileun, &
vort850tag, vort850units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, vort850, &
.false.)
CALL enswrtq (filenm, fileun, &
uwind850tag, uwind850units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, uwind850, &
.false.)
CALL enswrtq (filenm, fileun, &
vwind850tag, vwind850units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, vwind850, &
.false.)
CALL enswrtq (filenm, fileun, &
temp850tag, temp850units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, temp850, &
.false.)
CALL enswrtq (filenm, fileun, &
thcktag, thckunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, thck, &
.false.)
CALL enswrtq (filenm, fileun, &
rh700tag, rh700units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, rh700, &
.false.)
CALL enswrtq (filenm, fileun, &
omega700tag, omega700units, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, omega700, &
.false.)
CALL enswrtq (filenm, fileun, &
temp2mtag, temp2munits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, temp2m, &
.false.)
CALL enswrtq (filenm, fileun, &
dewp2mtag, dewp2munits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, dewp2m, &
.false.)
CALL enswrtq (filenm, fileun, &
accppttag, accpptunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, accppt, &
.false.)
CALL enswrtq (filenm, fileun, &
convppttag, convpptunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, convppt, &
.false.)
CALL enswrtq (filenm, fileun, &
srehtag, srehunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, sreh, &
.false.)
CALL enswrtq (filenm, fileun, &
litag, liunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, li, &
.false.)
CALL enswrtq (filenm, fileun, &
capetag, capeunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, cape, &
.false.)
CALL enswrtq (filenm, fileun, &
brntag, brnunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, brn, &
.false.)
CALL enswrtq (filenm, fileun, &
cintag, cinunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, cin, &
.false.)
CALL enswrtq (filenm, fileun, &
llmctag, llmcunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, llmc, &
.false.)
CALL enswrtq (filenm, fileun, &
pwtag, pwunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, pw, &
.false.)
CALL enswrtq (filenm, fileun, &
crftag, crfunits, &
qswcorn, qswcorw, qnecorn, qnecorw, &
qstag, &
qni, qnj, cmpref, &
.true.)
RETURN
END SUBROUTINE enswrtall
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SECTHRZA ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE secthrza(nx,ny,nz,s,z,ss1) 3
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Interpolate 3-D data to a given horizontal level.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Ming Xue & Hao Jin
! 12/18/92.
!
! MODIFICATION HISTORY:
!
! 9/10/94 (Weygandt & Y. Lu)
! Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! s 3-dimensional array of data to contour
! s1 2-dimensional array of data to contour
!
! 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 physical space (m)
!
! z01 value of x for first i grid point to plot
!
! OUTPUT:
! ss1 interpolated 3-D data to a given horizontal level
!
!-----------------------------------------------------------------------
!
! Parameters of output
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! Number of grid points in 3 directions
REAL :: s(nx,ny,nz) ! 3-dimensional array of data to contour
REAL :: z(nx,ny,nz) ! z coordinate of grid points
! in physical space (m)
REAL :: ss1(nx,ny) ! interpolated 3-D data to a
! given horizontal level
!
!-----------------------------------------------------------------------
!
! Common blocks for plotting control parameters
!
!-----------------------------------------------------------------------
!
REAL :: z01 ! the given height of the slice
COMMON /sliceh/z01
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Find index for interpolation
!
!-----------------------------------------------------------------------
!
DO j=1,ny-1
DO i=1,nx-1
IF(z01 <= z(i,j,1)) GO TO 11
IF(z01 >= z(i,j,nz-1)) GO TO 12
DO k=2,nz-2
IF(z01 >= z(i,j,k).AND.z01 < z(i,j,k+1)) GO TO 15
END DO
11 k=1
GO TO 15
12 k=nz-1
GO TO 15
15 ss1(i,j)=s(i,j,k)+(s(i,j,k+1)-s(i,j,k))* &
(z01-z(i,j,k))/(z(i,j,k+1)-z(i,j,k))
!-----------------------------------------------------------------------
!
! If the data point is below the ground level, set the
! data value to the missing value.
!
!-----------------------------------------------------------------------
! IF( z01.lt. z(i,j,2) ) ss1(i,j) = -9999.0
END DO
END DO
RETURN
END SUBROUTINE secthrza
SUBROUTINE ptqv2m(nx,ny,nz,nstyps,zp, & 1,5
u ,v ,w ,ptprt, pprt, qvprt, &
ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsfc,wetsfc,wetcanp, &
tem1,tem2, tem3,tem4, &
qvsfc,ptsfc, &
vke2,ptke2,qvke2)
IMPLICIT NONE
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'sfcphycst.inc'
! input
INTEGER :: nx,ny,nz,nstyps
REAL :: zp(nx,ny,nz)
REAL :: u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz)
REAL :: ptprt(nx,ny,nz),pprt(nx,ny,nz),qvprt(nx,ny,nz)
REAL :: ptbar(nx,ny,nz),pbar(nx,ny,nz),qvbar(nx,ny,nz)
REAL :: rhobar(nx,ny,nz)
INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type
REAL :: stypfrct(nx,ny,nstyps) ! Soil type fraction
INTEGER :: vegtyp (nx,ny) ! Vegetation type
REAL :: lai (nx,ny) ! Leaf Area Index
REAL :: roufns (nx,ny) ! Surface roughness
REAL :: veg (nx,ny) ! Vegetation fraction
REAL :: tsfc (nx,ny,0:nstyps) ! Soil Skin Temperature (K)
REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
! output
REAL :: tem1(nx,ny),tem2(nx,ny)
! temperary arrays
REAL :: tem3(nx,ny),tem4(nx,ny)
REAL :: ptsfc(nx,ny), qvsfc(nx,ny,0:nstyps)
REAL :: vke2(nx,ny),ptke2(nx,ny),qvke2(nx,ny)
! temperary variables
INTEGER :: i,j,k,is
REAL :: z1,z1drou,z1droup
REAL :: usuf,vsuf,wsuf,tsuf
REAL :: psfc,qvsatts, rhgs, pterm, pi
PARAMETER (pi=3.141592654)
REAL :: gumove,gvmove
DATA gumove /0.0/
DATA gvmove /0.0/
! field moisture capacity
REAL :: wfc(13)
DATA wfc / .135, .150, .195, .255, .240, .255, &
.322, .325, .310, .370, .367, 1.0E-25, 1.00/
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
PRINT *, 'starting PTQV2M'
IF (nstyp <= 1) THEN
! for nstyp=1, the DTAREAD subroutine returns tsfc etc. with (i,j,0)
! filled with data and (i,j,1) not assigned. This loop fills this Gap.
DO j = 1, ny
DO i = 1, nx
tsfc(i,j,1)=tsfc(i,j,0)
wetsfc(i,j,1)=wetsfc(i,j,0)
wetcanp(i,j,1)=wetcanp(i,j,0)
stypfrct(i,j,1)=1.0
END DO
END DO
END IF
DO j = 1, ny
DO i = 1, nx
qvsfc(i,j,0) = 0.0
tem1(i,j)=0.0
tem2(i,j)=0.0
tem3(i,j)=0.0
tem4(i,j)=0.0
END DO
END DO
! deduce qvsfc(effctive qv at surface) ,adapted from INISFC
PRINT *,'NSTYPS/NSTYP',nstyps,nstyp
DO is=1,nstyp
DO j = 1, ny-1
DO i = 1, nx-1
psfc = .5 * ( pbar(i,j,1) + pprt(i,j,1) &
+ pbar(i,j,2) + pprt(i,j,2) )
qvsatts = f_qvsat( psfc, tsfc(i,j,is) )
IF ( soiltyp(i,j,is) == 12 .OR. soiltyp(i,j,is) == 13) THEN
! Ice and water
rhgs = 1.0
ELSE
pterm = .5 + SIGN( .5, wetsfc(i,j,is) - wfc(soiltyp(i,j,is)) )
rhgs = pterm + (1.-pterm) &
* 0.5 * ( 1. - COS( wetsfc(i,j,is) * pi &
/ wfc(soiltyp(i,j,is)) ) )
END IF
qvsfc(i,j,is) = rhgs * qvsatts
END DO
END DO
END DO
DO is=1,nstyp
DO j = 1, ny-1
DO i = 1, nx-1
qvsfc(i,j,0) = qvsfc(i,j,0) + qvsfc(i,j,is)*stypfrct(i,j,is)
END DO
END DO
END DO
! calculate wind speed,pt and qv at k=2, pt at surface
DO j = 1, ny-1
DO i = 1, nx-1
usuf=gumove+0.5*(u(i,j,2)+u(i+1,j,2))
vsuf=gvmove+0.5*(v(i,j,2)+v(i,j+1,2))
wsuf= 0.5*(w(i,j,2)+w(i,j,3))
vke2(i,j)=MAX( SQRT(usuf*usuf+vsuf*vsuf+wsuf*wsuf),vsfcmin)
ptke2(i,j)=ptprt(i,j,2)+ptbar(i,j,2)
qvke2(i,j)=qvprt(i,j,2)+qvbar(i,j,2)
psfc=0.5*(pprt(i,j,1)+pbar(i,j,1) &
+pprt(i,j,2)+pbar(i,j,2))
ptsfc(i,j)=tsfc(i,j,0)*(100000.0/psfc)**rddcp
END DO
END DO
! calculate C_u and C_pt for nuetral condition stored in tem1 and
! tem2, respectively
DO j = 1, ny-1
DO i = 1, nx-1
z1=0.5*(zp(i,j,3)-zp(i,j,2))
z1=MIN(z1,z1limit)
z1drou=z1/roufns(i,j)
z1droup=(z1+roufns(i,j))/roufns(i,j)
IF (z1droup < 0.0) THEN
PRINT *,i,j,z1droup,z1,roufns(i,j)
STOP
END IF
tem1(i,j)=kv/LOG(z1droup)
tem2(i,j)=tem1(i,j)/prantl0l
IF (soiltyp(i,j,1) == 13.AND.stypfrct(i,j,1) >= 0.5) THEN
tem1(i,j)=SQRT((0.4+0.079*vke2(i,j))*1.e-3)
tem2(i,j)=1.14E-3/tem1(i,j)
END IF
END DO
END DO
! calculate C_pt according to land/water and stability
CALL cptca(nx,ny,nz,nx-1,ny-1,zp,roufns,vke2(1,1), &
ptsfc(1,1),ptke2(1,1),tem2(1,1),tem3(1,1),1)
CALL cptcwtra(nx,ny,nz,nx-1,ny-1,zp,soiltyp,vke2(1,1), &
ptsfc(1,1),ptke2(1,1),tem2(1,1),tem3(1,1),1)
CALL cptca(nx,ny,nz,nx-1,ny-1,zp,roufns,vke2(1,1), &
ptsfc(1,1),ptke2(1,1),tem2(1,1),tem4(1,1),2)
CALL cptcwtra(nx,ny,nz,nx-1,ny-1,zp,soiltyp,vke2(1,1), &
ptsfc(1,1),ptke2(1,1),tem2(1,1),tem4(1,1),2)
DO j = 1, ny-1
DO i = 1, nx-1
tem1(i,j)=ptsfc(i,j)+tem3(i,j)/tem4(i,j)*(ptke2(i,j)-ptsfc(i,j))
tem2(i,j)=qvsfc(i,j,0)+tem3(i,j)/tem4(i,j) &
*(qvke2(i,j)-qvsfc(i,j,0))
! tem2(i,j)=qvke2(i,j)
END DO
END DO
PRINT *, 'finishing PTQV2M'
RETURN
END SUBROUTINE ptqv2m
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CPTC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cptca(nx,ny,nz,iend,jend,zp,roufns,wspd,ptsfc,pt1, & 2
c_ptneu,c_pt,i1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute c_pt (a nondimensional temperature scale)
!
! The quantity c_pt is used by the subroutine SFCFLX to obtain
! surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
! AUTHORS: V. Wong, X. Song and N. Lin
! 9/10/1993
!
! For the unstable case (Bulk Richardson number bulkri < 0 ), the
! formulation is based on the paper "On the Analytical Solutions of
! Flux-Profile relationships for the Atmospheric Surface Layer" by
! D.W. Byun in J. of Applied Meteorlogy, July 1990, pp. 652-657.
!
! For stable case, the formulation is based on the paper "A Short
! History of the Operational PBL - Parameterization at ECMWF" by
! J.F.Louis, M. Tiedtke and J.F. Geleyn in "Workshop on Planetary
! Boundary Layer Parameterization", a publication by European Centre
! for Medium Range Weather Forecasts, 25-27 Nov. 1981.
!
! MODIFICATION HISTORY:
!
! 9/04/1994 (K. Brewster, V. Wong and X. Song)
! Facelift.
! 2/27/95 (V. Wong and X. Song)
!
! 02/07/96 (V.Wong and X.Song)
! Set an upper limiter, z1limit, for depth of the surface layer z1.
!
! 05/01/97 (V. Wong and X. Tan)
! Changed the computation of temperature scale
!
! 05/29/97 (V. Wong and X. Tan)
! Modified the formulation considering the height of the surface
! layer z1 may equal zero.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! iend i-index where evaluation ends.
! jend j-index where evaluation ends.
!
! zp The physical height coordinate defined at w-point of
! staggered grid.
! wspd Surface wind speed (m/s), defined as
! sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
! roufns Surface roughness
!
! ptsfc Potential temperature at the ground level (K)
! pt1 Potential temperature at the 1st scalar point above
! ground level, (K)
!
! c_ptneu Temperature scale (K) at neutral state, defined by
! surface heat flux / friction velocity
!
! OUTPUT:
!
! c_pt Temperature scale (K), defined by
! surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number grid points in 3 directions
INTEGER :: iend ! i-index where evaluation ends
INTEGER :: jend ! j-index where evaluation ends
REAL :: zp (nx,ny,nz) ! The physical height coordinate
! defined at w-point of staggered grid.
REAL :: roufns(nx,ny) ! Surface roughness length
REAL :: wspd (nx,ny) ! Surface wind speed (m/s)
REAL :: ptsfc(nx,ny) ! Potential temperature at the
! ground level (K)
REAL :: pt1 (nx,ny) ! Potential temperature at the
! 1st scalar
! point above ground level, (K)
REAL :: c_ptneu(nx,ny) ! Normalized temperature scale (K)
! at neutral state
REAL :: c_pt (nx,ny) ! Temperature scale (K)
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j
REAL :: z1 ! The height of 1st scalar point above
! ground level (m)
REAL :: dz0 ! z1-roufns, where roufns is defined as
! surface roughness length
REAL :: bulkri ! Bulk Richardson number
REAL :: stabp ! Monin-Obukhov STABility Parameter
! (zeta)
REAL :: y1,y0,psih ! Intermediate parameters needed
REAL :: z1drou,qb3pb2
REAL :: c7,c8
REAL :: sb,qb,pb,thetab,tb ! During computations
REAL :: a,b,c,d
REAL :: tempan,sqrtqb
REAL :: zt ! Thermal roughness length
REAL :: dzt,ztdrou
REAL :: z1droup,ztdroup
REAL :: dz00
INTEGER :: i1
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Following Byun (1990).
!
!-----------------------------------------------------------------------
!
PRINT *,'starting PTCP'
DO j=1,jend
DO i=1,iend
dz00 = 0.5*(zp(i,j,3)-zp(i,j,2))
dz00 = MIN(dz00,z1limit)
dz00 = dz00 - roufns(i,j)
bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dz00/ &
(wspd(i,j)*wspd(i,j))
IF (i1 == 1) THEN
z1 = 0.5*(zp(i,j,3)-zp(i,j,2))
ELSE
z1=2.0
END IF
z1 = MIN(z1,z1limit)
zt = ztz0*roufns(i,j)
dzt = z1-zt
ztdrou = z1/zt
ztdroup = (z1+zt)/zt
dz0 = z1-roufns(i,j)
z1drou = z1/roufns(i,j)
z1droup = (z1+roufns(i,j))/roufns(i,j)
IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
! Unstable case: See equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
bulkri = MAX (bulkri,-10.0)
sb =bulkri/prantl0l
qb=oned9*(c1l+c2l*sb*sb)
pb=oned54*(c3l+c4l*sb*sb)
qb3pb2=qb**3-pb*pb
c7 = (z1*dzt*LOG(z1droup)*LOG(z1droup))/(dz0*dz0*LOG(ztdroup))
IF( qb3pb2 >= 0.0 ) THEN
sqrtqb = SQRT(qb)
tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )
thetab=ACOS(tempan)
stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5l)
ELSE
tb =(SQRT(-qb3pb2)+ABS(pb))**oned3
stabp =c7*(-(tb+qb/tb)+c5l)
END IF
!
!-----------------------------------------------------------------------
!
! According to equation (15) in Byun (1990).
!
!-----------------------------------------------------------------------
!
c8=gammaml*stabp
y1=SQRT(1.0 - c8)
y0=SQRT(1.0 - c8/ztdrou)
psih=2.0*LOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
! Compute c_pt via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
c_pt(i,j)=kv / (prantl0l*(LOG(ztdroup)-psih))
ELSE
!
!-----------------------------------------------------------------------
!
! Stable case: See Louis et al (1981).
!
!-----------------------------------------------------------------------
!
a=kv/LOG(ztdroup)
b=5.0
d=5.0
c=SQRT(1.0+d*bulkri)
c_pt(i,j) = SQRT(1.0+2.0*b*bulkri/c)
c_pt(i,j) = a*c_pt(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))
END IF
END DO
END DO
RETURN
END SUBROUTINE cptca
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CPTCWTR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cptcwtra(nx,ny,nz,iend,jend,zp,soiltyp,wspd, & 2
ptsfc,pt1,c_ptwtrneu,c_ptwtr,i1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Compute c_ptwtr (the product of ustar and ptstar) at sea case.
! Note: a temperature scale defined as surface heat flux divided
! by the friction velocity.
!
! The quantity c_ptwtr is used by the subroutine SFCFLX to obtain
! surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
! AUTHORS: V. Wong and X. Song
! 8/04/1994
!
! MODIFICATION HISTORY:
!
! 2/27/95 (V.W. and X.S.)
!
! 1/12/96 (V.W. and X.S.)
! Changed the calculation related to zo over the sea.
! Added kvwtr to denote the Von Karman constant over the sea;
! Set a lower limiter for zo, zolimit, and an upper limiter for z1, z1limit.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! iend i-index where evaluation ends.
! jend j-index where evaluation ends.
!
! zp The physical height coordinate defined at w-point of
! staggered grid.
! wspd Surface wind speed (m/s), defined as
! sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
! ptsfc Potential temperature at the ground level (K)
! pt1 Potential temperature at the 1st scalar point above
! ground level, (K)
!
!
! OUTPUT:
!
! c_ptwtr The product of ustar and ptstar. ptstar is temperature
! scale (K), defined by surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz ! The number grid points in 3 directions
INTEGER :: iend ! i-index where evaluation ends
INTEGER :: jend ! j-index where evaluation ends
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of staggered grid.
INTEGER :: soiltyp(nx,ny) ! Soil type at each point.
REAL :: wspd (nx,ny) ! Surface wind speed (m/s)
REAL :: ptsfc(nx,ny) ! Potential temperature at the ground level
! (K)
REAL :: pt1 (nx,ny) ! Potential temperature at the 1st scalar
! point above ground level, (K)
REAL :: c_ptwtrneu(nx,ny) ! Product of ustar and ptstar
REAL :: c_ptwtr(nx,ny) ! Product of ustar and ptstar
!
!-----------------------------------------------------------------------
!
! Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j
REAL :: z1 ! The height of 1st scalar point above
! ground level (m)
REAL :: zo ! Defined as surface momentum roughness
! length
REAL :: zt ! Defined as surface heat transfer
! roughness length
REAL :: dzo ! z1-zo
REAL :: dzt ! z1-zt
REAL :: z1dzo ! z1/zo
REAL :: z1dzt ! z1/zt
REAL :: xcdn ! Cdn (sea)
REAL :: xcdh ! Hot-wired value of Cdh (sea)
REAL :: bulkri ! Bulk Richardson number
REAL :: stabp ! Monin-Obukhov STABility Parameter
! (zeta)
REAL :: y1,y0,psih ! Intermediate parameters needed
REAL :: qb3pb2
REAL :: c7,c8
REAL :: sb,qb,pb,thetab,tb ! during computations
REAL :: a,b,c,d
REAL :: tempan,sqrtqb
REAL :: z10,zo0,zt0,dzo0,dzt0
REAL :: i1
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
INCLUDE 'globcst.inc'
INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
PRINT *,'starting PTCPWTR'
xcdh = 1.1E-3
DO j=1,jend
DO i=1,iend
IF ( soiltyp(i,j) == 13) THEN
! IF (soiltyp(i,j,1).eq.13.and.stypfrct(i,j,1).ge.0.5) THEN
! calculate bulk Ri with the lowest level above the ground
xcdn = (0.4+0.07*wspd(i,j)) * 1.e-3
z10 = 0.5*(zp(i,j,3)-zp(i,j,2))
z10 = MIN(z10,z1limit)
!
! calculate zo and zt
!
zo0 =0.032*xcdn*wspd(i,j)*wspd(i,j)/9.8
zo0 = MAX(zo0,zolimit)
zt0 = z10 * EXP( -kv*SQRT(xcdn)/(prantl0w*xcdh) )
dzo0 = z10-zo0
dzt0 = z10-zt0
bulkri =(g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dzo0*dzo0/ &
(dzt0*wspd(i,j)*wspd(i,j))
! fined Ri bulk calcuation
xcdn = (0.4+0.07*wspd(i,j)) * 1.e-3
IF (i1 == 1) THEN
z1 = 0.5*(zp(i,j,3)-zp(i,j,2))
ELSE
z1=2.0
END IF
z1 = MIN(z1,z1limit)
!
! calculate zo and zt
!
zo =0.032*xcdn*wspd(i,j)*wspd(i,j)/9.8
zo = MAX(zo,zolimit)
zt = z1 * EXP( -kv*SQRT(xcdn)/(prantl0w*xcdh) )
dzo = z1-zo
z1dzo = z1/zo
dzt = z1-zt
z1dzt = z1/zt
IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
! Unstable case: A modified formulation, which is similar to
! equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
bulkri = MAX (bulkri,-10.0)
sb =bulkri/prantl0w
qb=oned9*(c1w+c2w*sb*sb)
pb=oned54*(c3w+c4w*sb*sb)
qb3pb2=qb**3-pb*pb
c7 = (z1*dzt/(dzo*dzo))*(ALOG(z1dzo)*ALOG(z1dzo)/ALOG(z1dzt))
IF( qb3pb2 >= 0.0 ) THEN
sqrtqb = SQRT(qb)
tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )
thetab=ACOS(tempan)
stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5w)
ELSE
tb =(SQRT(-qb3pb2)+ABS(pb))**oned3
stabp =c7*(-(tb+qb/tb)+c5w)
END IF
!
!-----------------------------------------------------------------------
!
! According to a modified equation, which is similar to equation (14)
! in Byun (1990).
!
!-----------------------------------------------------------------------
!
c8=gammamw*stabp
y1=SQRT(1.0 - c8)
y0=SQRT(1.0 - c8/z1dzt)
psih=2.0*ALOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
! Compute ptstar via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
c_ptwtr(i,j)=kv / (prantl0w*(ALOG(z1dzt)-psih))
ELSE
!
!-----------------------------------------------------------------------
!
! Stable case: With the modified formulation in Louis et al (1981).
!
!-----------------------------------------------------------------------
!
a=kv*kv/(prantl0w*LOG(z1dzo)*LOG(z1dzt))
b=5.0
d=5.0
c=SQRT(1.0+d*bulkri)
c_ptwtr(i,j) = SQRT(1.0+2.0*b*bulkri/c)
c_ptwtr(i,j) = a*c_ptwtr(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))
END IF
END IF
END DO
END DO
RETURN
END SUBROUTINE cptcwtra
SUBROUTINE psfsl(nx,ny,t2,qv2,p2,z2,z1,psf,psl) 1
IMPLICIT NONE
INTEGER :: nx,ny
REAL :: t2(nx,ny),qv2(nx,ny),p2(nx,ny)
REAL :: z2(nx,ny),z1(nx,ny)
REAL :: psf(nx,ny),psl(nx,ny) !output
REAL :: gammam,rd,rv,tsh,zsh,rog,fvirt
REAL :: tvsf,tvsl,tv2
INTEGER :: i,j
gammam=-6.5E-3
rd=2.8705E2
rv=4.615E2
rog=rd/9.8
fvirt=rv/rd-1.0
zsh=75.0
tsh=290.66
DO j=1,ny-1
DO i=1,nx-1
tv2=t2(i,j)*(1.+fvirt*qv2(i,j))
tvsf=tv2-gammam*(z2(i,j)-z1(i,j))
psf(i,j)=p2(i,j)*EXP(2.*(z2(i,j)-z1(i,j))/(rog*(tvsf+tv2)))
IF (z1(i,j) > zsh) THEN
tvsl=tvsf-gammam*z1(i,j)
IF (tvsl > tsh) THEN
IF (tvsf > tsh) THEN
tvsl=tsh-0.005*(tvsf-tsh)**2
ELSE
tvsl=tsh
END IF
END IF
ELSE
tvsl=tvsf
END IF
psl(i,j)=psf(i,j)*EXP(2.0*z1(i,j)/(rog*(tvsf+tvsl)))
psf(i,j)=-ALOG(psf(i,j))
END DO
END DO
RETURN
END SUBROUTINE psfsl
!
SUBROUTINE pintp(nx,ny,nz,s3d,nlgp3d,s2d,INDEX,nlgp01, & 15
t3d,qv3d,z2,zsf,nlgpsf)
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: i,j,k
REAL :: s3d(nx,ny,nz),nlgp3d(nx,ny,nz),s2d(nx,ny)
INTEGER :: INDEX
REAL :: nlgp01
REAL :: t3d(nx,ny,nz),qv3d(nx,ny,nz)
REAL :: z2(nx,ny),zsf(nx,ny)
REAL :: nlgpsf(nx,ny)
REAL :: gammam,rd,rv,tsh,zsh,rog,fvirt,gammas
REAL :: tvsf,tvsl,tv2
REAL :: tvkm,part
gammam=-6.5E-3
rd=2.8705E2
rv=4.615E2
rog=rd/9.8
fvirt=rv/rd-1.0
zsh=75.0
tsh=290.66
! print *,'starting PINTP'
DO j=1,ny-1
DO i=1,nx-1
IF(nlgp01 <= nlgp3d(i,j,2)) GO TO 11
IF(nlgp01 >= nlgp3d(i,j,nz-1)) GO TO 12
DO k=2,nz-2
IF(nlgp01 >= nlgp3d(i,j,k).AND.nlgp01 < nlgp3d(i,j,k+1)) GO TO 15
END DO
11 k=2
GO TO 15
12 k=nz-2
GO TO 15
15 s2d(i,j)=s3d(i,j,k)+(s3d(i,j,k+1)-s3d(i,j,k))* &
(nlgp01-nlgp3d(i,j,k))/(nlgp3d(i,j,k+1)-nlgp3d(i,j,k))
! if (i.eq.1.and.j.eq.1) then
! print *,i,j,k,nlgp01
! print *,s2d(i,j),s3d(i,j,k),s3d(i,j,k+1),s3d(i,j,k),
! : nlgp01,nlgp3d(i,j,k),nlgp3d(i,j,k+1),nlgp3d(i,j,k)
! stop
! endif
IF (nlgp01 <= nlgpsf(i,j)) THEN
IF (INDEX == 1.OR.INDEX == 2) THEN
s2d(i,j)=s3d(i,j,2)
ELSE IF (INDEX == 3) THEN
tvsf=t3d(i,j,2)*(1.+fvirt*qv3d(i,j,2))-gammam*(z2(i,j)-zsf(i,j))
IF (zsf(i,j) > zsh) THEN
tvsl=tvsf-gammam*zsf(i,j)
IF (tvsl > tsh) THEN
IF (tvsf > tsh) THEN
tvsl=tsh-0.005*(tvsf-tsh)**2
ELSE
tvsl=tsh
END IF
END IF
gammas=(tvsf-tvsl)/zsf(i,j)
ELSE
gammas=0.0
END IF
part=rog*(nlgp01-nlgpsf(i,j))
s2d(i,j)=zsf(i,j)+tvsf*part/(1-0.5*gammas*part)
END IF
ELSE IF (nlgp01 > nlgpsf(i,j).AND.nlgp01 <= nlgp3d(i,j,2)) THEN
IF (INDEX == 1) THEN
s2d(i,j)=s3d(i,j,2)
ELSE IF (INDEX == 2) THEN
s2d(i,j)=s2d(i,j)
ELSE IF (INDEX == 3) THEN
s2d(i,j)=s3d(i,j,2)+(zsf(i,j)-s3d(i,j,2))* &
(nlgp01-nlgp3d(i,j,2))/(nlgpsf(i,j)-nlgp3d(i,j,2))
END IF
ELSE IF (nlgp01 > nlgp3d(i,j,nz-1)) THEN
IF (INDEX == 1) THEN
s2d(i,j)=s3d(i,j,nz-1)
ELSE IF (INDEX == 2) THEN
s2d(i,j)=s3d(i,j,nz-1)
ELSE IF (INDEX == 3) THEN
tvkm=t3d(i,j,nz-1)*(1.0+fvirt*qv3d(i,j,nz-1))
s2d(i,j)=s3d(i,j,nz-1)+tvkm*rog*(nlgp01-nlgp3d(i,j,nz-1))
END IF
END IF
END DO
END DO
RETURN
END SUBROUTINE pintp
SUBROUTINE cllmc(nx,ny,nz,u,v,qv,algpsf,zp,rhobar,algpzs,llmc, & 1
x,y, tem1,tem2)
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: u(nx,ny,nz),v(nx,ny,nz),qv(nx,ny,nz)
REAL :: algpzs(nx,ny,nz),zp(nx,ny,nz),rhobar(nx,ny,nz)
REAL :: algpsf(nx,ny)
REAL :: llmc(nx-1,ny-1)
REAL :: tem1(nx,ny),tem2(nx,ny)
REAL :: x(nx,ny),y(nx,ny)
REAL :: ubar,vbar,qbar,mbar,dz
REAL :: algpst,pst,psf
INTEGER :: i,j,k,kt
PRINT *,'starting CLLMC'
DO j=1,ny-1
DO i=1,nx-1
llmc(i,j)=0
END DO
END DO
DO j=1,ny-1
DO i=1,nx-1
psf=EXP(-algpsf(i,j))
pst=psf-5000.0
algpst=-ALOG(pst)
ubar=0.0
vbar=0.0
qbar=0.0
mbar=0.0
DO k=2,nz-1
IF (algpst < algpzs(i,j,k+1)) THEN
kt=k
GO TO 300
END IF
END DO
PRINT *,'error, check CLLMC'
STOP
300 DO k=2,kt
dz=zp(i,j,k+1)-zp(i,j,k)
IF (k == kt) THEN
dz=dz*(algpst-0.5*(algpzs(i,j,k)+algpzs(i,j,k-1))) &
/(0.5*(algpzs(i,j,k+1)-algpzs(i,j,k-1)))
END IF
ubar=ubar+u(i,j,k)*rhobar(i,j,k)*dz
vbar=vbar+v(i,j,k)*rhobar(i,j,k)*dz
qbar=qbar+qv(i,j,k)*rhobar(i,j,k)*dz
mbar=mbar+rhobar(i,j,k)*dz
END DO
ubar=ubar/mbar
vbar=vbar/mbar
qbar=qbar/mbar
tem1(i,j)=ubar*qbar
tem2(i,j)=vbar*qbar
END DO
END DO
DO j=2,ny-2
DO i=2,nx-2
llmc(i,j)=-1000.0*((tem1(i+1,j)-tem1(i-1,j))/(x(i+1,j)-x(i-1,j)) &
+(tem2(i,j+1)-tem2(i,j-1))/(y(i,j+1)-y(i,j-1)))
END DO
END DO
PRINT *,'finishing CLLMC'
RETURN
END SUBROUTINE cllmc
SUBROUTINE cape_ncep(nx,ny,nz,tk1,qv1,p1,cape,cins,li, & 1,1
tk,qv,p,pint,plcl,peql,p1d,t1d,qv1d,lmh, &
item1,item2,item3,item4,item5,item6,item7, & ! work arrays
rtem1,rtem2,rtem3,rtem4,rtem5,rtem6,rtem7) ! work arrays
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: tk1(nx,ny,nz),qv1(nx,ny,nz),p1(nx,ny,nz)
REAL :: cape(nx,ny),cins(nx,ny),li(nx,ny)
REAL :: tk(nx,ny,nz),qv(nx,ny,nz),p(nx,ny,nz),pint(nx,ny,nz)
REAL :: plcl(nx,ny),peql(nx,ny)
REAL :: p1d(nx,ny),t1d(nx,ny),qv1d(nx,ny)
INTEGER :: lmh(nx,ny)
INTEGER :: i,j,l,lk,lmhij,itype
INTEGER :: item1(nx,ny),item2(nx,ny),item3(nx,ny),item4(nx,ny), &
item5(nx,ny),item6(nx,ny),item7(nx,ny)
REAL :: rtem1(nx,ny),rtem2(nx,ny),rtem3(nx,ny),rtem4(nx,ny), &
rtem5(nx,ny),rtem6(nx,ny),rtem7(nx,ny)
!----------------------------------------------------------------------
!***********************************************************************
!--------------PREPARATIONS---------------------------------------------
DO j=1,ny
DO i=1,nx
lmh(i,j)=nz-2
lmhij=nz-2
pint(i,j,1)=2500.
DO l=1,lmhij
lk=lmhij-l+2
p(i,j,l)=p1(i,j,lk)
tk(i,j,l)=tk1(i,j,lk)
qv(i,j,l)=qv1(i,j,lk)
pint(i,j,l+1)=(p1(i,j,lk)+p1(i,j,lk-1))*0.5
END DO
END DO
END DO
!
itype=1
CALL calcape(itype,p1d,t1d,qv1d, &
tk,qv,p,pint,lmh,nx,ny,nz-2,nz-1, &
cape,cins,li,plcl,peql, &
item1,item2,item3,item4,item5,item6,item7, & ! work arrays
rtem1,rtem2,rtem3,rtem4,rtem5,rtem6,rtem7) ! work arrays
!
! WRITE(6,*) ' CAPE= ',CAPE(I,J)
! WRITE(6,*) ' CIN= ',CINS(I,J)
! WRITE(6,*) ' LCL= ',PLCL(I,J)
! WRITE(6,*) ' EQ LVL= ',PEQL(I,J)
RETURN
END SUBROUTINE cape_ncep
SUBROUTINE calcape(itype,p1d,t1d,q1d, & 1,4
t,q,p,pint,lmh,im,jm,lm,lp1, &
cape,cins,li,plcl,peql, &
ieql,iptb,ithtb,ilres,jlres,ihres,jhres, & ! work arrays
lcl,tpar,pp,qq,psp,thesp,ape) ! work arrays
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . .
! SUBPROGRAM: CALCAPE COMPUTES CAPE AND CINS
! PRGRMMR: TREADON ORG: W/NMC2 DATE: 93-02-10
!
! ABSTRACT:
!
! THIS ROUTINE COMPUTES CAPE AND CINS GIVEN TEMPERATURE,
! PRESSURE, AND SPECIFIC HUMIDTY. IN "STORM AND CLOUD
! DYNAMICS" (1989, ACADEMIC PRESS) COTTON AND ANTHES DEFINE
! CAPE (EQUATION 9.16, P501) AS
!
! EL
! CAPE = SUM G * LN(THETAP/THETAA) DZ
! LCL
!
! WHERE,
! EL = EQUILIBRIUM LEVEL,
! LCL = LIFTING CONDENSTATION LEVEL,
! G = GRAVITATIONAL ACCELERATION,
! THETAP = LIFTED PARCEL POTENTIAL TEMPERATURE,
! THETAA = AMBIENT POTENTIAL TEMPERATURE.
!
! NOTE THAT THE INTEGRAND LN(THETAP/THETAA) APPROXIMATELY
! EQUALS (THETAP-THETAA)/THETAA. THIS RATIO IS OFTEN USED
! IN THE DEFINITION OF CAPE/CINS.
!
! TWO TYPES OF CAPE/CINS CAN BE COMPUTED BY THIS ROUTINE. THE
! SUMMATION PROCESS IS THE SAME FOR BOTH CASES. WHAT DIFFERS
! IS THE DEFINITION OF THE PARCEL TO LIFT. FOR ITYPE=1 THE
! PARCEL WITH THE WARMEST THETA-E IN A DPBND PASCAL LAYER ABOVE
! THE MODEL SURFACE IS LIFTED. THE ARRAYS P1D, T1D, AND Q1D
! ARE NOT USED. FOR ITYPE=2 THE ARRAYS P1D, T1D, AND Q1D
! DEFINE THE PARCEL TO LIFT IN EACH COLUMN. BOTH TYPES OF
! CAPE/CINS MAY BE COMPUTED IN A SINGLE EXECUTION OF THE POST
! PROCESSOR.
!
! THIS ALGORITHM PROCEEDS AS FOLLOWS.
! FOR EACH COLUMN,
! (1) INITIALIZE RUNNING CAPE AND CINS SUM TO 0.0
! (2) COMPUTE TEMPERATURE AND PRESSURE AT THE LCL USING
! LOOK UP TABLE (PTBL). USE EITHER PARCEL THAT GIVES
! MAX THETAE IN LOWEST DPBND ABOVE GROUND (ITYPE=1)
! OR GIVEN PARCEL FROM T1D,Q1D,...(ITYPE=2).
! (3) COMPUTE THE TEMP OF A PARCEL LIFTED FROM THE LCL.
! WE KNOW THAT THE PARCEL'S
! EQUIVALENT POTENTIAL TEMPERATURE (THESP) REMAINS
! CONSTANT THROUGH THIS PROCESS. WE CAN
! COMPUTE TPAR USING THIS KNOWLEDGE USING LOOK
! UP TABLE (SUBROUTINE TTBLEX).
! (4) FIND THE EQUILIBRIUM LEVEL. THIS IS DEFINED AS THE
! HIGHEST POSITIVELY BUOYANT LAYER.
! (IF THERE IS NO POSITIVELY BUOYANT LAYER, CAPE/CINS
! WILL BE ZERO)
! (5) COMPUTE CAPE/CINS.
! (A) COMPUTE THETAP. WE KNOW TPAR AND P.
! (B) COMPUTE THETAA. WE KNOW T AND P.
! (6) ADD G*(THETAP-THETAA)*DZ TO THE RUNNING CAPE OR CINS SUM.
! (A) IF THETAP > THETAA, ADD TO THE CAPE SUM.
! (B) IF THETAP < THETAA, ADD TO THE CINS SUM.
! (7) ARE WE AT EQUILIBRIUM LEVEL?
! (A) IF YES, STOP THE SUMMATION.
! (B) IF NO, CONTIUNUE THE SUMMATION.
! (8) ENFORCE LIMITS ON CAPE AND CINS (I.E. NO NEGATIVE CAPE)
!
! PROGRAM HISTORY LOG:
! 93-02-10 RUSS TREADON
! 93-06-19 RUSS TREADON - GENERALIZED ROUTINE TO ALLOW FOR
! TYPE 2 CAPE/CINS CALCULATIONS.
! 94-09-23 MIKE BALDWIN - MODIFIED TO USE LOOK UP TABLES
! INSTEAD OF COMPLICATED EQUATIONS.
! 94-10-13 MIKE BALDWIN - MODIFIED TO CONTINUE CAPE/CINS CALC
! UP TO AT HIGHEST BUOYANT LAYER.
! 98-03-25 MIKE BALDWIN - STAND ALONE VERSION FOR SAMEX
!
! USAGE: CALL CALCAPE(ITYPE,P1D,T1D,Q1D,
! & T,Q,P,PINT,LMH,IM,JM,LM,
! & CAPE,CINS)
! INPUT ARGUMENT LIST:
! ITYPE,P1D,T1D,Q1D,T,Q,P,PINT,LMH,IM,JM,LM
!
! CONTROL VARIABLES:
! ITYPE - INTEGER FLAG SPECIFYING HOW PARCEL TO LIFT IS
! IDENTIFIED. SEE COMMENTS ABOVE.
! IM - 1ST HORIZONTAL DIMENSION OF ARRAY
! JM - 2ND HORIZONTAL DIMENSION OF ARRAY
! LM - VERTICAL DIMENSION OF ARRAY
! LMH - ARRAY THAT CONTAINS THE NUMBER OF ABOVE GROUND LEVELS
! AT EACH GRID POINT DIMENSIONED (IM,JM)
!
! 3-D ARRAYS OF ATMOSPHERIC CONDITIONS:
! MID-LAYER VARIABLES:
! T - TEMP (K) DIMENSIONED TO (IM,JM,LM)
! Q - SPEC HUM (kg/kg) DIMENSIONED TO (IM,JM,LM)
! P - PRESSURE (Pa) DIMENSIONED TO (IM,JM,LM)
!
! LAYER-INTERFACE VARIABLES:
! PINT - PRESSURE (Pa) DIMENSIONED TO (IM,JM,LM+1)
!
! 2-D ARRAYS IF SENDING IN A SPECIFIC PARCEL TO LIFT (ITYPE=2):
! P1D(IM,JM) - ARRAY OF PRESSURE (Pa) OF PARCELS TO LIFT.
! T1D(IM,JM) - ARRAY OF TEMPERATURE (K) OF PARCELS TO LIFT.
! Q1D(IM,JM) - ARRAY OF SPECIFIC HUMIDITY (kg/kg) OF PARCELS TO LIFT.
!
! OUTPUT ARGUMENT LIST:
! CAPE(IM,JM) - CONVECTIVE AVAILABLE POTENTIAL ENERGY (J/KG)
! CINS(IM,JM) - CONVECTIVE INHIBITION (J/KG)
! PLCL(IM,JM) - PRESSURE OF LCL (Pa)
! PEQL(IM,JM) - PRESSURE OF EQ LEVEL (Pa)
!
! OUTPUT FILES:
! STDOUT - RUN TIME STANDARD OUT.
!
! SUBPROGRAMS CALLED:
! UTILITIES:
! TTBLEQ - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P
! TTBLE1 - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P
! TTBLEX - LOOKUP TABLE ROUTINE TO GET T FROM THETAE AND P
!
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN
! MACHINE : CRAY Y-MP
!$$$
!
!
!
! INCLUDE/SET PARAMETERS. CONSTANTS ARE FROM BOLTON (MWR, 1980).
!----------------------------------------------------------------------
! PARAMETER (IM=IMM,JM=JMM,LM=LMM,LP1=LM+1)
PARAMETER &
(h10e5=100000.e0 &
, epsq=2.e-12 &
, g=9.8E0,cp=1004.6E0,capa=0.28589641E0,rog=287.04/g)
!----------------------------------------------------------------------
PARAMETER &
(itb=076,jtb=134,itbq=152,jtbq=440)
PARAMETER (dpbnd=70.e2, &
elivw=2.72E6,elocp=elivw/cp)
!
! DECLARE VARIABLES.
!
INTEGER :: ieql(im,jm)
INTEGER :: lcl(im,jm)
!
REAL :: p1d(im,jm),t1d(im,jm),q1d(im,jm)
REAL :: cape(im,jm),cins(im,jm),li(im,jm)
REAL :: plcl(im,jm),peql(im,jm)
REAL :: tpar(im,jm,lm)
DIMENSION &
p (im,jm,lm),pint (im,jm,lm+1) &
,t(im,jm,lm),q(im,jm,lm) &
,lmh (im,jm)
DIMENSION &
qs0 (jtb),sqs (jtb),the0 (itb),sthe (itb) &
, the0q(itbq),stheq(itbq) &
,ptbl (itb,jtb),ttbl (jtb,itb),ttblq(jtbq,itbq)
!
DIMENSION &
iptb (im,jm),ithtb (im,jm) &
,pp (im,jm) &
,qq (im,jm) &
,psp (im,jm),thesp (im,jm) &
,ilres (im*jm),jlres (im*jm),ihres (im*jm),jhres (im*jm) &
,ape (im,jm,lm)
!
!
DATA itabl/0/,thl/210./,plq/70000./,ptt/2500./
!
SAVE ptbl,ttbl,rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0, &
ttblq,rdpq,rdtheq,plq,stheq,the0q
!-----------------------------------------------------------------------
!
! SET UP LOOKUP TABLE
!
IF (itabl == 0) THEN
CALL table1(ptbl,ttbl,ptt, &
rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
CALL tableq(ttblq,rdpq,rdtheq,plq,thl,stheq,the0q)
itabl=1
END IF
imjm=im*jm
!-----------------------------------------------------------------------
!
!
!**************************************************************
! START CALCAPE HERE.
!
!
! COMPUTE CAPE/CINS
!
! WHICH IS: THE SUM FROM THE LCL TO THE EQ LEVEL OF
! G * (LN(THETAP) - LN(THETAA)) * DZ
!
! (POSITIVE AREA FOR CAPE, NEGATIVE FOR CINS)
!
! WHERE:
! THETAP IS THE PARCEL THETA
! THETAA IS THE AMBIENT THETA
! DZ IS THE THICKNESS OF THE LAYER
!
! USING LCL AS LEVEL DIRECTLY BELOW SATURATION POINT
! AND EQ LEVEL IS THE HIGHEST POSITIVELY BUOYANT LEVEL.
!
! IEQL = EQ LEVEL
!
! INITIALIZE CAPE AND CINS ARRAYS
!
DO j = 1,jm
DO i = 1,im
cape(i,j) = 0.0
cins(i,j) = 0.0
thesp(i,j)= 0.0
ieql(i,j) = lm+1
END DO
END DO
DO l=1,lm
DO j = 1,jm
DO i = 1,im
tpar(i,j,l)= 0.0
IF (l <= lmh(i,j)) THEN
apests=p(i,j,l)
ape(i,j,l)=(h10e5/apests)**capa
IF( END IF
END DO
END DO
END DO
!
! NOTE THAT FOR TYPE 1 CAPE/CINS ARRAYS P1D, T1D, Q1D
! ARE DUMMY ARRAYS.
!
!-------FOR ITYPE=1-----------------------------------------------------
!---------FIND MAXIMUM THETA E LAYER IN LOWEST DPBND ABOVE GROUND-------
!-------FOR ITYPE=2-----------------------------------------------------
!---------FIND THETA E LAYER OF GIVEN T1D, Q1D, P1D---------------------
!--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
DO kb=1,lm
IF (itype == 1.OR.kb == 1) THEN
DO j=1,jm
DO i=1,im
lmhk =lmh(i,j)
psfck =p(i,j,lmhk)
pkl = p(i,j,kb)
IF (itype == 2.OR.(pkl >= psfck-dpbnd.AND.pkl <= psfck)) THEN
IF (itype == 1) THEN
tbtk =t(i,j,kb)
qbtk =q(i,j,kb)
apebtk =ape(i,j,kb)
!DHOU added the following statement to calculate the temperature of
! the parcil if it diabatically goes to 500hpa. This is stored in li.
li(i,j)=tbtk*(50000.0/p(i,j,kb))**capa
ELSE
pkl =p1d(i,j)
tbtk =t1d(i,j)
qbtk =q1d(i,j)
apebtk =(h10e5/pkl)**capa
END IF
!--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
tthbtk =tbtk*apebtk
tthk =(tthbtk-thl)*rdth
qq (i,j)=tthk-AINT(tthk)
ittbk =INT(tthk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
IF(ittbk < 1) THEN
ittbk =1
qq (i,j)=0.0
END IF
IF(ittbk >= jtb) THEN
ittbk =jtb-1
qq (i,j)=0.0
END IF
!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
bqs00k=qs0(ittbk)
sqs00k=sqs(ittbk)
bqs10k=qs0(ittbk+1)
sqs10k=sqs(ittbk+1)
!--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
bqk =(bqs10k-bqs00k)*qq(i,j)+bqs00k
sqk =(sqs10k-sqs00k)*qq(i,j)+sqs00k
tqk =(qbtk-bqk)/sqk*rdq
pp (i,j)=tqk-AINT(tqk)
iq =INT(tqk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
IF(iq < 1) THEN
iq =1
pp (i,j)=0.0
END IF
IF(iq >= itb) THEN
iq =itb-1
pp (i,j)=0.0
END IF
!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
it=ittbk
p00k=ptbl(iq ,it )
p10k=ptbl(iq+1,it )
p01k=ptbl(iq ,it+1)
p11k=ptbl(iq+1,it+1)
!--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
tpspk=p00k+(p10k-p00k)*pp(i,j)+(p01k-p00k)*qq(i,j) &
+(p00k-p10k-p01k+p11k)*pp(i,j)*qq(i,j)
apespk=(h10e5/tpspk)**capa
tthesk=tthbtk*EXP(elocp*qbtk*apespk/tthbtk)
!--------------CHECK FOR MAXIMUM THETA E--------------------------------
IF(tthesk > thesp(i,j)) THEN
psp (i,j)=tpspk
thesp(i,j)=tthesk
END IF
END IF
END DO
END DO
END IF
END DO
!-----CHOOSE LAYER DIRECTLY BELOW PSP AS LCL AND------------------------
!-----ENSURE THAT THE LCL IS ABOVE GROUND.------------------------------
!-------(IN SOME RARE CASES FOR ITYPE=2, IT IS NOT)---------------------
DO l=1,lm
DO j=1,jm
DO i=1,im
IF (l < lmh(i,j)) THEN
pkl = p(i,j,l)
IF (pkl < psp(i,j)) THEN
lcl(i,j)=l+1
plcl(i,j)=p(i,j,l+1)
END IF
END IF
END DO
END DO
END DO
!-----------------------------------------------------------------------
!---------FIND TEMP OF PARCEL LIFTED ALONG MOIST ADIABAT (TPAR)---------
!-----------------------------------------------------------------------
DO ivi=1,lm
l=lp1-ivi
!--------------SCALING PRESSURE & TT TABLE INDEX------------------------
knuml=0
knumh=0
DO j=1,jm
DO i=1,im
pkl = p(i,j,l)
IF(l <= lcl(i,j)) THEN
IF(pkl < plq)THEN
knuml=knuml+1
ilres(knuml)=i
jlres(knuml)=j
ELSE
knumh=knumh+1
ihres(knumh)=i
jhres(knumh)=j
END IF
END IF
END DO
END DO
!***
!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PLQ
!**
IF(knuml > 0)THEN
CALL ttblex(tpar(1,1,l),ttbl,im,jm,imjm,itb,jtb &
, knuml,ilres,jlres &
, p(1,1,l),pl,qq,pp,rdp,the0,sthe &
, rdthe,thesp,iptb,ithtb)
END IF
!***
!*** COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PLQ
!**
IF(knumh > 0)THEN
CALL ttblex(tpar(1,1,l),ttblq,im,jm,imjm,itbq,jtbq &
, knumh,ihres,jhres &
, p(1,1,l),plq,qq,pp,rdpq,the0q,stheq &
, rdtheq,thesp,iptb,ithtb)
END IF
!------------SEARCH FOR EQ LEVEL----------------------------------------
DO n=1,knumh
i=ihres(n)
j=jhres(n)
IF(tpar(i,j,l) > t(i,j,l)) THEN
ieql(i,j)=l
peql(i,j)=p(i,j,l)
END IF
END DO
DO n=1,knuml
i=ilres(n)
j=jlres(n)
IF(tpar(i,j,l) > t(i,j,l)) THEN
ieql(i,j)=l
peql(i,j)=p(i,j,l)
END IF
END DO
!-----------------------------------------------------------------------
END DO
!-----------------------------------------------------------------------
! DHOU added this section to calculate LI (Lifted Index=T500-Tp500)
DO j=1,jm
DO i=1,im
DO k=1,lm
IF (p(i,j,k+1) >= 50000.0.AND.p(i,j,k) < 50000.0) THEN
t500=t(i,j,k+1)-(ALOG( *(t(i,j,k+1)-t(i,j,k))/(ALOG( IF (psp(i,j) < 50000.0) THEN
tp500=li(i,j)
ELSE
tp500=tpar(i,j,k+1)-(ALOG( *(tpar(i,j,k+1)-tpar(i,j,k))/(ALOG( END IF
li(i,j)=t500-tp500
CYCLE
END IF
END DO
END DO
END DO
!------------COMPUTE CAPE AND CINS--------------------------------------
DO j=1,jm
DO i=1,im
lclk=lcl(i,j)
ieqk=ieql(i,j)
DO l=ieqk,lclk
presk=p(i,j,l)
dp=pint(i,j,l+1)-pint(i,j,l)
dzkl=t(i,j,l)*(q(i,j,l)*0.608+1.0)*rog*dp/presk
thetap=tpar(i,j,l)*(h10e5/presk)**capa
thetaa=t(i,j,l)*(h10e5/presk)**capa
IF (thetap < thetaa) &
cins(i,j)=cins(i,j)+g*(ALOG(thetap)-ALOG(thetaa))*dzkl
IF (thetap > thetaa) &
cape(i,j)=cape(i,j)+g*(ALOG(thetap)-ALOG(thetaa))*dzkl
END DO
END DO
END DO
!
! ENFORCE LOWER LIMIT OF 0.0 ON CAPE AND UPPER
! LIMIT OF 0.0 ON CINS.
!
DO j = 1,jm
DO i = 1,im
cape(i,j) = AMAX1(0.0,cape(i,j))
cins(i,j) = AMIN1(cins(i,j),0.0)
END DO
END DO
!
! END OF ROUTINE.
!
RETURN
END SUBROUTINE calcape
SUBROUTINE table1(ptbl,ttbl,pt & 1,2
, rdq,rdth,rdp,rdthe,pl,thl,qs0,sqs,sthe,the0)
! ******************************************************************
! * *
! * GENERATE VALUES FOR LOOK-UP TABLES USED IN CONVECTION *
! * *
! ******************************************************************
PARAMETER &
(itb=076,jtb=134)
PARAMETER &
(thh=365.,ph=105000. &
, pq0=379.90516,a1=610.78,a2=17.2693882,a3=273.16,a4=35.86 &
, r=287.04,cp=1004.6,eliwv=2.683E6,eps=1.e-9)
DIMENSION &
ptbl (itb,jtb),ttbl (jtb,itb),qsold (jtb),pold (jtb) &
, qs0 (jtb),sqs (jtb),qsnew (jtb) &
, y2p (jtb),app (jtb),aqp (jtb),pnew (jtb) &
, told (jtb),theold(jtb),the0 (itb),sthe (itb) &
, y2t (jtb),thenew(jtb),apt (jtb),aqt (jtb),tnew (jtb)
!--------------COARSE LOOK-UP TABLE FOR SATURATION POINT----------------
kthm=jtb
kpm=itb
kthm1=kthm-1
kpm1=kpm-1
!
pl=pt
!
dth=(thh-thl)/REAL(kthm-1)
dp =(ph -pl )/REAL(kpm -1)
!
rdth=1./dth
rdp=1./dp
rdq=kpm-1
!
th=thl-dth
!-----------------------------------------------------------------------
DO kth=1,kthm
th=th+dth
p=pl-dp
DO kp=1,kpm
p=p+dp
ape=(100000./p)**(r/cp)
qsold(kp)=pq0/p*EXP(a2*(th-a3*ape)/(th-a4*ape))
pold(kp)=p
END DO
!
qs0k=qsold(1)
sqsk=qsold(kpm)-qsold(1)
qsold(1 )=0.
qsold(kpm)=1.
!
DO kp=2,kpm1
qsold(kp)=(qsold(kp)-qs0k)/sqsk
!
IF((qsold(kp)-qsold(kp-1)) < eps) qsold(kp)=qsold(kp-1)+eps
!
END DO
!
qs0(kth)=qs0k
sqs(kth)=sqsk
!-----------------------------------------------------------------------
qsnew(1 )=0.
qsnew(kpm)=1.
dqs=1./REAL(kpm-1)
!
DO kp=2,kpm1
qsnew(kp)=qsnew(kp-1)+dqs
END DO
!
y2p(1 )=0.
y2p(kpm )=0.
!
CALL spline(jtb,kpm,qsold,pold,y2p,kpm,qsnew,pnew,app,aqp)
!
DO kp=1,kpm
ptbl(kp,kth)=pnew(kp)
END DO
!-----------------------------------------------------------------------
END DO
!--------------COARSE LOOK-UP TABLE FOR T(P) FROM CONSTANT THE----------
p=pl-dp
DO kp=1,kpm
p=p+dp
th=thl-dth
DO kth=1,kthm
th=th+dth
ape=(100000./p)**(r/cp)
qs=pq0/p*EXP(a2*(th-a3*ape)/(th-a4*ape))
told(kth)=th/ape
theold(kth)=th*EXP(eliwv*qs/(cp*told(kth)))
END DO
!
the0k=theold(1)
sthek=theold(kthm)-theold(1)
theold(1 )=0.
theold(kthm)=1.
!
DO kth=2,kthm1
theold(kth)=(theold(kth)-the0k)/sthek
!
IF((theold(kth)-theold(kth-1)) < eps) theold(kth)=theold(kth-1) + eps
!
END DO
!
the0(kp)=the0k
sthe(kp)=sthek
!-----------------------------------------------------------------------
thenew(1 )=0.
thenew(kthm)=1.
dthe=1./REAL(kthm-1)
rdthe=1./dthe
!
DO kth=2,kthm1
thenew(kth)=thenew(kth-1)+dthe
END DO
!
y2t(1 )=0.
y2t(kthm)=0.
!
CALL spline(jtb,kthm,theold,told,y2t,kthm,thenew,tnew,apt,aqt)
!
DO kth=1,kthm
ttbl(kth,kp)=tnew(kth)
END DO
!-----------------------------------------------------------------------
END DO
!
RETURN
END SUBROUTINE table1
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE tableq(ttblq & 1,1
, rdp,rdthe,pl,thl,sthe,the0)
! ******************************************************************
! * *
! * GENERATE VALUES FOR FINER LOOK-UP TABLES USED *
! * IN CONVECTION *
! * *
! ******************************************************************
PARAMETER &
(itb=152,jtb=440)
PARAMETER &
(thh=325.,ph=105000. &
, pq0=379.90516,a1=610.78,a2=17.2693882,a3=273.16,a4=35.86 &
, r=287.04,cp=1004.6,eliwv=2.683E6,eps=1.e-9)
DIMENSION &
ttblq (jtb,itb) &
, told (jtb),theold(jtb),the0 (itb),sthe (itb) &
, y2t (jtb),thenew(jtb),apt (jtb),aqt (jtb),tnew (jtb)
!--------------COARSE LOOK-UP TABLE FOR SATURATION POINT----------------
kthm=jtb
kpm=itb
kthm1=kthm-1
kpm1=kpm-1
!
dth=(thh-thl)/REAL(kthm-1)
dp =(ph -pl )/REAL(kpm -1)
!
rdp=1./dp
th=thl-dth
!--------------COARSE LOOK-UP TABLE FOR T(P) FROM CONSTANT THE----------
p=pl-dp
DO kp=1,kpm
p=p+dp
th=thl-dth
DO kth=1,kthm
th=th+dth
ape=(100000./p)**(r/cp)
qs=pq0/p*EXP(a2*(th-a3*ape)/(th-a4*ape))
told(kth)=th/ape
theold(kth)=th*EXP(eliwv*qs/(cp*told(kth)))
END DO
!
the0k=theold(1)
sthek=theold(kthm)-theold(1)
theold(1 )=0.
theold(kthm)=1.
!
DO kth=2,kthm1
theold(kth)=(theold(kth)-the0k)/sthek
!
IF((theold(kth)-theold(kth-1)) < eps) theold(kth)=theold(kth-1) + eps
!
END DO
!
the0(kp)=the0k
sthe(kp)=sthek
!-----------------------------------------------------------------------
thenew(1 )=0.
thenew(kthm)=1.
dthe=1./REAL(kthm-1)
rdthe=1./dthe
!
DO kth=2,kthm1
thenew(kth)=thenew(kth-1)+dthe
END DO
!
y2t(1 )=0.
y2t(kthm)=0.
!
CALL spline(jtb,kthm,theold,told,y2t,kthm,thenew,tnew,apt,aqt)
!
DO kth=1,kthm
ttblq(kth,kp)=tnew(kth)
END DO
!-----------------------------------------------------------------------
END DO
!
RETURN
END SUBROUTINE tableq
SUBROUTINE ttblex(tref,ttbl,im,jm,imxjm,itb,jtb & 4
, knum,iarr,jarr &
, pijl,pl,qq,pp,rdp,the0 &
, sthe,rdthe,thesp,iptb,ithtb)
! ******************************************************************
! * *
! * EXTRACT TEMPERATURE OF THE MOIST ADIABAT FROM *
! * THE APPROPRIATE TTBL *
! * *
! ******************************************************************
!----------------------------------------------------------------------
! PARAMETER(IM=1,JM=1,LM=50,IMXJM=IM*JM)
DIMENSION &
tref(im,jm),ttbl(jtb,itb),iarr(imxjm),jarr(imxjm) &
,pijl(im,jm),qq(im,jm),pp(im,jm),the0(itb) &
,sthe(itb),thesp(im,jm),iptb(im,jm),ithtb(im,jm)
!-----------------------------------------------------------------------
DO kk=1,knum
!--------------SCALING PRESSURE & TT TABLE INDEX------------------------
i=iarr(kk)
j=jarr(kk)
pk=pijl(i,j)
tpk=(pk-pl)*rdp
qq(i,j)=tpk-AINT(tpk)
iptb(i,j)=INT(tpk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
IF(iptb(i,j) < 1)THEN
iptb(i,j)=1
qq(i,j)=0.
END IF
IF(iptb(i,j) >= itb)THEN
iptb(i,j)=itb-1
qq(i,j)=0.
END IF
!--------------BASE AND SCALING FACTOR FOR THE--------------------------
iptbk=iptb(i,j)
bthe00k=the0(iptbk)
sthe00k=sthe(iptbk)
bthe10k=the0(iptbk+1)
sthe10k=sthe(iptbk+1)
!--------------SCALING THE & TT TABLE INDEX-----------------------------
bthk=(bthe10k-bthe00k)*qq(i,j)+bthe00k
sthk=(sthe10k-sthe00k)*qq(i,j)+sthe00k
tthk=(thesp(i,j)-bthk)/sthk*rdthe
pp(i,j)=tthk-AINT(tthk)
ithtb(i,j)=INT(tthk)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
IF(ithtb(i,j) < 1)THEN
ithtb(i,j)=1
pp(i,j)=0.
END IF
IF(ithtb(i,j) >= jtb)THEN
ithtb(i,j)=jtb-1
pp(i,j)=0.
END IF
!--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------
ith=ithtb(i,j)
ip=iptb(i,j)
t00k=ttbl(ith ,ip )
t10k=ttbl(ith+1,ip )
t01k=ttbl(ith ,ip+1)
t11k=ttbl(ith+1,ip+1)
!--------------PARCEL TEMPERATURE-------------------------------------
tref(i,j)=(t00k+(t10k-t00k)*pp(i,j)+(t01k-t00k)*qq(i,j) &
+(t00k-t10k-t01k+t11k)*pp(i,j)*qq(i,j))
END DO
!
RETURN
END SUBROUTINE ttblex
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
SUBROUTINE spline(jtb,nold,xold,yold,y2,nnew,xnew,ynew,p,q) 6
! ******************************************************************
! * *
! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
! * *
! * PROGRAMER Z. JANJIC, YUGOSLAV FED. HYDROMET. INST., BEOGRAD *
! * *
! * *
! * *
! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
! * SPECIFIED. *
! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
! * AND LE XOLD(NOLD). *
! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
! * *
! ******************************************************************
DIMENSION &
xold(jtb),yold(jtb),y2(jtb),p(jtb),q(jtb) &
,xnew(jtb),ynew(jtb)
!-----------------------------------------------------------------------
noldm1=nold-1
!
dxl=xold(2)-xold(1)
dxr=xold(3)-xold(2)
dydxl=(yold(2)-yold(1))/dxl
dydxr=(yold(3)-yold(2))/dxr
rtdxc=.5/(dxl+dxr)
!
p(1)= rtdxc*(6.*(dydxr-dydxl)-dxl*y2(1))
q(1)=-rtdxc*dxr
!
IF(nold == 3) GO TO 700
!-----------------------------------------------------------------------
k=3
!
100 dxl=dxr
dydxl=dydxr
dxr=xold(k+1)-xold(k)
dydxr=(yold(k+1)-yold(k))/dxr
dxc=dxl+dxr
den=1./(dxl*q(k-2)+dxc+dxc)
!
p(k-1)= den*(6.*(dydxr-dydxl)-dxl*p(k-2))
q(k-1)=-den*dxr
!
k=k+1
IF(k < nold) GO TO 100
!-----------------------------------------------------------------------
700 k=noldm1
!
200 y2(k)=p(k-1)+q(k-1)*y2(k+1)
!
k=k-1
IF(k > 1) GO TO 200
!-----------------------------------------------------------------------
k1=1
!
300 xk=xnew(k1)
!
DO k2=2,nold
IF(xold(k2) <= xk) CYCLE
kold=k2-1
GO TO 450
END DO
ynew(k1)=yold(nold)
GO TO 600
!
450 IF(k1 == 1) GO TO 500
IF(k == kold) GO TO 550
!
500 k=kold
!
y2k=y2(k)
y2kp1=y2(k+1)
dx=xold(k+1)-xold(k)
rdx=1./dx
!
!VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
! WRITE(6,5000) K,Y2K,Y2KP1,DX,RDX,YOLD(K),YOLD(K+1)
!5000 FORMAT(' K=',I4,' Y2K=',E12.4,' Y2KP1=',E12.4,' DX=',E12.4,' RDX='
! 2,E12.4,' YOK=',E12.4,' YOP1=',E12.4)
!AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
ak=.1666667*rdx*(y2kp1-y2k)
bk=.5*y2k
ck=rdx*(yold(k+1)-yold(k))-.1666667*dx*(y2kp1+y2k+y2k)
!
550 x=xk-xold(k)
xsq=x*x
!
ynew(k1)=ak*xsq*x+bk*xsq+ck*x+yold(k)
!
600 CONTINUE
k1=k1+1
IF(k1 <= nnew) GO TO 300
!-----------------------------------------------------------------------
RETURN
END SUBROUTINE spline
SUBROUTINE newmap(nqi,nqj,latn,lonn,dx, & 1,4
qswcorn,qswcorw,trulat,trulon,mapproj)
IMPLICIT NONE
INTEGER :: nqi,nqj,i,j
REAL :: latn(nqi,nqj),lonn(nqi,nqj)
REAL :: qswcorn,qswcorw,dx
REAL :: trulat(2),trulon
INTEGER :: mapproj
!
REAL :: x,y,x1,y1
REAL :: lat,lon
!
CALL setmapr(mapproj,1.0,trulat,trulon)
lat=qswcorn
lon=qswcorw
CALL lltoxy(1,1,lat,lon,x,y)
CALL setorig(1,x,y)
DO j=1,nqj
DO i=1,nqi
x1=(i-1)*dx*1000.0
y1=(j-1)*dx*1000.0
CALL xytoll(1,1,x1,y1,lat,lon)
latn(i,j)=lat
lonn(i,j)=lon
END DO
END DO
RETURN
END SUBROUTINE newmap
SUBROUTINE d2intrpa(xold,yold,nxold,nyold, & 1
xnew,ynew,nxnew,nynew,intrpl, &
aold1,anew1, &
aold2,anew2, &
aold3,anew3, &
aold4,anew4, &
aold5,anew5, &
aold6,anew6, &
aold7,anew7, &
aold8,anew8, &
aold9,anew9, &
aold10,anew10, &
aold11,anew11, &
aold12,anew12, &
aold13,anew13, &
aold14,anew14, &
aold15,anew15, &
aold16,anew16, &
aold17,anew17, &
aold18,anew18, &
aold19,anew19, &
aold20,anew20, &
aold21,anew21, &
aold22,anew22, &
aold23,anew23, &
aold24,anew24, &
aold25,anew25, &
aold26,anew26, &
aold27,anew27, &
aold28,anew28, &
aold29,anew29)
IMPLICIT NONE
INTEGER :: nxold,nyold, nxnew,nynew,intrpl
REAL :: xold(nxold,nyold),xnew(nxnew,nynew)
REAL :: yold(nxold,nyold),ynew(nxnew,nynew)
REAL :: aold1(nxold,nyold),anew1(nxnew,nynew)
REAL :: aold2(nxold,nyold),anew2(nxnew,nynew)
REAL :: aold3(nxold,nyold),anew3(nxnew,nynew)
REAL :: aold4(nxold,nyold),anew4(nxnew,nynew)
REAL :: aold5(nxold,nyold),anew5(nxnew,nynew)
REAL :: aold6(nxold,nyold),anew6(nxnew,nynew)
REAL :: aold7(nxold,nyold),anew7(nxnew,nynew)
REAL :: aold8(nxold,nyold),anew8(nxnew,nynew)
REAL :: aold9(nxold,nyold),anew9(nxnew,nynew)
REAL :: aold10(nxold,nyold),anew10(nxnew,nynew)
REAL :: aold11(nxold,nyold),anew11(nxnew,nynew)
REAL :: aold12(nxold,nyold),anew12(nxnew,nynew)
REAL :: aold13(nxold,nyold),anew13(nxnew,nynew)
REAL :: aold14(nxold,nyold),anew14(nxnew,nynew)
REAL :: aold15(nxold,nyold),anew15(nxnew,nynew)
REAL :: aold16(nxold,nyold),anew16(nxnew,nynew)
REAL :: aold17(nxold,nyold),anew17(nxnew,nynew)
REAL :: aold18(nxold,nyold),anew18(nxnew,nynew)
REAL :: aold19(nxold,nyold),anew19(nxnew,nynew)
REAL :: aold20(nxold,nyold),anew20(nxnew,nynew)
REAL :: aold21(nxold,nyold),anew21(nxnew,nynew)
REAL :: aold22(nxold,nyold),anew22(nxnew,nynew)
REAL :: aold23(nxold,nyold),anew23(nxnew,nynew)
REAL :: aold24(nxold,nyold),anew24(nxnew,nynew)
REAL :: aold25(nxold,nyold),anew25(nxnew,nynew)
REAL :: aold26(nxold,nyold),anew26(nxnew,nynew)
REAL :: aold27(nxold,nyold),anew27(nxnew,nynew)
REAL :: aold28(nxold,nyold),anew28(nxnew,nynew)
REAL :: aold29(nxold,nyold),anew29(nxnew,nynew)
INTEGER :: i,j,i1,i2,j1,j2
INTEGER :: io,jo
REAL :: xs,ys,s1,s2,s3,s4,sgrid
PRINT *,'starting D2INTRPA'
PRINT *,nxold,nyold
PRINT *,nxnew,nynew
PRINT *,xold(1,1),yold(1,1),xnew(1,1),ynew(1,1)
PRINT *,xold(1,nyold),yold(1,nyold),xnew(1,nynew),ynew(1,nynew)
PRINT *,xold(nxold,nyold),yold(nxold,nyold), &
xnew(nxnew,nynew),ynew(nxnew,nynew)
PRINT *,xold(nxold,1),yold(nxold,1),xnew(nxnew,1),ynew(nxnew,1)
DO j=1,nynew
DO i=1,nxnew
! Direct assignment in the case of intrpl=0 (No intepolation)
IF (intrpl == 0) THEN
io=i+1
jo=j+1
anew1(i,j)=aold1(io,jo)
anew2(i,j)=aold2(io,jo)
anew3(i,j)=aold3(io,jo)
anew4(i,j)=aold4(io,jo)
anew5(i,j)=aold5(io,jo)
anew6(i,j)=aold6(io,jo)
anew7(i,j)=aold7(io,jo)
anew8(i,j)=aold8(io,jo)
anew9(i,j)=aold9(io,jo)
anew10(i,j)=aold10(io,jo)
anew11(i,j)=aold11(io,jo)
anew12(i,j)=aold12(io,jo)
anew13(i,j)=aold13(io,jo)
anew14(i,j)=aold14(io,jo)
anew15(i,j)=aold15(io,jo)
anew16(i,j)=aold16(io,jo)
anew17(i,j)=aold17(io,jo)
anew18(i,j)=aold18(io,jo)
anew19(i,j)=aold19(io,jo)
anew20(i,j)=aold20(io,jo)
anew21(i,j)=aold21(io,jo)
anew22(i,j)=aold22(io,jo)
anew23(i,j)=aold23(io,jo)
anew24(i,j)=aold24(io,jo)
anew25(i,j)=aold25(io,jo)
anew26(i,j)=aold26(io,jo)
anew27(i,j)=aold27(io,jo)
anew28(i,j)=aold28(io,jo)
anew29(i,j)=aold29(io,jo)
CYCLE
END IF
! Interpolation, first finding the old points surroding the new point
xs=xnew(i,j)
ys=ynew(i,j)
IF (xs < xold(1,1)) THEN
i1=1
i2=2
ELSE IF (xs >= xold(nxold,1)) THEN
i1=nxold-1
i2=nxold
ELSE
DO io=1,nxold-1
IF (xs >= xold(io,1).AND.xs < xold(io+1,1)) THEN
i1=io
i2=io+1
EXIT
END IF
END DO
205 CONTINUE
END IF
IF (ys < yold(1,1)) THEN
j1=1
j2=2
ELSE IF (ys >= yold(1,nyold)) THEN
j1=nyold-1
j2=nyold
ELSE
DO jo=1,nyold-1
IF (ys >= yold(1,jo).AND.ys < yold(1,jo+1)) THEN
j1=jo
j2=jo+1
EXIT
END IF
END DO
305 CONTINUE
END IF
DO jo=1,nyold
DO io=1,nxold
! dirtect assignment for a new point, very close to an old point.
IF (ABS(xs-xold(io,jo)) < 1.0E2.AND. ABS(ys-yold(io,jo)) < 1.0E2) THEN
anew1(i,j)=aold1(io,jo)
anew2(i,j)=aold2(io,jo)
anew3(i,j)=aold3(io,jo)
anew4(i,j)=aold4(io,jo)
anew5(i,j)=aold5(io,jo)
anew6(i,j)=aold6(io,jo)
anew7(i,j)=aold7(io,jo)
anew8(i,j)=aold8(io,jo)
anew9(i,j)=aold9(io,jo)
anew10(i,j)=aold10(io,jo)
anew11(i,j)=aold11(io,jo)
anew12(i,j)=aold12(io,jo)
anew13(i,j)=aold13(io,jo)
anew14(i,j)=aold14(io,jo)
anew15(i,j)=aold15(io,jo)
anew16(i,j)=aold16(io,jo)
anew17(i,j)=aold17(io,jo)
anew18(i,j)=aold18(io,jo)
anew19(i,j)=aold19(io,jo)
anew20(i,j)=aold20(io,jo)
anew21(i,j)=aold21(io,jo)
anew22(i,j)=aold22(io,jo)
anew23(i,j)=aold23(io,jo)
anew24(i,j)=aold24(io,jo)
anew25(i,j)=aold25(io,jo)
anew26(i,j)=aold26(io,jo)
anew27(i,j)=aold27(io,jo)
anew28(i,j)=aold28(io,jo)
anew29(i,j)=aold29(io,jo)
GO TO 405
END IF
END DO
END DO
! interpolation
s1=SQRT((xs-xold(i1,j1))**2+(ys-yold(i1,j1))**2)
s2=SQRT((xs-xold(i2,j1))**2+(ys-yold(i2,j1))**2)
s3=SQRT((xs-xold(i2,j2))**2+(ys-yold(i2,j2))**2)
s4=SQRT((xs-xold(i1,j2))**2+(ys-yold(i1,j2))**2)
s1=1.0/s1
s2=1.0/s2
s3=1.0/s3
s4=1.0/s4
sgrid=1.0/(s1+s2+s3+s4)
anew1(i,j)=(aold1(i1,j1)*s1+aold1(i2,j1)*s2 &
+aold1(i2,j2)*s3+aold1(i1,j2)*s4)*sgrid
anew2(i,j)=(aold2(i1,j1)*s1+aold2(i2,j1)*s2 &
+aold2(i2,j2)*s3+aold2(i1,j2)*s4)*sgrid
anew3(i,j)=(aold3(i1,j1)*s1+aold3(i2,j1)*s2 &
+aold3(i2,j2)*s3+aold3(i1,j2)*s4)*sgrid
anew4(i,j)=(aold4(i1,j1)*s1+aold4(i2,j1)*s2 &
+aold4(i2,j2)*s3+aold4(i1,j2)*s4)*sgrid
anew5(i,j)=(aold5(i1,j1)*s1+aold5(i2,j1)*s2 &
+aold5(i2,j2)*s3+aold5(i1,j2)*s4)*sgrid
anew6(i,j)=(aold6(i1,j1)*s1+aold6(i2,j1)*s2 &
+aold6(i2,j2)*s3+aold6(i1,j2)*s4)*sgrid
anew7(i,j)=(aold7(i1,j1)*s1+aold7(i2,j1)*s2 &
+aold7(i2,j2)*s3+aold7(i1,j2)*s4)*sgrid
anew8(i,j)=(aold8(i1,j1)*s1+aold8(i2,j1)*s2 &
+aold8(i2,j2)*s3+aold8(i1,j2)*s4)*sgrid
anew9(i,j)=(aold9(i1,j1)*s1+aold9(i2,j1)*s2 &
+aold9(i2,j2)*s3+aold9(i1,j2)*s4)*sgrid
anew10(i,j)=(aold10(i1,j1)*s1+aold10(i2,j1)*s2 &
+aold10(i2,j2)*s3+aold10(i1,j2)*s4)*sgrid
anew11(i,j)=(aold11(i1,j1)*s1+aold11(i2,j1)*s2 &
+aold11(i2,j2)*s3+aold11(i1,j2)*s4)*sgrid
anew12(i,j)=(aold12(i1,j1)*s1+aold12(i2,j1)*s2 &
+aold12(i2,j2)*s3+aold12(i1,j2)*s4)*sgrid
anew13(i,j)=(aold13(i1,j1)*s1+aold13(i2,j1)*s2 &
+aold13(i2,j2)*s3+aold13(i1,j2)*s4)*sgrid
anew14(i,j)=(aold14(i1,j1)*s1+aold14(i2,j1)*s2 &
+aold14(i2,j2)*s3+aold14(i1,j2)*s4)*sgrid
anew15(i,j)=(aold15(i1,j1)*s1+aold15(i2,j1)*s2 &
+aold15(i2,j2)*s3+aold15(i1,j2)*s4)*sgrid
anew16(i,j)=(aold16(i1,j1)*s1+aold16(i2,j1)*s2 &
+aold16(i2,j2)*s3+aold16(i1,j2)*s4)*sgrid
anew17(i,j)=(aold17(i1,j1)*s1+aold17(i2,j1)*s2 &
+aold17(i2,j2)*s3+aold17(i1,j2)*s4)*sgrid
anew18(i,j)=(aold18(i1,j1)*s1+aold18(i2,j1)*s2 &
+aold18(i2,j2)*s3+aold18(i1,j2)*s4)*sgrid
anew19(i,j)=(aold19(i1,j1)*s1+aold19(i2,j1)*s2 &
+aold19(i2,j2)*s3+aold19(i1,j2)*s4)*sgrid
anew20(i,j)=(aold20(i1,j1)*s1+aold20(i2,j1)*s2 &
+aold20(i2,j2)*s3+aold20(i1,j2)*s4)*sgrid
anew21(i,j)=(aold21(i1,j1)*s1+aold21(i2,j1)*s2 &
+aold21(i2,j2)*s3+aold21(i1,j2)*s4)*sgrid
anew22(i,j)=(aold22(i1,j1)*s1+aold22(i2,j1)*s2 &
+aold22(i2,j2)*s3+aold22(i1,j2)*s4)*sgrid
anew23(i,j)=(aold23(i1,j1)*s1+aold23(i2,j1)*s2 &
+aold23(i2,j2)*s3+aold23(i1,j2)*s4)*sgrid
anew24(i,j)=(aold24(i1,j1)*s1+aold24(i2,j1)*s2 &
+aold24(i2,j2)*s3+aold24(i1,j2)*s4)*sgrid
anew25(i,j)=(aold25(i1,j1)*s1+aold25(i2,j1)*s2 &
+aold25(i2,j2)*s3+aold25(i1,j2)*s4)*sgrid
anew26(i,j)=(aold26(i1,j1)*s1+aold26(i2,j1)*s2 &
+aold26(i2,j2)*s3+aold26(i1,j2)*s4)*sgrid
anew27(i,j)=(aold27(i1,j1)*s1+aold27(i2,j1)*s2 &
+aold27(i2,j2)*s3+aold27(i1,j2)*s4)*sgrid
anew28(i,j)=(aold28(i1,j1)*s1+aold28(i2,j1)*s2 &
+aold28(i2,j2)*s3+aold28(i1,j2)*s4)*sgrid
anew29(i,j)=(aold29(i1,j1)*s1+aold29(i2,j1)*s2 &
+aold29(i2,j2)*s3+aold29(i1,j2)*s4)*sgrid
! if (i.eq.1.and.j.eq.1.or.i.eq.nxnew.and.j.eq.nynew) then
! print *,i1,j1,i2,j2
! print *,anew1(i,j),aold1(i1,j1),aold1(i2,j1),
! : aold1(i2,j2),aold1(i1,j2)
! print *,anew2(i,j),aold2(i1,j1),aold2(i2,j1),
! : aold2(i2,j2),aold2(i1,j2)
! endif
405 CONTINUE
END DO
END DO
PRINT *,'finishing D2INTRPA'
RETURN
END SUBROUTINE d2intrpa
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE FILZERO ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE filzero( a, n) 5
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Fill vector a with zeros.
!
!-----------------------------------------------------------------------
!
! AUTHOR: St. Paul, Dead Sea Scrolls
!
! MODIFICATIONS:
! 6/09/92 Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: n
REAL :: a(n)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO i=1,n
a(i)=0.0
END DO
RETURN
END SUBROUTINE filzero
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE IFILZERO ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE ifilzero( ia, n )
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Fill vector a with zeros.
!
!-----------------------------------------------------------------------
!
! AUTHOR: St. Paul, Dead Sea Scrolls
!
! MODIFICATIONS:
! 6/09/92 Added full documentation (K. Brewster)
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: n
INTEGER :: ia(n)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO i=1,n
ia(i)=0
END DO
RETURN
END SUBROUTINE ifilzero
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TEMPER ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE temper ( nx,ny,nz,theta, ppert, pbar, t ) 6
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Using a version of Poisson's formula, calculate temperature.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Joe Bradley
! 12/05/91
!
! MODIFICATIONS:
! Modified by Ming Xue so that arrays are only defined at
! one time level.
! 6/09/92 Added full documentation and phycst include file for
! rddcp=Rd/Cp (K. Brewster)
!
!-----------------------------------------------------------------------
!
! INPUT:
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
!
! theta Potential temperature (degrees Kelvin)
! ppert Perturbation pressure (Pascals)
! pbar Base state pressure (Pascals)
!
! OUTPUT:
!
! t Temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER :: nx,ny,nz
!
REAL :: theta(nx,ny,nz) ! potential temperature (degrees Kelvin)
REAL :: ppert(nx,ny,nz) ! perturbation pressure (Pascals)
REAL :: pbar (nx,ny,nz) ! base state pressure (Pascals)
!
REAL :: t (nx,ny,nz) ! temperature (degrees Kelvin)
!
!-----------------------------------------------------------------------
!
! Include file
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Calculate the temperature using Poisson's formula.
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
t(i,j,k) = theta(i,j,k) * &
(((ppert(i,j,k) + pbar(i,j,k)) / p0) ** rddcp)
END DO
END DO
END DO
RETURN
END SUBROUTINE temper
SUBROUTINE calender(year,mon,day,hour,ds,wkd) 1
! find the calender date for forecast initiated at hourZ,dayth of mon, Year
! with a lead time of ds seconds. the day of week is also found.
! effective only for date later than 1st Jan 1998, not exceeding 28 Feb. 2100
IMPLICIT NONE
INTEGER :: year,mon,day,hour,MIN,sec,ds,wkd
INTEGER :: mon1,dday
INTEGER :: y0,wkd0 !on Jan. 1st,the year of y0,the weekday index is wkd0
INTEGER :: i,dindex
y0=1998
wkd0=5
! cross minute hour and day
sec=0
MIN=0
sec=sec+ds
MIN=MIN+(sec-MOD(sec,60))/60
sec=MOD(sec,60)
hour=hour+(MIN-MOD(MIN,60))/60
MIN=MOD(MIN,60)
dday=day+(hour-MOD(hour,24))/24
hour=MOD(hour,24)
ds=MIN
! cross the month
IF ( (mon == 1.OR.mon == 3.OR.mon == 5.OR.mon == 7 &
.OR.mon == 8.OR.mon == 10.OR.mon == 12).AND.dday > 31) THEN
dday=dday-31
mon1=mon+1
ELSE IF ( (mon == 4.OR.mon == 6.OR.mon == 9.OR.mon == 11) &
.AND.dday > 30) THEN
dday=dday-30
mon1=mon+1
ELSE IF (mon == 2) THEN
IF (MOD(year,4) == 0.AND.dday > 29) THEN
dday=dday-29
mon1=mon+1
ELSE IF (MOD(year,4) /= 0.AND.dday > 28) THEN
dday=dday-28
mon1=mon+1
ELSE
mon1=mon
END IF
ELSE
mon1=mon
END IF
day=dday
mon=mon1
!* cross the year
IF (mon > 12) THEN
mon=mon-12
year=year+1
END IF
!* day of week
dindex=-1
DO i=y0,year-1
dindex=dindex+365
IF (MOD(i,4) == 0) dindex=dindex+1
END DO
DO i=1,mon-1
IF (i == 1.OR.i == 3.OR.i == 5.OR.i == 7 &
.OR.i == 8.OR.i == 10.OR.i == 12) THEN
dindex=dindex+31
ELSE IF (i == 4.OR.i == 6.OR.i == 9.OR.i == 11) THEN
dindex=dindex+30
ELSE IF (i == 2.AND.MOD(year,4) == 0) THEN
dindex=dindex+29
ELSE IF (i == 2.AND.MOD(year,4) /= 0) THEN
dindex=dindex+28
END IF
END DO
dindex=dindex+day
wkd=wkd0+dindex
PRINT *,wkd
wkd=MOD(wkd,7)
IF (wkd == 0) THEN
wkd=7
END IF
RETURN
END SUBROUTINE calender