Skip to content

Commit

Permalink
Add new scripts to generate tvd.prop for hybrid WENO/ELM
Browse files Browse the repository at this point in the history
  • Loading branch information
@josephzhang8
josephzhang8 committed May 19, 2023
1 parent 7e83f6f commit 9d7230e
Show file tree
Hide file tree
Showing 3 changed files with 343 additions and 24 deletions.
132 changes: 132 additions & 0 deletions src/Utility/Post-Processing-Fortran/combine_TVD.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
! Combine 2 tvd.prop.[1,2] based a specified region (final flag is
! from '2' inside region2.rgn). The 2 TVD files may be generated
! from gen_tvd_WENO.f90.
! Inputs:
! hgrid.gr3 (in any projection or lon/lat; b.c. part not needed)
! tvd.prop.[1,2]
! region2.rgn (gredit region format)
! Output: tvd.prop

! ifort -O2 -o combine_TVD UtilLib/pt_in_poly_test.f90 combine_TVD.f90

use pt_in_poly_test
implicit real*8(a-h,o-z)
integer :: nwild(3)
integer, allocatable :: i34(:),elnode(:,:),nne(:),indel(:,:),nnp(:), &
&indnd(:,:),isbnd(:),itvd(:),itvd1(:),itvd2(:)
real*8, allocatable :: xnd(:),ynd(:),dp(:),area(:),xpoly(:),ypoly(:)

open(14,file='hgrid.gr3',status='old')
read(14,*)
read(14,*) ne,np
allocate(xnd(np),ynd(np),dp(np),area(ne),i34(ne),elnode(4,ne),nne(np), &
&isbnd(np),itvd(ne),itvd1(ne),itvd2(ne))

do i=1,np
read(14,*) j,xnd(i),ynd(i),dp(i)
enddo !i

do i=1,ne
read(14,*) j,i34(i),elnode(1:i34(i),i)

n1=elnode(1,i)
n2=elnode(2,i)
n3=elnode(3,i)
area(i)=signa(xnd(n1),xnd(n2),xnd(n3),ynd(n1),ynd(n2),ynd(n3))
if(area(i)<=0) then
write(*,*)'Negative area at elem:',i
stop
endif

if(i34(i)==4) then
n4=elnode(4,i)
tmp=signa(xnd(n1),xnd(n3),xnd(n4),ynd(n1),ynd(n3),ynd(n4))
if(tmp<=0) then
write(*,*)'Negative area at elem:',i
stop
endif
area(i)=area(i)+tmp
endif
enddo !i=1,ne
close(14)

!Read in .rgn
open(10,file='region2.rgn',status='old')
read(10,*); read(10,*)npoly
if(npoly/=1) stop 'can only have 1 poly in .rgn'
read(10,*)nvertices
nvertices=nvertices+1 !repeat 1st vertex
allocate(xpoly(nvertices),ypoly(nvertices))
do i=1,nvertices-1
read(10,*)xpoly(i),ypoly(i)
enddo !i
xpoly(nvertices)=xpoly(1)
ypoly(nvertices)=ypoly(1)
close(10)

! Neighborhood
nne=0
do i=1,ne
do j=1,i34(i)
nd=elnode(j,i)
nne(nd)=nne(nd)+1
enddo
enddo !i
mnei=maxval(nne)

allocate(indel(mnei,np))
nne=0
! nnp=0 !estimate 1st
do i=1,ne
do j=1,i34(i)
nd=elnode(j,i)
nne(nd)=nne(nd)+1
if(nne(nd)>mnei) stop 'impossible'
indel(nne(nd),nd)=i
! nnp(nd)=nnp(nd)+i34(i)-1
enddo
enddo !i

! Output
open(8,file='tvd.prop.1',status='old')
open(9,file='tvd.prop.2',status='old')
do i=1,ne
read(8,*)j,tmp1
read(9,*)j,tmp2
itvd1(i)=nint(tmp1)
itvd2(i)=nint(tmp2)
enddo !i

