Revert some changes due to merge (b/c the reverting inside wwm_lr)

sample_inputs/sediment.nml

@@ -149,7 +149,7 @@ acrit = 0.2d0
 !- > tau_max: maximum shear stress value authorised for current and waves (in Pa) 
 !------------------------------------------------------------------------------
 tau_option = 1
-tau_max = 10.0d0 [Pa]
+tau_max = 10.0d0 ![Pa]
 zstress = 0.2d0 ![m]; only used if tau_option/=1
 
 !- Erosional formulations

src/Hydro/schism_init.F90

@@ -430,7 +430,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
 
       !Init defaults
       indvel=0; iupwind_mom=0; ihorcon=0; hvis_coef0=0.025_rkind; ishapiro=1; shapiro0=0.5_rkind; niter_shap=1
-      gen_wsett=real(1.d-4,rkind); flag_fib=1; ics=1; rearth_pole=6378206.4_rkind; rearth_eq=6378206.4_rkind; 
+      gen_wsett=real(0.d0,rkind); flag_fib=1; ics=1; rearth_pole=6378206.4_rkind; rearth_eq=6378206.4_rkind; 
       imm=0; ibdef=10; ihot=0; ihydraulics=0; izonal5=0; slam0=-124._rkind; sfea0=45._rkind; 
       ihdif=0; thetai=0.6_rkind; nrampbc=0; drampbc=1._rkind;  
       nramp=1; dramp=1._rkind; nadv=1; dtb_min=10._rkind; dtb_max=30._rkind; h0=0.01_rkind; nchi=0; dzb_min=0.5_rkind; dzb_decay=0._rkind;  
@@ -4487,7 +4487,8 @@ subroutine schism_init(iorder,indir,iths,ntime)
 
 #ifdef USE_AGE
       !Tracer age
-      !Method: all i.c. =0; conc=1 at relevant bnd(s), and itrtype=0 at ocean bnd 
+      !Method: all i.c. =0 except 1 in specific regions (e.g. near each bnd) for
+      !1st half of tracers (=0 for 2nd half). Generally set itrtype=0 at all bnd's
       if(myrank==0) write(16,*)'tracer age calculation evoked'
       nelem_age(:)=0 !init count for # of non-0 age tracers
       flag_ic(4)=1 !AGE must have type 1 i.c.

src/Hydro/schism_step.F90

@@ -37,6 +37,7 @@ subroutine schism_step(it)
 #endif
 
 #ifdef USE_FABM
+#include "fabm_version.h"
       USE fabm_schism, only: fabm_schism_do, fs, fabm_istart => istart
       USE fabm_schism, only: fabm_schism_write_output_netcdf
 #endif
@@ -8183,7 +8184,12 @@ subroutine schism_step(it)
 
 #ifdef USE_FABM
         do i=1,ntrs(11)
+!          call writeout_nc(id_out_var(noutput+i+4),trim(fs%model%state_variables(i)%name),2,nvrt,npa,tr_nd(i+fabm_istart-1,:,:))
+#if _FABM_API_VERSION_ < 1
           call writeout_nc(id_out_var(noutput+i+4),trim(fs%model%state_variables(i)%name),2,nvrt,npa,tr_nd(i+fabm_istart-1,:,:))
+#else
+          call writeout_nc(id_out_var(noutput+i+4),trim(fs%model%interior_state_variables(i)%name),2,nvrt,npa,tr_nd(i+fabm_istart-1,:,:))
+#endif
         end do
         noutput=noutput+ntrs(11)
 

src/ICM/icm.F90

@@ -1995,6 +1995,7 @@ subroutine calkwq(id,nv,ure,it)
 
     !NO3
     a=0.0
