Skip to content

Commit

Permalink
Part 2 of merging with lrdev_updates: hydro part.
Browse files Browse the repository at this point in the history 
TODO: add additional WWM outputs in scribe
  • Loading branch information
@josephzhang8
josephzhang8 committed Apr 27, 2022
1 parent 40baa99 commit 6524e19
Show file tree
Hide file tree
Showing 6 changed files with 181 additions and 64 deletions.
36 changes: 25 additions & 11 deletions sample_inputs/param.nml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,12 +188,18 @@
hmin_radstress = 1. !min. total water depth used only in radiation stress calculation [m]
! nrampwafo = 0 !ramp-up option for the wave forces (1: on; 0: off)
drampwafo = 0. !ramp-up period in days for the wave forces (no ramp-up if <=0)
turbinj = 0.15 !% of depth-induced wave breaking energy injected in turbulence (default: 0.15 (15%), as proposed by Feddersen, 2012)
turbinj = 0.15 !% of depth-induced wave breaking energy injected in turbulence
!(default: 0.15 (15%), as proposed by Feddersen, 2012)
turbinjds = 1.0 !% of wave energy dissipated through whitecapping injected in turbulence
!(default: 1 (100%), as proposed by Paskyabi et al. 2012)
alphaw = 0.5 !for itur=4 : scaling parameter for the surface roughness z0s = alphaw*Hm0.
!If negative z0s = abs(alphaw) e.g. z0s=0.2 m (Feddersen and Trowbridge, 2005)
! Vortex Force terms (off/on:0/1)
fwvor_advxy_stokes = 1 ! --> Stokes drift advection (xy), Coriolis
fwvor_advz_stokes = 1 ! --> Stokes drift advection (z) , Coriolis
fwvor_gradpress = 1 ! --> Pressure term
fwvor_breaking = 1 ! --> Wave breaking
fwvor_streaming = 1 ! --> Wave streaming
cur_wwm = 0 ! Coupling current in WWM
! 0: surface layer current
! 1: depth-averaged current
Expand Down Expand Up @@ -558,7 +564,7 @@
!-----------------------------------------------------------------------
! Sponge layer for elevation and vel.
! If inu_elev=0, no relaxation is applied to elev.
! If inu_elev=1, relax. constants (in 1/sec, i.e. relax*dt<=1) are specified in elev_nudge.gr3
! If inu_elev=1, relax. constants are specified in elev_nudge.gr3
! and applied to eta=0 (thus a depth=0 means no relaxation).
! Similarly for inu_uv (with input uv_nudge.gr3)
!-----------------------------------------------------------------------
Expand Down Expand Up @@ -792,14 +798,13 @@
iof_hydro(22) = 0 !eddy viscosity [m^2/s] {viscosity} 3D
iof_hydro(23) = 0 !turbulent kinetic energy {turbulentKineticEner} 3D
iof_hydro(24) = 0 !turbulent mixing length [m] {mixingLength} 3D
! iof_hydro(25) = 1 !z-coord {zCoordinates} 3D - this flag will be reset by code to be always on
iof_hydro(26) = 1 !horizontal velocity vector [m/s] {horizontalVelX,Y} 3D vector

iof_hydro(27) = 0 !horizontal velocity vector defined @side [m/s] {horizontalSideVelX,Y} 3D vector
iof_hydro(28) = 0 !vertical vel. @elem [m/s] {verticalVelAtElement} 3D
iof_hydro(29) = 0 !T @prism centers [C] {temperatureAtElement} 3D
iof_hydro(30) = 0 !S @prism centers [PSU] {salinityAtElement} 3D
iof_hydro(31) = 0 !Barotropic pressure gradient force vector (m.s-2) @side centers {pressure_gradient} -not implemented yet
! iof_hydro(25) = 1 !z-coord {zCoordinates} 3D - this flag will be reset by code to be always on
iof_hydro(26) = 0 !horizontal vel vector [m/s] {horizontalVelX,Y} 3D vector
iof_hydro(27) = 0 !horizontal vel vector defined @side [m/s] {horizontalSideVelX,Y} 3D vector
iof_hydro(28) = 0 !vertical vel. @elem [m/s] {verticalVelAtElement} 3D
iof_hydro(29) = 0 !T @prism centers [C] {temperatureAtElement} 3D
iof_hydro(30) = 0 !S @prism centers [PSU] {salinityAtElement} 3D
iof_hydro(31) = 0 !Barotropic pressure gradient force vector (m.s-2) @side centers {pressure_gradient} - not implemented yet

