Finished revising iout_sta=2; tested

Author: josephzhang8

src/Core/schism_glbl.F90

@@ -604,8 +604,8 @@ module schism_glbl
   real(rkind),save,allocatable :: xsta(:),ysta(:),zstal(:),zsta(:),arco_sta(:,:), &
                                   xsta_bmi(:),ysta_bmi(:),zsta_bmi(:), &
                                   &sta_out(:,:),sta_out_gb(:,:),sta_out3d(:,:,:), &
-                                  &zta_out3d(:,:,:),sta_out3d_gb(:,:,:),zta_out3d_gb(:,:,:)
-  integer,save,allocatable :: iep_sta(:),iep_flag(:),iof_sta(:),indx_out(:,:),indx_wwm_out(:)
+                                  &zta_out3d(:,:),sta_out3d_gb(:,:,:),zta_out3d_gb(:,:)
+  integer,save,allocatable :: idry_sta(:),iep_sta(:),iep_flag(:),iof_sta(:),indx_out(:,:),indx_wwm_out(:)
 
   ! Message passing arrays used in main
   integer,save,allocatable :: srqst(:),sstat(:,:),rrqst(:),rstat(:,:)

src/Hydro/schism_init.F90

@@ -4191,11 +4191,12 @@ subroutine schism_init(iorder,indir,iths,ntime)
 
 !       Allocate: zstal is vertical up; xsta, ysta, zsta are global coord. if ics=2
         if(iorder==0) then
-          allocate(xsta(nout_sta),ysta(nout_sta),zstal(nout_sta),zsta(nout_sta),iep_sta(nout_sta),iep_flag(nout_sta), &
+          allocate(xsta(nout_sta),ysta(nout_sta),zstal(nout_sta),zsta(nout_sta), &
+     &idry_sta(nout_sta),iep_sta(nout_sta),iep_flag(nout_sta), &
      &arco_sta(nout_sta,4),sta_out(nout_sta,nvar_sta),sta_out_gb(nout_sta,nvar_sta), &
      &sta_out3d(nvrt,nout_sta,nvar_sta),sta_out3d_gb(nvrt,nout_sta,nvar_sta), &
      &xsta_bmi(nout_sta),ysta_bmi(nout_sta),zsta_bmi(nout_sta), &
-     &zta_out3d(nvrt,nout_sta,nvar_sta),zta_out3d_gb(nvrt,nout_sta,nvar_sta),stat=istat)
+     &zta_out3d(nvrt,nout_sta),zta_out3d_gb(nvrt,nout_sta),stat=istat)
           if(istat/=0) call parallel_abort('MAIN: sta. allocation failure')
           iep_flag=0
           sta_out=0.0_rkind

src/Hydro/schism_step.F90

@@ -10330,15 +10330,17 @@ subroutine schism_step(it)
         do j=1,nvar_sta !excluding zcor
           if(iof_sta(j)==0.or.mod(it,nspool_sta)/=0) cycle
 
+          idry_sta(:)=0 !init
           do i=1,nout_sta
             ie=iep_sta(i)
             if(ie==0) then !no parent in this rank
               iep_flag(i)=0 !for comm. later
               sta_out(i,j)=0.d0
               sta_out3d(:,i,j)=0.d0
-              zta_out3d(:,i,j)=0.d0
+              zta_out3d(:,i)=0.d0 !same for all j
             else !is parent
               iep_flag(i)=1
+              idry_sta(i)=idry_e(ie) !save for comm later
               sta_out(i,j)=0.d0 !initialize
               if(j==1) then !elev.
                 swild2(1,1:i34(ie))=eta2(elnode(1:i34(ie),ie))
@@ -10371,7 +10373,7 @@ subroutine schism_step(it)
                 if(idry_e(ie)==1) then !dry
                   sta_out(i,j)=-1.d7 !-999.d0
                   sta_out3d(:,i,j)=-1.d7 !-999.d0
-                  zta_out3d(:,i,j)=-1.d7 !-999.d0
+                  zta_out3d(:,i)=-1.d7 !same for all j
                 else !wet
                   do m=1,i34(ie) !wet nodes
                     nd=elnode(m,ie)
@@ -10404,7 +10406,7 @@ subroutine schism_step(it)
                   itmp=minval(kbp(elnode(1:i34(ie),ie)))
                   do k=1,nvrt
                     if(k<itmp) then
-                      zta_out3d(k,i,j)=-1.d7
+                      zta_out3d(k,i)=-1.d7
                       sta_out3d(k,i,j)=-1.d7
                     else !at least 1 node has valid value
                       do m=1,i34(ie)
@@ -10419,7 +10421,7 @@ subroutine schism_step(it)
 !                          swild4(2,m)=swild2(k,m)
 !                        endif
                       enddo !m
-                      zta_out3d(k,i,j)=sum(arco_sta(i,1:i34(ie))*swild4(1,1:i34(ie)))
+                      zta_out3d(k,i)=sum(arco_sta(i,1:i34(ie))*swild4(1,1:i34(ie)))
                       sta_out3d(k,i,j)=sum(arco_sta(i,1:i34(ie))*swild4(2,1:i34(ie)))
                     endif !k
                   enddo !k
@@ -10432,10 +10434,14 @@ subroutine schism_step(it)
         enddo !j=1,nvar_sta
 
 !       Output by rank 0
+        call mpi_reduce(idry_sta,nwild2,nout_sta,itype,MPI_SUM,0,comm,ierr)
+        !Save back to idry_sta to free up nwild2; /=0 => dry
+        if(myrank==0) idry_sta=nwild2(1:nout_sta)
+
         call mpi_reduce(iep_flag,nwild2,nout_sta,itype,MPI_SUM,0,comm,ierr)
         call mpi_reduce(sta_out,sta_out_gb,nout_sta*nvar_sta,rtype,MPI_SUM,0,comm,ierr)
         call mpi_reduce(sta_out3d,sta_out3d_gb,nvrt*nout_sta*nvar_sta,rtype,MPI_SUM,0,comm,ierr)