+    b=0.0
     rtmp=-ANC(1)*(1.0-PrefN(k,1))*GP(klev,id,1)*PB1(k,1)-ANC(2)*(1.0-PrefN(k,2))*GP(klev,id,2)*PB2(k,1)-ANC(3)*(1.0-PrefN(k,3))*GP(klev,id,3)*PB3(k,1)
     b=b+rtmp
     if(iof_icm(84)==1) absNO3(klev,id)=rtmp

src/Sediment/sed_bedload.F90

@@ -41,7 +41,7 @@ SUBROUTINE sed_bedload_sd(ised,inea)
 !          routines                                                  !
 !          2013/03 - F.Ganthy : implement wave effects on bedload    !
 !          2017/06 - A. de Bakker: couple to schism                  !
-!          2020/02 - B.Mengual: > use of ustress/vstress, wdir       !
+!          2020/02 - B.Mengual: > use of [uv]bott, wdir       !
 !                       > wave asymmetry effect: rewritten and       !
 !                         introduction of a limiter w_asym_max       !
 !                       > bed slope effects: now under option with   !
@@ -118,8 +118,8 @@ SUBROUTINE sed_bedload_sd(ised,inea)
       !cff2 = sum(vv2(kb1,elnode(1:i34(inea),inea)))/i34(inea)
       !udir = DATAN2(cff2,cff1)
 
-      ! BM: see ustress, vstress estimates in sediment.F90
-      udir = ATAN2(vstress(inea),ustress(inea))
+      ! BM: see [uv]bott estimates in sediment.F90
+      udir = ATAN2(vbott(inea),ubott(inea))
 
 !--------------------------------------------------------------------
 !- Angle between current and wave directions                         
@@ -537,7 +537,7 @@ SUBROUTINE sed_bedload_wl14(ised,inea)
 !
 ! History: - 2019/02 - B.Mengual : Implementation of the original code
 !                                  of T.Guerin in Sed2d
-!          - 2020/02 - B.Mengual: > use of ustress/vstress, wdir       
+!          - 2020/02 - B.Mengual: > use of [uv]bott, wdir       
 !                         > wave asymmetry effect: rewritten and       
 !                           introduction of a limiter w_asym_max       
 !                         > bed slope effects: now under option with   
@@ -607,9 +607,9 @@ SUBROUTINE sed_bedload_wl14(ised,inea)
       !cff1 = sum(uu2(kb1,elnode(1:i34(inea),inea)))/i34(inea)
       !cff2 = sum(vv2(kb1,elnode(1:i34(inea),inea)))/i34(inea)
 
-      ! BM: see ustress, vstress estimates in sediment.F90
-      cff1 = ustress(inea)
-      cff2 = vstress(inea)
+      ! BM: see [uv]bott estimates in sediment.F90
+      cff1 = ubott(inea)
+      cff2 = vbott(inea)
 
       Uc = dsqrt(cff1*cff1+cff2*cff2)
       udir = DATAN2(cff2,cff1)

src/Sediment/sed_init.F90

@@ -147,10 +147,10 @@ SUBROUTINE sed_alloc()
       !BM
       ALLOCATE(wdir(nea),stat=i)
         IF(i/=0) CALL parallel_abort('Sed: wdir allocation failure')
-      ALLOCATE(ustress(nea),stat=i)
-        IF(i/=0) CALL parallel_abort('Sed: ustress allocation failure')
-      ALLOCATE(vstress(nea),stat=i)
-        IF(i/=0) CALL parallel_abort('Sed: vstress allocation failure')
+      ALLOCATE(ubott(nea),stat=i)
+        IF(i/=0) CALL parallel_abort('Sed: ubott allocation failure')
+      ALLOCATE(vbott(nea),stat=i)
+        IF(i/=0) CALL parallel_abort('Sed: vbott allocation failure')
       ALLOCATE(U_crest(nea),stat=i)
         IF(i/=0) CALL parallel_abort('Sed: U_crest allocation failure')
       ALLOCATE(U_trough(nea),stat=i)

src/Sediment/sed_mod.F90

