c*********************************************************************************** c PROGRAM owi22 c c owi22 will produce ADCIRC wind file from basin-scale/regional wind/pressure field c provided by OWI. It will require a grid file, two wind field files and two c pressure field files. The output file will contain wind velocity in knots and c srface atmospheric pressure in hPa. sb 8/15/2006 c program owi22 ! Declare variables ---------------------------------------------------- implicit none character(100) :: title integer :: nnAtl,nnR,nnB,ne,numAtl,i integer :: j,k,xi,yi real*8, dimension(:,:), allocatable :: uR,vR,pR,uB,vB,pB,Atlonglat real*8, dimension(:), allocatable :: latR,longR,latB,longB integer, dimension(:), allocatable :: nAtl real*8 :: Along, Alat real*8 :: ramp,rampfrac real*8 :: uu,vv,PP real*8 :: Penv integer :: date1R,date2R,date1B,date2B integer :: date1w,date2w integer :: date1p,date2p integer :: iLatR,iLongR,iCYMDHR,iMinR real*8 :: dxR,dyR,swlatR,swlongR integer :: iLatB,iLongB,iCYMDHB,iMinB real*8 :: dxB,dyB,swlatB,swlongB integer :: iLatw,iLongw,iCYMDHw,iMinw real*8 :: dxw,dyw,swlatw,swlongw integer :: iLatp,iLongp,iCYMDHp,iMinp real*8 :: dxp,dyp,swlatp,swlongp integer :: isnapR,updateR integer :: isnapB,updateB integer :: regionExists integer,allocatable :: swpointsR(:,:) integer,allocatable :: swpointsB(:,:) real*8 :: w,w1,w2,w3,w4 real*8,allocatable :: wR(:,:) real*8,allocatable :: wB(:,:) real*8,allocatable :: umax(:),pmin(:) CHARACTER FNAME*60 logical found logical SkipRegion integer idummy ! Set a value for Penv Penv=1013d0 !mb ! Open Atlantic grid file ---------------------------------------------- 31 WRITE(*,*) 'Enter the name of the ADCIRC UNIT 14 file:' READ(*,'(A60)') FNAME INQUIRE(FILE=FNAME,EXIST=FOUND) IF(FOUND) GOTO 32 WRITE(*,1010) FNAME GOTO 31 32 WRITE(*,1011) FNAME 1010 FORMAT(' File ',A60,/,' WAS NOT FOUND! Try again',/) 1011 FORMAT(' File ',A60,/,' WAS FOUND! Opening & Processing file',/) open(1,file=fname) read(1,*) title read(1,*) ne, nnAtl allocate(Atlonglat(nnAtl,2)) allocate(nAtl(nnAtl)) allocate(umax(nnAtl),pmin(nnAtl)) allocate(swpointsB(nnAtl,2),wB(nnAtl,4)) allocate(swpointsR(nnAtl,2),wR(nnAtl,4)) do i=1,nnAtl read(1,*) numAtl, Along, Alat Atlonglat(i,1)=Along Atlonglat(i,2)=Alat nAtl(i)=numAtl umax(i)=0 pmin(i)=Penv end do ! Create ADCIRC output file -------------------------------------------- 41 WRITE(*,*) 'Enter the name of the output unit 22 file:' READ(*,'(A60)') FNAME open(2,file=fname) 1000 format(i8,3e13.5) ! Ask if there is region-scale files ---------------------------------- 46 WRITE(*,*) 'Do you have region-scale files? (YES:1,NO:0)' READ(*,*) idummy IF(idummy.eq.1) THEN SkipRegion = .FALSE. ELSE IF(idummy.eq.0) THEN SkipRegion = .TRUE. ELSE GOTO 46 ENDIF ! Oceanweather basin win file ------------------------------------------------- 51 WRITE(*,*) 'Enter the name of the OWI Basin .win file:' READ(*,'(A60)') FNAME INQUIRE(FILE=FNAME,EXIST=FOUND) IF(FOUND) GOTO 52 WRITE(*,1010) FNAME GOTO 51 52 WRITE(*,1011) FNAME open(31,file=fname) ! Read in begining/ending dates of win file read(31,10) title,date1w,title,date2w 10 format(a55,i10,a5,i10) ! Oceanweather basin pre file ------------------------------------------------ 61 WRITE(*,*) 'Enter the name of the OWI Basin .pre file:' READ(*,'(A60)') FNAME INQUIRE(FILE=FNAME,EXIST=FOUND) IF(FOUND) GOTO 62 WRITE(*,1010) FNAME GOTO 61 62 WRITE(*,1011) FNAME open(41,file=fname) ! Read in begining/ending dates of win file read(41,10) title,date1p,title,date2p if(date1w.ne.date1p.or.date2w.ne.date2p) then print *, 'HEADER INFO IN WIN AND PRE FILES DO NOT MATCH' print *, 'EXECUTION WILL BE TERMINATED' stop endif date1B = date1w date2B = date2w ! Oceanweather region win file ----------------------------------------------- IF(SkipRegion) GOTO 83 71 WRITE(*,*) 'Enter the name of the OWI Region .win file:' READ(*,'(A60)') FNAME INQUIRE(FILE=FNAME,EXIST=FOUND) IF(FOUND) GOTO 72 WRITE(*,1010) FNAME GOTO 71 72 WRITE(*,1011) FNAME open(32,file=fname) ! Read in begining/ending dates of .win file read(32,10) title,date1w,title,date2w ! Oceanweather region pre file ----------------------------------------------- 81 WRITE(*,*) 'Enter the name of the OWI Region .pre file:' READ(*,'(A60)') FNAME INQUIRE(FILE=FNAME,EXIST=FOUND) IF(FOUND) GOTO 82 WRITE(*,1010) FNAME GOTO 81 82 WRITE(*,1011) FNAME open(42,file=fname) ! Read in begining/ending dates of .pre file read(42,10) title,date1p,title,date2p if(date1w.ne.date1p.or.date2w.ne.date2p) then print *, 'HEADER INFO IN WIN AND PRE FILES DO NOT MATCH' print *, 'EXECUTION WILL BE TERMINATED' stop endif date1R = date1w date2R = date2w 83 CONTINUE ! Open first file ------------------------------------------------------------ isnapB = 0 isnapR = 0 updateB = 1 updateR = 1 do while(.true.) !Process Basin Info --------------------------------------------------------- ! Increment counter isnapB = isnapB+1 11 format(t6,i4,t16,i4,t23,f6.0,t32,f6.0, & t44,f8.0,t58,f8.0,t69,i10,i2) ! Read Grid Specifications/Date read (31,11,end=9999) & iLatw,iLongw,dxw,dyw,swlatw,swlongw,iCYMDHw,iMinw ! Read Grid Specifications/Date read (41,11,end=9999) & iLatp,iLongp,dxp,dyp,swlatp,swlongp,iCYMDHp,iMinp if(iLatw.ne.iLatp.or.iLongw.ne.iLongp.or.dxw.ne.dxp.or. & dyw.ne.dyp.or.swlatw.ne.swlatp.or.swlongw.ne.swlongp.or. & iCYMDHw.ne.iCYMDHp.or.iMinw.ne.iMinp) then print *, 'SNAPSHOT HEADER IN WIN AND PRE FILES DO NOT MATCH' print *, 'EXECUTION WILL BE TERMINATED' stop endif ! Check if header info has changed from the previous snapshot if(isnapB.gt.1) then if(iLatw.ne.iLatB.or.iLongw.ne.iLongB.or.dxw.ne.dxB.or. & dyw.ne.dyB.or.swlatw.ne.swlatB.or. & swlongw.ne.swlongB) then updateB = 1 else updateB = 0 endif endif iCYMDHB = iCYMDHw iMinB = iMinw if(updateB.eq.1) then print *, 'BASIN COORDINATE UPDATED' iLatB = iLatw iLongB = iLongw dxB = dxw dyB = dyw swlatB = swlatw swlongB = swlongw ! Number of grid points nnB=iLatB*iLongB ! Allocate and create matrices if(isnapB.ne.1) then deallocate(uB,vB,pB,longB,latB) endif allocate(uB(iLongB,iLatB),vB(iLongB,iLatB),pB(iLongB,iLatB)) allocate(longB(iLongB),latB(iLatB)) ! Generate long&lat on each grid point do i=1,iLatB latB(i) = swlatB+(i-1)*dyB enddo do i=1,iLongB longB(i) = swlongB+(i-1)*dxB enddo ! Generate interpolation info (south west point and weights) do i=1,nnAtl if (Atlonglat(i,1)>=longB(1) .and. & Atlonglat(i,1)=latB(1) .and. & Atlonglat(i,2)=longB(j) .and. & Atlonglat(i,1)=latB(k) .and. & Atlonglat(i,2)=longR(1) .and. & Atlonglat(i,1)=latR(1) .and. & Atlonglat(i,2)=longR(j) .and. & Atlonglat(i,1)=latR(k) .and. & Atlonglat(i,2)