! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE initgemio ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE initgemio(nxlg,nylg,mapproj,trulat1,trulat2,trulon, & 1 lat1,lon1,lat2,lon2,iret) ! IMPLICIT NONE !----------------------------------------------------------------------- ! ! PURPOSE: Initial IO interface for GEMPAK package ! !----------------------------------------------------------------------- ! ! Author: Yunheng Wang (04/02/2007) ! !----------------------------------------------------------------------- INTEGER, INTENT(IN) :: nxlg, nylg INTEGER, INTENT(IN) :: mapproj REAL, INTENT(IN) :: trulat1, trulat2, trulon REAL, INTENT(IN) :: lat1, lon1, lat2, lon2 INTEGER, INTENT(OUT) :: iret !----------------------------------------------------------------------- ! ! GEMPAK variables ! !----------------------------------------------------------------------- ! CHARACTER(LEN=3), PARAMETER :: cproj(0:3) = (/'CED','STR','LCC','MER'/) ! CHARACTER (LEN=72) :: gdarea CHARACTER (LEN=72) :: chproj ! REAL, PARAMETER :: deltan = 1.0 ! REAL, PARAMETER :: deltax = -9999., deltay = -9999. LOGICAL, PARAMETER :: angflg = .TRUE. REAL :: angle1,angle2,angle3 REAL :: gbnds(4) LOGICAL :: notopn ! should be deleted !----------------------------------------------------------------------- ! Include files: !----------------------------------------------------------------------- INCLUDE 'GEMPRM.PRM' INCLUDE 'mp.inc' REAL :: rnvblk (llnnav), anlblk (llnanl) COMMON /gemblk/rnvblk,anlblk,notop !fpp$ expand (f_qvsat) !dir$ inline always f_qvsat ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@! ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ iret = 0 CALL in_bdta(iret) ! initializing GEMPAK sans TAE ! !----------------------------------------------------------------------- ! ! Initializing GEMPAK ! !----------------------------------------------------------------------- ! chproj = cproj(mapproj) angle2 = trulon ! ! modify projection if polar stereographic ! IF (chproj == 'STR') THEN angle1=90.0 angle3=0.0 IF(myproc == 0) THEN WRITE (*,*) '**********************************' WRITE (*,*) '**********************************' WRITE (*,*) ' ' WRITE (*,*) 'SETTING ANGLE1=90. FOR GEMPAK FILE' WRITE (*,*) 'SETTING ANGLE3=0.0 FOR GEMPAK FILE' WRITE (*,*) '(GEMPAK origin is at North Pole)' WRITE (*,*) ' ' WRITE (*,*) '**********************************' WRITE (*,*) '**********************************' END IF ELSE angle1=trulat1 angle3=trulat2 IF(myproc == 0) THEN WRITE (*,*) '**********************************' WRITE (*,*) '**********************************' WRITE (*,*) ' ' WRITE (*,*) ' SETTING ANGLE1=TRULAT1 ' WRITE (*,*) ' SETTING ANGLE3=TRULAT2 ' WRITE (*,*) ' ' WRITE (*,*) '**********************************' WRITE (*,*) '**********************************' END IF END IF gbnds(1)=lat1 gbnds(2)=lon1 gbnds(3)=lat2 gbnds(4)=lon2 ! WRITE(gdarea,'(f8.4,a1,f9.4,a1,f8.4,a1,f9.4)') & ! lat1,';',lon1,';',lat2,';',lon2 CALL gr_mnav(chproj, nxlg,nylg,lat1,lon1,lat2,lon2, & angle1,angle2,angle3,angflg, & rnvblk, iret) RETURN END SUBROUTINE initgemio ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE enswrtgem ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE enswrtgem (nx,ny, filename, lens, & 2,30 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) INTEGER :: nx,ny CHARACTER (LEN=*) :: filename CHARACTER (LEN=72) :: gdfile INTEGER :: lens, nchout INTEGER :: imxgrd PARAMETER (imxgrd=500) INTEGER :: ivsfc,ivprs,ivtheta,ivhgt PARAMETER (ivsfc=0,ivprs=1,ivtheta=2,ivhgt=3) INTEGER :: level(2) INTEGER :: ighdr(2) INTEGER :: ipktyp, nbits PARAMETER (ipktyp=1,nbits=16) CHARACTER (LEN=20) :: gemtime(2) CHARACTER (LEN=12) :: parm INTEGER :: nprgem,nhtgem INTEGER :: iprgem(nprgem),ihtgem(nhtgem) REAL :: hgtp(nx,ny,nprgem),uwndp(nx,ny,nprgem),vwndp(nx,ny,nprgem) REAL :: wwndp(nx,ny,nprgem),tmpp(nx,ny,nprgem),sphp(nx,ny,nprgem) REAL :: uwndh(nx,ny,nhtgem),vwndh(nx,ny,nhtgem),wwndh(nx,ny,nhtgem) REAL :: psf(nx,ny),mslp(nx,ny),vort500(nx,ny) REAL :: temp2m(nx,ny),dewp2m(nx,ny),qv2m(nx,ny) REAL :: u10m(nx,ny),v10m(nx,ny),hgtsfc(nx,ny) REAL :: refl1km(nx,ny),refl4km(nx,ny),cmpref(nx,ny) REAL :: accppt(nx,ny),convppt(nx,ny) REAL :: cape(nx,ny),mcape(nx,ny),cin(nx,ny),mcin(nx,ny),lcl(nx,ny) REAL :: srh01(nx,ny),srh03(nx,ny),uh25(nx,ny),sh01(nx,ny),sh06(nx,ny) REAL :: thck(nx,ny),li(nx,ny),brn(nx,ny),pw(nx,ny) INTEGER :: i,j,k,klev,iret,igdfil INTEGER :: year,month,day,hour,minute INTEGER :: iyr,ifhr,ifmin,ifile,ihd REAL :: ctime LOGICAL :: wrtflg,rewrite,notopn,gemExist INCLUDE 'mp.inc' INCLUDE 'GEMPRM.PRM' REAL :: rnvblk (llnnav), anlblk (llnanl) COMMON /gemblk/ rnvblk,anlblk,notopn !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Begining of executable code ... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! !----------------------------------------------------------------------- ! ! Build the GEMPAK grid time string ! It has format yymodd/hhmnFHHH ! yy: year mo: month dd: GMT day ! hh: GMT hour mn: minute ! F: seperation charcter ! HHH: forecast hour (000 = analysis) ! example time(1)='950126/1200F000' ! !----------------------------------------------------------------------- ! curtim = ctime gemtime(1)=' ' gemtime(2)=' ' iyr=MOD(year,100) ifhr=INT(curtim/3600.) ifmin=nint((curtim-(ifhr*3600.))/60.) IF(curtim == 0) THEN WRITE(gemtime(1), & '(i2.2,i2.2,i2.2,a1,i2.2,i2.2,a4)') & iyr,month,day,'/',hour,minute,'F000' ELSE WRITE(gemtime(1), & '(i2.2,i2.2,i2.2,a1,i2.2,i2.2,a1,i3.3,i2.2)') & iyr,month,day,'/',hour,minute,'F',ifhr,ifmin END IF IF(myproc==0) WRITE(6,'(a,a)') ' GEMPAK time string ',gemtime(1) ! !----------------------------------------------------------------------- ! ! Initialize header, grid area coordinates and analysis block ! !----------------------------------------------------------------------- ! ! IF(notopn) THEN ihd = 2 ighdr(1)=0 ighdr(2)=0 ! DO i=1,llnanl anlblk(i) = 0. END DO ! END IF ! !----------------------------------------------------------------------- ! ! Open/Create the grid file ! !----------------------------------------------------------------------- ! gdfile = filename(1:lens) ! WRITE (6, *) 'Opening GEMPAK file ',gdfile(1:lens) ! CALL gd_opnf ( gdfile, wrtflg, igdfil, navsz, rnvblk, & DO i=0,nprocs-1,max_fopen IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN inquire (file=gdfile,exist=gemExist) IF ( gemExist ) THEN IF(myproc==0) THEN WRITE (6, *) 'GEMPAK file ',trim(gdfile),' already exists!' stop " STOP!!!" END IF ELSE IF(myproc==0) THEN WRITE (6, *) 'Creating file ',trim(gdfile) CALL gd_cref ( gdfile, llnnav, rnvblk, llnanl, & anlblk, ihd, imxgrd, igdfil, iret ) END IF IF( iret /= 0 ) GO TO 950 END IF ! !----------------------------------------------------------------------- ! ! Set write flag to true ! !----------------------------------------------------------------------- ! wrtflg=.true. ! CALL gd_swrt( igdfil, wrtflg, iret ) IF(myproc==0) CALL gd_swrt( igdfil, wrtflg, iret ) IF( iret /= 0 ) GO TO 950 notopn=.false. !----------------------------------------------------------------------- ! ! Output one level fields ! !----------------------------------------------------------------------- ! level(1)=0 level(2)=-1 IF(myproc==0) PRINT *, ' Writing surface pressure' parm='PRES' CALL gemdump2d( igdfil, psf, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing MSLP pressure' parm='PMSL' CALL gemdump2d( igdfil, mslp, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) ! ! IF(myproc==0) PRINT *,' Writing 1-h total accumulated rainfall' parm='P01M' CALL gemdump2d( igdfil, accppt, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) ! IF(myproc==0) PRINT *,' Writing 1-h convective rainfall' ! parm='RAIN_C' ! CALL gemdump2d( igdfil, convppt, & ! nx, ny, ighdr, & ! gemtime, level, ivsfc, parm, & ! rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *,' Writing precipitable water' parm='PWAT' CALL gemdump2d( igdfil, pw, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 2-m temperature (C)' parm='TMPF' CALL gemdump2d( igdfil, temp2m, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 2-m dew point temperature (C)' parm='DWPF' CALL gemdump2d( igdfil, dewp2m, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 2-m specific humidity' parm='MIXR' CALL gemdump2d( igdfil, qv2m, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) ! IF(myproc==0) PRINT *, ' Writing 10-m u velocity' parm='UREL' CALL gemdump2d( igdfil, u10m, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) ! IF(myproc==0) PRINT *, ' Writing 10-m v velocity' parm='VREL' CALL gemdump2d( igdfil, v10m, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) ! IF(myproc==0) PRINT *, ' Writing 500 hPa vorticity' ! parm='VORT' ! CALL gemdump2d( igdfil, vort500, & ! nx, ny, ighdr, & ! gemtime, level, ivsfc, parm, & ! rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing surface geopotential height' parm='HGHT' CALL gemdump2d( igdfil, hgtsfc, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 1-km AGL reflectivity' parm='REFL1KM' CALL gemdump2d( igdfil, refl1km, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 4-km AGL reflectivity' parm='REFL4KM' CALL gemdump2d( igdfil, refl4km, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing composite reflectivity' parm='REFLCMP' CALL gemdump2d( igdfil, cmpref, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing surface-based CAPE' parm='CAPE' CALL gemdump2d( igdfil, cape, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing moist unstable CAPE' parm='MUCAPE' CALL gemdump2d( igdfil, mcape, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing surface-based CIN' parm='CINS' CALL gemdump2d( igdfil, cin, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing moist unstable CIN' parm='MUCIN' CALL gemdump2d( igdfil, mcin, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing surface-based LCL' parm='HLCL' CALL gemdump2d( igdfil, lcl, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 0-1 km AGL storm-relative helicity' parm='SRH01' CALL gemdump2d( igdfil, srh01, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 0-3 km AGL storm-relative helicity' parm='SRH03' CALL gemdump2d( igdfil, srh03, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing updraft helicity' parm='VHEL' CALL gemdump2d( igdfil, uh25, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 0-1 km AGL wind shear' parm='SHR01' CALL gemdump2d( igdfil, sh01, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) IF(myproc==0) PRINT *, ' Writing 0-6 km AGL wind shear' parm='SHR06' CALL gemdump2d( igdfil, sh06, & nx, ny, ighdr, & gemtime, level, ivsfc, parm, & rewrite, ipktyp, nbits, iret) ! !----------------------------------------------------------------------- ! ! Output constant pressure level data ! !----------------------------------------------------------------------- ! IF( nprgem > 0 ) THEN rewrite=.false. DO klev=1,nprgem level(1)=iprgem(klev) level(2)=-1 IF(myproc==0) PRINT *, ' Writing GEMPAK data at pr= ',iprgem(klev),' hPa' parm='HGHT' CALL gemdump2d( igdfil, hgtp(1,1,klev), & nx, ny, ighdr, & gemtime, level, ivprs, parm, & rewrite, ipktyp, nbits, iret) parm='UREL' CALL gemdump2d( igdfil, uwndp(1,1,klev), & nx, ny, ighdr, & gemtime, level, ivprs, parm, & rewrite, ipktyp, nbits, iret) parm='VREL' CALL gemdump2d( igdfil, vwndp(1,1,klev), & nx, ny, ighdr, & gemtime, level, ivprs, parm, & rewrite, ipktyp, nbits, iret) parm='VVEL' CALL gemdump2d( igdfil, wwndp(1,1,klev), & nx, ny, ighdr, & gemtime, level, ivprs, parm, & rewrite, ipktyp, nbits, iret) parm='TMPC' CALL gemdump2d( igdfil, tmpp(1,1,klev), & nx, ny, ighdr, & gemtime, level, ivprs, parm, & rewrite, ipktyp, nbits, iret) parm='MIXR' CALL gemdump2d( igdfil, sphp(1,1,klev), & nx, ny, ighdr, & gemtime, level, ivprs, parm, & rewrite, ipktyp, nbits, iret) END DO ! nprgem-loop END IF ! !----------------------------------------------------------------------- ! ! Output AGL height level data ! !----------------------------------------------------------------------- ! ! IF( nhtgem > 0 ) THEN ! rewrite=.false. ! ! DO klev=1,nhtgem ! level(1)=ihtgem(klev) ! level(2)=-1 ! IF(myproc==0) PRINT *, ' Writing GEMPAK data at z= ',ihtgem(klev),' km AGL' ! parm='UREL' ! CALL gemdump2d( igdfil, uwndh(1,1,klev), & ! nx, ny, ighdr, & ! gemtime, level, ivhgt, parm, & ! rewrite, ipktyp, nbits, iret) ! parm='VREL' ! CALL gemdump2d( igdfil, vwndh(1,1,klev), & ! nx, ny, ighdr, & ! gemtime, level, ivhgt, parm, & ! rewrite, ipktyp, nbits, iret) ! parm='VVEL' ! CALL gemdump2d( igdfil, wwndh(1,1,klev), & ! nx, ny, ighdr, & ! gemtime, level, ivhgt, parm, & ! rewrite, ipktyp, nbits, iret) ! END DO ! nhtgem-loop ! END IF ! IF(myproc==0) THEN ! print *,'igdfil,iret:',igdfil,iret IF (myproc == 0) CALL gd_clos(igdfil,iret) ! print *,'iret:',ire ! Generate a ready file ! CALL getunit( nchout ) ! OPEN (UNIT=nchout,FILE=trim(gdfile)//"_ready") ! WRITE (nchout,'(a)') trim(gdfile)//"_ready" ! CLOSE (nchout) ! CALL retunit (nchout ) ! END IF END IF IF (mp_opt > 0) CALL mpbarrier END DO RETURN 950 CONTINUE IF(myproc == 0) WRITE(6,'(a)') ' Error setting up GEMPAK file' STOP RETURN END SUBROUTINE enswrtgem ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE GEMDUMP2D ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE gemdump2d(igdfil,qfield, & 29,3 nx, ny, ighdr, & gemtime, level, ivflag, parm, & rewrite, ipktyp, nbits, iret) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Write out one single 2D field in gempak format. This interface is ! specially written for MPI dump to a joined GEMPAK file ! !----------------------------------------------------------------------- ! ! AUTHOR: Fanyou Kong ! 02/20/2007 !----------------------------------------------------------------------- ! IMPLICIT NONE INCLUDE 'mp.inc' INTEGER :: igdfil, iret, istat INTEGER :: nx,ny INTEGER :: ivflag INTEGER :: level(2) INTEGER :: ighdr(2) CHARACTER (LEN=20) :: gemtime(2) CHARACTER (LEN=12) :: parm REAL :: qfield(nx,ny) LOGICAL :: rewrite integer nbits, ipktyp INTEGER :: nxlg, nylg REAL, ALLOCATABLE :: out2d(:,:) !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Begining of executable code ... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ nxlg = nproc_x*(nx-3)+3 nylg = nproc_y*(ny-3)+3 ALLOCATE (out2d( nxlg,nylg ),stat=istat) CALL edgfill(qfield,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1) ! IF (mp_opt > 0 .AND. joindmp > 0) THEN IF (mp_opt > 0 ) THEN ! Write joined single file for the 2D field CALL mpimerge2d(qfield,nx,ny,out2d) ! CALL edgfill(qfield,1,nxlg,1,nxlg-1,1,nylg,1,nylg-1,1,1,1,1) IF (myproc == 0) THEN CALL gd_wpgd( igdfil, out2d, & nxlg,nylg, ighdr, & gemtime, level, ivflag, parm, & rewrite, ipktyp, nbits, iret) END IF CALL mpupdatei(iret,1) ELSE ! Write non-MPI file for the 3D field CALL gd_wpgd( igdfil, qfield, & nx,ny, ighdr, & gemtime, level, ivflag, parm, & rewrite, ipktyp, nbits, iret) END IF DEALLOCATE(out2d) RETURN END SUBROUTINE gemdump2d