SUBROUTINE postcore(nx,ny,nz,nzsoil, nstyps,ctime, & 2,60
hisfile,lenfil,outheader,lenbin,gemoutheader,lengem, &
icape,iaccu,iascii,ibinary,igempak, &
x,y,z,zp,uprt ,vprt ,wprt ,ptprt, pprt , &
qvprt, qc, qr, qi, qs, qh, &
ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, &
soiltyp,stypfrct,vegtyp,lai,roufns,veg, &
tsoil,qsoil,wetcanp, &
raing,rainc, &
tem1,tem2,tem3,tem4,tem5)
IMPLICIT NONE
!-----------------------------------------------------------------------
! Include files:
!-----------------------------------------------------------------------
! INCLUDE 'indtflg.inc'
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.inc'
! INCLUDE 'enscv.inc'
INCLUDE 'mp.inc'
INTEGER :: nx, ny, nz, nzsoil
INTEGER :: nstyps ! Maximum number of soil types in each
! grid box
!-----------------------------------------------------------------------
!
! Arrays to be read in:
!
!-----------------------------------------------------------------------
!
REAL :: x (nx) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL :: y (ny) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL :: z (nz) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL :: zp (nx,ny,nz) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s)
REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s)
REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s)
REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature
REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal)
REAL :: qvprt (nx,ny,nz)
REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg)
REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg)
REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg)
REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg)
REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg)
REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2)
REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s)
REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s)
REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s)
REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K)
REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal)
REAL :: rhobar(nx,ny,nz) ! Base state density rhobar
REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity
! (kg/kg)
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 :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3)
REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
REAL :: raing(nx,ny) ! Grid supersaturation rain
REAL :: rainc(nx,ny) ! Cumulus convective rain
!
!-----------------------------------------------------------------------
!
! Arrays derived from the read-in arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: xs (:) ! x-coor of sacalar point in computational
REAL, ALLOCATABLE :: ys (:) ! y-coor of sacalar point in computational
REAL, ALLOCATABLE :: zs (:,:,:) ! z-coor of sacalar point in computational
! space (m)
REAL, ALLOCATABLE :: zps (:,:,:) ! z-coor of sacalar point in physical
! space (m)
REAL, ALLOCATABLE :: zpc (:,:,:) ! z-coor of sacalar point in physical
! space (km)
REAL, ALLOCATABLE :: zagl (:,:,:) ! AGL height of sacalar point (km)
REAL, ALLOCATABLE :: u (:,:,:) ! Total u-velocity (m/s)
REAL, ALLOCATABLE :: v (:,:,:) ! Total v-velocity (m/s)
REAL, ALLOCATABLE :: w (:,:,:) ! Total w-velocity (m/s)
! from that of base state atmosphere (K)
REAL, ALLOCATABLE :: pt (:,:,:) !
REAL, ALLOCATABLE :: qv (:,:,:) !
REAL, ALLOCATABLE :: raint(:,:) ! Total rain (rainc+raing)
REAL, ALLOCATABLE :: hterain(:,:) ! The height of the terrain
!-----------------------------------------------------------------------
!
! Arrays for plots on constant pressure levels
!
!-----------------------------------------------------------------------
!
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 :: 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 :: thet(:,:) ! theta_E (K)
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 convergenc
REAL :: lat1,lon1,lat2,lon2
REAL :: pw1d,totw,rho
! REAL :: pi
!-----------------------------------------------------------------------
! OUTPUT arrays
!-----------------------------------------------------------------------
INTEGER :: houroday, minohour, dayoweek, dateomon, month1, year1
REAL, ALLOCATABLE :: hgtp(:,:,:),uwndp(:,:,:),vwndp(:,:,:),wwndp(:,:,:)
REAL, ALLOCATABLE :: tmpp(:,:,:),sphp(:,:,:)
REAL, ALLOCATABLE :: uwndh(:,:,:),vwndh(:,:,:),wwndh(:,:,:)
REAL, ALLOCATABLE :: psf(:,:),mslp(:,:),vort500(:,:)
REAL, ALLOCATABLE :: temp2m(:,:),dewp2m(:,:),qv2m(:,:)
REAL, ALLOCATABLE :: u10m(:,:),v10m(:,:),hgtsfc(:,:)
REAL, ALLOCATABLE :: refl1km(:,:),refl4km(:,:),cmpref(:,:)
REAL, ALLOCATABLE :: accppt(:,:),convppt(:,:)
REAL, ALLOCATABLE :: tpptold(:,:),cpptold(:,:)
REAL, ALLOCATABLE :: cape(:,:),mcape(:,:),cin(:,:),mcin(:,:),lcl(:,:)
REAL, ALLOCATABLE :: srh01(:,:),srh03(:,:),uh25(:,:),sh01(:,:),sh06(:,:)
REAL, ALLOCATABLE :: thck(:,:),li(:,:),brn(:,:),pw(:,:),f2(:,:)
REAL :: rlevel
INTEGER, PARAMETER :: destination = 0
INTEGER :: source
!
!-----------------------------------------------------------------------
! nprgem: number of pressure levels for GEMPAK file.
!
! iprgem: actual pressure levels defined by the user (mb)
!
INTEGER :: nprgem
!
PARAMETER (nprgem = 5)
!
INTEGER :: iprgem(nprgem)
!
DATA iprgem / 850, 700, 600, 500, 250/
!
!-----------------------------------------------------------------------
! nhtgem: number of height levels for GEMPAK file.
!
! ihtgem: actual height levels defined by the user (km)
!
INTEGER :: nhtgem
!
PARAMETER (nhtgem = 6)
!
INTEGER :: ihtgem(nhtgem)
!
DATA ihtgem / 1, 2, 3, 4, 5, 6/
!
!-----------------------------------------------------------------------
!
! Temporary work arrays for general use
!
!-----------------------------------------------------------------------
!
REAL :: tem1(nx,ny,nz)
REAL :: tem2(nx,ny,nz)
REAL :: tem3(nx,ny,nz)
REAL :: tem4(nx,ny,nz)
REAL :: tem5(nx,ny,nz)
! REAL :: vtwo1(nx,ny),vtwo2(nx,ny),vtwo3(nx,ny)
! REAL :: vtwo4(nx,ny),vtwo5(nx,ny),vtwo6(nx,ny)
! REAL :: qvsfc (nx,ny,nz) ! Soil Skin effective qv
REAL, ALLOCATABLE :: pt2m(:,:),t700(:,:)
!
!-----------------------------------------------------------------------
!
! 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,sumh
!
CHARACTER (LEN=256) :: filename
CHARACTER (LEN=256) :: hisfile
INTEGER :: lenfil,lenbin,lengem
INTEGER :: ihr
CHARACTER (LEN=4) :: ihrc
CHARACTER (LEN=40) :: tppname,cppname,varname,varunit
INTEGER :: istatus, iret, ierr
INTEGER :: i,j,k, itags,itagr
REAL :: ctime
REAL :: truelat(2)
! INTEGER :: nf
INTEGER :: lenstag, linstit, lmodel
REAL :: amax, amin
REAL :: ave1, ave2, ave3
!
!-----------------------------------------------------------------------
INTEGER :: icape,iaccu,iascii,ibinary,igempak
CHARACTER (LEN=*) :: outheader,gemoutheader
INTEGER :: dmp_out_joined
INTEGER :: klev
LOGICAL :: tppExist,cppExist
REAL, PARAMETER :: pi = 3.14159265
REAL, PARAMETER :: gamma=.0065 ! std lapse rate per meter
REAL, PARAMETER :: ex2=5.2558774 ! g/R/gamma
REAL, PARAMETER :: ex1=0.1903643 ! R*gamma/g
REAL :: t00, p00
!
!-----------------------------------------------------------------------
!
! Variables for mpi jobs
!
!-----------------------------------------------------------------------
INTEGER, PARAMETER :: fzone = 3
INTEGER :: nxlg, nylg ! global domain
INTEGER :: ii,jj,ia,ja
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directives for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(hgtp(nx,ny,nprgem),STAT=istatus)
hgtp=0
ALLOCATE(uwndp(nx,ny,nprgem),STAT=istatus)
uwndp=0
ALLOCATE(vwndp(nx,ny,nprgem),STAT=istatus)
vwndp=0
ALLOCATE(wwndp(nx,ny,nprgem),STAT=istatus)
wwndp=0
ALLOCATE(tmpp(nx,ny,nprgem),STAT=istatus)
tmpp=0
ALLOCATE(sphp(nx,ny,nprgem),STAT=istatus)
sphp=0
ALLOCATE(uwndh(nx,ny,nhtgem),STAT=istatus)
uwndh=0
ALLOCATE(vwndh(nx,ny,nhtgem),STAT=istatus)
vwndh=0
ALLOCATE(wwndh(nx,ny,nhtgem),STAT=istatus)
wwndh=0
ALLOCATE(psf(nx,ny),STAT=istatus)
psf=0
ALLOCATE(mslp(nx,ny),STAT=istatus)
mslp=0
ALLOCATE(vort500(nx,ny),STAT=istatus)
vort500=0
ALLOCATE(temp2m(nx,ny),STAT=istatus)
temp2m=0
ALLOCATE(dewp2m(nx,ny),STAT=istatus)
dewp2m=0
ALLOCATE(qv2m(nx,ny),STAT=istatus)
qv2m=0
ALLOCATE(u10m(nx,ny),STAT=istatus)
u10m=0
ALLOCATE(v10m(nx,ny),STAT=istatus)
v10m=0
ALLOCATE(hgtsfc(nx,ny),STAT=istatus)
hgtsfc=0
ALLOCATE(refl1km(nx,ny),STAT=istatus)
refl1km=0
ALLOCATE(refl4km(nx,ny),STAT=istatus)
refl4km=0
ALLOCATE(cmpref(nx,ny),STAT=istatus)
cmpref=0
ALLOCATE(accppt(nx,ny),STAT=istatus)
accppt=0
ALLOCATE(convppt(nx,ny),STAT=istatus)
convppt=0
ALLOCATE(tpptold(nx,ny),STAT=istatus)
tpptold=0
ALLOCATE(cpptold(nx,ny),STAT=istatus)
cpptold=0
ALLOCATE(cape(nx,ny),STAT=istatus)
cape=0
ALLOCATE(mcape(nx,ny),STAT=istatus)
mcape=0
ALLOCATE(cin(nx,ny),STAT=istatus)
cin=0
ALLOCATE(mcin(nx,ny),STAT=istatus)
mcin=0
ALLOCATE(lcl(nx,ny),STAT=istatus)
lcl=0
ALLOCATE(srh01(nx,ny),STAT=istatus)
srh01=0
ALLOCATE(srh03(nx,ny),STAT=istatus)
srh03=0
ALLOCATE(uh25(nx,ny),STAT=istatus)
uh25=0
ALLOCATE(sh01(nx,ny),STAT=istatus)
sh01=0
ALLOCATE(sh06(nx,ny),STAT=istatus)
sh06=0
ALLOCATE(thck(nx,ny),STAT=istatus)
thck=0
ALLOCATE(li(nx,ny),STAT=istatus)
li=0
ALLOCATE(brn(nx,ny),STAT=istatus)
brn=0
ALLOCATE(pw(nx,ny),STAT=istatus)
pw=0
ALLOCATE(f2(nx,ny),STAT=istatus)
f2=0
ALLOCATE(xs(nx),STAT=istatus)
xs=0
ALLOCATE(ys(ny),STAT=istatus)
ys=0
ALLOCATE(zs(nx,ny,nz),STAT=istatus)
zs=0
ALLOCATE(zps(nx,ny,nz),STAT=istatus)
zps=0
ALLOCATE(zpc(nx,ny,nz),STAT=istatus)
zpc=0
ALLOCATE(zagl(nx,ny,nz),STAT=istatus)
zagl=0
ALLOCATE(u(nx,ny,nz),STAT=istatus)
u=0
ALLOCATE(v(nx,ny,nz),STAT=istatus)
v=0
ALLOCATE(w(nx,ny,nz),STAT=istatus)
w=0
ALLOCATE(pt(nx,ny,nz),STAT=istatus)
pt=0
ALLOCATE(qv(nx,ny,nz),STAT=istatus)
qv=0
ALLOCATE(raint(nx,ny),STAT=istatus)
raint=0
ALLOCATE(hterain(nx,ny),STAT=istatus)
hterain=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(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(thet(nx,ny),STAT=istatus)
thet=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(pt2m(nx,ny),STAT=istatus)
pt2m=0
ALLOCATE(t700(nx,ny),STAT=istatus)
t700=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
nxlg = (nx-fzone)*nproc_x + fzone
nylg = (ny-fzone)*nproc_y + fzone
year1=year
month1=month
dateomon=day
houroday=hour
minohour=ctime
IF (myproc == 0) THEN
PRINT *, year1,month1,dateomon,houroday,ctime
PRINT *,'mapproj,trulat1,trulat2,trulon,sclfct:'
PRINT *, mapproj,trulat1,trulat2,trulon,sclfct
! Write time parameters
IF(hisfile(lenfil-4:lenfil-4) == '_' ) THEN
print *,'Write time parameter file:',trim(outheader)//'.time' &
//hisfile(lenfil-10:lenfil-5)
OPEN(4,file=trim(outheader)//'.time'//hisfile(lenfil-10:lenfil-5), &
form='formatted')
ELSE
print *,'Write time parameter file:',trim(outheader)//'.time' &
//hisfile(lenfil-5:lenfil)
OPEN(4,file=trim(outheader)//'.time'//hisfile(lenfil-5:lenfil), &
form='formatted')
END IF
WRITE(4,*) year,month,day,hour,minute,ctime
CLOSE(4)
END IF
! Convert the rainfal arries and store the old ones.
DO j=1,ny-1
DO i=1,nx-1
accppt(i,j)=raing(i,j)+rainc(i,j)
END DO
END DO
convppt = rainc
!!!!! This block is moved outside postcore
! 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
!
! tppname='tmp_accppt'
! cppname='tmp_conppt'
! varname='RAIN'
! varunit='mm'
!
! IF (int(ctime) == 0) THEN ! to clear the 1st time level values
! tpptold = accppt
! cpptold = convppt
! ELSE
! inquire (file=tppname,exist=tppExist)
! inquire (file=cppname,exist=cppExist)
! IF(tppExist) THEN
! CALL binread2d(nx,ny,trim(tppname),tpptold)
! ELSE
! tpptold = accppt
! END IF
! IF(cppExist) THEN
! CALL binread2d(nx,ny,trim(cppname),cpptold)
! ELSE
! cpptold = convppt
! END IF
! END IF
! CALL bindump2d(nx,ny,trim(tppname),varname,varunit,accppt)
! CALL bindump2d(nx,ny,trim(cppname),varname,varunit,convppt)
! CALL pptsto(nx,ny,accppt,tpptold)
! CALL pptsto(nx,ny,convppt,cpptold)
! END IF
!
!-----------------------------------------------------------------------
!
! Set coordinates at the grid volume center
!
!-----------------------------------------------------------------------
!
! if(nf == 1) then ! ONLY need to be done the first time level
DO k=1,nz-1
DO j=1,ny
DO i=1,nx
zs(i,j,k) = (z(k)+z(k+1))*0.5 ! in meter
zps(i,j,k)= (zp(i,j,k)+zp(i,j,k+1))*0.5 ! in meter
zpc(i,j,k)=zps(i,j,k)/1000.0 ! in km
zagl(i,j,k)=(zps(i,j,k) - zp(i,j,2))/1000.0 ! AGL height (km)
END DO
END DO
END DO
DO i=1,nx-1
xs(i) = (x(i)+x(i+1))*0.5
END DO
! xs(nx)=2.*(xs(nx-1)-xs(nx-2))
DO j=1,ny-1
ys(j) = (y(j)+y(j+1))*0.5
END DO
! ys(ny)=2.*(ys(ny-1)-ys(ny-2))
!
!-----------------------------------------------------------------------
! Set the map parameters and put the origin the same as used in defining
! x and y etc. The origin is actually point (2,2)
dx = x(3)-x(2)
dy = y(3)-y(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 )
swx = ctrx- (nxlg-3)*dx*0.5
swy = ctry- (nylg-3)*dy*0.5
CALL setorig( 1, swx, swy)
!-----------------------------------------------------------------------
! Calculate the Coriolis para. f=2wsin(lat) at the data's scalar points
DO j=1,ny-1
DO i=1,nx-1
CALL xytoll( 1,1, xs(i),ys(j),lat1,lon1)
f2(i,j)=SIN(lat1*pi/180.0)*2*omega
END DO
END DO
!
!-----------------------------------------------------------------------
! 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, x(1),y(1),lat1,lon1)
CALL xytoll( 1,1, x(nx),y(ny),lat2,lon2)
! IF(myproc == 0) PRINT *,lat1,lon1,lat2,lon2,' Data Map corners'
! PRINT *,lat1,lon1,lat2,lon2,' Data Map corners',myproc
! PRINT *,qswcorn,qswcorw,qnecorn,qnecorw,'Output Map corners'
IF(mp_opt > 0) THEN
source = nprocs -1
CALL inctag
IF(myproc == source) THEN
itags = gentag
CALL mpsendr(lat2,1,destination,itags,ierr)
itags = gentag + 1
CALL mpsendr(lon2,1,destination,itags,ierr)
END IF
IF(myproc == 0) THEN
itagr = gentag
CALL mprecvr(lat2,1,source,itagr,ierr)
itagr = gentag + 1
CALL mprecvr(lon2,1,source,itagr,ierr)
END IF
END IF
IF(myproc == 0) THEN
WRITE(6,'(/a,4(a,F7.2),a/)') ' ARPS domain ', &
'SW: (',lat1,',',lon1, ') NE: (',lat2,',',lon2,')'
print *,'Write map parameter file:',trim(outheader)//'.map'
OPEN(4,file=trim(outheader)//'.map',form='formatted')
WRITE(4,*) mapproj,trulon,trulat1,trulat2,lat1,lon1,lat2,lon2
CLOSE(4)
END IF
! END IF
!
!-----------------------------------------------------------------------
!
! Initializing GEMPAK (if igempak = 1)
!
!-----------------------------------------------------------------------
!
IF (igempak ==1) THEN
CALL initgemio(nxlg,nylg,mapproj,trulat1,trulat2,trulon, &
lat1,lon1,lat2,lon2,iret)
IF( iret /= 0 ) GO TO 950
END IF
li=0.0
cape=0.0
brn=0.0
cin=0.0
pw=0.0
! combine prt and bar
DO k=1,nz-1
DO j=1,ny-1
DO i=1,nx-1
u(i,j,k)=0.5*(uprt(i,j,k)+ubar(i,j,k)+ &
uprt(i+1,j,k)+ubar(i+1,j,k))
v(i,j,k)=0.5*(vprt(i,j,k)+vbar(i,j,k)+ &
vprt(i,j+1,k)+vbar(i,j+1,k))
w(i,j,k)=0.5*(wprt(i,j,k)+wbar(i,j,k)+ &
wprt(i,j,k+1)+wbar(i,j,k+1))
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)
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(w,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(qv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL edgfill(pt,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)
! convert u,v,w to scalar points.
!
!----------------------------------------------------------------------
!
! 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)
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)
CALL edgfill(tem4,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
CALL a3dmax0(qv,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, amax,amin)
IF(myproc==0) print *,'QV -- amax, amin: ',amax, amin
CALL a3dmax0(tem2,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, amax,amin)
IF(myproc==0) print *,'T -- amax, amin: ',amax, amin
CALL a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, amax,amin)
IF(myproc==0) print *,'P -- amax, amin: ',amax, amin
CALL a3dmax0(tem3,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1, amax,amin)
IF(myproc==0) print *,'Td -- amax, amin: ',amax, amin
ice=1
CALL reflec(nx,ny,nz,rhobar,qr,qs,qh,tem5)
! CALL reflec_ferrier(nx,ny,nz,rhobar,qr,qs,qh,tem2,tem5)
CALL hintrp1(nx,ny,nz,2,nz-2,tem5,zagl,1.0,refl1km) ! interpolate to 1-km AGL
CALL hintrp1(nx,ny,nz,2,nz-2,tem5,zagl,4.0,refl4km) ! interpolate to 4-km AGL
DO j=1,ny-1
DO i=1,nx-1
u10m(i,j)=u(i,j,2)
v10m(i,j)=v(i,j,2)
hgtsfc(i,j)=zp(i,j,2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! 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), mslp
! CALL psfsl(nx,ny,tem2(1,1,2),qv(1,1,2),tem1(1,1,2), &
! zps(1,1,2),zp(1,1,2),psf,mslp)
rlevel=-ALOG(70000.0)
CALL hintrp1(nx,ny,nz,2,nz-2,tem2,algpzc,rlevel,t700)
DO j=1,ny
DO i=1,nx
p00 = 0.01*(tem1(i,j,2))
IF(p00 <= 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*zps(i,j,2)/tem1(i,j,1))**ex2
! mslp(i,j)=p00*exp(9.8*zps(i,j,2)/(286.5*(tem2(i,j,2)+0.5*gamma*zps(i,j,2))))
!write(0,*) 'before mslp'
!
mslp(i,j)=p00*((t00+gamma*zps(i,j,2))/t00)**ex2
!write(0,*) 'after mslp'
psf(i,j)=p00 *exp(9.8*(zps(i,j,2)-zp(i,j,2))/(286.5*tem2(i,j,2)))
END DO
END DO
!print *,zps(5,5,2)-zp(5,5,2),zps(5,5,2),tem2(5,5,2),tem2(5,5,2)+0.5*gamma*zps(5,5,2)
!print *,tem1(5,5,2),psf(5,5),mslp(5,5)
! 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)/(287.0*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 ! unit: g/m2
pw(i,j)=pw1d ! unit: mm or kg/m2
END DO
END DO
!CALL a3dmax0(pw,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
!print *,'PWAT -- amax, amin: ',amax, amin
!if(sfcphy > 0) then
! 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)
!CALL a3dmax0(qv2m,1,nx,1,nx-1,1,ny,1,ny,1,1,1,1, amax,amin)
!print *,'Qv -- amax, amin: ',amax, amin
!endif
!-----------------------------------------------------------------------
! Generate pressure level fields
IF( nprgem > 0 ) THEN
DO klev=1,nprgem
rlevel = -ALOG(float(iprgem(klev))*100.0)
IF(myproc == 0) &
PRINT *, ' Fields at pressure level ',iprgem(klev),' hPa'
CALL hintrp1(nx,ny,nz,2,nz-2,zps,algpzc,rlevel,hgtp(1,1,klev))
CALL hintrp1(nx,ny,nz,2,nz-2,u,algpzc,rlevel,uwndp(1,1,klev))
CALL hintrp1(nx,ny,nz,2,nz-2,v,algpzc,rlevel,vwndp(1,1,klev))
CALL hintrp1(nx,ny,nz,2,nz-2,w,algpzc,rlevel,wwndp(1,1,klev))
CALL hintrp1(nx,ny,nz,2,nz-2,tem2-273.15,algpzc,rlevel,tmpp(1,1,klev))
CALL hintrp1(nx,ny,nz,2,nz-2,1000.*qv,algpzc,rlevel,sphp(1,1,klev))
ENDDO
END IF
! absolute vorticity at 500 hPa
CALL vrtcalc(uwndp(1,1,4),vwndp(1,1,4),nx,ny,vort500,f2,dx,dy)
rlevel=-ALOG(100000.0)
CALL hintrp1(nx,ny,nz,2,nz-2,zps,algpzc,rlevel,thck)
DO j=1,ny-1
DO i=1,nx-1
thck(i,j) = hgtp(1,1,4)-thck(i,j)
! temp850(i,j)=temp850(i,j)-273.15
END DO
END DO
!-----------------------------------------------------------------------
! Generate constant height level fields
IF( nhtgem > 0 ) THEN
DO klev=1,nhtgem
rlevel = ihtgem(klev)
IF(myproc == 0) &
PRINT *, ' Fields at AGL height ',ihtgem(klev),' km'
CALL hintrp1(nx,ny,nz,2,nz-2,u,zagl,rlevel,uwndh(1,1,klev))
CALL hintrp1(nx,ny,nz,2,nz-2,v,zagl,rlevel,vwndh(1,1,klev))
CALL hintrp1(nx,ny,nz,2,nz-2,w,zagl,rlevel,wwndh(1,1,klev))
ENDDO
END IF
! Wind shear: 0-1 km AGL; 0-6km AGL
DO j=1,ny-1
DO i=1,nx-1
sh01(i,j) = (SQRT( (uwndh(i,j,1)-u10m(i,j))**2 + &
( vwndh(i,j,1)-v10m(i,j))**2 ))/1000.
sh06(i,j) = (SQRT( (uwndh(i,j,6)-u10m(i,j))**2 + &
(vwndh(i,j,6)-v10m(i,j))**2 ))/6000.
END DO
END DO
! temp2m, dewp2m scheme 1, direct interpolation from grid values
z01=2.0
! CALL secthrz(nx,ny,nz,tem3,zs,dewp2m)
z01=2.0
! CALL secthrz(nx,ny,nz,tem2,zs,temp2m)
! 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
dewp2m(i,j)=(tem3(i,j,2) - 273.15)*9.0/5.0 + 32.0 ! F unit
temp2m(i,j)=(tem2(i,j,2) - 273.15)*9.0/5.0 + 32.0 ! F unit
qv2m(i,j)=1000.*qv(i,j,2)
IF (dewp2m(i,j) > temp2m(i,j)) dewp2m(i,j)=temp2m(i,j)
END DO
END DO
!CALL a3dmax0(qv2m,1,nx,1,nx-1,1,ny,1,ny,1,1,1,1, amax,amin)
!print *,'Qv -- amax, amin: ',amax, amin
!-----------------------------------------------------------------------
!
! 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 == 1) THEN
IF (myproc == 0) PRINT *, ' calling arps_be'
CALL arps_be1(nx,ny,nz, &
tem1,zpc,tem2,tem4, &
lcl,lfc,el,twdf,li,cape,mcape,cin,mcin,capst, &
wrk1,wrk2,wrk3,wrk4,wrk5,wrk6, &
wrk7,wrk8,wrk9,wrk10,wrk11,wrk12,t1)
IF (myproc == 0) PRINT *, ' back from arps_be'
ELSE
IF (myproc == 0) PRINT *,'icape option is incorrect, chose 1 only'
STOP
END IF
! Temperary assigning
! mcape = cape
! mcin = cin
!
!-----------------------------------------------------------------------
!
! Calculate helicity and other shear-related paramaters
!
!-----------------------------------------------------------------------
!
IF (myproc == 0) PRINT *, ' calling calcshr'
CALL xytomf(nx,ny,x,y,b1)
CALL calcshr(nx,ny,nz,x,y,zp,b1, &
tem1,tem2,u,v,cape, &
shr37,ustrm,vstrm,srlfl,srmfl,srh03,brn,brnu,blcon, &
t1,t2,t3,t4,tem4)
IF (myproc == 0) PRINT *, ' back from calcshr'
! Calculate 0-1km AGL storm-relative helicity
DO j=1,ny-1
DO i=1,nx-1
sumh=0.
DO k=3,nz-2
IF(zps(i,j,k) > zp(i,j,2)+1000.) EXIT
sumh=sumh + &
( u(i,j,k)*v(i,j,k-1) ) - &
( v(i,j,k)*u(i,j,k-1) )
ENDDO
sumh=sumh + &
( uwndh(i,j,1)*v(i,j,k-1) ) - &
( vwndh(i,j,1)*u(i,j,k-1) )
srh01(i,j)=sumh + (vwndh(i,j,1)-v(i,j,2))*ustrm(i,j) - &
(uwndh(i,j,1)-u(i,j,2))*vstrm(i,j)
ENDDO
ENDDO
! Updraft Helicity (2 - 5 km layer)
tem5 = 0.0
uh25 = 0.0
CALL vrtcalc(uwndh(1,1,2),vwndh(1,1,2),nx,ny,tem5(1,1,2),tem5(1,1,1),dx,dy)
CALL vrtcalc(uwndh(1,1,3),vwndh(1,1,3),nx,ny,tem5(1,1,3),tem5(1,1,1),dx,dy)
CALL vrtcalc(uwndh(1,1,4),vwndh(1,1,4),nx,ny,tem5(1,1,4),tem5(1,1,1),dx,dy)
CALL vrtcalc(uwndh(1,1,5),vwndh(1,1,5),nx,ny,tem5(1,1,5),tem5(1,1,1),dx,dy)
! tem5(1,1,2) - tem5(1,1,5) return relative vorticity for the levels
DO j=1,ny-1
DO i=1,nx-1
ave1 = (wwndh(i,j,2)*tem5(i,j,2)+wwndh(i,j,3)*tem5(i,j,3))*0.5*1000.0
ave2 = (wwndh(i,j,3)*tem5(i,j,3)+wwndh(i,j,4)*tem5(i,j,4))*0.5*1000.0
ave3 = (wwndh(i,j,4)*tem5(i,j,4)+wwndh(i,j,5)*tem5(i,j,5))*0.5*1000.0
! require upward motion throughout the depth
IF(wwndh(i,j,2)>0.0 .AND. wwndh(i,j,3)>0.0 .AND. wwndh(i,j,4)>0.0 &
.AND. wwndh(i,j,5)>0.0) THEN
uh25(i,j) = ave1+ave2+ave3 ! updarf helicity
END IF
ENDDO
ENDDO
!CALL a3dmax0(wwndh,1,nx,1,nx-1,1,ny,1,ny-1,1,6,1,6, amax,amin)
!IF(myproc==0) print *,'wwndh -- amax, amin: ',amax, amin
CALL a3dmax0(uh25,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin)
IF(myproc==0) print *,'uh25 -- amax, amin: ',amax, amin
!CALL a3dmax(uh25,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1, amax,amin, &
! imax,jmax,kmax, imin,jmin,kmin)
!IF(myproc==0) print *,'uh25 -- amax, amin: ',amax, amin,imax,jmax
!-----------------------------------------------------------------------
! Interpolate the 2_d variables
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! 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
! Find forecast hour for name construction
ihr = int(ctime/3600.)
WRITE(ihrc,'(a1,i3.3)') 'f',ihr
if(iascii /= 0) then ! will add in the future
! filename=trim(outheader)//'.asc'//hisfile(nf)(lenfil-5:lenfil)
! lenfil = len_trim(filename)
! IF(myproc == 0) WRITE(6,'(a,a/)') &
! ' ASCII output file: ', filename(1:lenfil)
!
! CALL enswrtasc (filename(1:lenfil), 2, dxnew, &
! houroday,minohour,dayoweek,dateomon,month1,year1, &
! 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, pw,cmpref, &
! refl,u10m,v10m)
endif
if(ibinary /= 0) then
IF(hisfile(lenfil-4:lenfil-4) == '_' ) THEN
filename=outheader(1:lenbin)//'.bin'//hisfile(lenfil-10:lenfil-5)
ELSE
filename=outheader(1:lenbin)//'.bin'//hisfile(lenfil-5:lenfil)
END IF
lenfil = len_trim(filename)
IF(myproc == 0) WRITE(6,'(a,a/)') &
' Binary output file: ', filename(1:lenfil)
CALL enswrtbin (nx,ny,filename(1:lenfil), lenfil, &
ctime,year,month,day,hour,minute, &
iprgem,nprgem,ihtgem,nhtgem, &
hgtp,uwndp,vwndp,wwndp,tmpp,sphp, &
uwndh,vwndh,wwndh, &
psf,mslp,vort500,temp2m,dewp2m,qv2m, &
u10m,v10m,hgtsfc,refl1km,refl4km,cmpref, &
accppt,convppt,cape,mcape,cin,mcin,lcl, &
srh01,srh03,uh25,sh01,sh06,thck,li,brn,pw)
endif
if(igempak /= 0) then
! filename=trim(gemoutheader)//'_'//hisfile(nf)(lenfil-5:lenfil)//'.gem'
! filename=trim(gemoutheader)//ihrc//'.gem'
filename=gemoutheader(1:lengem)//ihrc
lenfil = len_trim(filename)
IF(myproc == 0) WRITE(6,'(a,a/)') &
' GEMPAK output file: ', filename(1:lenfil)
CALL enswrtgem (nx,ny,filename(1:lenfil), lenfil, &
ctime,year,month,day,hour,minute, &
iprgem,nprgem,ihtgem,nhtgem, &
hgtp,uwndp,vwndp,wwndp,tmpp,sphp, &
uwndh,vwndh,wwndh, &
psf,mslp,vort500,temp2m,dewp2m,qv2m, &
u10m,v10m,hgtsfc,refl1km,refl4km,cmpref, &
accppt,convppt,cape,mcape,cin,mcin,lcl, &
srh01,srh03,uh25,sh01,sh06,thck,li,brn,pw)
! Generate a ready file
IF(myproc == 0) THEN
! CALL getunit( nchout )
OPEN (UNIT=4,FILE=trim(filename)//"_ready")
WRITE (4,'(a)') trim(filename)//"_ready"
CLOSE (4)
! CALL retunit (nchout )
END IF
endif
!-----------------------------------------------------------------------
DEALLOCATE(hgtp)
DEALLOCATE(uwndp)
DEALLOCATE(vwndp)
DEALLOCATE(wwndp)
DEALLOCATE(tmpp)
DEALLOCATE(sphp)
DEALLOCATE(uwndh)
DEALLOCATE(vwndh)
DEALLOCATE(wwndh)
DEALLOCATE(psf)
DEALLOCATE(mslp)
DEALLOCATE(vort500)
DEALLOCATE(temp2m)
DEALLOCATE(dewp2m)
DEALLOCATE(qv2m)
DEALLOCATE(u10m)
DEALLOCATE(v10m)
DEALLOCATE(hgtsfc)
DEALLOCATE(refl1km)
DEALLOCATE(refl4km)
DEALLOCATE(cmpref)
DEALLOCATE(accppt)
DEALLOCATE(convppt)
DEALLOCATE(tpptold)
DEALLOCATE(cpptold)
DEALLOCATE(cape)
DEALLOCATE(mcape)
DEALLOCATE(cin)
DEALLOCATE(mcin)
DEALLOCATE(lcl)
DEALLOCATE(srh01)
DEALLOCATE(srh03)
DEALLOCATE(uh25)
DEALLOCATE(sh01)
DEALLOCATE(sh06)
DEALLOCATE(thck)
DEALLOCATE(li)
DEALLOCATE(brn)
DEALLOCATE(pw)
DEALLOCATE(f2)
DEALLOCATE(xs)
DEALLOCATE(ys)
DEALLOCATE(zs)
DEALLOCATE(zps)
DEALLOCATE(zpc)
DEALLOCATE(zagl)
DEALLOCATE(u)
DEALLOCATE(v)
DEALLOCATE(w)
DEALLOCATE(pt)
DEALLOCATE(qv)
DEALLOCATE(raint)
DEALLOCATE(hterain)
DEALLOCATE(algpzc)
DEALLOCATE(wrk1)
DEALLOCATE(wrk2)
DEALLOCATE(wrk3)
DEALLOCATE(wrk4)
DEALLOCATE(wrk5)
DEALLOCATE(wrk6)
DEALLOCATE(wrk7)
DEALLOCATE(wrk8)
DEALLOCATE(wrk9)
DEALLOCATE(wrk10)
DEALLOCATE(wrk11)
DEALLOCATE(wrk12)
DEALLOCATE(lfc)
DEALLOCATE(el)
DEALLOCATE(twdf)
DEALLOCATE(mbe)
DEALLOCATE(thet)
DEALLOCATE(brnu)
DEALLOCATE(srlfl)
DEALLOCATE(srmfl)
DEALLOCATE(shr37)
DEALLOCATE(ustrm)
DEALLOCATE(vstrm)
DEALLOCATE(capst)
DEALLOCATE(blcon)
DEALLOCATE(pt2m)
DEALLOCATE(t700)
DEALLOCATE(b1)
DEALLOCATE(t1)
DEALLOCATE(t2)
DEALLOCATE(t3)
DEALLOCATE(t4)
RETURN
950 CONTINUE
IF(myproc == 0) WRITE(6,'(a)') ' Error setting up GEMPAK file'
RETURN
END SUBROUTINE postcore
SUBROUTINE pptsto(nx,ny,a,b) 3
INTEGER :: nx,ny,i,j
REAL :: a(nx,ny),b(nx,ny),c
DO j=1,ny
DO i=1,nx
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 vrtcalc(u,v,nx,ny,b,f,dx,dy) 8,2
IMPLICIT NONE
INTEGER :: nx,ny
REAL :: u(nx,ny),v(nx,ny),b(nx,ny),f(nx,ny)
REAL :: dx,dy
INTEGER :: i,j
DO j=1,ny
DO i=1,nx
b(i,j)=0.0
END DO
END DO
DO j=2,ny-1
DO i=2,nx-1
b(i,j)=0.5*((v(i+1,j)-v(i-1,j))/dx &
-(u(i,j+1)-u(i,j-1))/dy) &
+f(i,j)
END DO
END DO
CALL edgfill(b,1,nx,2,nx-1,1,ny,2,ny-1,1,1,1,1)
RETURN
END SUBROUTINE vrtcalc
SUBROUTINE ptqv2m(nx,ny,nz,nstyps,zp, & 1,10
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 cptc(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 cptcwtr(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 cptc(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 cptcwtr(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 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 calender(year,mon,day,hour,ds,wkd) 2
! 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