! 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 (combined or uncombined) nc outputs and compute average field ! at a particular Z-level (use above surface/below bottom to get surface/bottom). ! Skip dry times for 3D variables. ! Works for mixed quad/tri on NODE based variables only. ! Input: schout*.nc (combined or uncombined); local_to_global*; vgrid.in; screen ! Output: average.out (gredit for scalar or xmgr5 format for vector) ! ! ifort -O2 -mcmodel=medium -assume byterecl -CB -o compute_average3.exe ../UtilLib/extract_mod.f90 ../UtilLib/compute_zcor.f90 compute_average3.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 parameter(nbyte=4) character(len=30) :: file63,varname character(len=12) :: it_char character(len=48) :: data_format ! integer,allocatable :: i34(:),elnode(:,:) integer :: iday1, iday2, iskipst allocatable :: sigma(:),cs(:),ztot(:) allocatable :: outvar(:,:,:,:),icum(:,:,:),eta2(:,:),ztmp(:,:) allocatable :: xcj(:),ycj(:),dps(:),kbs(:),xctr(:),yctr(:),dpe(:),kbe(:) allocatable :: idry(:),outs(:,:,:),residual(:,:),icounter(:) allocatable :: ic3(:,:),isdel(:,:),isidenode(:,:),kbp(:),sigma_lcl(:,:) real*8,allocatable :: timeout(:) pi=3.1415926 !... Set max array size for system memory print*, 'Input max array size (e.g., <=2.e9):' read(*,*)max_array_size print*, 'max_array_size read in=',max_array_size print*, 'Do you work on uncombined (0) or combined (1) nc?' read(*,*)icomb if(icomb/=0.and.icomb/=1) stop 'Unknown icomb' print*, 'Input variable name (e.g. salt):' read(*,'(a30)')varname varname=adjustl(varname); len_var=len_trim(varname) print*, 'Input start and end stack #s to read:' read(*,*) iday1,iday2 print*, 'Input stride in stacks (1 to include all):' read(*,*) iskipst print*, 'Input start and end record #s in start|end stack respectively:' !' read(*,*) irec_start,irec_end print*, 'Input z-coord. (<=0 below MSL):' read(*,*) z00 ! if(mod(iday2-iday1,iskipst)/=0) then ! write(*,*)'should be n skips over stack1 and stack2:',iday1,iday2,iskipst ! stop ! endif ! file63=adjustl(file63) ! len_file63=len_trim(file63) !... 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(iday1) endif print*, 'After header:',ne,np,nrec,i34(ne), & &elnode(1:i34(ne),ne),nvrt,h0,x(np),y(np),dp(np) !,start_time !... Read in time records in segments for mem last_dim=np nrec3=max_array_size/(2*nvrt*last_dim)-1 nrec3=min(nrec,max(nrec3,1)) print*, '# of records in each stack is: ',nrec, & &'; actual # of records that can be read in each time is:',nrec3 allocate(timeout(nrec),ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),kbp(np), & &outvar(2,nvrt,last_dim,nrec3),eta2(np,nrec3),ztmp(nvrt,np),residual(np,2),icounter(np),idry(np)) outvar=-huge(1.0) !test mem call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp) print*, 'last element',elnode(1:i34(ne),ne) !... Time iteration !... out=-99 !init. ztmp=-99 residual=0 ! ndays=iday2-iday1+1 icounter=0 !counter for each node (wet/dry) do iday=iday1,iday2,iskipst !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! write(it_char,'(i12)')iday ! it_char=adjustl(it_char) ! leng=len_trim(it_char) ! file63='schout_'//it_char(1:leng)//'.nc' ! iret=nf90_open(trim(adjustl(file63)),NF90_NOWRITE,ncid) ! !time is double ! iret=nf90_inq_varid(ncid,'time',itime_id) ! iret=nf90_get_var(ncid,itime_id,timeout,(/1/),(/nrec/)) if(icomb==0) then !uncombined call get_timeout(iday,nrec,timeout,icomb) else call get_timeout(iday,nrec,timeout) endif ! print*, 'time=',timeout irec1=1 !start record loop1: do irec2=min(nrec,irec1+nrec3-1) if(icomb==0) then !uncombined do irank=0,nproc-1 call get_outvar_multirecord(iday,varname,irec1,irec2,np,last_dim,nvrt,nrec3,outvar,i23d,ivs,eta2,irank) enddo !irank else call get_outvar_multirecord(iday,varname,irec1,irec2,np,last_dim,nvrt,nrec3,outvar,i23d,ivs,eta2) endif ! print*, 'done reading records:',irec1,irec2 !Available now: outvar(2,nvrt,np|ne,irec2-irec1+1),i23d,ivs,eta2(np,irec2-irec1+1) do irec=1,irec2-irec1+1 !offeset record # !---------------------------------------------------------------------------- irec_real=irec1+irec-1 !actual record # if(mod(i23d-1,3)==0) then !2D !print*,'irec=',irec,iday1,irec_start,iday2,irec_end if(.not.(iday==iday1.and.irec_realirec_end)) then do i=1,np icounter(i)=icounter(i)+1 residual(i,1:ivs)=residual(i,1:ivs)+outvar(1:ivs,1,i,irec) enddo !i endif else !3D !Compute z coordinates do i=1,np if(ivcor==1) then !localized if(dp(i)+eta2(i,irec)<=h0) then idry(i)=1 else !wet idry(i)=0 do k=kbp(i),nvrt ztmp(k,i)=(eta2(i,irec)+dp(i))*sigma_lcl(k,i)+eta2(i,irec) enddo !k endif !dp else if(ivcor==2) then !SZ call zcor_SZ_single(dp(i),eta2(i,irec),h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot, & &sigma,ztmp(:,i),idry(i),kbpl) kbp(i)=kbpl endif !ivcor if(idry(i)==0) then !wet !Interplate in vertical if(z00>=ztmp(nvrt,i)) then !above F.S. k0=nvrt-1; rat=1 else if(z00<=ztmp(kbp(i),i)) then !below bottom; extrapolate k0=kbp(i); rat=0 else !above bottom; cannot be above F.S. k0=0 do k=kbp(i),nvrt-1 if(z00>=ztmp(k,i).and.z00<=ztmp(k+1,i)) then k0=k rat=(z00-ztmp(k,i))/(ztmp(k+1,i)-ztmp(k,i)) exit endif enddo !k endif !ztmp if(k0==0) then write(*,*)'failed to find a vertical level:',irec_real,i,ifs,z2,ztmp(:,i) stop endif if(.not.(iday==iday1.and.irec_realirec_end)) then icounter(i)=icounter(i)+1 do m=1,ivs tmp=outvar(m,k0,i,irec)*(1-rat)+outvar(m,k0+1,i,irec)*rat residual(i,m)=residual(i,m)+tmp enddo !m endif !not endif !idry enddo !i=1,np endif !2/3D !---------------------------------------------------------------------------- enddo !irec if(irec2==nrec) exit loop1 irec1=irec1+nrec3 end do loop1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ enddo !iday ! Output do i=1,np if(icounter(i)==0) then residual(i,:)=0 else residual(i,:)=residual(i,:)/icounter(i) endif !icounter enddo !i open(65,file='average.out') !For wave dir if(varname(1:len_var).eq.'WWM_18') then !Read in sub-samples open(20,file='nodeflags.bp',status='old') read(20,*); read(20,*) eta2(:,1)=residual(:,1) !temp save do i=1,np read(20,*)j,tmp,tmp,tmp2 residual(i,1)=-sin(eta2(i,1)/180*pi) residual(i,2)=-cos(eta2(i,1)/180*pi) if(nint(tmp2)==0) residual(i,1:2)=0 enddo !i ivs=2 !xyuv format endif !file63 if(ivs==1) then write(65,*)iday1,iday2 write(65,*)ne,np do i=1,np write(65,*)i,x(i),y(i),residual(i,1) enddo !i do i=1,ne write(65,*)i,i34(i),elnode(1:i34(i),i) enddo !i else !vectors do i=1,np write(65,*)x(i),y(i),residual(i,1:2) enddo !i endif close(65) print*, 'Finished!' stop end