Add new option for source/sink: if_source=-1 (netcdf input source.nc)

src/Core/schism_glbl.F90

@@ -154,7 +154,7 @@ module schism_glbl
   character(len= 8),save :: a_8
   character(len= 4),save :: a_4
   integer,save :: ncid_nu(natrm),ncid_tr3D(natrm),ncid_elev2D,ncid_uv3D,ncid_schout, &
- &nstride_schout,nrec2_schout,istack0_schout
+ &nstride_schout,nrec2_schout,istack0_schout,ncid_source
 
   ! ADT for global-to-local linked-lists
   type :: llist_type
@@ -367,8 +367,9 @@ module schism_glbl
                                   &vmo(:,:,:),vfa(:,:,:),eth(:,:), &
                                   &qthcon(:),uthnd(:,:,:),vthnd(:,:,:), &
                                   &ath(:,:,:,:),carea(:),clen(:),eta_mean(:),q_block(:),vnth_block(:,:), &
-                                  &dir_block(:,:),q_block_lcl(:),ath3(:,:,:,:)
+                                  &dir_block(:,:),q_block_lcl(:)
   real(4),save,allocatable :: ath2(:,:,:,:,:) !used to read *.nc for b.c. time series
+  real(4),save,allocatable :: ath3(:,:,:,:) !used to read source/sink inputs
 
   ! Land boundary segment data
   integer,save :: nland_global                 ! Global number of land bndry segments

src/Hydro/misc_subs.F90

@@ -558,7 +558,7 @@ subroutine other_hot_init(time)
         endif !ntrs
       enddo !i
 
