#include "cppdefs.h" MODULE dissip_inw_mod #if defined INWAVE_MODEL # if defined WDISS_ROELVINK || defined WDISS_GAMMA ! !======================================================================= ! ! ! This routine computes the energy dissipation ! ! ! !======================================================================= ! implicit none PRIVATE PUBLIC :: dissip_inw_tile CONTAINS ! !*********************************************************************** SUBROUTINE dissip_inw (ng, tile, nout) !*********************************************************************** ! USE mod_param USE mod_grid USE mod_inwave_params USE mod_inwave_vars USE mod_ocean USE mod_coupling USE mod_stepping USE mod_forces ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, nout ! ! Local variable declarations. ! # include "tile.h" ! !# ifdef PROFILE ! CALL wclock_on (ng, iNLM, 35) !# endif CALL dissip_inw_tile(ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nout, & # ifdef MASKING & GRID(ng) % rmask, & # endif # ifdef WET_DRY & GRID(ng) % rmask_wet, & & GRID(ng) % umask_wet, & & GRID(ng) % vmask_wet, & # endif & FORCES(ng) % Dissip_break, & & FORCES(ng) % Dissip_wcap, & & WAVEP(ng) % h_tot, & & WAVEP(ng) % AC, WAVEP(ng) % Tr, WAVEP(ng) % kwc) !# ifdef PROFILE ! CALL wclock_off (ng, iNLM, 35) !# endif RETURN END SUBROUTINE dissip_inw ! !*********************************************************************** SUBROUTINE dissip_inw_tile(ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nout, & # ifdef MASKING & rmask, & # endif # ifdef WET_DRY & rmask_wet, umask_wet, vmask_wet, & # endif & Dissip_break, & & Dissip_wcap, & & h_tot, & & AC, Tr, kwc) !*********************************************************************** ! USE mod_param USE mod_inwave_params USE mod_boundary USE mod_grid USE mod_scalars USE exchange_2d_mod # ifdef DISTRIBUTE USE mp_exchange_mod, ONLY : mp_exchange2d # endif USE bc_2d_mod ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nout ! # ifdef ASSUMED_SHAPE # ifdef MASKING real(r8), intent(in) :: rmask(LBi:,LBj:) # endif # ifdef WET_DRY real(r8), intent(in) :: rmask_wet(LBi:,LBj:) real(r8), intent(in) :: umask_wet(LBi:,LBj:) real(r8), intent(in) :: vmask_wet(LBi:,LBj:) # endif real(r8), intent(inout) :: Dissip_break(LBi:,LBj:) real(r8), intent(inout) :: Dissip_wcap(LBi:,LBj:) real(r8), intent(in) :: h_tot(LBi:,LBj:) real(r8), intent(inout) :: AC(LBi:,LBj:,:,:) real(r8), intent(in) :: Tr(LBi:,LBj:,:) real(r8), intent(in) :: kwc(LBi:,LBj:,:) # else # ifdef MASKING real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) # endif # ifdef WET_DRY real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj) # endif real(r8), intent(inout) :: Dissip_break(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: Dissip_wcap(LBi:UBi,LBj:UBj) real(r8), intent(in) :: h_tot(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: AC(LBi:UBi,LBj:UBj,ND,3) real(r8), intent(in) :: Tr(LBi:UBi,LBj:UBj,ND) real(r8), intent(in) :: kwc(LBi:UBi,LBj:UBj,ND) # endif ! ! Local variable declarations. ! integer :: i, j, d real(r8) :: EW, oEW, TRM, H, Qb, Hmax_r, diff, Emax_r real(r8) :: twopi, otwopi, ogrho0, cff real(r8), parameter :: Trmin=1.0_r8 real(r8), parameter :: alfa=1.0_r8 real(r8), parameter :: break=0.45_r8 real(r8), parameter :: n_r=15.0_r8 real(r8), parameter :: EWlim=0.00001_r8 real(r8):: EN(ND) # include "set_bounds.h" twopi=2.0_r8*pi otwopi=1.0_r8/twopi ogrho0=1.0_r8/(g*rho0) ! DO j=Jstr,Jend DO i=Istr,Iend EW=0.0_r8 # ifdef WDISS_ROELVINK TRM=0.0_r8 # endif DO d=1,ND !======================================================================= ! Compute the energy from action balance and wave heigth !======================================================================= EN(d)=AC(i,j,d,nout)*twopi/(MAX(Trmin,Tr(i,j,d))) !======================================================================= ! Compute the total energy !======================================================================= EW=EW+EN(d) !======================================================================= ! Compute the mean wave number and intrinsic periods ! What we do is give more importance to those wave ! numbers with more energy !======================================================================= # ifdef WDISS_ROELVINK TRM=TRM+Tr(i,j,d)*EN(d) # endif ENDDO # ifdef WDISS_ROELVINK cff=1.0_r8/(max(EW,EWlim)) TRM=TRM*cff # endif ! EW=MAX(EW,EWlim) EW=MAX(EW,0.0_r8) !this was needed # ifdef WDISS_ROELVINK !======================================================================= ! Compute the wave height. This is based on Hrms. !======================================================================= H=(8.0_r8*EW*ogrho0)**0.5_r8 # endif !======================================================================= ! Compute the energy dissipation !======================================================================= IF (h_tot(i,j).ge.Dcrit(ng)) THEN Hmax_r=break*(MAX(h_tot(i,j),0.0_r8)) # ifdef WDISS_ROELVINK Qb=MIN(1.0_r8,1.0_r8-EXP(-(H/Hmax_r)**n_r)) IF (TRM.gt.0.0001_r8) THEN Dissip_break(i,j)=2.0_r8*alfa/TRM*EW*Qb*dt(ng) ELSE Dissip_break(i,j)=0.0_r8 END IF # elif defined WDISS_GAMMA Emax_r=0.125_r8*g*rho0*Hmax_r**2.0_r8 diff=EW-Emax_r Dissip_break(i,j)=MAX(0.0_r8,diff) # endif ELSE Dissip_break(i,j)=0.0_r8 END IF # ifdef MASKING Dissip_break(i,j)=Dissip_break(i,j)*rmask(i,j) # endif # ifdef WET_DRY Dissip_break(i,j)=Dissip_break(i,j)*rmask_wet(i,j) # endif !======================================================================= ! Distribute dissipation over directions and recompute Ac !======================================================================= oEW=1.0_r8/MAX(EW,EWlim) DO d=1,ND IF ((h_tot(i,j).ge.Dcrit(ng)).and.(EW.gt.EWlim)) THEN EN(d)=MAX(0.0_r8,EN(d)-Dissip_break(i,j)*EN(d)*oEW) AC(i,j,d,nout)=EN(d)*Tr(i,j,d)*otwopi ELSE AC(i,j,d,nout)=0.0_r8 ENDIF # ifdef MASKING AC(i,j,d,nout)=AC(i,j,d,nout)*rmask(i,j) # endif # ifdef WET_DRY AC(i,j,d,nout)=AC(i,j,d,nout)*rmask_wet(i,j) # endif ENDDO Dissip_wcap(i,j)=0.0_r8 Dissip_break(i,j)=Dissip_break(i,j)/(dt(ng)*rho0) ENDDO ENDDO ! ! Apply boundary conditions. ! CALL bc_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & Dissip_break) CALL bc_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & Dissip_wcap) # ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, iNLM, 2, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & Dissip_break, Dissip_wcap) # endif RETURN END SUBROUTINE dissip_inw_tile # endif #endif END MODULE dissip_inw_mod