!-----------------------------------------------------------------------
! Outputs from optional modules. Only uncomment these if the USE_* is on
Expand Down Expand Up @@ -839,7 +844,16 @@
! iof_wwm(25) = 0 !Charnock coefficient {CharnockCoeff} 2D
! iof_wwm(26) = 0 !Rougness length {rougnessLength} 2D
! iof_wwm(27) = 0 !WWM_energy vector {waveEnergyDirX,Y} 2D vector
! iof_wwm(28) = 0 !Wave force vector (m.s-2) computed by wwm @side centers and whole levels {waveForceX,Y} 3D vector
! iof_wwm(28) = 0 !Wave force vector (m.s-2) computed by wwm @side centers and whole levels {waveForceX,Y} 2D vector

! iof_wwm(29) = 0 !Horizontal Stokes velocity (m.s-1) @nodes and whole levels {stokes_hvel}
! iof_wwm(30) = 0 !Vertical Stokes velocity (m.s-1) @nodes and whole levels {stokes_wvel}
! iof_wwm(31) = 0 !Roller contribution to horizontal Stokes velocity (m.s-1) @nodes and whole levels {roller_stokes_hvel}
! iof_wwm(32) = 0 !Roller energy dissipation rate (W/m²) @nodes {Drol}
! iof_wwm(33) = 0 !Total wave energy dissipation rate by depth-induced breaking (W/m²) @nodes {wave_sbrtot}
! iof_wwm(34) = 0 !Total wave energy dissipation rate by bottom friction (W/m²) @nodes {wave_sbftot}
! iof_wwm(35) = 0 !Total wave energy dissipation rate by whitecapping (W/m²) @nodes {wave_sdstot}
! iof_wwm(36) = 0 !Total wave energy input rate from atmospheric forcing (W/m²) @nodes {wave_sintot}

!-----------------------------------------------------------------------
! Tracer module outputs. In most cases, actual # of outputs depends on # of tracers used
Expand Down
4 changes: 1 addition & 3 deletions sample_inputs/wwminput.nml.WW3
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@
!no wind (0). Try MESIN=1, LSOURCESWAM=F or MESIN=2 and LSOURCESWAM=T
IFRIC = 1 ! Formulation for atmospheric boundary layer, (IFRIC = 1 for MESIN = 1, IFRIC = 4 for MESIN=5); only used if cycle 3 is switched on
MESBF = 1 ! Bottom friction: 1 - JONSWAP (Default); 2 - Madsen et al. (1989); 3 - SHOWEX
FRICC = 0.067 ! if MESBF=1: JONSWAP bottom friction coefficient [0.038,0.067]. If MESBF=2: physical bottom roughness (ignored if given in rough.gr3). If MESBF=3: D50
FRICC = 0.067 ! if MESBF=1: JONSWAP bottom friction coefficient [0.038,0.067]. If MESBF=2: physical bottom roughness (ignored if given in rough.gr3). If MESBF=3: D50 (if negative read from SHOWEX_D50.gr3)
MESBR = 1 ! Shallow water wave breaking; 0: off; 1: on
IBREAK = 1 ! Wave breaking formulation: 1 - Battjes and Janssen (1978)
! 2 - Thornton and Guza (1983) - Constant weighting function
Expand Down Expand Up @@ -207,8 +207,6 @@
! 3: R-K3 (if ICOMP=0 or 1) - slow; 4: Dynamic Splitting (experimental)
! 6: Sub-time steps for breaking term integration (subroutine INT_SHALLOW_SOURCETERMS)

ROLMETHOD = 2 ! Roller source term integration (1 or 2)

DMETHOD = 2
! This switch controls the numerical method in directional space.
! DMETHOD = 0
Expand Down
12 changes: 9 additions & 3 deletions src/Core/schism_glbl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ module schism_glbl
integer,save :: ipre,ipre2,indvel,imm,ihot,ics,iwbl,iharind,nws,impose_net_flux,iwindoff, &
&ibc,ibdef,ihorcon,nstep_wwm,icou_elfe_wwm, &
&fwvor_advxy_stokes,fwvor_advz_stokes,fwvor_gradpress,fwvor_breaking, &
&cur_wwm,wafo_obcramp,iwind_form,irec_nu,itur,ihhat,inu_elev, &
&fwvor_streaming,cur_wwm,wafo_obcramp,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, &
Expand All @@ -103,7 +103,7 @@ module schism_glbl
&vnh1,vnh2,vnf1,vnf2,rnday,btrack_nudge,hmin_man, &
&prmsl_ref,hmin_radstress,eos_a,eos_b,eps1_tvd_imp,eps2_tvd_imp, &
&xlsc0,rearth_pole,rearth_eq,hvis_coef0,disch_coef(10),hw_depth,hw_ratio, &
&slr_rate,rho0,shw,gen_wsett,turbinj,h1_bcc,h2_bcc,vclose_surf_frac, &
&slr_rate,rho0,shw,gen_wsett,turbinj,turbinjds,alphaw,h1_bcc,h2_bcc,vclose_surf_frac, &
&hmin_airsea_ex,hmin_salt_ex,shapiro0,loadtide_coef,h_massconsv,rinflation_icm

