! !**************************************************************************************** ! * ! Read binary format v5.0 (hybrid S-Z models) of *hvel.67 (sidecenter vel.) and ! compute element volumes and fluxes between elements. ! WARNING: does not work for lat/lon! ! ! Usage: ! %proc MODEL_DATA_PATH OUTPUT_PATH FIRST_STACK LAST_STACK ! ! Inputs: *_elev.61, *hvel.67, *temp.70, *salt.70, flux_regions.gr3 (depth=0,1,...M; ! element flags based on max. of 3 node depths); screen inputs ! Outputs: fluxes.out ! volume.out ! History: !**************************************************************************************** ! COMPILATION ! ifort -Bstatic -O2 -assume byterecl -o selfe-delwaq selfe-delwaq.f90 selfe_geometry.f90 waq_hyd.f90 ! wrselafin.f90 wrwaqfil.f90 wrwaqhyd.f90 juldat.f90 timdat.f90 ! gfortran -Bstatic -Wall -O2 -o selfe-delwaq selfe-delwaq.f90 selfe_geometry.f90 waq_hyd.f90 wrselafin.f90 ! wrwaqfil.f90 wrwaqhyd.f90 juldat.f90 timdat.f90 program read_out use timers implicit none integer, parameter :: nbyte=4 integer(4), allocatable :: idry (:) !< dry-flag for nodes integer(4), allocatable :: idry_s (:) !< dry-flag for sides (edges) integer(4), allocatable :: idry_e (:) !< dry flag for elements integer(4), allocatable :: kbp00 (:) !< actual bottom located at index kbp00 (or kbp or kbp2) integer(4), allocatable :: kbs (:) !< actual bottom located at index for sides integer(4), allocatable :: kbe (:) !< actual bottom located at index for elements integer(4), allocatable :: kbp2 (:) !< actual bottom located at index kbp00 (or kbp or kbp2) integer(4), allocatable :: kbs2 (:) !< actual bottom located at index for sides integer(4), allocatable :: kbe2 (:) !< actual bottom located at index for elements integer(4), allocatable :: nm (:,:) !< (ne,4) node nrs. (connectivity table) integer(4), allocatable :: openbndnodes(:) !< open boundary nodes integer(4), allocatable :: isopenbnd (:) !< back pointer for nodes real (4), allocatable :: zs (:,:) !< real (4), allocatable :: ze (:,:) !< real (4), allocatable :: ztot (:) !< z-coordinates of each Z-level real (4), allocatable :: ztmp (:,:) !< real (4), allocatable :: zeOld (:,:) !< integer(4), allocatable :: ic3 (:,:) !< integer(4), allocatable :: js (:,:) !< integer(4), allocatable :: is (:,:) !< integer(4), allocatable :: isidenode (:,:) !< real (4), allocatable :: sigma (:) !< sigma-coordinates of each S-level real (4), allocatable :: cs (:) !< real (4), allocatable :: dp (:) !< depths of the nodes below reference level real (4), allocatable :: dps (:) !< real (4), allocatable :: dpe (:) !< real (4), allocatable :: x (:) !< x-values of the nodes real (4), allocatable :: y (:) !< y-values of the nodes real (4), allocatable :: xcj (:) !< real (4), allocatable :: ycj (:) !< real (4), allocatable :: eta2 (:) !< real (4), allocatable :: suv21 (:,:) !< real (4), allocatable :: suv22 (:,:) !< real (4), allocatable :: area_el (:) !< real (4), allocatable :: area_e (:) !< real (4), allocatable :: sideuv (:,:) !< real (4), allocatable :: area_cell (:) !< real (4), allocatable :: sne (:,:) !< real (4), allocatable :: dihv (:,:) !< real (4), allocatable :: centeru (:) !< real (4), allocatable :: centerv (:) !< real (8), allocatable :: hfluxes (:,:) !< horizontal Selfe fluxes along the sides real (8), allocatable :: vfluxes (:,:) !< vertical fluxes between the layers real (8), allocatable :: nodal_flux (:,:) !< real , allocatable :: flx_adv (:,:,:) !< original horizontal flux !< (1: the local x-driection) and vertical flux (2: positive upward) !< Transformation tensor for side frame: sframe(i,j,isd) real (4), allocatable :: sframe (:,:,:) real (4), allocatable :: distj (:) !< Side lengths real (4), allocatable :: volume (:,:) !< real (4), allocatable :: volumeOld (:,:) !< real (4), allocatable :: volumeError(:,:) !< real (4), allocatable :: cum_influx (:,:) !< integer(4), allocatable :: aggre_point (:) !< Delwaq horizontal aggregation pointer integer(4), allocatable :: verti_pnt (:) !< Delwaq Selfe => Delwaq layer interface pointer integer(4), allocatable :: faces_pnt (:) !< Delwaq Selfe => Delwaq horizontal exchanges pointer integer(4), allocatable :: layer_pnt (:) !< Delwaq Selfe => Delwaq layer number pointer integer(4), allocatable :: attribute (:,:) !< Delwaq attribute array for activity of Z-layers real (4), allocatable :: bfluxes (:,:) !< Delwaq boundary fluxes real (4), allocatable :: wfluxes (:,:) !< Delwaq horizontal fluxes along the sides real (4), allocatable :: uvelos (:,:) !< Delwaq u-velocities in the nodes real (4), allocatable :: vvelos (:,:) !< Delwaq v-velocities in the nodes character(256) path !< path of the hydrodynamic result files character(256) outpath !< path of the Delwaq output files character( 48) data_format !< data format should be 'DataFormat v5.0' character( 48) version !< version string character( 48) start_time !< start time string character( 48) variable_nm !< variable name string character( 48) variable_dim !< variable dimension string integer ( 4) nrec !< nr. of output time steps real ( 4) dtout !< output time step integer ( 4) nspool !< skip integer ( 4) ivs !< ivs (=1 or 2 for scalar or vector) integer ( 4) i23d !< i23d (=2 or 3 for 2D or 3D) integer ( 4) nvrt !< total # of vertical levels integer ( 4) kz !< nr. of Z-levels real ( 4) h0 !< dry threshold real ( 4) h_s !< Z-sigma transition depth real ( 4) h_c !< ~thickness of finely resolved surf layer real ( 4) theta_b !< sigma zoom param: 1: resolve surf+bottom, 0: only surf real ( 4) theta_f !< sigma zoom param: ->0: uniform sigma, >>1: coarser in middle real ( 4) w1, w2, w3, w4 !< help variables for area-weighting of horizontal velocites around nodes integer ( 4) np !< nr. of nodes integer ( 4) ne !< nr. of elements integer ( 4) i34 !< element type (currently must be 3) integer ( 4) nolay !< Delwaq number of layers integer nopenbnds, nopenbndnodes, nopn integer first_stack, last_stack, nstacks real :: xel(3), yel(3) real snx,sny,midx,midy real el_depth, uint, vint real hmod2 real thetan, time integer i, ie,ie1,ie2, ics, istat, itmp, isd integer irec, irec_tmp, itsflux integer irec0, irec0_u, irec0_dihv integer iblock, istack, it integer j,k, kbpl, kin integer n1, n2, n3, n4 integer ne2, ne_e, np2, ns, ns2, ns_e integer l, m, mm character(60) out_prefix character(12) it_char character( 3) stack_char,stack_char2 real (4) fact real (4), allocatable :: dist_dw (:) , area_dw(:,:) integer(4), allocatable :: openbnds(:,:) integer save_time , idummy integer(4) ithndl(20) ! handle to time this subroutine data ithndl / 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & & 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 / call timini ( ) timon = .true. call timstrt( "selfe-delwaq", ithndl(1) ) save_time = -999 call timstrt( "initialisation", ithndl(2) ) open ( 10, file='selfe-delwaq.inp' ) open ( 20, file='selfe-delwaq.log' ) read ( 10, * ) path read ( 10, * ) outpath read ( 10, * ) first_stack read ( 10, * ) last_stack write ( 20, * ) '' write ( 20, * ) 'Path: ', trim(path) write ( 20, * ) 'Output Dir: ', trim(outpath) write ( 20, * ) 'First stack: ', first_stack write ( 20, * ) 'Last stack: ', last_stack nstacks = last_stack-first_stack+1 itsflux = 0 ! 1 -> compute T,S ics = 1 ! Cartesian coordinates write(stack_char,'(i3)') first_stack stack_char = adjustl(stack_char) write(stack_char2,'(i3)') last_stack stack_char2 = adjustl(stack_char2) ! Prefix for output files write(out_prefix,'(a60)') trim(OUTPATH)//'/stacks_'//trim(stack_char)//'-'//trim(stack_char2)//'_' out_prefix= adjustl(out_prefix) ! Open input files write(it_char,'(i12)') first_stack it_char = adjustl(it_char) open ( 63, file=trim(PATH)//trim(it_char)//'_elev.61', status='old', form='binary' ) open ( 67, file=trim(PATH)//trim(it_char)//'_dihv.66', status='old', form='binary' ) open ( 69, file=trim(PATH)//trim(it_char)//'_uvel.67', status='old', form='binary' ) open ( 70, file=trim(PATH)//trim(it_char)//'_vvel.67', status='old', form='binary' ) ! open(68,file=trim(PATH)//'1_vert.69',status='old',access='direct',recl=nbyte) if(itsflux==1) then open(65,file=trim(PATH)//trim(it_char)//'_temp.70',status='old',access='direct',recl=nbyte) open(66,file=trim(PATH)//trim(it_char)//'_salt.70',status='old',access='direct',recl=nbyte) endif !itsflux !..... Elevation file ! - Header irec=0 read ( 63 ) data_format if ( data_format.ne.'DataFormat v5.0' ) then write( 20, * ) 'This code reads only v5.0: ',data_format stop endif irec=irec+48/nbyte read ( 63 ) version irec=irec+48/nbyte read ( 63 ) start_time irec=irec+48/nbyte read ( 63 ) variable_nm irec=irec+48/nbyte read ( 63 ) variable_dim irec=irec+48/nbyte write ( 20, '(a48)' ) data_format write ( 20, '(a48)' ) version write ( 20, '(a48)' ) start_time write ( 20, '(a48)' ) variable_nm write ( 20, '(a48)' ) variable_dim read ( 63 ) nrec ! # of output time steps (int) read ( 63 ) dtout ! output time step (real) read ( 63 ) nspool ! skip (int) read ( 63 ) ivs ! ivs (=1 or 2 for scalar or vector) read ( 63 ) i23d ! i23d (=2 or 3 for 2D or 3D) read ( 63 ) nvrt ! total # of vertical levels read ( 63 ) kz ! # of Z-levels read ( 63 ) h0 ! dry threshold read ( 63 ) h_s ! Z-sigma transition depth read ( 63 ) h_c ! ~thickness of finely resolved surf layer read ( 63 ) theta_b ! sigma zoom param: 1: resolve surf+bottom, 0: only surf read ( 63 ) theta_f ! sigma zoom param: ->0: uniform sigma, >>1: coarser in middle write ( 20, * ) 'nrec ',nrec write ( 20, * ) 'dtout ',dtout write ( 20, * ) 'nspool ',nspool write ( 20, * ) 'ivs ',ivs write ( 20, * ) 'i23d ',i23d write ( 20, * ) 'nvrt ',nvrt write ( 20, * ) 'kz ',kz write ( 20, * ) 'h0 ',h0 write ( 20, * ) 'h_s ',h_s write ( 20, * ) 'h_c ',h_c write ( 20, * ) 'theta_b ',theta_b write ( 20, * ) 'theta_f ',theta_f irec=irec+12 ! kz Z-levels means kz-1 Z-layers for Delwaq (kz = 1 for a no-Zlayers model) ! ztot(kz) is not read for Z-layers but equals h_s. ! is read for S-layers and equals -1.0 ! nvrt-kz+1 S-values are read (inclusive of -1.0 and 0.0) so there are nvrt-kz sigma layers ! The total is nvrt - 1 water layers in between nvrt water levels allocate ( ztot(nvrt), sigma(nvrt), cs(nvrt) ) ztot = 0.0 ; sigma = 0.0 ; cs = 0.0 ! this is not automatically done with some compilers / modes read ( 63 ) ( ztot(k), k = 1, kz-1 ) ! z-coordinates of each z-level do k = 1, nvrt-kz+1 read ( 63 ) sigma(k) ! sigma-coordinates of each S-level cs(k) = (1-theta_b)* sinh(theta_f* sigma(k) )/sinh(theta_f ) + & & theta_b *(tanh(theta_f*(sigma(k)+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5) enddo irec=irec+nvrt irec_tmp=irec !save for non-standard outputs ! Horizontal grid read ( 63 ) np ! # of nodes read ( 63 ) ne ! # of elements irec=irec+2 write ( 20, * ) ' np, ne: ',np,ne allocate ( x (np), y (np), dp (np), kbp00(np), nm (ne,4), & & idry (np), idry_e (ne), eta2 (np), dpe (ne), kbe (ne) , & & ztmp(nvrt,np), ze (nvrt,ne), zeOld(nvrt,ne), kbp2 (np), kbe2(ne) , & & area_el (ne), area_cell (np) ) allocate( aggre_point(np), verti_pnt(nvrt), layer_pnt(nvrt) ) open ( 11, file='aggrepoint.dat' ) read ( 11, * ) read ( 11, * ) aggre_point close ( 11 ) verti_pnt = 0 ; layer_pnt = 0 ; i = 0 do k = 1, nvrt read ( 10, * ) verti_pnt(nvrt-k+1) do mm = 1, verti_pnt(nvrt-k+1) i = i+1 if ( i .eq. nvrt ) then write ( 20, * ) ' Too many layers specified!' stop endif layer_pnt(nvrt-i) = k ! from Selfe layer number to Delwaq layer number enddo if ( i .eq. nvrt-1 ) then nolay = k write ( 20, * ) ' Nr. of Delwaq layers: ', nolay exit endif enddo verti_pnt(1) = 1 ; i = 1 do k = nolay, 1, -1 i = i + 1 verti_pnt(i) = verti_pnt(i-1) + verti_pnt(nvrt-k+1) ! Selfe interface administration of enddo ! the Delwaq layers (so bottom-up) if ( verti_pnt(nolay+1) .ne. nvrt ) then write ( 20, * ) ' Not enough layers specified!' stop endif write ( 238, * ) layer_pnt write ( 238, * ) verti_pnt read ( 63 ) ( x(m), y(m), dp(m), kbp00(m), m = 1, np ) irec=irec+4*np ! Make the Delwaq active(1) - inactive(0) indicator for the z-layer administration ! + second indicator: 1 = surface, 2 is mid, 3 is bed layer. allocate ( attribute(np,nolay) ) ; attribute = 0 do i = 1, np do k = kz, nvrt-1 ! these are the sigma layers attribute(i,layer_pnt(k)) = 21 ! middle layers enddo if ( nvrt-1 .gt. kz ) attribute(i, 1 ) = 11 ! surface layer if ( kz .eq. 1 ) attribute(i, nolay ) = 31 ! bed layer if only sigma if ( nvrt-1 .eq. 1 ) attribute(i, 1 ) = 01 ! one mixed layer if ( kz .eq. 1 ) cycle ! no z-layers present if ( dp(i) .lt. h_s ) then ! there are z-layers attribute(i, layer_pnt(kz) ) = 31 ! but here below the bed cycle endif do k = kz-1, 1, -1 ! these are the Z-layers attribute(i,layer_pnt(k)) = 21 ! middle layers if ( dp(i) .lt. -ztot(k) ) then if ( layer_pnt(k) .gt. layer_pnt(kz-1) ) then ! more z-layers if ( layer_pnt(kz-1) .eq. 1 ) & & attribute(i, 1 ) = 11 ! surface layer if only z-layers attribute(i,layer_pnt(k)) = 31 ! bed layer else ! only 1 z-layer attribute(i,layer_pnt(k)) = 31 ! this is the bed layer or .... if ( layer_pnt(kz-1) .eq. 1 ) & & attribute(i, 1 ) = 01 ! mixed layer if only 1 z-layer endif exit endif enddo enddo do m = 1, ne read ( 63 ) i34, ( nm(m,mm), mm = 1, i34 ) ! node # (connectivity table) irec = irec + i34 + 1 enddo !m irec0 = irec ! read file 14 hgrid.gr3 to get open boundary nodes open ( 14, file=trim(PATH)//'../hgrid.gr3', status='old' ) read ( 14, * ) ; read(14,*) ne2, np2 if ( ne /= ne2 ) then write ( 20, * ) np , ne write ( 20, * ) np2, ne2 stop 'mismatch in hgrid nodes' endif do i = 1, np2 ; read ( 14, * ) ; enddo ! skip node coordinates do i = 1, ne2 ; read ( 14, * ) ; enddo ! skip connectivity table read ( 14, * ) nopenbnds read ( 14, * ) nopenbndnodes write( 20, * ) nopenbnds, nopenbndnodes allocate ( openbnds(2,nopenbnds) ) allocate ( openbndnodes(nopenbndnodes), isopenbnd(np), stat=istat ) if ( istat /=0 ) stop 'Failed to allocate open bnd arrays' isopenbnd = 0 l = 1 do i = 1, nopenbnds read ( 14, * ) nopn do j = 1, nopn read ( 14, * ) openbndnodes(l) isopenbnd( openbndnodes(l) ) = 1 l = l + 1 if ( j .eq. 1 ) openbnds(1,i) = j if ( j .eq. nopn ) openbnds(2,i) = j enddo enddo close(14) ! dihv.66 read ( 67 ) ( idummy , k = 1, irec_tmp ) read ( 67 ) ne2, ne_e read ( 67 ) ( idummy , k = 1, 4*ne + 4*ne_e ) if ( ne /= ne2 ) stop 'mismatch in dihv.66 nodes' allocate ( dihv(2,ne), stat=istat ) if ( istat /= 0 ) stop 'Failed to allocate hvel arrays' irec0_dihv = irec_tmp + 2 + 4*ne + 4*ne_e ! hvel.67 read ( 69 ) ( idummy , k = 1, irec_tmp ) read ( 69 ) ns, ns_e write ( 20, * ) ' ns, nse:',ns,ns_e allocate ( dps (ns) , kbs (ns) , kbs2(ns) , zs(nvrt,ns), idry_s(ns), & & suv21(nvrt,ns), suv22(nvrt,ns), stat=istat ) if ( istat /= 0 ) stop 'Failed to allocate hvel arrays' allocate ( sideuv(nvrt,2), centeru(nvrt) , centerv(nvrt), stat=istat ) irec0_u = irec_tmp + 2 read ( 69 ) ( idummy, idummy, dps(m), kbs(m), m = 1, ns ) read ( 69 ) ( idummy, m = 1, 4*ns_e ) irec0_u=irec0_u+4*ns irec0_u=irec0_u+4*ns_e ! allocate fluxes arrays allocate(volume(np,nvrt),volumeOld(np,nvrt),volumeError(np,nvrt),cum_influx(np,nvrt),stat=istat) ! allocate(volume(nvrt,ne),stat=istat) if(istat/=0) stop 'Failed to allocate volumes' allocate(flx_adv(2,nvrt,ns),stat=istat) if(istat/=0) stop 'Transport: fail to allocate' allocate(hfluxes(ns,nvrt),vfluxes(np,nvrt),nodal_flux(3,nvrt),stat=istat) allocate ( wfluxes(ns,nolay) ) if(istat/=0) stop 'Transport: fail to allocate (2)' ! Compute geometry call compute_nside(np,ne,nm,ns2) if(ns/=ns2) stop 'mismatch in sides' write( 20, * ) 'done compute_nside' allocate(ic3(ne,3),js(ne,3),is(ns2,2),isidenode(ns2,2),xcj(ns2),ycj(ns2),stat=istat) if(istat/=0) stop 'Failed to allocate side arrays' ! ns: # of sides ! ic3(ne,3): 3 neighboring elements of an element ! js(ne,3): 3 sides of an element ! is(ns,2): 2 adjacent elements of a side ! isidenode(ns,2): 2 end nodes of a side ! xcj(ns),ycj(ns): x,y of each side center call selfe_geometry(np,ne,ns2,x,y,nm,ic3,js,is,isidenode,xcj,ycj) write ( 20, * ) 'done selfe_geometry' write ( 20, * ) 'last element',(nm(ne,j),j=1,3) ! Allocate and compute side frame tensor for ics=1 or 2 ! sframe(i,j,isd): where j is the axis id, i is the component id, isd is the local side id ! For ics=1, only sframe(1:2,1:2,isd) are used allocate(sframe(3,3,ns),stat=istat); if(istat/=0) stop 'AQUIRE_HGRID: sframe allocation failure' allocate(distj(ns),sne(nvrt,3),area_e(ne),stat=istat); if(istat/=0) stop 'AQUIRE_HGRID: distj allocation failure' ! Build side data for augmented subdomain do i=1,ns ! for all sides !Internal wet sides ie1=is(i,1); ie2=is(i,2) ! the elements sharing the side n1=isidenode(i,1); n2=isidenode(i,2) ! side end nodes distj(i)=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2) ! side length ! znd zero for cartesian ! distj(jsj)=sqrt((x(n2)-x(n1))**2+(y(n2)-y(n1))**2+(znd(n2)-znd(n1))**2) enddo thetan=-1.e10 !max. deviation between ze and zs axes sframe=0 !for ics=1 do j=1,ns n1=isidenode(j,1) n2=isidenode(j,2) if(ics==1) then ! cartesian thetan=atan2(x(n1)-x(n2),y(n2)-y(n1)) sframe(1,1,j)=cos(thetan) !old snx sframe(2,1,j)=sin(thetan) !old sny sframe(1,2,j)=-sframe(2,1,j) sframe(2,2,j)=sframe(1,1,j) else !lat/lon ! omitted endif !ics enddo !j=1,ns !... Element area area_el=0 area_cell=0 do i=1,ne xel(1:3)=x(nm(i,1:3)) yel(1:3)=y(nm(i,1:3)) area_el(i)=((xel(1)-xel(3))*(yel(2)-yel(3))-(xel(2)-xel(3))*(yel(1)-yel(3)))/2d0 do j=1,3 ! node k = nm(i,j) area_cell(k) = area_cell(k) + area_el(i)/3 enddo enddo !i=1,nea allocate ( uvelos(nolay,np), vvelos(nolay,np) ) !*********delwaq********** allocate ( dist_dw(ns), area_dw(ns,nvrt-1), bfluxes(nopenbndnodes,nvrt-1) ) call waq_hyd ( start_time , dtout , nrec , nolay , sigma , & & out_prefix , ns , np , version , nopenbndnodes ) call wrselafin ( out_prefix , version , start_time , ne , np , & & nm , x , y , dp , nopenbnds , & & openbnds , nopenbndnodes , openbndnodes ) call wrwaqfil ( out_prefix , np , ns , ne , nolay , & & dp , area_cell , isidenode , distj , is , & & x , y , xcj , ycj , nm , & & dist_dw , nopenbndnodes , openbndnodes ) write ( 545 ) attribute if ( timon ) call timstop( ithndl(2) ) !*********delwaq******* *** !... Time iteration !... if ( timon ) call timstrt ( "time loop", ithndl(11)) do istack = first_stack , last_stack if ( timon ) call timstrt ( "stack file opening", ithndl(3) ) write ( it_char, '(i12)' ) istack it_char = adjustl(it_char) if ( istack .ne. first_stack ) then open ( 63, file=trim(PATH)//trim(it_char)//'_elev.61', status='old', form='binary' ) read ( 63 ) ( idummy, k=1,irec0 ) open ( 67, file=trim(PATH)//trim(it_char)//'_dihv.66', status='old', form='binary' ) read ( 67 ) ( idummy, k=1,irec0_dihv ) open ( 69, file=trim(PATH)//trim(it_char)//'_uvel.67', status='old', form='binary' ) read ( 69 ) ( idummy, k=1,irec0_u ) endif open ( 70, file=trim(PATH)//trim(it_char)//'_vvel.67', status='old', form='binary' ) read ( 70 ) ( idummy, k=1,irec0_u ) if ( timon ) call timstop( ithndl(3) ) do iblock=1,nrec if ( timon ) call timstrt ( "read 2D info", ithndl(4) ) read ( 63 ) time, it if (iblock == 1 .or. iblock == nrec) then write ( 20, * ) 'i=',trim(it_char),' iblock=',iblock,' time (days) =',time/86400 endif read ( 63 ) eta2 read ( 63 ) ( idummy, k=1,np ) read ( 69 ) ( idummy, k=1,ns+2 ) read ( 70 ) ( idummy, k=1,ns+2 ) read ( 67 ) ( idummy, k=1,ne+2 ) ! Compute z coordinates ztmp=-1.e22 do i=1,np if(eta2(i)+dp(i)=ztot(k).and.-dp(i)1.e18) then write ( 20, * )'Dry side (1):',i,k,n1,n2,ztmp(k,n1),ztmp(k,n2) stop endif enddo !k enddo !i=1,ns zeOld=ze do i=1,ne if(idry_e(i)==1) cycle !Compute zcor for wet elem. n1=nm(i,1); n2=nm(i,2); n3=nm(i,3) do k=kbe2(i),nvrt ze(k,i)=(ztmp(k,n1)+ztmp(k,n2)+ztmp(k,n3))/3 if(abs(ze(k,i))>1.e18) then write ( 20, * )'Dry elem:',i,ztmp(k,n1),ztmp(k,n2),ztmp(k,n3) stop endif enddo !k enddo !i=1,ne if ( timon ) call timstop ( ithndl(5) ) if ( timon ) call timstrt ( "compute velocities", ithndl(12) ) uvelos = 0.0 ; vvelos = 0.0 do isd = 1, ns ! loop edges if ( idry_s(isd) == 1 ) cycle n1 = isidenode(isd,1) n2 = isidenode(isd,2) w1 = area_el(is(ns,1)) + area_el(is(ns,2)) ! areas of elements at both sides of the edge n3 = 0 do k = kbs2(isd), nvrt-1 ! This is the Selfe nodal layer number (separation plane) if ( n3 .ne. layer_pnt(k) ) then n3 = layer_pnt(k) ! This is the Delwaq layer number, condensed, reversed order w2 = w1 / (zs(verti_pnt(n3+1),isd)-zs(verti_pnt(n3),isd)) endif w3 = w2 * (zs(k+1,isd)-zs(k,isd)) * ( suv21(k+1,isd) + suv21(k,isd) ) w4 = w2 * (zs(k+1,isd)-zs(k,isd)) * ( suv22(k+1,isd) + suv22(k,isd) ) uvelos(n3,n1) = uvelos(n3,n1) + w3 uvelos(n3,n2) = uvelos(n3,n2) + w3 vvelos(n3,n1) = vvelos(n3,n1) + w4 vvelos(n3,n2) = vvelos(n3,n2) + w4 enddo enddo do i=1,np if ( idry(i) == 1 ) cycle w1 = 1.0 / ( area_cell(i) * 2.0 ) uvelos(:,i) = uvelos(:,i) * w1 vvelos(:,i) = vvelos(:,i) * w1 enddo if ( timon ) call timstop ( ithndl(12) ) if ( timon ) call timstrt ( "compute h-fluxes", ithndl(6) ) hfluxes=0 vfluxes=0 ! loop over elems do ie=1,ne if ( idry_e(ie) == 1 ) cycle el_depth = ze(nvrt,ie) - ze(kbe2(ie),ie) ! compute hvel at triangle center centeru = 0 ; centerv = 0 ! vertical values uint = 0 ; vint = 0 ! vertical integral do j = 1, 3 ! sides isd = js(ie,j) do k = kbe2(ie), nvrt-1 centeru(k) = centeru(k) + suv21(k+1,isd) + suv21(k,isd) centerv(k) = centerv(k) + suv22(k+1,isd) + suv22(k,isd) enddo enddo do k = kbe2(ie), nvrt-1 uint = uint + centeru(k)*(ze(k+1,ie)-ze(k,ie)) vint = vint + centerv(k)*(ze(k+1,ie)-ze(k,ie)) enddo ! correct 3D velocity field fact = dihv(1,ie)/uint if( abs(fact) < 1.6 ) then ! correct by multiplying the profile do k=kbe2(ie),nvrt-1 centeru(k) = centeru(k)*fact enddo else ! correct by substracting the error uint = ( uint - dihv(1,ie) ) / el_depth ! correction do k=kbe2(ie),nvrt-1 centeru(k) = centeru(k) - uint enddo endif fact = dihv(2,ie)/vint if( abs(fact) < 1.6 ) then ! correct by multiplying the profile do k=kbe2(ie),nvrt-1 centerv(k) = centerv(k)*fact enddo else ! correct by substracting the error vint = ( vint - dihv(2,ie) ) / el_depth ! correction do k=kbe2(ie),nvrt-1 centerv(k) = centerv(k) - vint enddo endif ! compute the 3 nodal fluxes ( flux towards node nm(ie,i) ) nodal_flux = 0.0 do i=1,3 ! nodes n1 = nm(ie,i) ! other nodes n2 = nm(ie,modulo(i+0,3)+1) n3 = nm(ie,modulo(i+1,3)+1) ! opposite edge mid point midx = 0.5*(x(n2)+x(n3)) midy = 0.5*(y(n2)+y(n3)) ! vector normal to opposite edge (towards node n1) snx = -( y(n2)-y(n3) ) ! normal sny = +( x(n2)-x(n3) ) if ( (x(n1)-midx)*snx+(y(n1)-midy)*sny < 0 ) then snx=-snx sny=-sny endif ! flux towards node n1 from element ie do k=kbe2(ie),nvrt-1 nodal_flux(i,k) = 0.5*( centeru(k)*snx + centerv(k)*sny )*(ze(k+1,ie)-ze(k,ie)) enddo enddo ! assemble fluxes on sides do j=1,3 ! sides isd = js(ie,j) ! end nodes n1 = isidenode(isd,1) n2 = isidenode(isd,2) do i=1,3 ! nodes if(nm(ie,i)==n2) then ! end node is target node do k=kbe2(ie),nvrt-1 hfluxes(isd,k) = hfluxes(isd,k) + nodal_flux(i,k)/3 enddo elseif (nm(ie,i)==n1) then ! start node is target node do k=kbe2(ie),nvrt-1 hfluxes(isd,k) = hfluxes(isd,k) - nodal_flux(i,k)/3 enddo endif enddo enddo enddo if ( timon ) call timstop ( ithndl(6) ) if ( timon ) call timstrt ( "compute cum-h-fluxes", ithndl(7) ) ! compute volume change based on horizontal fluxes volumeOld = volume volume = 0.0 wfluxes = 0.0 cum_influx = 0.0 do isd=1,ns if ( idry_s(isd) == 1 ) cycle ! choose the end node n1 = isidenode(isd,1) n2 = isidenode(isd,2) do k = kbs2(isd), nvrt-1 n3 = layer_pnt(k) cum_influx(n1,n3) = cum_influx(n1,n3) - hfluxes(isd,k) cum_influx(n2,n3) = cum_influx(n2,n3) + hfluxes(isd,k) wfluxes (isd,n3) = wfluxes (isd,n3) + hfluxes(isd,k) enddo enddo if ( timon ) call timstop ( ithndl(7) ) ! compute new volume from z coords ! compute vertical flux from geometry and h. fluxes if ( timon ) call timstrt ( "compute volums & v-fluxes", ithndl(8) ) do i=1,np if(idry(i)==1) cycle ! bottom to top if ( isopenbnd(i) .eq. 1 ) cycle do k = nolay-layer_pnt(kbp2(i))+1, nolay n1 = max(kbp2(i)+1,verti_pnt(k+1)) n2 = max(kbp2(i) ,verti_pnt(k )) volume(i,k) = area_cell(i)*(ztmp(n1,i) - ztmp(n2,i)) ! recursion on vertical fluxes, vfluxes(i,kbp2(i)) = 0 vfluxes(i,k+1) = -(volume(i,k)-volumeOld(i,k))/dtout +cum_influx(i,nolay-k+1) + vfluxes(i,k) enddo enddo if ( timon ) call timstop ( ithndl(8) ) if ( timon ) call timstrt ( "compute boundary-fluxes", ithndl(9) ) bfluxes = 0.0d+00 do i = 1, nopenbndnodes n3 = openbndnodes(i) if ( idry(n3) == 1 ) cycle ! exclude dry bnd nodes ? (should never happen...) do k = nolay-layer_pnt(kbp2(n3))+1, nolay n1 = max(kbp2(n3)+1,verti_pnt(k+1)) n2 = max(kbp2(n3) ,verti_pnt(k )) volume(n3,k) = area_cell(n3)*(ztmp(n1,n3) - ztmp(n2,n3)) bfluxes(i,k) = (volume(n3,k)-volumeOld(n3,k))/dtout - cum_influx(n3,nolay-k+1) enddo enddo if ( timon ) call timstop ( ithndl(9) ) if ( timon ) call timstrt ( "write Delwaq files", ithndl(10) ) write ( 540 ) int(time), ( volume ( : , k ) , k = nolay, 1, -1 ) if ( save_time .ne. -999 ) then write ( 541 ) save_time do k = 1, nolay write( 541 ) bfluxes(:,nolay-k+1) ! open boundary fluxes write( 541 ) wfluxes(:,k) ! normal horizontal fluxes enddo write ( 541 ) ( -1.0*sngl(vfluxes(:,k)) , k = nolay, 2, -1 ) ! vertical fluxes endif do isd = 1, ns if ( idry_s(isd) == 1 ) cycle ! for dry sides kbs2 == -1 which fails n1 = isidenode(isd,1) n2 = isidenode(isd,2) if ( n1 .gt. 0 .and. n2 .gt. 0 ) then do k = nolay-layer_pnt(kbs2(isd))+1, nolay n3 = verti_pnt(k+1) ! max(kbs2(isd) ,verti_pnt(k+1)) n4 = verti_pnt(k ) ! max(kbs2(isd) ,verti_pnt(k )) area_dw(isd,k) = dist_dw(isd)*(ztmp(n3,n1)-ztmp(n4,n1)+ztmp(n3,n2)-ztmp(n4,n2))*0.5 enddo elseif ( n1 .gt. 0 ) then do k = nolay-layer_pnt(kbs2(isd))+1, nolay n3 = verti_pnt(k+1) ! max(kbs2(isd) ,verti_pnt(k+1)) n4 = verti_pnt(k ) ! max(kbs2(isd) ,verti_pnt(k )) area_dw(isd,k) = dist_dw(isd)*(ztmp(n3,n1)-ztmp(n4,n1)) enddo elseif ( n2 .gt. 0 ) then do k = nolay-layer_pnt(kbs2(isd))+1, nolay n3 = verti_pnt(k+1) ! max(kbs2(isd) ,verti_pnt(k+1)) n4 = verti_pnt(k ) ! max(kbs2(isd) ,verti_pnt(k )) area_dw(isd,k) = dist_dw(isd)*(ztmp(n3,n2)-ztmp(n4,n2)) enddo endif do k = nolay-layer_pnt(kbs2(isd))+1, nolay area_dw(isd,k) = max( 1.0, area_dw(isd,k) ) enddo enddo if ( save_time .ne. -999 ) then write ( 542 ) save_time do k = nolay, 1, -1 write ( 542 ) ( 100.0 , i = 1, nopenbndnodes ) ! open boundary areas write ( 542 ) area_dw(:,k) ! normal horizontal crosssectional areas enddo write ( 542 ) ( area_cell(:) , k=nolay, 2, -1 ) ! vertical areas crosssectional areas write ( 543 ) save_time, (uvelos(k,:),k=1,nolay) write ( 544 ) save_time, (vvelos(k,:),k=1,nolay) endif save_time = int(time) if ( timon ) call timstop ( ithndl(10) ) enddo !iblock=1,nrec enddo !istack=1,nstacks if ( timon ) call timstop( ithndl(11)) if ( timon ) call timstrt ( "write Delwaq files", ithndl(10) ) bfluxes = 0.0 wfluxes = 0.0 vfluxes = 0.0 uvelos = 0.0 vvelos = 0.0 write ( 541 ) int(time), ( bfluxes(:,k), wfluxes(:,k), k=1,nolay ), (sngl(vfluxes(:,k)),k=nolay,2,-1) write ( 542 ) int(time), ( bfluxes(:,k), wfluxes(:,k), k=1,nolay ), (sngl(vfluxes(:,k)),k=nolay,2,-1) write ( 543 ) int(time), uvelos write ( 544 ) int(time), vvelos if ( timon ) call timstop ( ithndl(10) ) if ( timon ) then call timstop ( ithndl(1) ) call timdump ( 'selfe-delwaq-timers.out' ) endif close ( 63 ) close ( 69 ) close ( 70 ) stop end