@@ -288,8 +288,8 @@ MODULE sed_mod
        REAL(rkind), ALLOCATABLE :: poron(:)   ! (npa)
        REAL(rkind), ALLOCATABLE :: Qaccun(:), Qaccvn(:) ! (npa) 
 
-       !BM Velocity components used in bottom shear stress estimates
-       REAL(rkind), ALLOCATABLE :: ustress(:), vstress(:) ! (nea)
+       !Near bottom velocity components used in bottom shear stress estimates
+       REAL(rkind), ALLOCATABLE :: ubott(:), vbott(:) ! (nea)
 
        !BM Wave-induced bedload transport caused by acceleration-skewness
        REAL(rkind), ALLOCATABLE :: Qaccu(:) ! (nea)

src/Sediment/sediment.F90

@@ -130,7 +130,7 @@ subroutine sediment(it,moitn,mxitn,rtol,dave,tot_bedmass)
       USE sed_mod
       USE schism_glbl, ONLY : rkind,nvrt,ne,nea,npa,np,irange_tr,ntrs,idry_e,   &
      & idry,area,znl,dt,i34,elnode,xctr,yctr,   &
-     & kbe,ze,pi,nne,indel,tr_el,flx_bt,    &
+     & kbp,kbe,ze,pi,nne,indel,tr_el,flx_bt,    &
      & errmsg,ielg,iplg,nond_global,iond_global,&
      & ipgl,nope_global,np_global,dp,h0,dpe,    &
      & iegl,out_wwm,pi,eta2,dp00,dldxy,we,time_stamp, &
@@ -293,9 +293,9 @@ subroutine sediment(it,moitn,mxitn,rtol,dave,tot_bedmass)
       T_trough = 0.d0
       Uorbi_elfrink = 0.0d0
 
-      ! Bottom shear stress
-      ustress=0.0d0
-      vstress=0.0d0
+      ! Near bottom vel
+!      ubott=0.0d0
+!      vbott=0.0d0
 
       ! RUNTIME
       time=dt*it
@@ -388,60 +388,6 @@ subroutine sediment(it,moitn,mxitn,rtol,dave,tot_bedmass)
       ENDDO ! End loop nea
 #endif
 
-!---------------------------------------------------------------------
-! - Compute bottom stress
-!   Test of whether WWM is used or not is performed within the subroutine
-!---------------------------------------------------------------------
-
-!BM: Several methods to compute current-induced bottom shear stress 
-!    (tau_option, see sediment.in).
-!    Current components (ustress, vstress) are then used in sed_current_stress 
-!    to compute BSS, or in sed_bedload routines to derive current direction.
-
-      DO i=1,nea
-        IF (idry_e(i)==1) CYCLE  
-        kb0 = kbe(i)
-        kb1 = kbe(i)+1
-        dzb=ze(kb1,i)-ze(kb0,i)
-
-        IF ((tau_option .EQ. 1) .OR. (tau_option .EQ. 3 .AND. dzb .GE. zstress)) THEN
-            ustress(i) = sum(uu2(kb1,elnode(1:i34(i),i)))/i34(i)
-            vstress(i) = sum(vv2(kb1,elnode(1:i34(i),i)))/i34(i)
-
-        ELSE IF ((tau_option .EQ. 2) .OR. (tau_option .EQ. 3 .AND. dzb .LT. zstress)) THEN
-          DO j=1,i34(i)
-            nd=elnode(j,i)
-            kk=nvrt
-            DO k=kb1,nvrt
-              IF (znl(k,nd)-znl(kb0,nd) .GE. zstress) THEN
-                kk=k
-                EXIT
-              END IF
-            END DO
-
-            IF (kk .LT. nvrt) THEN
-              wkk=(zstress+znl(kb0,nd)-znl(kk-1,nd))/(znl(kk,nd)-znl(kk-1,nd))
-              wkkm1=(znl(kk,nd)-zstress-znl(kb0,nd))/(znl(kk,nd)-znl(kk-1,nd))
-              ustress(i)=ustress(i)+(wkkm1*uu2(kk-1,nd) + wkk*uu2(kk,nd))
-              vstress(i)=vstress(i)+(wkkm1*vv2(kk-1,nd) + wkk*vv2(kk,nd))
-
-            ELSE
-              ustress(i)=ustress(i)+dav(1,nd)
-              vstress(i)=vstress(i)+dav(2,nd)
-
-            END IF
-          END DO
-          ustress(i) = ustress(i)/i34(i)
-          vstress(i) = vstress(i)/i34(i)
-
-        ELSE
-          call parallel_abort('SED: unknown tau_option')
-        END IF ! test on tau_option and dzb
-      END DO
-
-      CALL exchange_e2d(ustress)
-      CALL exchange_e2d(vstress)
-
       ! Compute current/wave-induced and total shear stresses
       CALL sed_wavecurrent_stress()
 