! Misc. variables shared between routines
Expand Down Expand Up @@ -585,14 +585,20 @@ module schism_glbl
! WWM
!#ifdef USE_WWM
integer,save :: msc2,mdc2
real(rkind),save,allocatable :: wwave_force(:,:,:), jpress(:), sbr(:,:), sbf(:,:), srol(:,:)
real(rkind),save,allocatable :: wwave_force(:,:,:), jpress(:), sbr(:,:), sbf(:,:), srol(:,:),sds(:,:), eps_w(:), eps_r(:),eps_br(:)
real(rkind),save,allocatable :: stokes_hvel(:,:,:), stokes_wvel(:,:), stokes_hvel_side(:,:,:), stokes_wvel_side(:,:)
real(rkind),save,allocatable :: roller_stokes_hvel(:,:,:), roller_stokes_hvel_side(:,:,:)
!real(rkind),save,allocatable :: stokes_vel(:,:,:), stokes_w_nd(:,:), stokes_vel_sd(:,:,:)
real,save,allocatable :: out_wwm(:,:), out_wwm_rol(:,:), out_wwm_windpar(:,:)
real(rkind),save,allocatable :: tanbeta_x(:), tanbeta_y(:) ! MP from KM: bottom slope
real(rkind),save,allocatable :: curx_wwm(:),cury_wwm(:) !BM: coupling current
real(rkind),save,allocatable :: wafo_opbnd_ramp(:) !BM: ramp on wave forces at open boundary
real(rkind),save,allocatable :: taub_wc(:)
real(rkind),save,allocatable :: delta_wbl(:)
real(rkind),save,allocatable :: wave_sbrtot(:)
real(rkind),save,allocatable :: wave_sbftot(:)
real(rkind),save,allocatable :: wave_sintot(:)
real(rkind),save,allocatable :: wave_sdstot(:)
!#endif

! TIMOR
Expand Down
22 changes: 15 additions & 7 deletions src/Hydro/misc_subs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3864,7 +3864,7 @@ subroutine wbl_GM(taubx,tauby,z0,ubm,wfr,wdir,z0b,fw,delta_wc,icount,iabnormal)
!integer MadsenFlag !0 - Madsen2004, 1 - Madsen79
!Local
real(rkind) :: rkappa,rkn,taub,phi_c,phi_cw,rmu,rmu2,c_mu,tmp,tau_wm, &
&cm_ubm,aa
&cm_ubm,aa,wdir_math

! sanity check
if(z0<0._rkind.or.ubm<0._rkind.or.wfr<0._rkind) then
Expand All @@ -3888,7 +3888,14 @@ subroutine wbl_GM(taubx,tauby,z0,ubm,wfr,wdir,z0b,fw,delta_wc,icount,iabnormal)
! Ubm = Wheight*wr/Sinh(Wnum*Depth) !orbital vel.
taub=sqrt(taubx*taubx+tauby*tauby)
phi_c=atan2(tauby,taubx) !current dir
phi_cw=phi_c+wdir/180._rkind*pi+pi/2._rkind !convert to math convention
!convert to math convention
wdir_math = 180._rkind + 90._rkind - wdir
if (wdir_math .GE. 360.0_rkind) then
wdir_math = MOD (wdir_math, 360.0_rkind)
else if (wdir_math .LT. 0.) then
wdir_math = MOD (wdir_math, 360.0_rkind) + 360.0_rkind
endif
phi_cw=phi_c+wdir_math/180._rkind*pi