open(12,file='tvd.prop',status='replace')
small1=1.d-4
ray_angle=3.13192d0
itvd(:)=1 !init
do i=1,ne
!Beware lon jumps
xctr=sum(xnd(elnode(1:i34(i),i)))/i34(i)
yctr=sum(ynd(elnode(1:i34(i),i)))/i34(i)
call pt_in_poly_ray_method_double(nvertices,small1,ray_angle,xpoly,ypoly,xctr,yctr,in_out,inters,npoly)
if(in_out==1) then
itvd(i)=itvd2(i)
else
itvd(i)=itvd1(i)
endif !in_out
enddo !i

do i=1,ne
write(12,*)i,itvd(i)
enddo !i
close(12)

stop
end

function signa(x1,x2,x3,y1,y2,y3)
!... Compute signed area formed by pts 1,2,3
implicit real*8(a-h,o-z)

signa=((x1-x3)*(y2-y3)-(x2-x3)*(y1-y3))/2

return
end

71 changes: 47 additions & 24 deletions src/Utility/Post-Processing-Fortran/compute_average5.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
! Skip dry times for 3D variables.
! Works for mixed quad/tri on NODE based variables only.
! Input: out2d* and corresponding nc file if the variable is 3D; vgrid.in; screen
! Output: average.out (gredit fromat)
! Output: average.out (gredit or xyuv format)
!
! ifort -O2 -mcmodel=medium -assume byterecl -CB -o compute_average5.exe ../UtilLib/extract_mod2.f90 ../UtilLib/compute_zcor.f90 compute_average5.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff
!********************************************************************************
Expand All @@ -31,28 +31,32 @@ program read_out
use extract_mod2
use compute_zcor
parameter(nbyte=4)
character(len=30) :: file62,file63,varname
character(len=30) :: file62,file63,file64,varname,varname2
character(len=12) :: it_char
character(len=48) :: data_format

integer,allocatable :: i34(:),elnode(:,:)
integer :: iday1, iday2, iskipst
allocatable :: sigma(:),cs(:),ztot(:)
allocatable :: outvar(:,:),icum(:,:,:),eta2(:),ztmp(:,:)
allocatable :: outvar(:,:,:),icum(:,:,:),eta2(:),ztmp(:,:)
allocatable :: xcj(:),ycj(:),dp(:),dps(:),kbs(:),xctr(:),yctr(:),dpe(:),kbe(:)
allocatable :: idry(:),outs(:,:,:),residual(:,:),icounter(:)
allocatable :: ic3(:,:),isdel(:,:),isidenode(:,:),kbp(:),sigma_lcl(:,:)
real :: wild(2)
real*8,allocatable :: timeout(:),xnd(:),ynd(:)
integer :: char_len,start_2d(2),start_3d(3),start_4d(4), &
&count_2d(2),count_3d(3),count_4d(4),dimids(100),idims(100)

pi=3.1415926
ivs=1 !for vectors, use own script to combine components

print*, 'Input variable name (e.g. salinity):'
print*, 'Input variable name (e.g. horizontalVelX):'
read(*,'(a30)')varname
varname=adjustl(varname); len_var=len_trim(varname)

print*, 'Is the var scalar (1) or vector (2)?'
read(*,*) ivs
if(ivs/=1.and.ivs/=2) stop 'check ivs'

print*, 'Input start and end stack #s to read:'
read(*,*) iday1,iday2

Expand Down Expand Up @@ -87,14 +91,13 @@ program read_out
!... Read in time records in segments for mem
last_dim=np
allocate(ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),timeout(nrec), &
&outvar(nvrt,last_dim),eta2(np),ztmp(nvrt,np),residual(np,2),icounter(np),idry(np))
outvar=-huge(1.0) !test mem
&outvar(nvrt,last_dim,ivs),eta2(np),ztmp(nvrt,np),residual(np,2),icounter(np),idry(np))
outvar=-99. !test mem

call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp)

!... Time iteration
!...
outvar=-99
ztmp=-99
residual=0
icounter=0 !counter for each node (wet/dry)
Expand All @@ -112,12 +115,23 @@ program read_out

!Find nc file
file63=varname(1:len_var)//'_'//it_char(1:leng)//'.nc'
if(ivs==2) varname2=varname(1:len_var-1)//'Y'
inquire(file=file63,exist=lexist)
if(lexist) then
i23d=2 !3D var
else
i23d=1 !2D
if(lexist) then !3D var
i23d=2
if(ivs==2) then
file64=varname2(1:len_var)//'_'//it_char(1:leng)//'.nc'
inquire(file=file64,exist=lexist2)
if(.not.lexist2) then
print*, 'Missing y-component:',file64
print*, file63
stop
endif
endif
else !2D
i23d=1
file63=file62
file64=file62
endif

