PROGRAM bin2gem,32 ! !################################################################## !################################################################## !###### ###### !###### PROGRAM BIN2GEM ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! ! !----------------------------------------------------------------------- ! PURPOSE: ! Read in binary dump produced by arpspost and write out in GEMPAK ! ! Author : Fanyou Kong ! History: March 31, 2007: First code implementation ! !----------------------------------------------------------------------- ! ! MODIFIED: ! !----------------------------------------------------------------------- ! IMPLICIT NONE INCLUDE 'GEMPRM.PRM' INCLUDE 'mp.inc' INTEGER :: nx,ny,nz INTEGER :: ireturn, iret, ierr, istatus CHARACTER (LEN=256) :: filename,filnam CHARACTER (LEN=256) :: binheader,gemoutheader INTEGER :: lengbf,lenfil,lenbin,lengem INTEGER :: ihr CHARACTER (LEN=4) :: ihrc CHARACTER (LEN=6) :: tlevel(100) CHARACTER (LEN=6) :: varname,varunit INTEGER :: i,j,k, itags,itagr INTEGER :: nf NAMELIST /bin2gem_data/nf,tlevel,binheader,gemoutheader REAL :: lat1,lon1,lat2,lon2 REAL :: ctime INTEGER :: year,month,day,hour,minute 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 :: cape(:,:),mcape(:,:),cin(:,:),mcin(:,:),lcl(:,:) REAL, ALLOCATABLE :: srh01(:,:),srh03(:,:),uh25(:,:),sh01(:,:),sh06(:,:) REAL, ALLOCATABLE :: thck(:,:),li(:,:),brn(:,:),pw(:,:) !----------------------------------------------------------------------- ! 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/ CHARACTER (LEN=6) :: hgtptag(5),uwndptag(5) CHARACTER (LEN=6) :: vwndptag(5),wwndptag(5) CHARACTER (LEN=6) :: tmpptag(5),qvptag(5) CHARACTER (LEN=6) :: uwndhtag(6),vwndhtag(6),wwndhtag(6) DATA hgtptag / 'hgt850','hgt700','hgt600','hgt500','hgt250'/ DATA uwndptag / 'u850__','u700__','u600__','u500__','u250__'/ DATA vwndptag / 'v850__','v700__','v600__','v500__','v250__'/ DATA wwndptag / 'w850__','w700__','w600__','w500__','w250__'/ DATA tmpptag / 'tmp850','tmp700','tmp600','tmp500','tmp250'/ DATA qvptag / 'sph850','sph700','sph600','sph500','sph250'/ ! DATA uwndhtag /'u1km__','u2km__','u3km__','u4km__','u5km__','u6km__'/ DATA vwndhtag /'v1km__','v2km__','v3km__','v4km__','v5km__','v6km__'/ DATA wwndhtag /'w1km__','w2km__','w3km__','w4km__','w5km__','w6km__'/ !----------------------------------------------------------------------- ! ! GEMPAK variables ! !----------------------------------------------------------------------- ! CHARACTER (LEN=72) :: gdfile,gdarea REAL :: rnvblk (llnnav), anlblk (llnanl) CHARACTER (LEN=3) :: cproj(0:3) DATA cproj /'CED','STR','LCC','MER'/ CHARACTER (LEN=72) :: chproj REAL :: deltan REAL :: deltax,deltay REAL :: angle1,angle2,angle3 LOGICAL :: angflg PARAMETER ( deltan= 1., & deltax= -9999., & deltay= -9999., & angflg=.true.) REAL :: gbnds(4) LOGICAL :: wrtflg,rewrite,notopn,gemExist,tppExist,cppExist COMMON /gemblk/rnvblk,anlblk,notopn nproc_x = 1 nproc_y = 1 readsplit = 0 joindmp = 1 nprocs = 1 max_fopen = 1 READ(5,bin2gem_data,END=200) write(6,bin2gem_data) GO TO 20 200 WRITE(6,'(a)') & 'Error reading NAMELIST block bin2gem_data. ', & 'Default values used.' 20 CONTINUE !----------------------------------------------------------------------- ! ! Get dimension nx,ny ! !----------------------------------------------------------------------- OPEN(4,FILE=trim(binheader)//'.sfpres'//tlevel(1),STATUS='old', & FORM='unformatted',ERR=9000, IOSTAT=istat) READ(4, ERR=9000, END=9000, IOSTAT=istat) nx,ny,nz CLOSE(4) print *,'nx,ny,nz:',nx,ny,nz ! !----------------------------------------------------------------------- ! ! Allocate arrays ! !----------------------------------------------------------------------- ! 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(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 ! Read in map parameters OPEN(4,file=trim(binheader)//'.map') READ(4,*) mapproj,trulon,trulat1,trulat2,lat1,lon1,lat2,lon2 CLOSE(4) WRITE(6,'(/a,4(a,F7.2),a/)') ' ARPS domain ', & 'SW: (',lat1,',',lon1, ') NE: (',lat2,',',lon2,')' ! !----------------------------------------------------------------------- ! ! Initializing GEMPAK ! !----------------------------------------------------------------------- ! ! notopn = .true. DO ntime=1,nf CALL in_bdta(iret) ! initializing GEMPAK sans TAE chproj = cproj(mapproj) angle2 = trulon ! ! modify projection if polar stereographic ! IF (chproj == 'STR') THEN angle1=90.0 angle3=0.0 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 (*,*) '**********************************' ELSE angle1=trulat1 angle3=trulat2 WRITE (*,*) '**********************************' WRITE (*,*) '**********************************' WRITE (*,*) ' ' WRITE (*,*) ' SETTING ANGLE1=TRULAT1 ' WRITE (*,*) ' SETTING ANGLE3=TRULAT2 ' WRITE (*,*) ' ' WRITE (*,*) '**********************************' WRITE (*,*) '**********************************' 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, nx,ny,lat1,lon1,lat2,lon2, & angle1,angle2,angle3,angflg, & rnvblk, iret) IF( iret /= 0 ) GO TO 950 ! !----------------------------------------------------------------------- ! ! Build analysis block ! !----------------------------------------------------------------------- ! ! CALL gr_mban ( deltan, deltax, deltay, & ! gbnds, gbnds, gbnds, anlblk, iret ) ! ! IF( iret /= 0 ) GO TO 950 ! ! Read 2D Fields ... filnam = trim(binheader)//'.sfpres'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),psf) print *,trim(filnam) filnam = trim(binheader)//'.mspres'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),mslp) print *,trim(filnam) filnam = trim(binheader)//'.accppt'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),accppt) print *,trim(filnam) filnam = trim(binheader)//'.pwat__'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),pw) print *,trim(filnam) CALL a3dmax0(pw,1,nx,1,nx-1,1,ny,1,ny,1,1,1,1, amax,amin) print *,'PWAT -- amax, amin: ',amax, amin filnam = trim(binheader)//'.temp2m'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),temp2m) print *,trim(filnam) filnam = trim(binheader)//'.dewp2m'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),dewp2m) print *,trim(filnam) filnam = trim(binheader)//'.qv2m__'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),qv2m) print *,trim(filnam) filnam = trim(binheader)//'.u10m__'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),u10m) print *,trim(filnam) filnam = trim(binheader)//'.v10m__'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),v10m) print *,trim(filnam) filnam = trim(binheader)//'.hgtsfc'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),hgtsfc) print *,trim(filnam) ! filnam = trim(binheader)//'.vrt500'//tlevel(ntime) ! CALL binread2d(nx,ny,trim(filnam),vort500) ! print *,trim(filnam) filnam = trim(binheader)//'.ref1km'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),refl1km) print *,trim(filnam) filnam = trim(binheader)//'.ref4km'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),refl4km) print *,trim(filnam) filnam = trim(binheader)//'.cmpref'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),cmpref) print *,trim(filnam) filnam = trim(binheader)//'.sbcape'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),cape) print *,trim(filnam) filnam = trim(binheader)//'.mucape'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),mcape) print *,trim(filnam) filnam = trim(binheader)//'.sbcins'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),cin) print *,trim(filnam) filnam = trim(binheader)//'.mucins'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),mcin) print *,trim(filnam) filnam = trim(binheader)//'.sblcl_'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),lcl) print *,trim(filnam) filnam = trim(binheader)//'.srh01_'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),srh01) print *,trim(filnam) filnam = trim(binheader)//'.srh03_'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),srh03) print *,trim(filnam) filnam = trim(binheader)//'.vhel__'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),uh25) print *,trim(filnam) filnam = trim(binheader)//'.shr01_'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),sh01) print *,trim(filnam) filnam = trim(binheader)//'.shr06_'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),sh06) print *,trim(filnam) filnam = trim(binheader)//'.shr06_'//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),sh06) print *,trim(filnam) IF( nprgem > 0 ) THEN DO klev=1,nprgem filnam = trim(binheader)//'.'//hgtptag(klev)//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),hgtp(1,1,klev)) print *,trim(filnam) filnam = trim(binheader)//'.'//uwndptag(klev)//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),uwndp(1,1,klev)) print *,trim(filnam) filnam = trim(binheader)//'.'//vwndptag(klev)//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),vwndp(1,1,klev)) print *,trim(filnam) filnam = trim(binheader)//'.'//wwndptag(klev)//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),wwndp(1,1,klev)) print *,trim(filnam) filnam = trim(binheader)//'.'//tmpptag(klev)//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),tmpp(1,1,klev)) print *,trim(filnam) filnam = trim(binheader)//'.'//qvptag(klev)//tlevel(ntime) CALL binread2d(nx,ny,trim(filnam),sphp(1,1,klev)) print *,trim(filnam) END DO ! nprgem-loop END IF !CALL a3dmax0(tmpp(1,1,1),1,nx,1,nx-1,1,ny,1,ny,1,1,1,1, amax,amin) !print *,'T 850 -- amax, amin: ',amax, amin ! Read in time parameters OPEN(4,file=trim(binheader)//'.time'//tlevel(ntime)) READ(4,*) year,month,day,hour,minute,ctime CLOSE(4) print *,'year,month,day,hour,minute,ctime:' print *,year,month,day,hour,minute,ctime ! Find forecast hour for name construction ihr = int(ctime/3600.) WRITE(ihrc,'(a1,i3.3)') 'f',ihr filename=trim(gemoutheader)//ihrc lenfil = len_trim(filename) WRITE(6,'(a,a/)') & ' GEMPAK output file: ', filename(1:lenfil) !if(0>1) then 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) !end if ! CALL getunit( nchout ) OPEN (UNIT=4,FILE=trim(filename)//"_ready") WRITE (4,'(a)') trim(filename)//"_ready" CLOSE (4) ! CALL retunit (nchout ) ENDDO ! ntime loop STOP 9000 CONTINUE ! I/O errors CLOSE(4) PRINT *, 'IO ERRORS' STOP 950 CONTINUE WRITE(6,'(a)') ' Error setting up GEMPAK file' STOP END PROGRAM bin2gem