rmu=0._rkind !init. guess
c_mu=1._rkind
Expand Down Expand Up @@ -3951,8 +3958,9 @@ end subroutine wbl_GM
! Compute the bottom shear stress due to both wave and current following
! Soulsby (Ch5, Dynamics of Marine Sand, 1997)
! Authors: Kévin Martins, Xavier Bertin, Joseph Zhang
! March 2022, LRU team : correction of a mistake in tau_bot formula
!===============================================================================
subroutine wbl_Soulsby97(Uc_x,Uc_y,z0,sigma,uorb,bthick,Cdp)
subroutine wbl_Soulsby97(Uc_x,Uc_y,z0,sigma,uorb,bthick,Cdp,tau_bot)
! Inputs:
! (Uc_x,Uc_y) - components of the current velocity at the top of the bottom cell;
! z0 - bottom roughness (no waves; m);
Expand All @@ -3969,10 +3977,10 @@ subroutine wbl_Soulsby97(Uc_x,Uc_y,z0,sigma,uorb,bthick,Cdp)
use schism_msgp, only : parallel_abort
implicit none
real(rkind), intent(in) :: Uc_x, Uc_y, z0, sigma, uorb, bthick
real(rkind), intent(inout) :: Cdp
real(rkind), intent(inout) :: Cdp,tau_bot

! Local
real(rkind) :: epsi, Uc, tau_c, tau_w, fw, tau_bot
real(rkind) :: epsi, Uc, tau_c, tau_w, fw

! Some constant
epsi = 0.000001_rkind
Expand All @@ -3991,11 +3999,11 @@ subroutine wbl_Soulsby97(Uc_x,Uc_y,z0,sigma,uorb,bthick,Cdp)
tau_c = Cdp*Uc*Uc ! Norm of the the current-induced shear stress (skin friction) [m^2/s/s]

! Compute wave-induced bottom stress
fw = 1.39_rkind*(sigma*z0/uorb)**0.52_rkind ! Friction factor
fw = min(0.3_rkind,1.39_rkind*(sigma*z0/uorb)**0.52) ! Friction factor
tau_w = 0.5_rkind*fw*uorb*uorb ! Norm of the the wave-induced shear stress

! Compute the combination of both
tau_bot = tau_c*(1._rkind+1.2_rkind*(tau_w/max(epsi,tau_w+tau_c)))
tau_bot = tau_c*(1._rkind+1.2_rkind*(tau_w/max(epsi,tau_w+tau_c))**3.2)

