! Copyright 2014 College of William and Mary
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
! http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!****************************************************************************************
!
! Extract binary results in a part of hgrid.gr3
! Works for mixed tri/quads outputs.
! Inputs: binary file and hgrid.gr3; subflag.gr3 (hgrid.gr3 based)
! Outputs: ?_sub_<file63> (new binary file); fort.13 (output grid
! for mlab)
!
! ifort -Bstatic -assume byterecl -O3 -o extract_subregion2 extract_subregion2.f90
! pgf90 -O2 -mcmodel=medium -Mbounds -o extract_subregion2 extract_subregion2.f90
!
! History: use subflag.gr3 instead of subgrid.gr3 so much faster than extract_subregion
!
program read_out
parameter(nbyte=4)
character*30 file63
character*12 it_char
character*48 start_time,version,variable_nm,variable_dim
character*48 data_format
allocatable :: sigma(:),cs(:),ztot(:),x(:),y(:),dp(:),kbp00(:)
integer, allocatable :: i34(:),elnode(:,:),i34out(:),nmout(:,:)
allocatable :: icum(:,:,:),eta2(:),imapelem2(:)
allocatable :: xout(:),yout(:),dpout(:),imapnode(:),imapnode2(:),ifl(:)
!long int for large files
integer(kind=8) :: irec,irecout
print*, 'Input file to read from (without *_):'
read(*,'(a30)')file63
print*, 'Input start and end file #:'
read(*,*) istart,iend
! Read hgrid.gr3
open(14,file='hgrid.gr3',status='old')
read(14,*)
read(14,*)ne,np
allocate(x(np),y(np),dp(np),kbp00(np),i34(ne),elnode(4,ne),i34out(ne),nmout(4,ne), &
&eta2(np),xout(np),yout(np),dpout(np),imapnode(np),imapnode2(np),ifl(np),imapelem2(ne))
do i=1,np
read(14,*)j,x(i),y(i),dp(i)
enddo !i
do i=1,ne
read(14,*)j,i34(i),elnode(1:i34(i),i)
enddo !i
close(14)
! Read subflag.gr3
open(14,file='subflag.gr3',status='old')
read(14,*)
read(14,*)ne2,np2
if(np2/=np) stop 'subflag.gr3 should be based on hgrid.gr3'
!'
do i=1,np
read(14,*)j,xtmp,ytmp,tmp
ifl(i)=tmp
enddo !i
close(14)
!... Find nodes in the original grid; create mapping
imapnode=0 !from new to original grid
imapnode2=0 !from original to new grid
imapelem2=0 !from original to new grid
neout=0
npout=0
do i=1,ne
itmp=minval(ifl(elnode(1:i34(i),i)))
if(itmp/=0) then
neout=neout+1
imapelem2(i)=neout
i34out(neout)=i34(i)
do j=1,i34(i)
nd=elnode(j,i)
if(imapnode2(nd)==0) then !new node
npout=npout+1
imapnode2(nd)=npout
imapnode(npout)=nd
xout(npout)=x(nd)
yout(npout)=y(nd)
dpout(npout)=dp(nd)
endif !imapnode2
enddo !j
nmout(1:i34out(neout),neout)=imapnode2(elnode(1:i34(i),i))
endif !itmp
enddo !i=1,ne
! Output grid
write(13,*)
write(13,*)neout,npout
do i=1,npout
write(13,*)i,xout(i),yout(i),dpout(i)
enddo !i
do i=1,neout
write(13,*)i,i34out(i),nmout(1:i34out(i),i)
enddo !i
close(13)
print*, '# of nodes in output grid= ',npout
print*, 'Starting time iteration...'
!... Time iteration
!...
it_tot=(istart-1)*nrec
do iday=istart,iend
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!... Header
!...
write(it_char,'(i12)')iday
open(63,file=it_char//'_'//file63,status='old',access='direct',recl=nbyte)
open(65,file=it_char//'_sub_'//file63,access='direct',recl=nbyte)
print*, 'Doing ',it_char//'_'//file63
irec=0
irecout=0
do m=1,48/nbyte
read(63,rec=irec+m) data_format(nbyte*(m-1)+1:nbyte*m)
write(65,rec=irecout+m) data_format(nbyte*(m-1)+1:nbyte*m)
enddo
if(data_format.ne.'DataFormat v5.0') then
print*, 'This code reads only v5.0: ',data_format
stop
endif
irec=irec+48/nbyte
irecout=irecout+48/nbyte
do m=1,48/nbyte
read(63,rec=irec+m) version(nbyte*(m-1)+1:nbyte*m)
write(65,rec=irecout+m) version(nbyte*(m-1)+1:nbyte*m)
enddo
irec=irec+48/nbyte
irecout=irecout+48/nbyte
do m=1,48/nbyte
read(63,rec=irec+m) start_time(nbyte*(m-1)+1:nbyte*m)
write(65,rec=irecout+m) start_time(nbyte*(m-1)+1:nbyte*m)
enddo
irec=irec+48/nbyte
irecout=irecout+48/nbyte
do m=1,48/nbyte
read(63,rec=irec+m) variable_nm(nbyte*(m-1)+1:nbyte*m)
write(65,rec=irecout+m) variable_nm(nbyte*(m-1)+1:nbyte*m)
enddo
irec=irec+48/nbyte
irecout=irecout+48/nbyte
do m=1,48/nbyte
read(63,rec=irec+m) variable_dim(nbyte*(m-1)+1:nbyte*m)
write(65,rec=irecout+m) variable_dim(nbyte*(m-1)+1:nbyte*m)
enddo
irec=irec+48/nbyte
irecout=irecout+48/nbyte
! write(*,'(a48)')data_format
! write(*,'(a48)')version
! write(*,'(a48)')start_time
! write(*,'(a48)')variable_nm
! write(*,'(a48)')variable_dim
read(63,rec=irec+1) nrec
read(63,rec=irec+2) dtout
read(63,rec=irec+3) nspool
read(63,rec=irec+4) ivs
read(63,rec=irec+5) i23d
irec=irec+5
write(65,rec=irecout+1) nrec
write(65,rec=irecout+2) dtout
write(65,rec=irecout+3) nspool
write(65,rec=irecout+4) ivs
write(65,rec=irecout+5) i23d
irecout=irecout+5
! print*, 'i23d=',i23d,' nrec= ',nrec
! Vertical grid
read(63,rec=irec+1) nvrt
read(63,rec=irec+2) kz
read(63,rec=irec+3) h0
read(63,rec=irec+4) h_s
read(63,rec=irec+5) h_c
read(63,rec=irec+6) theta_b
read(63,rec=irec+7) theta_f
irec=irec+7
write(65,rec=irecout+1) nvrt
write(65,rec=irecout+2) kz
write(65,rec=irecout+3) h0
write(65,rec=irecout+4) h_s
write(65,rec=irecout+5) h_c
write(65,rec=irecout+6) theta_b
write(65,rec=irecout+7) theta_f
irecout=irecout+7
if(iday==istart) then
allocate(sigma(nvrt),cs(nvrt),ztot(nvrt),icum(np,nvrt,2))
endif
! print*, 'nvrt=',nvrt,' theta_f= ',theta_f
do k=1,kz-1
read(63,rec=irec+k) ztot(k)
write(65,rec=irecout+k) ztot(k)
enddo
do k=kz,nvrt
kin=k-kz+1
read(63,rec=irec+k) sigma(kin)
write(65,rec=irecout+k) sigma(kin)
! cs(kin)=(1-theta_b)*sinh(theta_f*sigma(kin))/sinh(theta_f)+ &
! &theta_b*(tanh(theta_f*(sigma(kin)+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5)
enddo
irec=irec+nvrt
irecout=irecout+nvrt
! print*, 'Done vgrid...'
! Horizontal grid
read(63,rec=irec+1) np
read(63,rec=irec+2) ne
irec=irec+2
! if(np.gt.mnp.or.ne.gt.mne) then
! write(*,*)'Too many nodes/elements',np,ne
! stop
! endif
do m=1,np
! read(63,rec=irec+1)x(m)
! read(63,rec=irec+2)y(m)
! read(63,rec=irec+3)dp(m)
read(63,rec=irec+4)kbp00(m)
irec=irec+4
enddo !m=1,np
do m=1,ne
! read(63,rec=irec+1)i34
irec=irec+1
do mm=1,i34(m)
! read(63,rec=irec+1)elnode(mm,m)
irec=irec+1
enddo !mm
enddo !m
write(65,rec=irecout+1) npout
write(65,rec=irecout+2) neout
irecout=irecout+2
do m=1,npout
write(65,rec=irecout+1)xout(m)
write(65,rec=irecout+2)yout(m)
write(65,rec=irecout+3)dpout(m)
write(65,rec=irecout+4)kbp00(imapnode(m))
irecout=irecout+4
enddo !m=1,npout
do m=1,neout
write(65,rec=irecout+1)i34out(m)
irecout=irecout+1
do mm=1,i34out(m)
write(65,rec=irecout+1)nmout(mm,m)
irecout=irecout+1
enddo !mm
enddo !m
! print*, 'last element: ',nmout(1:i34(neout),neout)
!... Compute relative record # for a node and level for 3D outputs
!...
if(iday==istart) then
icount=0
do i=1,np
do k=max0(1,kbp00(i)),nvrt
do m=1,ivs
icount=icount+1
icum(i,k,m)=icount
enddo !m
enddo !k
enddo !i=1,np
endif !iday==istart
do it1=1,nrec
read(63,rec=irec+1) time
read(63,rec=irec+2) it
irec=irec+2
it_tot=it_tot+1
time=it_tot*dtout
write(65,rec=irecout+1) time
write(65,rec=irecout+2) it_tot
irecout=irecout+2
! print*, 'time=',time/86400
do i=1,np
read(63,rec=irec+i) eta2(i)
enddo !i
irec=irec+np
do i=1,npout
write(65,rec=irecout+i) eta2(imapnode(i))
enddo !i
irecout=irecout+npout
if(i23d==2) then
do i=1,npout
do m=1,ivs
read(63,rec=irec+(imapnode(i)-1)*ivs+m) tmp
write(65,rec=irecout+1) tmp
irecout=irecout+1
enddo !m
enddo !i
irec=irec+np*ivs
else !i23d=3
do i=1,npout
nd=imapnode(i)
do k=max0(1,kbp00(nd)),nvrt
do m=1,ivs
read(63,rec=irec+icum(nd,k,m)) tmp
write(65,rec=irecout+1) tmp
irecout=irecout+1
enddo !m
enddo !k
enddo !i
irec=irec+icum(np,nvrt,ivs)
endif !i23d
enddo !it1=1,nrec
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
enddo !iday
stop
end