iret=nf90_open(trim(adjustl(file63)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid)
Expand All @@ -133,6 +147,13 @@ program read_out
!'
if(idims(ndims)/=nrec) stop 'last dim is not time'

if(ivs==2) then !vector
iret=nf90_open(trim(adjustl(file64)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid2)
if(iret/=nf90_NoErr) stop 'Failed to open file64'
iret=nf90_inq_varid(ncid2,varname2(1:len_var),ivarid2)
if(iret/=nf90_NoErr) stop 'Var2 not found'
endif !ivs

do irec=1,nrec
!----------------------------------------------------------------------------
!Get elev
Expand All @@ -144,21 +165,23 @@ program read_out
if(i23d==1) then !2D
start_2d(1)=1; start_2d(2)=irec
count_2d(1)=npes; count_2d(2)=1
iret=nf90_get_var(ncid,ivarid1,outvar(1,1:npes),start_2d,count_2d)
iret=nf90_get_var(ncid,ivarid1,outvar(1,1:npes,1),start_2d,count_2d)
if(ivs==2) iret=nf90_get_var(ncid2,ivarid2,outvar(1,1:npes,2),start_2d,count_2d)
else !3D
start_3d(1:2)=1; start_3d(3)=irec
count_3d(1)=nvrt; count_3d(2)=npes; count_3d(3)=1
iret=nf90_get_var(ncid,ivarid1,outvar(:,1:npes),start_3d,count_3d)
iret=nf90_get_var(ncid,ivarid1,outvar(:,1:npes,1),start_3d,count_3d)
if(ivs==2) iret=nf90_get_var(ncid2,ivarid2,outvar(:,1:npes,2),start_3d,count_3d)
endif

!Available now: outvar(nvrt,np|ne), eta2(np)
!Available now: outvar(nvrt,np,ivs), eta2(np)

if(i23d==1) then !2D
!print*,'irec=',irec,iday1,irec_start,iday2,irec_end
if(.not.(iday==iday1.and.irec<irec_start.or.iday==iday2.and.irec>irec_end)) then
do i=1,np
icounter(i)=icounter(i)+1
residual(i,1)=residual(i,1)+outvar(1,i)
residual(i,:)=residual(i,:)+outvar(1,i,:)
enddo !i
endif
else !3D
Expand Down Expand Up @@ -204,19 +227,19 @@ program read_out
icounter(i)=icounter(i)+1

if(z00<=-1.e8) then !depth average
av=0.
wild=0.
do k=kbp(i),nvrt-1
av=av+(outvar(k,i)+outvar(k+1,i))*0.5*(ztmp(k+1,i)-ztmp(k,i))
wild(1:2)=wild(1:2)+(outvar(k,i,:)+outvar(k+1,i,:))*0.5*(ztmp(k+1,i)-ztmp(k,i))
enddo !k
toth=ztmp(nvrt,i)-ztmp(kbp(i),i)
if(toth<=0.) then
write(*,*)'Negative depth at node:',i,toth
stop
endif
residual(i,1)=residual(i,1)+av/toth
residual(i,:)=residual(i,:)+wild(1:2)/toth
else !not depth average
do m=1,ivs
tmp=outvar(k0,i)*(1-rat)+outvar(k0+1,i)*rat
tmp=outvar(k0,i,m)*(1-rat)+outvar(k0+1,i,m)*rat
residual(i,m)=residual(i,m)+tmp
enddo !m
endif !z00
Expand Down Expand Up @@ -265,10 +288,10 @@ program read_out
do i=1,ne
write(65,*)i,i34(i),elnode(1:i34(i),i)
enddo !i
! else !vectors
! do i=1,np
! write(65,*)x(i),y(i),residual(i,1:2)
! enddo !i
else !vectors
do i=1,np
write(65,*)xnd(i),ynd(i),residual(i,1:2)
enddo !i
endif !ivs
close(65)

Expand Down
Loading

0 comments on commit 9d7230e

Please sign in to comment.