-        call mpi_reduce(zta_out3d,zta_out3d_gb,nvrt*nout_sta*nvar_sta,rtype,MPI_SUM,0,comm,ierr)
+        call mpi_reduce(zta_out3d,zta_out3d_gb,nvrt*nout_sta,rtype,MPI_SUM,0,comm,ierr)
 
         if(myrank==0) then
 !          write(290,*)nwild2(1:nout_sta)
@@ -10448,30 +10454,30 @@ subroutine schism_step(it)
                 sta_out_gb(j,i)=1.d7 !-9999.d0
                 if(i>4) then !3D only
                   sta_out3d_gb(:,j,i)=1.d7 !-9999.d0
-                  zta_out3d_gb(:,j,i)=1.d7 !-9999.d0
+                  zta_out3d_gb(:,j)=1.d7 !-9999.d0
                 endif
               else
                 sta_out_gb(j,i)=sta_out_gb(j,i)/dble(nwild2(j))
                 if(i>4) then !3D only
                   sta_out3d_gb(:,j,i)=sta_out3d_gb(:,j,i)/dble(nwild2(j))
-                  zta_out3d_gb(:,j,i)=zta_out3d_gb(:,j,i)/dble(nwild2(j))
+                  zta_out3d_gb(:,j)=zta_out3d_gb(:,j)/dble(nwild2(j))
                 endif
               endif
             enddo !j
-            write(250+i,'(e24.16,6000(1x,e14.6))')time,sta_out_gb(:,i)
+            write(250+i,'(e24.16,6000(1x,e15.6e3))')time,sta_out_gb(:,i)
             if(iout_sta==2.and.i>4) then
-              write(250+i,'(e24.16,300000(1x,e14.6))')time,sta_out3d_gb(:,:,i) !,zta_out3d_gb(:,:,i)
+              write(250+i,'(e24.16,300000(1x,e15.6e3))')time,sta_out3d_gb(:,:,i) !,zta_out3d_gb(:,:,i)
               !Add zcor output: do it only once
               if(.not.ltmp) then
                 ltmp=.true.
-                write(250+nvar_sta+1,*) !empty line to be in consistent form as other 3D 
-                write(250+nvar_sta+1,'(e24.16,300000(1x,e14.6))')time,zta_out3d_gb(:,:,i)
+                write(250+nvar_sta+1,'(300000(1x,i5))')idry_sta(:)
+                write(250+nvar_sta+1,'(e24.16,300000(1x,e15.6e3))')time,zta_out3d_gb(:,:)
               endif !
             endif !iout_sta
           enddo !i=1,nvar_sta
 
           write(16,*)'done station outputs...'
-        endif !myrank
+        endif !myrank==0
       endif !iout_sta/=0
 
 #ifdef USE_HA

src/Utility/Post-Processing-Fortran/read_staout.f90

@@ -9,6 +9,7 @@
   character(len=30) :: fname
   integer :: time_dim,ndims(100),nwild(100)
   logical :: ltmp,lfirst
+  integer, allocatable :: idry_sta(:)
   real, allocatable :: wild0(:),wild(:,:),zcor(:,:),xsta(:),ysta(:),zsta(:)
   real :: a1d(1)
   integer :: iempty(nfiles)
@@ -20,7 +21,7 @@
 
   open(9,file='station.in',status='old')
   read(9,*); read(9,*)nsta
-  allocate(wild0(nsta),wild(nvrt,nsta),zcor(nvrt,nsta),xsta(nsta),ysta(nsta),zsta(nsta))
+  allocate(wild0(nsta),wild(nvrt,nsta),zcor(nvrt,nsta),xsta(nsta),ysta(nsta),zsta(nsta),idry_sta(nsta))
   do i=1,nsta
     read(9,*)j,xsta(i),ysta(i),zsta(i)
   enddo !i
@@ -78,6 +79,11 @@
   j=nf90_def_var(ncid,'station_location_z',NF90_FLOAT,ndims(1),nwild(19))
   j=nf90_put_att(ncid,nwild(19),'name','Z coorinates at each station in meters')
 
+  !Dry flags
+  ndims(1)=nsta_dim; ndims(2)=time_dim
+  j=nf90_def_var(ncid,'Station_dry_flags',NF90_INT,ndims(1:2),nwild(20))
+  j=nf90_put_att(ncid,nwild(20),'name','Non-zero if dry at station')
+
   do i=1,19 !nfiles+1+6
     j=nf90_def_var_deflate(ncid,nwild(i),0,1,4)
   enddo !i
@@ -130,11 +136,11 @@
       nt0=0
       do it=1,ntime
         nt0=nt0+1
-        !odd lines are empty for zcor
+        !odd lines are dry flags for zcor
         if(i<=nfiles) then
           read(10,*,end=99,err=99)time,wild0(:)
         else !empty line for zcor
-          read(10,*,end=99,err=99)
+          read(10,*,end=99,err=99)idry_sta(:)
         endif
         if(i>4) then !3D (incl zcor)
           read(10,*,end=99,err=99)time,wild(:,:) !,zcor(:,:)
@@ -155,6 +161,7 @@
           endif !i
         else !zcor
           j=nf90_put_var(ncid,nwild(16),wild(:,:),(/1,1,nt0/),(/nvrt,nsta,1/))
+          j=nf90_put_var(ncid,nwild(20),idry_sta,(/1,nt0/),(/nsta,1/))
         endif !i<=nfiles
       enddo !it=1,ntime