More work on WW3 coupler. Only issue left is the unit of BHD (wave

pressure).

src/Core/schism_glbl.F90

@@ -610,9 +610,10 @@ module schism_glbl
   real(rkind),save,allocatable :: wave_sintot(:)
   real(rkind),save,allocatable :: wave_sdstot(:)
   real(rkind),save,allocatable :: wave_svegtot(:)
-  real(rkind),save,allocatable :: wave_hs(:),wave_wnm(:),wave_pres(:),wave_stokes_x(:), &
- &wave_stokes_y(:),wave_ocean_flux_x(:),wave_ocean_flux_y(:),wave_flux_friction_x(:), &
- &wave_flux_friction_y(:)
+  real(rkind),save,allocatable :: wave_hs(:),wave_dir(:),wave_tm1(:),wave_wnm(:), &
+ &wave_pres(:),wave_stokes_x(:),wave_stokes_y(:),wave_ocean_flux_x(:), &
+ &wave_ocean_flux_y(:),wave_flux_friction_x(:),wave_flux_friction_y(:), &
+ &wave_orbu(:),wave_orbv(:)
 
 ! TIMOR
 !#ifdef USE_TIMOR

src/Hydro/misc_subs.F90

@@ -6376,33 +6376,40 @@ end subroutine compute_wave_force_lon
 !=========================================================================
 #ifdef USE_WW3
 !     Grab necessary arrays from WW3 and save into SCHISM global arrays
-      SUBROUTINE get_WW3_arrays(WW3__OHS,WW3__WNM,WW3__BHD,WW3_USSX,WW3_USSY, &
-     &WW3_TWOX,WW3_TWOY,WW3_TBBX,WW3_TBBY)
-        USE schism_glbl, ONLY: rkind,errmsg,np,npa,wave_hs,wave_wnm,wave_pres,wave_stokes_x, &
-     &wave_stokes_y,wave_ocean_flux_x,wave_ocean_flux_y,wave_flux_friction_x, &
-     &wave_flux_friction_y
+      SUBROUTINE get_WW3_arrays(WW3__OHS,WW3__DIR,WW3_T0M1,WW3__WNM,WW3__BHD,WW3_USSX,WW3_USSY, &
+     &WW3_TWOX,WW3_TWOY,WW3_TBBX,WW3_TBBY,WW3_UBRX,WW3_UBRY)
+        USE schism_glbl, ONLY: rkind,errmsg,np,npa,wave_hs,wave_dir,wave_tm1, &
+     &wave_wnm,wave_pres,wave_stokes_x,wave_stokes_y,wave_ocean_flux_x, &
+     &wave_ocean_flux_y,wave_flux_friction_x,wave_flux_friction_y, &
+     &wave_orbu,wave_orbv
         USE schism_msgp
         IMPLICIT NONE
 
         !No ghost
-        REAL(rkind),intent(in) :: WW3__OHS(np),WW3__WNM(np),WW3__BHD(np), &
-     &WW3_USSX(np),WW3_USSY(np),WW3_TWOX(np),WW3_TWOY(np),WW3_TBBX(np), &
-     &WW3_TBBY(np)
+        REAL(rkind),intent(in) :: WW3__OHS(np),WW3__DIR(np),WW3_T0M1(np), &
+     &WW3__WNM(np),WW3__BHD(np),WW3_USSX(np),WW3_USSY(np),WW3_TWOX(np), &
+     &WW3_TWOY(np),WW3_TBBX(np),WW3_TBBY(np),WW3_UBRX,WW3_UBRY
 
         REAL(rkind) :: tmp
 
         wave_hs(1:np)=WW3__OHS !Sig wave height [m]
+        wave_dir(1:np)=WW3__DIR !mean wave dir [deg]
+        wave_tm1(1:np)=WW3_T0M1 !mean wave period [s]
         wave_wnm(1:np)=WW3__WNM !mean wave number [1/m]
-        wave_pres(1:np)=WW3__BHD !wave-induced Bernoulli head pressure [N/m]
+        wave_pres(1:np)=WW3__BHD !wave-induced Bernoulli head pressure [N/m or Pa?]
         wave_stokes_x(1:np)=WW3_USSX !Stokes drift [m/s]
         wave_stokes_y(1:np)=WW3_USSY 
         wave_ocean_flux_x(1:np)=WW3_TWOX !wave-ocean mom flux [m2/s2]
         wave_ocean_flux_y(1:np)=WW3_TWOY
         wave_flux_friction_x(1:np)=WW3_TBBX !Momentum flux due to bottom friction [m2/s2]
         wave_flux_friction_y(1:np)=WW3_TBBY