src/Utility/Pre-Processing/NWM/Elev_IC/gen_elev.f90

@@ -4,20 +4,20 @@
 !     ifort -O2 -mcmodel=medium -CB -Bstatic -o gen_elev gen_elev.f90
       implicit real*8(a-h,o-z)
       allocatable :: x(:),y(:),dp(:),area(:),vso(:),r_rough(:),nlayers(:)
-      integer, allocatable :: i34(:),elnode(:,:),i_rain(:),isource(:),iDB(:)
+      integer, allocatable :: i34(:),elnode(:,:),i_rain(:),isource(:),i_region(:)
 
-      h0=1e-6
+      h0=1e-4
 
       open(10,file='include.gr3',status='old')
       read(10,*); read(10,*)
 
       open(14,file='hgrid.gr3')
       read(14,*)
       read(14,*)ne,np
-      allocate(x(np),y(np),dp(np),r_rough(np),nlayers(np),iDB(np),i_rain(np),isource(ne),i34(ne),elnode(4,ne),area(ne),vso(ne))
+      allocate(x(np),y(np),dp(np),r_rough(np),nlayers(np),i_region(np),i_rain(np),isource(ne),i34(ne),elnode(4,ne),area(ne),vso(ne))
       do i=1,np
         read(14,*)j,x(i),y(i),dp(i)
-        read(10,*)dummy,dummy,dummy,iDB(i)
+        read(10,*)dummy,dummy,dummy,i_region(i)
       enddo !i
       close(10)
 
@@ -29,10 +29,10 @@
       open(9,file='elev.ic',status='replace')
       write(9,*); write(9,*)ne,np
       do i=1,np
-        if (iDB(i) ==0 ) then
+        if (i_region(i) ==0 ) then
           write(9,'(i8,3(1x,f15.6))')i,x(i),y(i),max(0.d0,-dp(i))
         else
-          write(9,'(i8,3(1x,f15.6))')i,x(i),y(i),max(0.d0,-dp(i)-1e-6)
+          write(9,'(i8,3(1x,f15.6))')i,x(i),y(i),max(0.d0,-dp(i)-h0)
         endif
       enddo 
       do i=1,ne

src/Utility/Pre-Processing/NWM/NWM_coupling/README

@@ -4,8 +4,8 @@ The main purpose of coupling_nwm.f90 is to look for the intersections between SC
 
 To run this script, 3 inputs are required:
 1) hgrid.ll (SCHISM horizontal grid)
-2) NWM_shp_ll.nc (shape file containing NWM segments)
-3) NWM outputs
+2) NWM_shp_ll.nc (nc format file containing information of NWM segments. It has 4 variables: featureID, lat, lon, and ORIG_FID. It is originated from a geo-database downloaded from https://water.noaa.gov/about/nwm. The process of making this file is described in NWM_nc_note.txt. The NWM_shp_ll.nc provided under this folder only contains East Coast and Gulf of Mexico region.)
+3) NWM outputs: a directory or symbolic link named "NWM_DATA" containing NWM outputs (e.g. 200508080900.CHRTOUT_DOMAIN1.comp) convering the simulation period.
 Some parameters are also required:
 1) Input nudging ratio (suggest 1.e-3)
 2) Input number of days