if(Uc==0._rkind) then
!keep original
Expand Down
26 changes: 15 additions & 11 deletions src/Hydro/schism_init.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,8 +197,8 @@ subroutine schism_init(iorder,indir,iths,ntime)
&moitn0,mxitn0,rtol0,iflux,inter_mom,h_bcc1,inu_elev,inu_uv, &
&ihhat,kr_co,rmaxvel,velmin_btrack,btrack_nudge,ibtrack_test,irouse_test, &
&inunfl,shorewafo,ic_elev,nramp_elev,inv_atm_bnd,prmsl_ref,s1_mxnbt,s2_mxnbt, &
&iharind,icou_elfe_wwm,drampwafo,nstep_wwm,hmin_radstress,turbinj, &
&fwvor_advxy_stokes,fwvor_advz_stokes,fwvor_gradpress,fwvor_breaking,wafo_obcramp, &
&iharind,icou_elfe_wwm,drampwafo,nstep_wwm,hmin_radstress,turbinj,turbinjds,alphaw, &
&fwvor_advxy_stokes,fwvor_advz_stokes,fwvor_gradpress,fwvor_breaking,fwvor_streaming,wafo_obcramp, &
&iwbl,cur_wwm,if_source,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, &
&level_age,vclose_surf_frac,iadjust_mass_consv0,ipre2, &
Expand Down Expand Up @@ -431,7 +431,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
!just add the output statements in _step and flags in param.nml (same
!order). Flags for modules other than hydro are only used inside USE_*
if(iorder==0) then
allocate(iof_hydro(40),iof_wwm(30),iof_gen(max(1,ntracer_gen)),iof_age(max(1,ntracer_age)),level_age(ntracer_age/2), &
allocate(iof_hydro(40),iof_wwm(40),iof_gen(max(1,ntracer_gen)),iof_age(max(1,ntracer_age)),level_age(ntracer_age/2), &
&iof_sed(3*sed_class+20),iof_eco(max(1,eco_class)),iof_icm(44),iof_cos(20),iof_fib(5), &
&iof_sed2d(14),iof_ice(10),iof_ana(20),iof_marsh(2),iof_dvd(max(1,ntrs(12))), &
!dim of srqst7 increased to account for 2D elem/side etc
Expand Down Expand Up @@ -470,6 +470,8 @@ subroutine schism_init(iorder,indir,iths,ntime)
inunfl=0; shorewafo=0; ic_elev=0; nramp_elev=0; inv_atm_bnd=0; prmsl_ref=101325._rkind;
s1_mxnbt=0.5_rkind; s2_mxnbt=3.5_rkind;
iharind=0; icou_elfe_wwm=0; drampwafo=0.d0; nstep_wwm=1; hmin_radstress=1._rkind; turbinj=0.15_rkind;
alphaw=1.0_rkind; fwvor_advxy_stokes=1; fwvor_advz_stokes=1;
fwvor_gradpress=1; fwvor_breaking=1; fwvor_streaming=1; wafo_obcramp=0;
fwvor_advxy_stokes=1; fwvor_advz_stokes=1; fwvor_gradpress=1; fwvor_breaking=1; wafo_obcramp=0;
iwbl=0; cur_wwm=0; if_source=0; dramp_ss=2._rkind; ieos_type=0; ieos_pres=0; eos_a=-0.1_rkind; eos_b=1001._rkind;
slr_rate=120._rkind; rho0=1000._rkind; shw=4184._rkind; isav=0; nstep_ice=1; h1_bcc=50._rkind; h2_bcc=100._rkind
Expand Down Expand Up @@ -1483,17 +1485,19 @@ subroutine schism_init(iorder,indir,iths,ntime)
#ifdef USE_WWM
if(iorder==0) then
allocate(out_wwm(npa,35),out_wwm_windpar(npa,10), &
& out_wwm_rol(npa,35), &
& out_wwm_rol(npa,35),taub_wc(npa), &
& stokes_hvel(2,nvrt,npa), stokes_wvel(nvrt,npa),stokes_hvel_side(2,nvrt,nsa), stokes_wvel_side(nvrt,nsa), &
& roller_stokes_hvel(2,nvrt,npa), roller_stokes_hvel_side(2,nvrt,nsa), &
& jpress(npa), sbr(2,npa), sbf(2,npa), srol(2,npa), &
& nne_wwm(np), stat=istat)
& jpress(npa), sbr(2,npa), sbf(2,npa), srol(2,npa), sds(2,npa), &
& nne_wwm(np), eps_w(npa), eps_r(npa), eps_br(npa), delta_wbl(npa), &
& wave_sbrtot(npa),wave_sbftot(npa),wave_sdstot(npa),wave_sintot(npa), stat=istat)
if(istat/=0) call parallel_abort('MAIN: WWM allocation failure')
endif !iorder
out_wwm=0.d0; out_wwm_windpar=0.d0; out_wwm_rol=0.d0
jpress=0.d0; sbr=0.d0; sbf=0.d0; srol=0.d0
out_wwm=0.d0; out_wwm_windpar=0.d0; out_wwm_rol=0.d0; eps_w=0.d0; eps_r=0.d0; eps_br=0.d0
jpress=0.d0; sbr=0.d0; sbf=0.d0; srol=0.d0; sds=0.d0; taub_wc=0.d0
stokes_hvel=0.d0; stokes_wvel=0.d0; stokes_hvel_side=0.d0; stokes_wvel_side=0.d0
roller_stokes_hvel=0.d0; roller_stokes_hvel_side=0.d0
roller_stokes_hvel=0.d0; roller_stokes_hvel_side=0.d0; delta_wbl=0.D0
wave_sbrtot=0.0D0; wave_sbftot=0.0D0; wave_sintot=0.0D0; wave_sdstot=0.0D0
!BM: coupling current for WWM
allocate(curx_wwm(npa),cury_wwm(npa),stat=istat)
if(istat/=0) call parallel_abort('MAIN: (2) WWM alloc failure')
Expand Down Expand Up @@ -5167,7 +5171,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
!... otherwise they will be assigned values below.
if(itur==4) then
#ifdef USE_GOTM
call init_turbulence(8,'gotmturb.inp',nvrt-1) !GOTM starts from level 0
call init_turbulence(8,'gotmturb.nml',nvrt-1) !GOTM starts from level 0
call init_tridiagonal(nvrt-1)
#endif
endif
Expand Down Expand Up @@ -6819,7 +6823,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
call mpi_send(iths,1,itype,nproc_schism-i,110,comm_schism,ierr)
call mpi_send(ntime,1,itype,nproc_schism-i,111,comm_schism,ierr)
call mpi_send(iof_hydro,40,itype,nproc_schism-i,112,comm_schism,ierr)
call mpi_send(iof_wwm,30,itype,nproc_schism-i,113,comm_schism,ierr)
call mpi_send(iof_wwm,40,itype,nproc_schism-i,113,comm_schism,ierr)
!Make sure char len is 20 in scribe_io also!
call mpi_send(out_name,counter_out_name*20,MPI_CHAR,nproc_schism-i,114,comm_schism,ierr)
call mpi_send(ncount_2delem,1,itype,nproc_schism-i,115,comm_schism,ierr)
Expand Down
Loading

0 comments on commit 6524e19

Please sign in to comment.