+        wave_orbu(1:np)=WW3_UBRX !near bed orbital vel [m/s]
+        wave_orbv(1:np)=WW3_UBRY 
 
         !Exchange
         call exchange_p2d(wave_hs)
+        call exchange_p2d(wave_dir)
+        call exchange_p2d(wave_tm1)
         call exchange_p2d(wave_wnm)
         call exchange_p2d(wave_pres)
         call exchange_p2d(wave_stokes_x)
@@ -6411,9 +6418,12 @@ SUBROUTINE get_WW3_arrays(WW3__OHS,WW3__WNM,WW3__BHD,WW3_USSX,WW3_USSY, &
         call exchange_p2d(wave_ocean_flux_y)
         call exchange_p2d(wave_flux_friction_x)
         call exchange_p2d(wave_flux_friction_y)
+        call exchange_p2d(wave_orbu)
+        call exchange_p2d(wave_orbv)
 
-        tmp=sum(wave_hs+wave_wnm+wave_pres+wave_stokes_x+wave_stokes_y+ &
-     &wave_ocean_flux_x+wave_ocean_flux_y+wave_flux_friction_x+wave_flux_friction_y)
+        tmp=sum(wave_hs+wave_dir+wave_tm1+wave_wnm+wave_pres+wave_stokes_x+wave_stokes_y+ &
+     &wave_ocean_flux_x+wave_ocean_flux_y+wave_flux_friction_x+wave_flux_friction_y+ &
+     &wave_orbu+wave_orbv)
         if(tmp/=tmp) then
           write(errmsg,*)'WW3 input has nan:',tmp
           call parallel_abort(errmsg)
@@ -6454,7 +6464,7 @@ SUBROUTINE STOKES_STRESS_INTEGRAL_SCHISM
           D_loc = max(znl(nvrt,ip)-znl(kbp(ip),ip),hmin_radstress) !>0
 
           !new40
-          jpress(ip)=wave_pres(ip)/rho0 ![m2/s2]
+          jpress(ip)=wave_pres(ip)/rho0 !needs to be [m2/s2]
 
           k_loc=wave_wnm(ip) !MIN(KDMAX/DEP(IP),WK(IS,IP))
           kD_loc=k_loc*D_loc !MIN(KDMAX,WK(IS,IP)*D_loc)
@@ -6699,7 +6709,7 @@ SUBROUTINE COMPUTE_CONSERVATIVE_VF_TERMS_SCHISM
             STCOR_x_loc = cori(is)*Vst_loc
             STCOR_y_loc = -cori(is)*Ust_loc
 
-            ! Saving wave forces
+            ! Saving wave forces [m/s/s]
             wwave_force(1,k,is)=wwave_force(1,k,is)+VF_x_loc+STCOR_x_loc-dJ_dx_loc
             wwave_force(2,k,is)=wwave_force(2,k,is)+VF_y_loc+STCOR_y_loc-dJ_dy_loc
           END DO !k
@@ -6869,7 +6879,7 @@ SUBROUTINE COMPUTE_STREAMING_VF_TERMS_SCHISM
 
           IF(kbs(is)+1 == nvrt) THEN !2D
             ! N.B. average between the two adjacent nodes
-            !new40: WW3_TBBX [m2/s2] devided by delta_wbl to get m/s/s
+            !new40: WW3_TBBX [m2/s2] devided by delta_wbl or htot to get m/s/s
             Fws_x_loc=-(wave_flux_friction_x(n1)+wave_flux_friction_x(n2))/2.d0/htot !m/s/s
             Fws_y_loc=-(wave_flux_friction_y(n1)+wave_flux_friction_y(n2))/2.d0/htot
             ! Saving wave streaming

src/Hydro/schism_init.F90

@@ -1630,13 +1630,14 @@ subroutine schism_init(iorder,indir,iths,ntime)
 !Additional WW3 arrays
 #ifdef  USE_WW3
       if(iorder==0) then
-        allocate(wave_hs(npa),wave_wnm(npa),wave_pres(npa),wave_stokes_x(npa), &
+        allocate(wave_hs(npa),wave_dir(npa),wave_tm1(npa),wave_wnm(npa),wave_pres(npa),wave_stokes_x(npa), &
      &wave_stokes_y(npa),wave_ocean_flux_x(npa),wave_ocean_flux_y(npa), &
-     &wave_flux_friction_x(npa),wave_flux_friction_y(npa),stat=istat)
+     &wave_flux_friction_x(npa),wave_flux_friction_y(npa),wave_orbu(npa), &
+     &wave_orbv(npa),stat=istat)
         if(istat/=0) call parallel_abort('INIT: alloc WW3')
-        wave_hs=0.d0; wave_wnm=0.d0; wave_pres=0.d0; wave_stokes_x=0.d0
-        wave_stokes_y=0.d0; wave_ocean_flux_x=0.d0; wave_ocean_flux_y=0.d0
-        wave_flux_friction_x=0.d0; wave_flux_friction_y=0.d0
+        wave_hs=0.d0; wave_dir=0.d0; wave_tm1=0.d0; wave_wnm=0.d0; wave_pres=0.d0
+        wave_stokes_x=0.d0; wave_stokes_y=0.d0; wave_ocean_flux_x=0.d0; wave_ocean_flux_y=0.d0
+        wave_flux_friction_x=0.d0; wave_flux_friction_y=0.d0; wave_orbu=0.d0; wave_orbv=0.d0
       endif !iorder=0
 #endif /*USE_WW3*/
 

src/Hydro/schism_step.F90

@@ -2387,13 +2387,28 @@ subroutine schism_step(it)
 !              endif
               !WBL
 #if defined USE_WWM || defined USE_WW3
-              ! Quantities used in both formulations for the WBL
-              ubm = out_wwm(i,22)    ! orbital vel.
+#ifdef USE_WWM
+              ubm = out_wwm(i,22)  !orbital vel. [m/s]
+              aorb = out_wwm(i,23) !orbital amp=U_orb/ <angular freq> [m]
+              wdir = out_wwm(i,18) !wave direction [deg]
               if(out_wwm(i,12)==0.d0) then
                 wfr=0.d0
               else
                 wfr = 2.d0*pi/out_wwm(i,12)  ! angular freq.; out_wwm is not real*8
               endif
+#endif /*USE_WWM*/
+
+#ifdef USE_WW3
+              ubm=sqrt(wave_orbu(i)**2.d0+wave_orbv(i)**2.d0) !orbital vel. [m/s]
+              aorb=ubm*wave_tm1(i)/2.d0/pi !orbital amp [m]
+              wdir=wave_dir(i)
+              if(wave_tm1(i)==0.d0) then
+                wfr=0.d0
+              else
+                wfr=2.d0*pi/wave_tm1(i) !angular freq
+              endif
+#endif /*USE_WW3*/
+
               vmag = sqrt(uu2(kbp(i)+1,i)**2.d0+vv2(kbp(i)+1,i)**2.d0) !current magnitude
 
               ! Wave boundary layer
@@ -2402,7 +2417,6 @@ subroutine schism_step(it)
               else if(iwbl == 1) then ! Grant and Madsen type of WBL
                 taubx = Cdp(i)*vmag*uu2(kbp(i)+1,i)
                 tauby = Cdp(i)*vmag*vv2(kbp(i)+1,i)
-                wdir = out_wwm(i,18) !wave direction
                 !call wbl_GM(taubx,tauby,rough_p(i),ubm,wfr,wdir,z0b,fw,delta_wc,iter,ifl)
                 call wbl_GM(taubx,tauby,rough_p(i),ubm,wfr,wdir,z0b,fw,delta_wbl(i),iter,ifl)
                 !z0b_save(i) = z0b ! (z0b, T. Guérin) 
@@ -2414,7 +2428,6 @@ subroutine schism_step(it)
                 taub_wc(i) = Cdp(i)*vmag*vmag !(uu2(kbp(i)+1,i)**2.d0+vv2(kbp(i)+1,i)**2.d0) 
               else if(iwbl == 2) then! Soulsby (1997) type of WBL
                 call wbl_Soulsby97(uu2(kbp(i)+1,i),vv2(kbp(i)+1,i),rough_p(i),wfr,ubm,bthick,Cdp(i),taub_wc(i))
-                aorb = out_wwm(i,23) ! orbital amp.
                 delta_wbl(i) = 0.09D0*30.D0*rough_p(i)*(aorb/(30.D0*rough_p(i)))**0.82D0
               endif !iwbl