@@ -29,7 +29,8 @@ Inputs files:
       dist: distance in meters if inbr=0
       [inbr=0 and dist=4000 is recommended]
 Output files
-      source_sink.in.1, msource.th.1, vsource.th.1, vsink.th.1
+      vsource.th.1, vsink.th.1
+msource.th and source_sink.in will not change.
 
 To run those scripts, after compile. simply by
 ./coupling_nwm

src/Utility/Pre-Processing/NWM/NWM_coupling/coupling_nwm.f90

@@ -4,7 +4,8 @@
 !It assumes the land boundary is the 1st land boundary and there is only
 !one land boundary.
 !Inputs files: 
-!      NWM_shp_ll.nc : which is converted from the shape file using the NWM streams. 
+!      NWM_shp_ll.nc: which is converted from the shape file using the NWM streams. 
+!      NWM_DATA: A folder or symlink containing NWM model outputs 
 !      hgrid.ll: the gr3 format SCHISM grid; please only keep land
 !      bnd segments needed for coupling (islands are OK). The # of land
 !      bnd nodes is only used for dimensioning so give it a large number
@@ -32,7 +33,7 @@
 !
 ! 
 !
-!ifort -O2 -CB -g -traceback -o coupling_nwm  ~/git/schism/src/Utility/UtilLib/julian_date.f90 ~/git/schism/src/Utility/UtilLib/schism_geometry.f90 ~/git/schism/src/Utility/UtilLib/pt_in_poly_test.f90 coupling_nwm.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -L$NETCDF/lib -lnetcdf -lnetcdff 
+!ifort -O2 -CB -traceback -o coupling_nwm ../../../UtilLib/julian_date.f90 ../../../UtilLib/schism_geometry.f90 ../../../UtilLib/pt_in_poly_test.f90 coupling_nwm.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -L$NETCDF/lib -lnetcdf -lnetcdff
 
        program coupling_nwm
        use netcdf
@@ -41,8 +42,7 @@ program coupling_nwm
        use ifport
        implicit real*8(a-h,o-z)
 
-       character(len=*),parameter::REPODir='/sciclone/home20/whuang07/&
-                                             data10/NWM/CHRTOUT/'
+       character(len=*),parameter::REPODir='./NWM_DATA/'
 
        character(len=*),parameter::FILE_NAME='NWM_shp_ll.nc'
        integer :: ncid
@@ -82,16 +82,19 @@ program coupling_nwm
        !of NWM)
        print*, 'Input nudging ratio (suggest 1.e-3):'
        read*, epsilon
+       write(*,*) "Nudging ratio: ", epsilon
 
        !print*, 'Input rain_rate:'
        !rain_rate=0.03 !m/hour
 
        print*, 'Input number of days'
        read*,nday
+       write(*,*) "Number of days: ", nday
        ntime=nday*24
 
        print*, 'Enter start time - dd mm yyyy (e.g. 24 07 1988)'
        read(*,*) idd,imm,iyy
+       write(*,*) "start time: ", idd, imm, iyy
 
 
        open(14,file='hgrid.ll',status='old')!lambert projection, NWM has shorter precision
@@ -501,6 +504,7 @@ program coupling_nwm
 
         NWMfile=REPODir//yy//mm1//dd1&
                      //hh1//hh2//'00.CHRTOUT_DOMAIN1.comp'
+        write(*,*)'reading ', NWMfile
         write(98,*)'Inside time iteration:',i,NWMfile 
 
 

src/Utility/Pre-Processing/NWM/Vgrid/gen_vqs.f90

