! 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 nc outputs from scribe I/O versions for multiple files at all nodes ! Works for mixed tri/quad outputs on NODE/ELEMENT based vars. ! Inputs: screen; nc file and out2d*.nc; vgrid.in (in this dir or ../) ! Outputs: extract.out (ascii); optional: max/min/avg 2D output (*_[max|min|avg].gr3 and *_violators.bp) !**************************************************************************************** ! ifort -O2 -assume byterecl -o read_output10_allnodes.exe ../UtilLib/extract_mod2.f90 ../UtilLib/schism_geometry.f90 ../UtilLib/compute_zcor.f90 read_output10_allnodes.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff program read_out use netcdf use extract_mod2 use schism_geometry_mod use compute_zcor ! parameter(nbyte=4) character(len=30) :: file63,varname,outname(3),file62 character(len=12) :: it_char ! character*48 data_format logical::lexist integer :: nx(4,4,3),dimids(100),idims(100) allocatable :: sigma(:),cs(:),ztot(:) allocatable :: outvar(:,:),out(:,:,:),icum(:,:,:),eta2(:),ztmp(:,:) allocatable :: xcj(:),ycj(:),dp(:),dps(:),kbs(:),xctr(:),yctr(:),dpe(:),kbe(:) allocatable :: zs(:,:),ze(:,:),idry(:),outs(:,:,:),oute(:,:,:),rstat2d(:,:),rviolator(:) allocatable :: ic3(:,:),isdel(:,:),isidenode(:,:),kbp(:),sigma_lcl(:,:),kbp00(:) integer, allocatable :: elside(:,:),idry_e(:),nne(:),indel(:,:),iviolator(:) integer,allocatable :: i34(:),elnode(:,:) integer :: char_len,start_2d(2),start_3d(3),start_4d(4), & &count_2d(2),count_3d(3),count_4d(4) real*8,allocatable :: timeout(:),xnd(:),ynd(:) integer :: icomp_stats(3) print*, 'Input NODE-based variable name to read from nc (e.g. elevation):' read(*,'(a30)')varname !!' print*, '<<<<<var name: ',varname print*, 'Is the var node-based or ele-based? 0: node based; 1: element based' read(*,*)inode_ele !!' print*, '<<<<<inode_ele: ',inode_ele varname=adjustl(varname); len_var=len_trim(varname) print*, 'Input start and end stack # to read:' read(*,*) iday1,iday2 print*, '<<<<<start and end stack #: ',iday1,iday2 print*, 'Do you want to compute stats for 2D var? (0/1:min; 0/1:max; 0/1:time avg)' !!' read(*,*) icomp_stats(1),icomp_stats(2),icomp_stats(3) print*, '<<<<<icomp_stats: ',icomp_stats outname(1)='min';outname(2)='max';outname(3)='avg' print*, 'Input a threshold: values below the threshold will be output as a bp file.' !!' read(*,*) thres print*, '<<<<<threshold: ',thres print*, 'How do you want to initialize the variable values for min/max?' !!' print*, '0: do nothing; 1: intialized to -dp (useful for maxelev)' !!' read(*,*) iinitial !! print*, '<<<<<initialization flag: ',iinitial if (iinitial.ne.0 .and. iinitial.ne.1) stop 'wrong initialization flag' !!' if(varname(1:len_var).eq.'elev' .and. iinitial==1) then is_elev=1 print*, '<<<<<special treatment will be implemented for elev' else is_elev=0 endif open(65,file='extract.out') !... Get basic info from out2d*.nc !Returned vars: ne,np,ns,nrec,[xnd ynd dp](np), !elnode,i34,nvrt,h0,dtout,kbp call get_dims(iday1,np,ne,ns,nvrt,h0) allocate(xnd(np),ynd(np),dp(np),kbp(np),i34(ne),elnode(4,ne),stat=istat) if(istat/=0) stop 'alloc (1)' call readheader(iday1,np,ne,ns,kbp,i34,elnode,nrec,xnd,ynd,dp,dtout) print* print*, 'After header:',ne,np,nrec,i34(ne),elnode(1:i34(ne),ne),nvrt,h0,xnd(np),ynd(np),dp(np) !,start_time ! Compute geometry allocate(nne(np)) do k=3,4 do i=1,k do j=1,k-1 nx(k,i,j)=i+j if(nx(k,i,j)>k) nx(k,i,j)=nx(k,i,j)-k if(nx(k,i,j)<1.or.nx(k,i,j)>k) then write(*,*)'nx wrong',i,j,k,nx(k,i,j) stop endif enddo !j enddo !i enddo !k nne=0 do i=1,ne do j=1,i34(i) nd=elnode(j,i) nne(nd)=nne(nd)+1 enddo enddo mnei=maxval(nne) allocate(indel(mnei,np),stat=istat) if(istat/=0) stop 'Failed to alloc. indel' nne=0 do i=1,ne do j=1,i34(i) nd=elnode(j,i) nne(nd)=nne(nd)+1 if(nne(nd)>mnei) then write(*,*)'Too many neighbors',nd stop endif indel(nne(nd),nd)=i enddo enddo !i call compute_nside(np,ne,i34,elnode(1:4,1:ne),ns2) if(ns2/=ns) then write(*,*)'Mismatch in side:',ns2,ns stop endif allocate(ic3(4,ne),elside(4,ne),isdel(2,ns2),isidenode(2,ns2),xcj(ns2),ycj(ns2),stat=istat) if(istat/=0) stop 'Allocation error: side(0)' call schism_geometry_single(np,ne,ns2,real(xnd),real(ynd),i34,elnode(1:4,1:ne),ic3(1:4,1:ne), & &elside(1:4,1:ne),isdel,isidenode,xcj,ycj) !For dimensioning purpose if(np>ns.or.ne>ns) stop 'ns is not largest' if (inode_ele==0) then last_dim=np elseif (inode_ele==1) then last_dim=ne endif print*, '<<<<<last_dim: ',last_dim allocate(idry_e(ne),rstat2d(3,ne),ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np), & &outvar(nvrt,last_dim),eta2(np),idry(np),ztmp(nvrt,np),timeout(nrec), & &rviolator(ne),iviolator(ne),kbp00(np)) outvar=-huge(1.0) !test mem ! Read vgrid.in call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp) ! Calculate kbp00 if(ivcor==1) then kbp00=kbp else do i=1,np call zcor_SZ_single(dp(i),1.e8,h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot,sigma, & &ztmp(:,i),idry2,kbp00(i)) enddo !i endif !ivcor !... Time iteration !... outvar=-99 !init. ztmp=-99 if(iinitial==1) then print*,'setting initial elev = -dp for all (dry and wet) nodes' rstat2d(1,1:np)=-dp !for elev, min is -h rstat2d(2,1:np)=-dp !for elev, min is -h rstat2d(3,1:np)=-dp !for elev, min is -h else if (icomp_stats(1)==1) then !min rstat2d(1,:)=huge(1.0) elseif (icomp_stats(2)==1) then !max rstat2d(2,:)=-huge(1.0) elseif (icomp_stats(3)==1) then !avg rstat2d(3,:)=0.0 endif endif do iday=iday1,iday2 !1,ndays !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ write(it_char,'(i12)')iday it_char=adjustl(it_char) leng=len_trim(it_char) file62='out2d_'//it_char(1:leng)//'.nc' iret=nf90_open(trim(adjustl(file62)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid4) !time is double iret=nf90_inq_varid(ncid4,'time',itime_id) iret=nf90_get_var(ncid4,itime_id,timeout,(/1/),(/nrec/)) print*, 'time=',timeout !,trim(adjustl(file63)) !Find nc file file63=varname(1:len_var)//'_'//it_char(1:leng)//'.nc' inquire(file=file63,exist=lexist) if(lexist) then i23d=2 !3D var else i23d=1 !2D file63=file62 endif iret=nf90_open(trim(adjustl(file63)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid) iret=nf90_inq_varid(ncid,varname(1:len_var),ivarid1) if(iret/=nf90_NoErr) stop 'Var not found' iret=nf90_Inquire_Variable(ncid,ivarid1,ndims=ndims,dimids=dimids) if(ndims>100) stop 'increase dimension of dimids & idims' do i=1,ndims iret=nf90_Inquire_Dimension(ncid,dimids(i),len=idims(i)) enddo !i npes=idims(ndims-1) !np|ne|ns if(npes/=np.and.npes/=ne) stop 'can only handle node- or elem-based' !' if(idims(ndims)/=nrec) stop 'last dim is not time' ! iret=nf90_get_att(ncid,ivarid1,'i23d',i23d) ! if(i23d<=0.or.i23d>6) stop 'wrong i23d' ! if(i23d>3.and.ics==2) stop 'ics=2 with elem-based var' ! iret=nf90_get_att(ncid,ivarid1,'ivs',ivs) ! !print*, 'i23d:',i23d,ivs,idims(1:ndims) ! do irec=1,nrec !Get elev iret=nf90_inq_varid(ncid4,'elevation',itmp) start_2d(1)=1; start_2d(2)=irec count_2d(1)=np; count_2d(2)=1 iret=nf90_get_var(ncid4,itmp,eta2,start_2d,count_2d) 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) 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) endif !Available now: outvar(nvrt,np|ne), eta2(np) if(i23d==1) then !2D ! Output: time, 2D variable at all nodes if(maxval(icomp_stats(1:3))/=0) & &write(65,'(e14.6,1000000(1x,e14.4))')timeout(irec)/86400,(outvar(1,i),i=1,np) !Compute stats (magnitude for vectors) if(sum(icomp_stats(:))/=0) then if(is_elev==1 .and. icomp_stats(2)==1) then !maxelev !Mark wet elem do i=1,ne if(minval(outvar(1,elnode(1:i34(i),i))+dp(elnode(1:i34(i),i)))>0) then idry_e(i)=0 else idry_e(i)=1 endif enddo !i do i=1,ne if(idry_e(i)==0) then !wet ifl=1 !isolated wet loop2: do j=1,i34(i) nd=elnode(j,i) do m=1,nne(nd) ie=indel(m,nd) if(i/=ne.and.idry_e(ie)==0) then ifl=0 exit loop2 endif enddo !m enddo loop2 !j if(ifl==0) then ! not isolated wet do j=1,i34(i) nd=elnode(j,i) tmp=outvar(1,nd) if(tmp+dp(nd)>0) rstat2d(2,nd)=max(rstat2d(2,nd),tmp) enddo endif endif !idry_e enddo !i=1,ne else !other variables, or min_ele or avg_ele do i=1,last_dim tmp=outvar(1,i) if (icomp_stats(1)==1) then !min rstat2d(1,i)=min(rstat2d(1,i),tmp) endif if(icomp_stats(2)==1) then !max rstat2d(2,i)=max(rstat2d(2,i),tmp) endif if(icomp_stats(3)==1) then !min rstat2d(3,i)=tmp+rstat2d(3,i) endif enddo !i endif !is_elev endif ! icomp_stats else !if(i23d==3) then !3D !Compute z coordinates do i=1,last_dim 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) endif do k=max0(1,kbp00(i)),nvrt ! Output: time, node #, level #, z-coordinate, 3D variable if(idry(i)==1) then write(65,*)timeout(irec)/86400,i,k,-1.e6,outvar(k,i) else write(65,*)timeout(irec)/86400,i,k,ztmp(k,i),outvar(k,i) endif enddo !k !Debug !write(65,*)x(i),y(i),out(i,1,1:ivs) enddo !i endif !i23d enddo !irec !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ enddo !iday nrec_total=nrec*(iday2-iday1+1) if(icomp_stats(3)==1) then rstat2d(3,:)=rstat2d(3,:)/real(nrec_total) endif print*, '<<<<<total rec: ', nrec_total do k=1,3 if(icomp_stats(k)==0) cycle if (inode_ele==0) then !node-based, write gr3 open(12,file=trim(adjustl(varname))//'_'//trim(adjustl(outname(k)))//'.gr3',status='replace') write(12,*)iday1,iday2 write(12,*)ne,np do i=1,np write(12,*)i,xnd(i),ynd(i),rstat2d(k,i) enddo !i do i=1,ne write(12,*)i,i34(i),elnode(1:i34(i),i) enddo !i close(12) elseif(inode_ele==1) then !ele based, write prop open(12,file=trim(adjustl(varname))//'_'//trim(adjustl(outname(k)))//'.prop',status='replace') do i=1,ne write(12,*)i,rstat2d(k,i) enddo endif close(12) enddo !k=1,3 !record violators do k=1,3 nrec=0; rviolator=0.0; iviolator=0 if (icomp_stats(k)==0) cycle do i=1,last_dim if (rstat2d(k,i)<thres) then nrec=nrec+1 rviolator(nrec)=rstat2d(k,i) !val iviolator(nrec)=i !id endif enddo open(13,file=trim(adjustl(varname))//'_'//trim(adjustl(outname(k)))//'_violators.bp',status='replace') write(13,*) 'violators <', thres write(13,*) nrec if (inode_ele==0) then !node-based do i=1,nrec ip=iviolator(i) write(13,*) i,xnd(ip),ynd(ip),rviolator(i),ip enddo elseif (inode_ele==1) then !ele based do i=1,nrec ie=iviolator(i) write(13,*) i,sum(xnd(elnode(1:i34(ie),ie)))/i34(ie),& & sum(ynd(elnode(1:i34(ie),ie)))/i34(ie), rviolator(i),ie enddo endif close(13) enddo !k=1,3 print*, 'Finished!' stop end