! 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 in (x,y,time) from station.xyt (time in sec); e.g., casts; ! for 3D variables (surface values for 2D variables) DEFINED AT NODES or ELEM. ! Interpolation in time, and ! add extra times before and after to examine phase errors. ! Works for mixed tri/quad outputs from scribe I/O versions. ! Inputs: (1) nc files and out2d*.nc; ! (2) station.xyt (bp format): make sure all times (in sec) are after 1st record (to ensure interpolation in time); ! pad extra days before and after if necessary. ! (3) read_output_xyt.in: ! 1st line: variable name (e.g. elev)\n ! 2nd line: invalid value (for out of domain, dry etc)\n ! 3rd line: window (hours) b4 and after the cast, stride (hrs) - used to ! examine the phase error. If window=0, each cast is repeated twice ! there are no other extra casts. \n ! 4th line: inode_elem (1: node based; 2: elem based) ! (4) vgrid.in (in this dir or ../) ! Outputs: fort.18; fort.11 (fatal errors); fort.12: nonfatal errors. ! The total # of 'virtual' casts for each actual cast is 2*window/stride+2 ! ! ifort -mcmodel=medium -assume byterecl -CB -O2 -o read_output10_xyt.exe ../UtilLib/extract_mod2.f90 ../UtilLib/compute_zcor.f90 ../UtilLib/pt_in_poly_test.f90 read_output10_xyt.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 compute_zcor use pt_in_poly_test character(len=30) :: file63,varname,file62 character(len=12) :: it_char logical :: lexist integer,allocatable :: kbp(:),kbp00(:),node3(:,:),iep(:),iday(:,:),irecord(:,:),irank_read(:), & &i34(:),elnode(:,:) real,allocatable :: sigma(:),cs(:),ztot(:),out2(:,:),eta2(:),arco(:,:),ztmp(:),x00(:),y00(:), & &out4(:),t00(:),times(:,:),sigma_lcl(:,:),ztmp2(:,:),outvar(:,:),dp(:) real*8,allocatable :: timeout(:),xnd(:),ynd(:) integer :: nodel(3),dimids(100),idims(100),char_len,start_2d(2),start_3d(3),start_4d(4), & &count_2d(2),count_3d(3),count_4d(4) open(10,file='read_output_xyt.in',status='old') read(10,'(a30)')varname ! Junk used for 3D variables: below bottom; dry spot; no parents read(10,*)rjunk read(10,*)window,wstride !in hours read(10,*)inode_elem close(10) if(icomb/=0.and.icomb/=1) stop 'Unknown icomb' varname=adjustl(varname); len_var=len_trim(varname) if(wstride==0) stop 'wstride=0' nextra=2*window/wstride+1 !extra casts in addition to each cast (for phase error) !... 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(1,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(1,np,ne,ns,kbp,i34,elnode,nrec,xnd,ynd,dp,dtout) print*, 'After header:',ne,np,ns,nvrt,nrec,i34(ne), & &elnode(1:i34(ne),ne),h0,xnd(np),ynd(np),dp(np),dtout !,start_time ! Read in station.xyt open(10,file='station.xyt',status='old') read(10,*) read(10,*) nxy nxy2=nxy*(1+nextra) allocate(x00(nxy2),y00(nxy2),t00(nxy2),eta2(np),stat=istat) if(istat/=0) stop 'Falied to allocate (1)' do i=1,nxy read(10,*)k,xtmp,ytmp,ttmp ! Check if time is before first record if(ttmpnrec.or.irecord(2,i)>nrec) then write(11,*)'Record # overflow: ',i,irecord(:,i) stop endif if(t00(i)times(2,i)) then write(11,*)'Wrong time bounds:',i,t00(i),times(:,i),iday(:,i),irecord(:,i) stop endif enddo !i=1,nxy !... Time iteration !... do i=1,nxy loop1: do l=1,2 !2 times print*, 'reading stack ',iday(l,i),' for point ',i write(it_char,'(i12)')iday(l,i) 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 ii=1,ndims iret=nf90_Inquire_Dimension(ncid,dimids(ii),len=idims(ii)) enddo !ii 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' irec=irecord(l,i) !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) out2(l,:)=0 if(i23d==1) then !2D do j=1,3 !nodes nd=node3(i,j) if(inode_elem==1) then !node out2(l,1)=out2(l,1)+arco(i,j)*outvar(1,nd) else !elem out2(l,1)=outvar(1,iep(i)) endif enddo !j else !3D ! Do interpolation etal=0; dep=0; idry=0 do j=1,3 nd=node3(i,j) if(eta2(nd)+dp(nd)