@@ -7,7 +7,7 @@
 !            (3) screen; (4) transect.bp (optional transect bp file)
 !    Outputs: vgrid.in; vgrid_master.out;  transect*.out; debug outputs (fort*)
 !    Use plot_VQS.m to viz vgrid_master.out; transect*.out
-!    pgf90 -O2 -mcmodel=medium -Mbounds -Bstatic -o gen_vqs gen_vqs.f90 ~/SELFE/svn/trunk/src/Utility/UtilLib/schism_geometry.f90
+!    Sample compilation command:
 !    ifort -Bstatic -O2 -CB -o gen_vqs gen_vqs.f90 schism_geometry.f90
 
       implicit real*8(a-h,o-z)
@@ -20,11 +20,18 @@
 
 !     Read  vgrid.in
       !m_vqs: # of master grids
-      m_vqs=13; n_sd=9
-      dz_bot_min=0.5 !min. bottom layer thickness [m]
+      m_vqs=22; n_sd=18
+      dz_bot_min=0.1 !min. bottom layer thickness [m]
       allocate(hsm(m_vqs),nv_vqs(m_vqs),a_vqs(m_vqs))
-      hsm=(/0.5,2.,3.,4.,5.,10.,50.,100.,350.,1050.,2000.,5000.,8500./)
-      nv_vqs(1:n_sd)=(/2,4,5,6,7,8,10,12,21/) !# of levels for each master grid (increasing with depth)
+      hsm=(/1., 2., 3., &
+          & 4., 6., 8., 12., 18., &
+          & 25., 33., 42., 52., 67., &
+          & 83., 100., 150., 230., 350., &
+          & 1050., 2000., 5000., 8500./)
+      nv_vqs(1:n_sd)=(/2,3,5, &
+                     & 7,9,11,13,15, &
+                     & 17,17,18,18,19, &
+                     & 20,21,22,24,26/) !# of levels for each master grid (increasing with depth) 
       nv_vqs(n_sd+1)=nv_vqs(n_sd)+9
       nv_vqs(n_sd+2)=nv_vqs(n_sd+1)+5
       nv_vqs(n_sd+3)=nv_vqs(n_sd+2)+5
@@ -59,27 +66,24 @@
       print*, 'nvrt in master vgrid=',nvrt_m
       allocate(z_mas(nvrt_m,m_vqs))
       z_mas=-1.e5
+      theta_b=0.0
       do m=1,n_sd !m_vqs
+        if (m<=7) then
+          theta_f=0.0001
+        elseif (m<=17) then
+          theta_f=min(1.0,max(0.0001, (m-4)/10.0)) * 3.0
+        else
+          theta_f=4.4
+        endif
+        if (m==14) theta_f=theta_f-0.1
+        if (m==15) theta_f=theta_f+0.1
+        if (m==16) theta_f=theta_f+0.55
+        if (m==17) theta_f=theta_f+0.97
         do k=1,nv_vqs(m)
           sigma=(k-1.0)/(1.0-nv_vqs(m))
-
-!         Alternative transformations below
-!         Option 1: quadratic 
-!          a_vqs(m)=max(-1.d0,a_vqs0-(m-1)*0.03)
-!          tmp=a_vqs(m)*sigma*sigma+(1+a_vqs(m))*sigma !transformed sigma
-!          z_mas(k,m)=tmp*(etal+hsm(m))+etal
-
-!          Option 2: S
-          theta_b=0
-          if(m<=4) then
-            theta_f=2.5
-          else if(m==5) then
-            theta_f=2.65
-          endif
           cs=(1-theta_b)*sinh(theta_f*sigma)/sinh(theta_f)+ &
-     &theta_b*(tanh(theta_f*(sigma+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5)
+            &theta_b*(tanh(theta_f*(sigma+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5)
           z_mas(k,m)=etal*(1+sigma)+hsm(1)*sigma+(hsm(m)-hsm(1))*cs
-
         enddo !k
       enddo !m