Merge hybrid_ELM into master

Squashed commit of the following:

commit feef498
Author: Fei Ye <feiye@vims.edu>
Date:   Mon Apr 6 15:03:24 2020 -0400

    Added hybrid ELM for itr_met==4. Passed preliminary tests (itr_met=3 and 4) Mistral.

commit e72e16b
Author: Joseph Zhang <yjzhang@vims.edu>
Date:   Sun Apr 5 14:49:43 2020 -0400

    Edit notes

commit e9d2348
Author: Joseph Zhang <yjzhang@vims.edu>
Date:   Sun Apr 5 14:47:28 2020 -0400

    SMall fixes in make files

commit 2bd700a
Author: Joseph Zhang <yjzhang@vims.edu>
Date:   Sat Apr 4 19:18:34 2020 -0400

    Add hybridized ELM in transport to bypass small transport time step.
    Tested with CORIE, lock exchange.

commit ccf5d97
Author: Joseph Zhang <yjzhang@vims.edu>
Date:   Fri Apr 3 12:37:24 2020 -0400

    minor edit

Author: josephzhang8

CONTRIBUTING.md

@@ -100,6 +100,9 @@ nothing to commit (working directory clean)
 
 `$ git push --set-upstream origin new_branch`
 
+(Get updates from other with pull)
+`$ git pull origin new_branch`
+
 When ready to merge your branch with master, there are two approaches: via pull request (PR) or local merge+push
 
 2.1 Local merge + push

cmake/SCHISM.local.cyclops

