Add SAL

Author: josephzhang8

sample_inputs/param.nml

@@ -58,6 +58,16 @@
 !-----------------------------------------------------------------------
   itransport_only = 0
 
+!-----------------------------------------------------------------------
+! Option to add self-attrection loading tide (SAL) into tidal potential 
+! (usually for global model). If /=1, needs inputs: loadtide_[FREQ].gr3, 
+! where [FREQ] are freq names (shared with tidal potential) and the _two_
+! 'depths' inside are amplitude (m) and phases (degrees behind GMT), 
+! interpolated from global tide model (e.g. FES2014); 
+! also shares some parameters with tidal potential in bctides.in (cut-off depth, frequencies).
+!-----------------------------------------------------------------------
+  iloadtide = 0
+
 ! Starting time
   start_year = 2000 !int
   start_month = 1 !int

src/Core/schism_glbl.F90

@@ -92,7 +92,7 @@ module schism_glbl
                   &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,ielm_transport,max_subcyc, &
-                  &itransport_only,meth_sink
+                  &itransport_only,meth_sink,iloadtide
   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
@@ -362,7 +362,7 @@ module schism_glbl
                              &jspc(:)          
   integer,save,allocatable :: isblock_nd(:,:),isblock_el(:),iq_block_lcl(:),iq_block(:)
   real(rkind),save,allocatable :: trobc(:,:),vobc1(:),vobc2(:),tamp(:), &
-                                  &tnf(:),tfreq(:),tear(:),amig(:),ff(:),face(:), &
+                                  &tnf(:),tfreq(:),tear(:),amig(:),ff(:),face(:),rloadtide(:,:,:), &
                                   &emo(:,:,:),efa(:,:,:),umo(:,:,:),ufa(:,:,:), &
                                   &vmo(:,:,:),vfa(:,:,:),eth(:,:), &
                                   &qthcon(:),uthnd(:,:,:),vthnd(:,:,:), &

src/Hydro/schism_init.F90

@@ -113,6 +113,8 @@ subroutine schism_init(iorder,indir,iths,ntime)
       character(len=10) :: timestamp
 !      character(len=72) :: it_char
       character(len=72) :: fgb  ! Processor specific global output file name
+      character(len=6) :: char6
+      character(len=6), allocatable :: tp_name(:)
       integer :: lfgb       ! Length of processor specific global output file name
       real(4) :: floatout,floatout2
       real(rkind) :: double1 !for hotstart.in
@@ -184,7 +186,8 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &iwbl,cur_wwm,if_source,nramp_ss,dramp_ss,ieos_type,ieos_pres,eos_a,eos_b,slr_rate, &
      &rho0,shw,isav,nstep_ice,iunder_deep,h1_bcc,h2_bcc,hw_depth,hw_ratio, &
      &ibtrack_openbnd,level_age,vclose_surf_frac,iadjust_mass_consv0,ipre2, &
-     &ielm_transport,max_subcyc,i_hmin_airsea_ex,hmin_airsea_ex,itransport_only,meth_sink
+     &ielm_transport,max_subcyc,i_hmin_airsea_ex,hmin_airsea_ex,itransport_only,meth_sink, &
+     &iloadtide
 
      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, &
@@ -464,6 +467,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
       ielm_transport=0; max_subcyc=10
       hmin_airsea_ex=0.2_rkind
       itransport_only=0; meth_sink=0
+      iloadtide=0
 
       !Output elev, hvel by detault
       iof_hydro=0; iof_wwm=0; iof_gen=0; iof_age=0; iof_sed=0; iof_eco=0; iof_dvd=0
@@ -2466,8 +2470,12 @@ subroutine schism_init(iorder,indir,iths,ntime)
       read(31,*) ntip,tip_dp !cut-off depth for applying tidal potential
       if(ntip>0) then
         if(iorder==0) then
-          allocate(tamp(ntip),tnf(ntip),tfreq(ntip),jspc(ntip),tear(ntip),stat=istat)
-          if(istat/=0) call parallel_abort('MAIN: allocation failure for tamp etc')
+          allocate(tp_name(ntip),tamp(ntip),tnf(ntip),tfreq(ntip),jspc(ntip),tear(ntip),stat=istat)
+          if(istat/=0) call parallel_abort('INIT: allocation failure for tamp etc')
+          if(iloadtide/=0) then !loading tide (SAL)
+            allocate(rloadtide(2,ntip,npa),stat=istat)
+            if(istat/=0) call parallel_abort('INIT: alloc failure for SAL')
+          endif !iloadtide/
         endif !iorder
 !'
         open(32,file=in_dir(1:len_in_dir)//'hgrid.ll',status='old')
@@ -2489,14 +2497,36 @@ subroutine schism_init(iorder,indir,iths,ntime)
         lreadll=.true.
 
         do i=1,ntip
-          read(31,*) !tag
+          read(31,'(a6)')tp_name(i) !tag
           read(31,*) jspc(i),tamp(i),tfreq(i),tnf(i),tear(i)
           if(jspc(i)<0.or.jspc(i)>2) then
             write(errmsg,*)'Illegal tidal species #',jspc(i)
             call parallel_abort(errmsg)
           endif
           tear(i)=tear(i)*pi/180.d0
         enddo !i
+
+        if(iloadtide/=0) then !loading tide
+          do i=1,ntip
+            char6=adjustl(tp_name(i))
+            itmp=len_trim(char6)
+!Debug
+!            write(12,*)'SAL grid name:',in_dir(1:len_in_dir)//'loadtide_'//char6(1:itmp)//'.gr3'
+ 
+            !.gr3 has both amp, phase
+            open(32,file=in_dir(1:len_in_dir)//'loadtide_'//char6(1:itmp)//'.gr3',status='old')
+            read(32,*); read(32,*)
+            do j=1,np_global
+              read(32,*)k,xtmp,ytmp,tmp1,tmp2
+              if(ipgl(j)%rank==myrank) then
+                ii=ipgl(j)%id
+                rloadtide(1,i,ii)=tmp1 !amp [m]
+                rloadtide(2,i,ii)=tmp2*pi/180.d0 !phase [radian]
+              endif
+            enddo !j
+            close(32)
+          enddo !i
+        endif !iloadtide
       endif !ntip>0
 
 !...  Boundary forcing freqs.

src/Hydro/schism_step.F90

@@ -410,15 +410,19 @@ subroutine schism_step(it)
       enddo !i
 !$OMP end do
 
-!...  Earth tidal potential at nodes: pre-compute to save time
-!...
+!...  Earth tidal potential and loading tide at nodes: pre-compute to save time
+!... 
 !$OMP do
       do i=1,npa
         etp(i)=0.d0
         do j=1,ntip
           ncyc=int(tfreq(j)*time/2.d0/pi)
           arg=tfreq(j)*time-real(ncyc,rkind)*2.d0*pi+jspc(j)*xlon(i)+tear(j)
           etp(i)=etp(i)+0.69d0*ramp*tamp(j)*tnf(j)*fun_lat(jspc(j),i)*cos(arg)
+
+          if(iloadtide/=0) then !loading tide
+            etp(i)=etp(i)+rloadtide(1,j,i)*cos(tfreq(j)*time-rloadtide(2,j,i))
+          endif !iloadtide/
         enddo !j
       enddo !i
 !$OMP end do