!
!##################################################################
!##################################################################
!###### ######
!###### 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