@@ -55,5 +55,5 @@ set(NetCDF_FORTRAN_DIR "$ENV{NETCDF_FORTRAN}" CACHE PATH "Path to NetCDF Fortran
 
 ##set(MPI_VERSION CACHE STRING  "2")
 ###Compile flags
-set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -xCORE-AVX2 -fma -align array64byte -finline-functions" CACHE STRING "Fortran flags" FORCE)
+set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -xCORE-AVX2 -fma -align array64byte -finline-functions -mcmodel medium -shared-intel" CACHE STRING "Fortran flags" FORCE)
 

mk/Make.defs.cyclops

@@ -41,13 +41,13 @@ ENV = cyclops
 # Environment 
 ################################################################################
 
-FCP = mpiifort
+FCP = mpif90
 FLD = $(FCP)
 FCS = $(FCP)
 PPFLAGS := $(PPFLAGS) -DMPIVERSION=2 
 FCPFLAGS = $(PPFLAGS) -O2 -xCORE-AVX2 #-qopenmp
 #FCSFLAGS = $(FFCPFLAGS)
-FLDFLAGS = -O2 -xCORE-AVX2 -fma -align array64byte -finline-functions 
+FLDFLAGS = -O2 -xCORE-AVX2 -fma -align array64byte -finline-functions -mcmodel medium -shared-intel #-fPIC  
 #####Libraries
 MTSLIBS = -L./ParMetis-3.1-Sep2010/ -lparmetis -lmetis
 CDFLIBS = -L$(NETCDF_FORTRAN)/lib -L$(NETCDF)/lib -lnetcdf -lnetcdff

sample_inputs/param.nml

@@ -401,6 +401,10 @@
   eps1_tvd_imp = 1.e-4 !suggested value is 1.e-4, but for large suspended load, need to use a smaller value (e.g. 1.e-9)
   eps2_tvd_imp = 1.e-14  
 
+  !Optional hybridized ELM transport for efficiency
+  ielm_transport = 0 !1: turn on. Used only with itr_met>=3
+  max_subcyc = 10 !used only if ielm_transport/=0. Max # of subcycling per time step in transport allowed
+
   !if itr_met = 4, the following parameters are needed
   !if itr_met=4 and ipre=1, diagnostic outputs are generated for weno accuracy and stencil quality, 
   !  see subroutine weno_diag in src/Hydro/misc_subs.F90 for details

src/Core/schism_glbl.F90

@@ -81,14 +81,16 @@ module schism_glbl
                       a2_cm03,schk,schpsi
 
   integer,parameter :: natrm=12 !# of _available_ tracer models at the moment (including T,S)
+  integer,parameter :: mntracers=30 !max # of tracers, used only for dimensioning btrack arrays. Must >=ntracers
+
   !Parameters from param.nml
   integer,save :: ipre,ipre2,indvel,imm,ihot,ics,iwbl,iharind,nws,impose_net_flux,iwindoff, &
                   &ibc,nrampbc,nrampwind,nrampwafo,nramp,nramp_ss,ibdef,ihorcon,nstep_wwm,icou_elfe_wwm, &
                   &iwind_form,irec_nu,itur,ihhat,inu_elev, &
                   &inu_uv,ibcc_mean,iflux,iout_sta,nspool_sta,nhot,nhot_write, &
                   &moitn0,mxitn0,nchi,ibtrack_test,nramp_elev,islip,ibtp,inunfl,shorewafo, &
                   &inv_atm_bnd,ieos_type,ieos_pres,iupwind_mom,inter_mom,ishapiro,isav, &
-                  &nstep_ice,niter_shap,iunder_deep,ibtrack_openbnd,flag_fib
+                  &nstep_ice,niter_shap,iunder_deep,ibtrack_openbnd,flag_fib,ielm_transport,max_subcyc
   integer,save :: ntrs(natrm),nnu_pts(natrm),mnu_pts
   integer,save,dimension(:),allocatable :: iof_hydro,iof_wwm,iof_gen,iof_age,iof_sed,iof_eco, &
      &iof_icm,iof_cos,iof_fib,iof_sed2d,iof_ice,iof_ana,iof_marsh,iof_dvd,iadjust_mass_consv
@@ -114,7 +116,9 @@ module schism_glbl
   real(rkind),save :: q2min,tempmin,tempmax,saltmin,saltmax, &
                       &vis_coe1,vis_coe2,h_bcc1,velmin_btrack,h_tvd,rmaxvel1,rmaxvel2, &
                       &difnum_max_l2,wtime1,wtime2,fluxsu00,srad00,cmiu0, &
-                      &cpsi2,rpub,rmub,rnub,cpsi1,psimin,eps_min,tip_dp,sav_di0,sav_h0,sav_nv0
+                      &cpsi2,rpub,rmub,rnub,cpsi1,psimin,eps_min,tip_dp,sav_di0,sav_h0,sav_nv0, &
+                      &dtb_min_transport
+
 !  logical,save :: lm2d !2D or 3D model
   logical,save :: lhas_quad=.false. !existence of quads
   logical,save :: lflbc !flag to indicate existence of ifltype/=0
@@ -172,7 +176,7 @@ module schism_glbl
     real(rkind) :: rt2        ! time remaining from left-over from previous subdomain 
     real(rkind) :: ut,vt,wt  ! Current backtracking sub-step velocity
     real(rkind) :: xt,yt,zt  ! Current backtracking sub-step point
-    real(rkind) :: sclr(4)     ! Backtracked values for some tracers
+    real(rkind) :: sclr(4+mntracers)     ! Backtracked values for tracers etc
     real(rkind) :: gcor0(3)  ! global coord. of the starting pt (for ics=2)
     real(rkind) :: frame0(3,3) ! frame tensor at starting pt (for ics=2)
   end type bt_type

src/Hydro/bktrk_subs.F90

@@ -64,9 +64,9 @@ subroutine init_inter_btrack
   if(ierr/=MPI_SUCCESS) call parallel_abort('INIT_INTER_BTRACK: mpi_get_address',ierr)
   types(1)=itype
 
-  ! Second part of bt_type is block of 26 doubles (including arrays)
+  ! Second part of bt_type is block of 26+mntracers doubles (including arrays)
   ! (starting at bttmp%dtbk)
-  blockl(2)=26
+  blockl(2)=26+mntracers
 #if MPIVERSION==1
   call mpi_address(bttmp%dtbk,displ(2),ierr)
 #elif MPIVERSION==2
@@ -710,7 +710,7 @@ subroutine btrack(l_ns,ipsgb,ifl_bnd,j0,iadvf,gcor0,frame0,dtbk, &
 !             nnel,jlev: initial and final element and level;
 !
 !       Output:
-!             sclr(4): btrack'ed values of some variables;
+!             sclr(4+mntracers): btrack'ed values of some variables;
 !             iexit: logical flag indicating backtracking exits augmented subdomain. If
 !                    iexit=.true., nnel is inside the aug. domain and should also be inside
 !                    one of the neighboring processes. (xt,yt) is inside nnel.
@@ -727,7 +727,7 @@ subroutine btrack(l_ns,ipsgb,ifl_bnd,j0,iadvf,gcor0,frame0,dtbk, &
       real(rkind), intent(in) :: gcor0(3),frame0(3,3),dtbk,vis_coe
       real(rkind), intent(inout) :: time_rm,time_rm2,uuint,vvint,wwint,xt,yt,zt
       integer, intent(inout) :: nnel,jlev
-      real(rkind), intent(out) :: sclr(4)
+      real(rkind), intent(out) :: sclr(4+mntracers)
       logical, intent(out) :: iexit
 
       !Local
@@ -941,6 +941,11 @@ subroutine btrack(l_ns,ipsgb,ifl_bnd,j0,iadvf,gcor0,frame0,dtbk, &
 !     Error: Kriging for wvel as well?
       if(l_ns==3) return
 
+      if(zrat<0._rkind.or.zrat>1._rkind) then
+        write(errmsg,*)'BTRACK: zrat wrong:',jlev,zrat
+        call parallel_abort(errmsg)
+      endif
+
 !     Calc max/min for ELAD
 !     If inter_mom/=0, sclr() will be updated below
       if(ibtrack_test==1) then
@@ -970,6 +975,18 @@ subroutine btrack(l_ns,ipsgb,ifl_bnd,j0,iadvf,gcor0,frame0,dtbk, &
           sclr(3)=max(sclr(3),vv2(jlev,nd),vv2(jlev-1,nd))
           sclr(4)=min(sclr(4),vv2(jlev,nd),vv2(jlev-1,nd))
         enddo !j
+
+        !Interp tracers
+        sclr(5:4+ntracers)=0.d0 
+        do i=1,ntracers
+          do j=1,i34(nnel)
+            nd=elnode(j,nnel)
+            !tmp=tr_nd(i,jlev,nd)*(1._rkind-zrat)+tr_nd(i,jlev-1,nd)*zrat
+            !Transport uses split, so keep vertical level at original
+            tmp=tr_nd(i,j0,nd)
+            sclr(4+i)=sclr(4+i)+tmp*arco(j)
+          enddo !j
+        enddo !i=1,ntracers
       endif !ibtrack_test
 
 !     Kriging for vel. (excluding bnd sides)
@@ -1058,11 +1075,6 @@ subroutine btrack(l_ns,ipsgb,ifl_bnd,j0,iadvf,gcor0,frame0,dtbk, &
       endif !Kriging
 
 !     nnel wet
-      if(zrat<0._rkind.or.zrat>1._rkind) then
-        write(errmsg,*)'BTRACK: zrat wrong:',jlev,zrat
-        call parallel_abort(errmsg)
-      endif
-
       end subroutine btrack
 
 !===============================================================================

src/Hydro/schism_init.F90

@@ -178,7 +178,8 @@ subroutine schism_init(indir,iths,ntime)
      &iharind,icou_elfe_wwm,nrampwafo,drampwafo,nstep_wwm,hmin_radstress,turbinj, &
      &iwbl,if_source,nramp_ss,dramp_ss,ieos_type,ieos_pres,eos_a,eos_b,slr_rate, &
      &rho0,shw,isav,sav_cd,nstep_ice,iunder_deep,h1_bcc,h2_bcc,hw_depth,hw_ratio, &
-     &ibtrack_openbnd,level_age,vclose_surf_frac,iadjust_mass_consv0,ipre2
+     &ibtrack_openbnd,level_age,vclose_surf_frac,iadjust_mass_consv0,ipre2, &
+     &ielm_transport,max_subcyc
 
      namelist /SCHOUT/iof_hydro,iof_wwm,iof_gen,iof_age,iof_sed,iof_eco,iof_icm,iof_cos,iof_fib, &
      &iof_sed2d,iof_ice,iof_ana,iof_marsh,iof_dvd, &
@@ -387,6 +388,7 @@ subroutine schism_init(indir,iths,ntime)
       !The big tracer arrays are: tr_el(ntracers,nvrt,nea2),tr_nd0(ntracers,nvrt,npa)
       !The order of each tracer modules can be seen above
       ntracers=sum(ntrs(:)) !including T,S
+     if(mntracers<ntracers) call parallel_abort('INIT: mntracers<ntracers')
 
       !'Init. ranges for each model. These index into 1:ntracers
       do i=1,natrm !# of _available_ tracer models at the moment
@@ -446,6 +448,7 @@ subroutine schism_init(indir,iths,ntime)
       vclose_surf_frac=1.0
       iadjust_mass_consv0=0 !Enforce mass conservation for a tracer 
       ipre2=0
+      ielm_transport=0; max_subcyc=10
 
       !Output elev, hvel by detault
       iof_hydro(1)=1; iof_hydro(25)=1
@@ -1473,6 +1476,11 @@ subroutine schism_init(indir,iths,ntime)
       endif      
 !...  TWO-PHASE-MIXTURE TSINGHUA GROUP------------------
 
+      if(ielm_transport/=0) then
+        if(max_subcyc<=0) call parallel_abort('INIT: max_subcyc<=0')
+        dtb_min_transport=dt/max_subcyc !min dt for transport allowed
+      endif 
+
 !...  Check parameter read in from param.in
       if(myrank==0) write(16,*)'done reading param.nml'
 !'
@@ -1655,7 +1663,7 @@ subroutine schism_init(indir,iths,ntime)
 !'
 
 !     All other arrays
-      allocate(sdbt(4,nvrt,nsa), & !webt(nvrt,nea), bubt(2,nea), & 
+      allocate(sdbt(2+ntracers,nvrt,nsa), & !webt(nvrt,nea), bubt(2,nea), & 
          &  windx1(npa),windy1(npa),windx2(npa),windy2(npa),windx(npa),windy(npa), &
          &  tau(2,npa),tau_bot_node(3,npa),iadv(npa),windfactor(npa),pr1(npa),airt1(npa),shum1(npa), &
          &  pr2(npa),airt2(npa),shum2(npa),pr(npa),sflux(npa),srad(npa),tauxz(npa),tauyz(npa), &

src/Hydro/schism_step.F90

@@ -197,7 +197,7 @@ subroutine schism_step(it)
                      &xlmax(nvrt),cpsi3(2:nvrt),cpsi2p(2:nvrt),q2ha(2:nvrt),xlha(2:nvrt), &
                      &chi(nsa),chi2(nsa),vsource(nea),sav_c2(nsa),sav_beta(nsa)
       real(rkind) :: swild(max(100,nsa+nvrt+12+ntracers)),swild2(nvrt,12),swild10(max(4,nvrt),12), &
-     &swild3(20+ntracers),swild4(2,4)
+     &swild3(20+mntracers),swild4(2,4)
 !#ifdef USE_SED
       real(rkind) :: swild_m(6,ntracers),swild_w(3),q2fha(2:nvrt),q2fpha(2:nvrt),epsftmp(nvrt), &
                      &Tpzzntr(nvrt),Dpzzntr(nvrt)  
@@ -3137,7 +3137,7 @@ subroutine schism_step(it)
 !     swild96(1:2,nvrt,nsa): \epsilon (over/under-shoots in ELAD) for u,v (if ibtrack_test=1, 1->T and 2 is not used)
 !     swild97(1:2,nvrt,nsa): u,v in the next iteration (ELAD). If ibtrack_test=1, 1->T (2 not used)
 !     swild98(1:4,nvrt,nsa): 1:2 max/min for u; 3:4 for v (only 1:2 are used for T)
-      allocate(swild96(2,nvrt,nsa),swild97(2,nvrt,nsa),swild98(4,nvrt,nsa),stat=istat)
+      allocate(swild98(4,nvrt,nsa),stat=istat)
       if(istat/=0) call parallel_abort('MAIN: fail to allocate swild98 (3.2)')
 
 !'    Debug: test backtracking alone
@@ -3324,6 +3324,7 @@ subroutine schism_step(it)
           if(vmag<=velmin_btrack) then !No activity 
             sdbt(1,j,isd0)=su2(j,isd0)
             sdbt(2,j,isd0)=sv2(j,isd0)
+            sdbt(3:2+ntracers,j,isd0)=(tr_nd(1:ntracers,j,n1)+tr_nd(1:ntracers,j,n2))*0.5d0
             swild98(1:2,j,isd0)=su2(j,isd0) !max/min
             swild98(3:4,j,isd0)=sv2(j,isd0)
 
@@ -3426,6 +3427,7 @@ subroutine schism_step(it)
                 swild98(2,j,isd0)=swild3(3) !min
               else
                 swild98(1:4,j,isd0)=swild3(1:4) 
+                sdbt(3:2+ntracers,j,isd0)=swild3(5:4+ntracers) 
               endif
 
                   !Check for ics=2 and zonal flow
@@ -3557,6 +3559,7 @@ subroutine schism_step(it)
                 swild98(2,j,isd0)=btlist(ibt)%sclr(3) !min
               else
                 swild98(1:4,j,isd0)=btlist(ibt)%sclr(1:4)
+                sdbt(3:2+ntracers,j,isd0)=btlist(ibt)%sclr(5:4+ntracers)
               endif
 
 !              xyzs(isd0,j,1)=btlist(ibt)%xt; xyzs(isd0,j,2)=btlist(ibt)%yt; xyzs(isd0,j,3)=btlist(ibt)%zt;
@@ -3588,9 +3591,22 @@ subroutine schism_step(it)
 #ifdef INCLUDE_TIMING
       cwtmp=mpi_wtime()
 #endif
-      call exchange_s3d_4(sdbt)
+!      call exchange_s3d_4(sdbt)
       call exchange_s3d_4(swild98)
-!      call exchange_e3dw(webt)
+
+      allocate(swild96(2,nvrt,nsa),stat=istat)
+      if(istat/=0) call parallel_abort('MAIN: fail to allocate swild96 (3.2)')
+      swild96=sdbt(1:2,:,:)
+      call exchange_s3d_2(swild96)
+      sdbt(1:2,:,:)=swild96
+      deallocate(swild96)
+
+      allocate(swild96(ntracers,nvrt,nsa),stat=istat)
+      if(istat/=0) call parallel_abort('MAIN: fail to allocate swild96 (3.3)')
+      swild96=sdbt(3:2+ntracers,:,:)
+      call exchange_s3d_tr2(swild96)
+      sdbt(3:2+ntracers,:,:)=swild96
+      deallocate(swild96)
 
       if(ibtrack_test==1) call exchange_s3dw(tsd)
 
@@ -3600,6 +3616,9 @@ subroutine schism_step(it)
 
 !     ELAD for kriging
       if(inter_mom/=0) then
+        allocate(swild96(2,nvrt,nsa),swild97(2,nvrt,nsa),stat=istat)
+        if(istat/=0) call parallel_abort('MAIN: fail to allocate swild96 (2.2)')
+
 !$OMP parallel default(shared) private(iter,i,ie,k,suru,surv,ll,j,id,kin)
         do iter=1,15 !100
           !Calc epsilon
@@ -3689,6 +3708,7 @@ subroutine schism_step(it)
         enddo !iter
 
 !$OMP end parallel
+        deallocate(swild96,swild97)
       endif !inter_mom; ELAD
 
 !     Debug
@@ -3830,7 +3850,7 @@ subroutine schism_step(it)
       wtmp1=wtmp2
 #endif
 
-      deallocate(swild96,swild97,swild98)
+      deallocate(swild98)
 
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       else