-      if(if_source==1) then
+      if(if_source==1) then !ASCII
         if(nsources>0) then
           open(63,file=in_dir(1:len_in_dir)//'vsource.th',status='old') !values (>=0) in m^3/s
           read(63,*)tmp,ath3(1:nsources,1,1,1)
@@ -606,6 +606,44 @@ subroutine other_hot_init(time)
         endif !nsinks
       endif !if_source
 
+      if(if_source==-1) then !nc
+        if(nsources>0) then
+          ninv=time/th_dt3(1)
+          th_time3(1,1)=dble(ninv)*th_dt3(1)
+          th_time3(2,1)=th_time3(1,1)+th_dt3(1)
+          j=nf90_inq_varid(ncid_source, "vsource",mm)
+          if(j/=NF90_NOERR) call parallel_abort('MISC: vsource')
+!Error: consider bcast
+          j=nf90_get_var(ncid_source,mm,ath3(1:nsources,1,1,1),(/1,ninv+1/),(/nsources,1/))
+          if(j/=NF90_NOERR) call parallel_abort('MISC: vsource(2)')
+          j=nf90_get_var(ncid_source,mm,ath3(1:nsources,1,2,1),(/1,ninv+2/),(/nsources,1/))
+          if(j/=NF90_NOERR) call parallel_abort('MISC: vsource(3)')
+
+          !msource
+          ninv=time/th_dt3(3)
+          th_time3(1,3)=dble(ninv)*th_dt3(3)
+          th_time3(2,3)=th_time3(1,3)+th_dt3(3)
+          j=nf90_inq_varid(ncid_source, "msource",mm)
+          if(j/=NF90_NOERR) call parallel_abort('MISC: msource')
+          j=nf90_get_var(ncid_source,mm,ath3(1:nsources,1:ntracers,1,3),(/1,1,ninv+1/),(/nsources,ntracers,1/))
+          if(j/=NF90_NOERR) call parallel_abort('MISC: msource(2)')
+          j=nf90_get_var(ncid_source,mm,ath3(1:nsources,1:ntracers,2,3),(/1,1,ninv+2/),(/nsources,ntracers,1/))
+          if(j/=NF90_NOERR) call parallel_abort('MISC: msource(2)')
+        endif !nsources>0
+
+        if(nsinks>0) then
+          ninv=time/th_dt3(2)
+          th_time3(1,2)=dble(ninv)*th_dt3(2)
+          th_time3(2,2)=th_time3(1,2)+th_dt3(2)
+          j=nf90_inq_varid(ncid_source, "vsink",mm)
+          if(j/=NF90_NOERR) call parallel_abort('MISC: vsink')
+          j=nf90_get_var(ncid_source,mm,ath3(1:nsinks,1,1,2),(/1,ninv+1/),(/nsinks,1/))
+          if(j/=NF90_NOERR) call parallel_abort('MISC: vsink(2)')
+          j=nf90_get_var(ncid_source,mm,ath3(1:nsinks,1,2,2),(/1,ninv+2/),(/nsinks,1/))
+          if(j/=NF90_NOERR) call parallel_abort('MISC: vsink(3)')
+        endif !nsinks>0
+      endif !if_source=-1
+
 #ifdef USE_SED
 !...  Sediment model initialization
       call sed_init

src/Hydro/schism_init.F90

@@ -1428,13 +1428,10 @@ subroutine schism_init(iorder,indir,iths,ntime)
       endif
 
 
-!     Volume and mass sources/sinks option
-!      call get_param('param.in','if_source',1,if_source,tmp,stringvalue)
-      if(if_source/=0.and.if_source/=1) call parallel_abort('Wrong if_source')
+!     Volume and mass sources/sinks option (-1:nc; 1:ASCII)
+      if(iabs(if_source)>1) call parallel_abort('Wrong if_source')
 
-      if(if_source==1) then
-!        call get_param('param.in','nramp_ss',1,nramp_ss,tmp,stringvalue)
-!        call get_param('param.in','dramp_ss',2,itmp,dramp_ss,stringvalue)
+      if(if_source/=0) then
         if(dramp_ss<=0) call parallel_abort('INIT: wrong dramp_ss')
       endif
 
@@ -2934,8 +2931,8 @@ subroutine schism_init(iorder,indir,iths,ntime)
 
       endif !ihydraulics/=0
 
-!     Read in source_sink.in and open t.h. files
-      if(if_source==1) then
+!     Read in source/sink info 
+      if(if_source==1) then !ASCII
         open(31,file=in_dir(1:len_in_dir)//'source_sink.in',status='old')
         read(31,*)nsources
         if(iorder==0) then
@@ -2958,6 +2955,60 @@ subroutine schism_init(iorder,indir,iths,ntime)
         close(31)
       endif !if_source
 
+      if(if_source==-1) then !nc
+        j=nf90_open(in_dir(1:len_in_dir)//'source.nc',OR(NF90_NETCDF4,NF90_NOWRITE),ncid_source)
+        if(j/=NF90_NOERR) call parallel_abort('init: source.nc')
+        j=nf90_inq_dimid(ncid_source,'nsources',mm)
+        j=nf90_inquire_dimension(ncid_source,mm,len=nsources)
+        if(j/=NF90_NOERR) call parallel_abort('init: nsources')
+        j=nf90_inq_dimid(ncid_source,'nsinks',mm)
+        j=nf90_inquire_dimension(ncid_source,mm,len=nsinks)
+        if(j/=NF90_NOERR) call parallel_abort('init: nsinks')
+        j=nf90_inq_dimid(ncid_source,'ntracers',mm)
+        j=nf90_inquire_dimension(ncid_source,mm,len=itmp)
+        if(itmp/=ntracers) call parallel_abort('init: wrong ntracers in source.nc')
+
+        j=nf90_inq_varid(ncid_source, "time_step_vsource",mm)
+        if(j/=NF90_NOERR) call parallel_abort('init: time_step_vsource')
+        j=nf90_get_var(ncid_source,mm,floatout)
+        if(j/=NF90_NOERR) call parallel_abort('init: time_step_vsource(2)')
+        if(floatout<dt) call parallel_abort('INIT: dt_vsource wrong')
+        th_dt3(1)=dble(floatout)
+
+        j=nf90_inq_varid(ncid_source, "time_step_msource",mm)
+        if(j/=NF90_NOERR) call parallel_abort('init: time_step_msource')
+        j=nf90_get_var(ncid_source,mm,floatout)
+        if(j/=NF90_NOERR) call parallel_abort('init: time_step_msource(2)')
+        if(floatout<dt) call parallel_abort('INIT: dt_msource wrong')
+        th_dt3(3)=dble(floatout)
+
+        j=nf90_inq_varid(ncid_source, "time_step_vsink",mm)
+        if(j/=NF90_NOERR) call parallel_abort('init: time_step_vsink')
+        j=nf90_get_var(ncid_source,mm,floatout)
+        if(j/=NF90_NOERR) call parallel_abort('init: time_step_vsink(2)')
+        if(floatout<dt) call parallel_abort('INIT: dt_vsink wrong')
+        th_dt3(2)=dble(floatout)
+
+        if(iorder==0) then
+          allocate(ieg_source(max(1,nsources)),ieg_sink(max(1,nsinks)), &
+     &ath3(max(1,nsources,nsinks),ntracers,2,nthfiles3),stat=istat)
+          if(istat/=0) call parallel_abort('INIT: ieg_source failure(3)')
+        endif
+
+        if(nsources>0) then
+          j=nf90_inq_varid(ncid_source, "source_elem",mm)
+          if(j/=NF90_NOERR) call parallel_abort('init: source_elem')
+          j=nf90_get_var(ncid_source,mm,ieg_source(1:nsources),(/1/),(/nsources/))
+          if(j/=NF90_NOERR) call parallel_abort('init: source_elem(2)')
+        endif !nsources
+
+        if(nsinks>0) then
+          j=nf90_inq_varid(ncid_source, "sink_elem",mm)
+          if(j/=NF90_NOERR) call parallel_abort('init: sink_elem')
+          j=nf90_get_var(ncid_source,mm,ieg_sink(1:nsinks),(/1/),(/nsinks/))
+          if(j/=NF90_NOERR) call parallel_abort('init: sink_elem(2)')
+        endif !nsinks
+      endif !if_source=-1
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 !  Initialize model for cold and hot start

src/Hydro/schism_step.F90

@@ -1363,24 +1363,41 @@ subroutine schism_step(it)
 !     Other tracers: 0 (otherwise additional nutrients from rain will fall onto
 !     water) 
       vsource=0 !init; dimension [m^3/s]; includes sinks as well
-      if(if_source==1) then
+      if(if_source/=0) then
         !init all first; dimension same as concentration (psu etc)
         msource=0.d0 
         !Exceptions
         msource(1:2,:)=-9999.d0 !junk so ambient values will be used
         if(nsources>0) then
           if(time>th_time3(2,1)) then !not '>=' to avoid last step
             ath3(:,1,1,1)=ath3(:,1,2,1)
-            read(63,*)tmp,ath3(1:nsources,1,2,1)
             th_time3(1,1)=th_time3(2,1)
             th_time3(2,1)=th_time3(2,1)+th_dt3(1)
+
+            if(if_source==1) then
+              read(63,*)tmp,ath3(1:nsources,1,2,1)
+            else !nc
+              itmp2=time/th_dt3(1)+2
+              j=nf90_inq_varid(ncid_source, "vsource",mm)
+              j=nf90_get_var(ncid_source,mm,ath3(1:nsources,1,2,1),(/1,itmp2/),(/nsources,1/))
+              if(j/=NF90_NOERR) call parallel_abort('STEP: vsource')
+            endif !if_source
           endif !time
 
+          !msource
           if(time>th_time3(2,3)) then !not '>=' to avoid last step
             ath3(:,:,1,3)=ath3(:,:,2,3)  
-            read(65,*)tmp,ath3(1:nsources,1:ntracers,2,3)
             th_time3(1,3)=th_time3(2,3)
             th_time3(2,3)=th_time3(2,3)+th_dt3(3)
+
+            if(if_source==1) then
+              read(65,*)tmp,ath3(1:nsources,1:ntracers,2,3)
+            else !nc
+              itmp2=time/th_dt3(3)+2
+              j=nf90_inq_varid(ncid_source, "msource",mm)
+              j=nf90_get_var(ncid_source,mm,ath3(1:nsources,1:ntracers,2,3),(/1,1,itmp2/),(/nsources,ntracers,1/))
+              if(j/=NF90_NOERR) call parallel_abort('STEP: msource')
+            endif !if_source
           endif !time
 
           rat=(time-th_time3(1,1))/th_dt3(1)
@@ -1390,17 +1407,18 @@ subroutine schism_step(it)
           endif
 
           do i=1,nsources
-            if(ath3(i,1,1,1)<0.d0.or.ath3(i,1,2,1)<0.d0) then
+            if(ath3(i,1,1,1)<0..or.ath3(i,1,2,1)<0.) then
               write(errmsg,*)'STEP: wrong sign vsource',it,i,ath3(i,1,1:2,1)
               call parallel_abort(errmsg)
             endif
 
             if(iegl(ieg_source(i))%rank==myrank) then
               ie=iegl(ieg_source(i))%id
-              vsource(ie)=vsource(ie)+((1-rat)*ath3(i,1,1,1)+rat*ath3(i,1,2,1))*ramp_ss
+              vsource(ie)=vsource(ie)+((1.d0-rat)*ath3(i,1,1,1)+rat*ath3(i,1,2,1))*ramp_ss
             endif !ielg
           enddo !i
 
+          !msource
           rat=(time-th_time3(1,3))/th_dt3(3)
           if(rat<-small1.or.rat>1.d0+small1) then
             write(errmsg,*) 'STEP: rat out in msource.th:',rat,time,th_time3(1:2,3)
@@ -1415,14 +1433,22 @@ subroutine schism_step(it)
               endif !ielg
             enddo !i
           enddo !j
-        endif !nsources
+        endif !nsources>0
 
         if(nsinks>0) then
           if(time>th_time3(2,2)) then !not '>=' to avoid last step
             ath3(:,1,1,2)=ath3(:,1,2,2)
-            read(64,*)tmp,ath3(1:nsinks,1,2,2)
             th_time3(1,2)=th_time3(2,2)
             th_time3(2,2)=th_time3(2,2)+th_dt3(2)
+
+            if(if_source==1) then
+              read(64,*)tmp,ath3(1:nsinks,1,2,2)
+            else !nc
+              itmp2=time/th_dt3(2)+2
+              j=nf90_inq_varid(ncid_source, "vsink",mm)
+              j=nf90_get_var(ncid_source,mm,ath3(1:nsinks,1,2,2),(/1,itmp2/),(/nsinks,1/))
+              if(j/=NF90_NOERR) call parallel_abort('STEP: vsink')
+            endif !if_source
           endif !time
 
           rat=(time-th_time3(1,2))/th_dt3(2)
@@ -1432,7 +1458,7 @@ subroutine schism_step(it)
           endif
 
           do i=1,nsinks
-            if(ath3(i,1,1,2)>0.d0.or.ath3(i,1,2,2)>0.d0) then
+            if(ath3(i,1,1,2)>0..or.ath3(i,1,2,2)>0.) then
               write(errmsg,*)'STEP: wrong sign vsink',it,i,ath3(i,1,1:2,2)
               call parallel_abort(errmsg)
             endif
@@ -1442,8 +1468,8 @@ subroutine schism_step(it)
               vsource(ie)=vsource(ie)+((1.d0-rat)*ath3(i,1,1,2)+rat*ath3(i,1,2,2))*ramp_ss
             endif !ielg
           enddo !i
-        endif !nsinks
-      endif !if_source
+        endif !nsinks>0
+      endif !if_source/=0
 
 !...  Volume sources from evap and precip
       if(isconsv/=0) then
@@ -6814,7 +6840,7 @@ subroutine schism_step(it)
 !$OMP ta,ie,kin,swild_m,swild_w,tmp0,vnf)
 
 !       Point sources/sinks using operator splitting (that guarentees max.
-!       principle); at bottom layer
+!       principle); imposed at bottom layer.
 !       Do nothing for net sinks
 !Error: need to reconcile with ICM