! 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, combined or uncombined outputs. ! Inputs: (1) nc files; ! (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: ! variable name (e.g. elev)\n ! invalid value (for out of domain, dry etc)\n ! 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 ! icomb: 0 (uncombined nc); 1 (combined nc) ! (4) vgrid.in (in this dir or ../) ! Outputs: fort.1[89]; 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_output9_xyt.exe ../UtilLib/extract_mod.f90 \ !../UtilLib/compute_zcor.f90 ../UtilLib/pt_in_poly_test.f90 read_output9_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_mod use compute_zcor use pt_in_poly_test character(len=30) :: file63,varname character(len=12) :: it_char integer,allocatable :: kbp(:),node3(:,:),iep(:),iday(:,:),irecord(:,:),irank_read(:) dimension swild(3),out3(2,2) real,allocatable :: sigma(:),cs(:),ztot(:), & &out(:,:,:),out2(:,:,:),eta2(:),arco(:,:),ztmp(:),x00(:),y00(:), & &out4(:,:),t00(:),times(:,:),sigma_lcl(:,:),ztmp2(:,:),outvar(:,:,:) real*8,allocatable :: timeout(:) integer :: nodel(3) 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,*)icomb 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) !... Header !Returned vars: ne,np,ns,nrec,[x y dp](np), !elnode,i34,nvrt,h0,dtout !If icomb=0, additonal vars: !nproc,iegl_rank,iplg,ielg,islg,np_lcl(:),ne_lcl(:),ns_lcl(:) if(icomb==0) then !uncombined call get_global_geo else call readheader(1) endif print*, 'After header:',ne,np,ns,nrec,i34(ne), & &elnode(1:i34(ne),ne),nvrt,h0,x(np),y(np),dp(np) !,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 ! write(it_char,'(i12)')iday(l,i) ! it_char=adjustl(it_char); leng=len_trim(it_char) ! file63='schout_'//it_char(1:leng)//'.nc' ! iret=nf90_open(trim(adjustl(file63)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid) ! !time is double ! iret=nf90_inq_varid(ncid,'time',itime_id) ! iret=nf90_get_var(ncid,itime_id,timeout,(/1/),(/nrec/)) ! call get_timeout(iday(l,i),nrec,timeout) if(icomb==0) then !uncombined call get_timeout(iday(l,i),nrec,timeout,icomb) else call get_timeout(iday(l,i),nrec,timeout) endif ! print*, 'time=',timeout,trim(adjustl(file63)) irec=irecord(l,i) if(icomb==0) then !uncombined do irank=0,nproc-1 if(irank_read(irank)>0) then call get_outvar(iday(l,i),irec,varname,np,last_dim,nvrt,outvar,i23d,ivs,eta2,irank) endif enddo !irank else call get_outvar(iday(l,i),irec,varname,np,last_dim,nvrt,outvar,i23d,ivs,eta2) endif ! !Get elev ! start_2d(1)=1; start_2d(2)=irec ! count_2d(1)=np; count_2d(2)=1 ! iret=nf90_get_var(ncid,ielev_id,eta2,start_2d,count_2d) ! ! 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 !i ! npes=idims(ndims-1) ! 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' ! iret=nf90_get_att(ncid,ivarid1,'ivs',ivs) ! !print*, 'i23d:',i23d,ivs,idims(1:ndims) ! if(ivs==1) then !scalar ! if(mod(i23d-1,3)==0) 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,1:npes),start_2d,count_2d) ! else !3D ! start_3d(1:2)=1; start_3d(3)=irec ! count_3d(2)=npes; count_3d(1)=nvrt; count_3d(3)=1 ! iret=nf90_get_var(ncid,ivarid1,outvar(1,:,1:npes),start_3d,count_3d) ! endif ! else !vector ! if(mod(i23d-1,3)==0) then !2D ! start_3d(1:2)=1; start_3d(3)=irec ! count_3d(2)=npes; count_3d(1)=2; count_3d(3)=1 ! iret=nf90_get_var(ncid,ivarid1,outvar(1:2,1,1:npes),start_3d,count_3d) ! else if(ndims-1==3) then !3D vector ! start_4d(1:3)=1; start_4d(4)=irec ! count_4d(3)=npes; count_4d(2)=nvrt; count_4d(1)=2; count_4d(4)=1 ! iret=nf90_get_var(ncid,ivarid1,outvar(:,:,1:npes),start_4d,count_4d) ! else ! stop 'Unknown type(2)' ! endif ! endif !ivs !Available now: outvar(2,nvrt,np|ne) !However, for uncombined nc, values in untouched ranks are junk out2(l,:,:)=0 if(mod(i23d-1,3)==0) then !2D do j=1,3 !nodes nd=node3(i,j) do m=1,ivs if(i23d<=3) then !node out2(l,1,m)=out2(l,1,m)+arco(i,j)*outvar(m,1,nd) else if (i23d<=6) then !elem out2(l,1,m)=outvar(m,1,iep(i)) endif enddo !m 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)