! 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. ! !**************************************************************************************** ! * ! Read binary format v5.0 (hybrid S-Z models) of *hvel.67 (sidecenter vel.) and ! compute fluxes at the borders between multiple regions (excluding open/land boundaries). ! WARNING: does not work for lat/lon! ! Works with mixed grids and LSC2 ! ! Inputs: vgrid.in, *_elev.61, *hvel.67, fluxflag.prop (depth=-1,0,1,...); ! screen inputs. If itsflux/=0, *temp.70, *salt.70 (not completed) ! Outputs: fluxes.out (if the difference between region #s=1, and neith=-1) ! History: !**************************************************************************************** ! ifort -cpp -O2 -mcmodel=medium -CB -Bstatic -o compute_fluxes_ns.WW ../UtilLib/schism_geometry.f90 ../UtilLib/compute_zcor.f90 compute_fluxes_ns.f90 program read_out use compute_zcor use schism_geometry_mod parameter(nbyte=4) character*30 file63 character*12 it_char character*48 start_time,version,variable_nm,variable_dim character*48 data_format integer,allocatable :: elnode(:,:),elside(:,:),i34(:) allocatable :: sigma(:),cs(:),ztot(:),x(:),y(:),dp(:),kbp(:),kfp(:) allocatable :: icum(:,:,:),eta2(:),ztmp(:,:) allocatable :: xcj(:),ycj(:),dps(:),kbs(:),xctr(:),yctr(:),dpe(:),kbe(:) allocatable :: zs(:,:),ze(:,:),idry(:),idry_s(:),idry_e(:) allocatable :: ic3(:,:),isdel(:,:),isidenode(:,:),nwild(:),iflux_e(:) allocatable :: fluxes(:),suv2(:,:,:),tsel(:,:,:),kbp2(:),kbs2(:),kbe2(:),sigma_lcl(:,:) ! print*, 'Outputs in Cartesian (1) or lat/lon (2)?' ! read(*,*) ics ics=1 print*, 'Input start and end stack #:' read(*,*)iday1,iday2 print*, 'Do you want to compute T,S fluxes? 0: no; 1: yes' read(*,*)itsflux !... Header open(63,file='1_elev.61',status='old',access='direct',recl=nbyte) open(64,file='1_hvel.67',status='old',access='direct',recl=nbyte) if(itsflux==1) then open(65,file='1_temp.70',status='old',access='direct',recl=nbyte) open(66,file='1_salt.70',status='old',access='direct',recl=nbyte) endif irec=0 do m=1,48/nbyte read(63,rec=irec+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 do m=1,48/nbyte read(63,rec=irec+m) version(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+48/nbyte do m=1,48/nbyte read(63,rec=irec+m) start_time(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+48/nbyte do m=1,48/nbyte read(63,rec=irec+m) variable_nm(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+48/nbyte do m=1,48/nbyte read(63,rec=irec+m) variable_dim(nbyte*(m-1)+1:nbyte*m) enddo irec=irec+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 print*, 'ivs=',ivs,'; i23d=',i23d,' ;nrec= ',nrec ! Vertical grid (obsolete; will be overwritten later) 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 ! do k=1,kz-1 ! read(63,rec=irec+k) ztot(k) ! enddo ! do k=kz,nvrt ! kin=k-kz+1 ! read(63,rec=irec+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 irec_tmp=irec !save for non-standard outputs ! Horizontal grid read(63,rec=irec+1) np read(63,rec=irec+2) ne irec=irec+2 allocate(x(np),y(np),dp(np),kbp(np),kfp(np),elnode(4,ne),idry(np),idry_e(ne), & &eta2(np),xctr(ne),yctr(ne),dpe(ne),kbe(ne),ztmp(nvrt,np),ze(nvrt,ne), & &tsel(2,nvrt,ne),kbp2(np),kbe2(ne),i34(ne),stat=istat) if(istat/=0) stop 'Falied to allocate (2)' 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)kbp(m) irec=irec+4 enddo !m=1,np do m=1,ne read(63,rec=irec+1)i34(m) irec=irec+1 do mm=1,i34(m) read(63,rec=irec+1)elnode(mm,m) irec=irec+1 enddo !mm enddo !m irec0=irec ! hvel.67 read(64,rec=irec_tmp+1) ns read(64,rec=irec_tmp+2) ns_e print*, 'side grid:',ns,ns_e allocate(dps(ns),kbs(ns),kbs2(ns),zs(nvrt,ns),idry_s(ns),suv2(2,nvrt,ns),stat=istat) if(istat/=0) stop 'Allocation error: side(8)' irec0_u=irec_tmp+2 do m=1,ns read(64,rec=irec0_u+3)dps(m) read(64,rec=irec0_u+4)kbs(m) irec0_u=irec0_u+4 enddo !m ! do m=1,ns_e ! read(64,rec=irec0_u+1)i34 ! irec0_u=irec0_u+1 ! do mm=1,i34 ! read(64,rec=irec0_u+1)itmp ! irec0_u=irec0_u+1 ! if(m==ns_e) print*, 'last elem in side grid:',itmp ! enddo !mm ! enddo !m irec0_u=irec0_u+4*ns_e !pure tri for sidecenters.gr3 ! *.70 if(itsflux==1) then read(65,rec=irec_tmp+1) ne1 read(65,rec=irec_tmp+2) ne_e if(ne1/=ne) stop 'mismatch in elem.' print*, 'Elem. grid:',ne1,ne_e irec0_TS=irec_tmp+2 do m=1,ne read(65,rec=irec0_TS+3)dpe(m) read(65,rec=irec0_TS+4)kbe(m) irec0_TS=irec0_TS+4 enddo !m irec0_TS=irec0_TS+4*ne_e endif !itsflux ! Compute geometry call compute_nside(np,ne,i34,elnode,ns2) if(ns/=ns2) stop 'mismatch in sides' print*, 'done compute_nside' allocate(ic3(4,ne),elside(4,ne),isdel(2,ns),isidenode(2,ns),xcj(ns),ycj(ns),stat=istat) if(istat/=0) stop 'Allocation error: side(7)' call schism_geometry_single(np,ne,ns2,x,y,i34,elnode,ic3,elside,isdel,isidenode,xcj,ycj) print*, 'done schism_geometry' print*, 'last element',(elnode(j,ne),j=1,i34(ne)) ! Read vgrid.in allocate(ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),kbp(np),stat=istat) if(istat/=0) stop 'Falied to allocate (1)' call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp) ! Read in fluxflag.prop allocate(iflux_e(ne),stat=istat) open(32,file='fluxflag.prop',status='old') do i=1,ne read(32,*)j,tmp1 itmp=tmp1 if(itmp<-1) stop 'fluxflag.prop wrong' iflux_e(i)=itmp enddo close(32) ! Only compute if diff(region #)=1; fluxes(i) is from i to i-1 max_reg=maxval(iflux_e) allocate(fluxes(max_reg),stat=istat) if(istat/=0) stop 'Falied to allocate (6)' open(60,file='fluxes.out') !... Compute relative record # for a node and level for 3D outputs !... ! icount=0 ! do i=1,np ! do k=max0(1,kbp(i)),nvrt ! do m=1,ivs ! icount=icount+1 ! icum(i,k,m)=icount ! enddo !m ! enddo !k ! enddo !i=1,np !... Time iteration !... it_tot=0 do iday=iday1,iday2 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ write(it_char,'(i12)')iday open(63,file=it_char//'_elev.61',status='old',access='direct',recl=nbyte) open(64,file=it_char//'_hvel.67',status='old',access='direct',recl=nbyte) if(itsflux==1) then open(65,file=it_char//'_temp.70',status='old',access='direct',recl=nbyte) open(66,file=it_char//'_salt.70',status='old',access='direct',recl=nbyte) endif irec=irec0 irec_u=irec0_u irec_TS=irec0_TS 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 print*, 'time in days =',time/86400 do i=1,np read(63,rec=irec+i) eta2(i) enddo !i irec=irec+np+np !done with (63 for this step !Compute z coordinates do i=1,np if(ivcor==1) then !localized if(dp(i)+eta2(i)<=h0) then idry(i)=1 else !wet idry(i)=0 do k=kbp(i),nvrt ztmp(k,i)=(eta2(i)+dp(i))*sigma_lcl(k,i)+eta2(i) enddo !k endif !wet/dry else if(ivcor==2) then !SZ call zcor_SZ_single(dp(i),eta2(i),h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot, & &sigma,ztmp(:,i),idry(i),kbpl) if(kbpl/=kbp(i)) stop 'mismatch (2)' endif enddo !i=1,np ! idry_e, idry_s do i=1,ne idry_e(i)=maxval(idry(elnode(1:i34(i),i))) kbe2(i)=minval(kbp(elnode(1:i34(i),i))) enddo !i idry_s=1 !init. do i=1,ns kbs2(i)=minval(kbp(isidenode(1:2,i))) do j=1,2 ie=isdel(j,i) if(ie/=0.and.idry_e(max(1,ie))==0) idry_s(i)=0 enddo !j !Check if(idry_s(i)==0) then itmp=maxval(idry(isidenode(1:2,i))) if(itmp/=0) then write(*,*)'Wet side has dry nodes:',i,isidenode(1:2,i) stop endif endif enddo !i=1,ns ! hvel.67 irec_u=irec_u+ns+2 suv2=0 !for below bottom do i=1,ns do k=max0(1,kbs(i)),nvrt do m=1,2 read(64,rec=irec_u+1) suv2(m,k,i) irec_u=irec_u+1 enddo !m enddo !k !Debug !if(iday==ndays.and.it1==nrec) write(97,*)xcj(i),ycj(i),suv2(1:2,nvrt,i) enddo !i !Compute zcor for wet sides (deal with some inconsistency) do i=1,ns if(idry_s(i)==1) cycle n1=isidenode(1,i); n2=isidenode(2,i) do k=kbs2(i),nvrt k1=max(k,kbp(n1)) k2=max(k,kbp(n2)) zs(k,i)=(ztmp(k1,n1)+ztmp(k2,n2))/2 if(abs(zs(k,i))>1.e18) then write(*,*)'Dry side (1):',i,k,n1,n2,ztmp(k,n1),ztmp(k,n2) stop endif enddo !k enddo !i=1,ns ! T,S if(itsflux==1) then irec_TS=irec_TS+ne+2 tsel=-99 do i=1,ne do k=max0(1,kbe(i)),nvrt read(65,rec=irec_TS+1) tsel(1,k,i) read(66,rec=irec_TS+1) tsel(2,k,i) irec_TS=irec_TS+1 enddo !k !Debug !if(iday==ndays.and.it1==nrec) write(99,*)i,tsel(2,nvrt,i) enddo !i do i=1,ne if(idry_e(i)==1) cycle !Compute zcor for wet elem. do k=kbe2(i),nvrt ze(k,i)=0 do j=1,i34(i) nd=elnode(j,i) kl=max(k,kbp(nd)) ze(k,i)=ze(k,i)+ztmp(kl,nd)/i34(i) enddo !j if(abs(ze(k,i))>1.e18) then write(*,*)'Dry elem:',i stop endif enddo !k !Extend valid values to bottom for wet elem. !After this all wet elem. have valid values from 1:nvrt kk0=-1 do k=1,nvrt if(tsel(1,k,i)>=0) then kk0=k; exit endif enddo !k if(kk0<=0) then write(*,*)'Weird T:',i,kk0,tsel(1,:,i) stop else if(tsel(2,kk0,i)<0) then write(*,*)'Weird S:',i,kk0,tsel(2,:,i) stop endif do k=1,kk0-1 tsel(1:2,k,i)=tsel(1:2,kk0,i) enddo !k enddo !i=1,ne endif !itsflux ! Fluxes fluxes=0 do i=1,ns if(isdel(2,i)==0.or.idry_s(i)==1) cycle !Internal wet sides ie1=isdel(1,i); ie2=isdel(2,i) if(iflux_e(ie1)==-1.or.iflux_e(ie2)==-1) cycle if(iabs(iflux_e(ie1)-iflux_e(ie2))/=1) cycle n1=isidenode(1,i); n2=isidenode(2,i) distj=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2) thetan=atan2(x(n1)-x(n2),y(n2)-y(n1)) snx=cos(thetan); sny=sin(thetan) ireg_hi=max(iflux_e(ie1),iflux_e(ie2)) ireg_lo=min(iflux_e(ie1),iflux_e(ie2)) if(ireg_hi==iflux_e(ie1)) then isgn=1 !'hi' to 'lo' is same as elem. '1' to '2' else isgn=-1 endif do k=kbs2(i),nvrt-1 !Normal vel. from 1 to 2 vnn=(suv2(1,k+1,i)+suv2(1,k,i))/2*snx+(suv2(2,k+1,i)+suv2(2,k,i))/2*sny fluxes(ireg_hi)=fluxes(ireg_hi)+isgn*vnn*distj*(zs(k+1,i)-zs(k,i)) !m^3/s enddo !k enddo !i=1,ns ! Output. fluxes(i) is from region i to i-1 write(60,'(f16.6,60000(1x,e14.4))')time/86400,fluxes(:) enddo !it1=1,nrec !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ enddo !iday stop end