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