#include "cppdefs.h" MODULE celer_inw_mod #if defined INWAVE_MODEL ! !svn $Id: celer_inw.F 732 2008-09-07 01:55:51Z jcwarner $ !======================================================================! ! ! ! This routine computes the group celerities needed to solve the ! ! action density equations. ! ! ! !======================================================================! ! implicit none PRIVATE PUBLIC :: celer_inw CONTAINS ! !*********************************************************************** SUBROUTINE celer_inw (ng, tile) !*********************************************************************** ! USE mod_param USE mod_grid USE mod_ocean USE mod_stepping USE mod_inwave_vars USE mod_inwave_params USE mod_inwave_bound ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile ! ! Local variable declarations. ! # include "tile.h" ! !# ifdef PROFILE ! CALL wclock_on (ng, iNLM, 35) !# endif CALL celer_inw_tile(ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nrhs(ng), nstp(ng), nnew(ng), & # ifdef MASKING & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif # ifdef WET_DRY & GRID(ng) % rmask_wet, & & GRID(ng) % umask_wet, & & GRID(ng) % vmask_wet, & # endif & GRID(ng) % angler, & & GRID(ng) % pm, & & GRID(ng) % pn, & & WAVEP(ng) % h_tot, & & WAVEP(ng) % u_rho, & & WAVEP(ng) % v_rho, & & WAVEP(ng) % cx, & & WAVEP(ng) % cy, & & WAVEP(ng) % ct, & & WAVEP(ng) % Ta, & & WAVEP(ng) % Tr, & & WAVEP(ng) % kwc, & & WAVEP(ng) % cwc, & & WAVEG(ng) % Inwavecircle, & & WAVEG(ng) % wd) !# ifdef PROFILE ! CALL wclock_off (ng, iNLM, 35) !# endif RETURN END SUBROUTINE celer_inw ! !*********************************************************************** SUBROUTINE celer_inw_tile(ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nrhs, nstp, nnew, & # ifdef MASKING & rmask, umask, vmask, & # endif # ifdef WET_DRY & rmask_wet, umask_wet, vmask_wet, & # endif & angler, pm, pn, & & h_tot,u_rho,v_rho, & & cx, cy, ct, Ta, Tr, & & kwc, cwc, Inwavecircle, wd) !*********************************************************************** ! USE mod_param USE mod_scalars USE mod_inwave_params USE cx3dbc_mod USE cy3dbc_mod USE ct3dbc_mod USE ct3dbc_dir_mod USE Tr3dbc_mod USE exchange_3d_mod USE bc_3d_mod # ifdef DISTRIBUTE USE mp_exchange_mod, ONLY : mp_exchange3d # endif ! ! 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) :: nrhs, nstp, nnew, Inwavecircle # ifdef ASSUMED_SHAPE # ifdef MASKING real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(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(in) :: angler(LBi:,LBj:) real(r8), intent(in) :: pm(LBi:,LBj:) real(r8), intent(in) :: pn(LBi:,LBj:) real(r8), intent(inout) :: cx(LBi:,LBj:,:) real(r8), intent(inout) :: cy(LBi:,LBj:,:) real(r8), intent(inout) :: ct(LBi:,LBj:,:) real(r8), intent(in) :: kwc(LBi:,LBj:,:) real(r8), intent(in) :: cwc(LBi:,LBj:,:) real(r8), intent(in) :: Tr(LBi:,LBj:,:) real(r8), intent(inout) :: Ta(LBi:,LBj:,:) real(r8), intent(in) :: h_tot(LBi:,LBj:) real(r8), intent(in) :: u_rho(LBi:,LBj:) real(r8), intent(in) :: v_rho(LBi:,LBj:) real(r8), intent(in) :: wd(:) # else # ifdef MASKING real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(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(in) :: angler(LBi:,LBj:) real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj) real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: cx(LBi:UBi,LBj:UBj,ND) real(r8), intent(inout) :: cy(LBi:UBi,LBj:UBj,ND) real(r8), intent(inout) :: ct(LBi:UBi,LBj:UBj,ND+1) real(r8), intent(in) :: kwc(LBi:UBi,LBj:UBj,ND) real(r8), intent(in) :: cwc(LBi:UBi,LBj:UBj,ND) real(r8), intent(in) :: Tr(LBi:UBi,LBj:UBj,ND) real(r8), intent(inout) :: Ta(LBi:UBi,LBj:UBj,ND) real(r8), intent(in) :: h_tot(LBi:UBi,LBj:UBj) real(r8), intent(in) :: u_rho(LBi:UBi,LBj:UBj) real(r8), intent(in) :: v_rho(LBi:UBi,LBj:UBj) real(r8), intent(in) :: wd(ND) # endif ! ! Local variable declarations. integer :: i, j, k, d real(r8) :: twopi, otwopi, halfpi real(r8) :: alfa_wave, kh, wa, wr real(r8) :: theta_cur real(r8) :: cff1 real(r8) :: dir_edge real(r8) :: G1, cr, cgr, cgrx, cgry real(r8) :: dudx, dudy, dvdx, dvdy, dhdx, dhdy real(r8) :: cff, cosde, sinde real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: u_dir real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: u_block real(r8), dimension(IminS:ImaxS,JminS:JmaxS,ND) :: cx_rho real(r8), dimension(IminS:ImaxS,JminS:JmaxS,ND) :: cy_rho real(r8), dimension(IminS:ImaxS,JminS:JmaxS,ND) :: ct_r real(r8), parameter :: Tamin=1.0_r8 real(r8), parameter :: Tamax=5000.0_r8 real(r8), parameter :: kDmax = 100.0_r8 ! # include "set_bounds.h" ! twopi=2.0_r8*pi halfpi=0.5_r8*pi otwopi=1.0_r8/twopi DO d=1,ND DO j=Jstr-1,Jend+1 DO i=Istr-1,Iend+1 ! !======================================================================= ! Compute xi and etai components of the wave number !======================================================================= ! alfa_wave =(1.5_r8*pi-wd(d))-angler(i,j) # ifdef DOPPLER ! !======================================================================= ! Compute current component in the direction of waves u_dir !======================================================================= ! IF (u_rho(i,j).eq.0.0_r8) THEN theta_cur=0.5_r8*pi*SIGN(1.0_r8,v_rho(i,j)) ELSE theta_cur=ATAN2(v_rho(i,j),u_rho(i,j)) ENDIF u_dir(i,j)=SQRT(u_rho(i,j)**2+v_rho(i,j)**2)* & & COS(alfa_wave-theta_cur) u_block(i,j) = -(twopi/MAX(Tamin,Tr(i,j,d)))/(kwc(i,j,d)) # endif # ifdef DOPPLER IF((h_tot(i,j).le.Dcrit(ng)) & & .or.((u_dir(i,j).lt.0.0_r8) & & .and.(u_dir(i,j).le.u_block(i,j)))) THEN # else IF(h_tot(i,j).le.Dcrit(ng)) THEN # endif ! !======================================================================! ! If dry or waves are blocked set celerities to zero ! !======================================================================! ! cx_rho(i,j,d)=0.0_r8 cy_rho(i,j,d)=0.0_r8 ELSE ! !=======================================================================! ! Compute the relative group velocities ! !=======================================================================! ! kh=MIN(kwc(i,j,d)*h_tot(i,j),kDmax) cr=sqrt(g/kwc(i,j,d)*tanh(kh)) G1=2.0_r8*kh/sinh(2.0_r8*kh) cgr=0.5_r8*cr*(1.0_r8+G1) ! ! These angles are refered to the local grid ! cgrx=cgr*cos(alfa_wave) cgry=cgr*sin(alfa_wave) ! !======================================================================== ! Compute the absolute group velocities in space direction (Xi and ETAi) !======================================================================== ! # ifdef DOPPLER cx_rho(i,j,d)=cgrx+u_rho(i,j) cy_rho(i,j,d)=cgry+v_rho(i,j) # else cx_rho(i,j,d)=cgrx cy_rho(i,j,d)=cgry # endif ENDIF ENDDO ENDDO # ifdef ACT_ADVECTION DO j=Jstr,Jend DO i=Istr,Iend # ifdef DOPPLER IF((h_tot(i,j).le.Dcrit(ng)) & & .or.((u_dir(i,j).lt.0.0_r8) & & .and.(u_dir(i,j).le.u_block(i,j)))) THEN # else IF(h_tot(i,j).le.Dcrit(ng)) THEN # endif ! !======================================================================! ! If dry or waves are blocked set celerities to zero ! !======================================================================! ! ct_r(i,j,d)=0.0_r8 ELSE ! Note: change in angle convention to solve blah blah blah. dir_edge=-(wd(d)-halfpi) dhdx=(pm(i+1,j)+pm(i,j))*(h_tot(i+1,j)-h_tot(i-1,j)) dhdy=(pn(i,j+1)+pn(i,j))*(h_tot(i,j+1)-h_tot(i,j-1)) kh=MIN(kwc(i,j,d)*h_tot(i,j),kDmax) ct_r(i,j,d)=twopi/(MAX(Tamin,Tr(i,j,d))*sinh(2.0_r8*kh))* & & (sin(dir_edge)*dhdx-cos(dir_edge)*dhdy) # ifdef DOPPLER dudx=pm(i,j)*(u_rho(i+1,j)-u_rho(i,j)) dudy=0.5_r8* & & (0.5_r8*pn(i+1,j)* & & (u_rho(i+1,j+1)-u_rho(i+1,j-1))+ & & 0.5_r8*pn(i,j)* & & (u_rho(i,j+1)-u_rho(i,j-1))) dvdx=0.5_r8* & & (0.5_r8*pm(i,j+1)* & & (v_rho(i+1,j+1)-v_rho(i-1,j+1))+ & & 0.5_r8*pm(i,j)* & & (v_rho(i+1,j)-v_rho(i-1,j))) dvdy=pn(i,j)*(v_rho(i,j+1)-v_rho(i,j)) ! cosde=cos(dir_edge) sinde=sin(dir_edge) cff=cosde* & & (sinde*pm(i,j)*dudx- & & cosde*pn(i,j)*dudy)+ & & sinde* & & (sinde*pm(i,j)*dvdx- & & cosde*pn(i,j)*dvdy) ct_r(i,j,d)=ct_r(i,j,d)+cff # endif ENDIF ENDDO ENDDO # endif DO j=Jstr,Jend DO i=Istr,Iend wr=twopi/MAX(Tamin,Tr(i,j,d)) wa=wr # ifdef DOPPLER wa=wa+kwc(i,j,d)*u_dir(i,j) # endif Ta(i,j,d)=MAX(Tamin,twopi/wa) Ta(i,j,d)=MIN(Ta(i,j,d),Tamax) ENDDO ENDDO ENDDO ! Apply nonperiodic boundary conditions in xi and etai space. ! CALL bc_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, ND, & & Ta) # ifdef DISTRIBUTE CALL mp_exchange3d (ng, tile, iNLM, 1, & & LBi, UBi, LBj, UBj, 1, ND, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & Ta) # endif ! !======================================================================= ! Interpolate cx and cy celerities to cell faces. !======================================================================= ! # ifdef ACX_ADVECTION DO d=1,ND DO j=Jstr,Jend DO i=IstrU,Iend cx(i,j,d)=0.5_r8*(cx_rho(i-1,j,d)+cx_rho(i,j,d)) # ifdef MASKING cx(i,j,d)=cx(i,j,d)*umask(i,j) # endif # ifdef WET_DRY cx(i,j,d)=cx(i,j,d)*umask_wet(i,j) # endif END DO END DO ENDDO # endif # ifdef ACY_ADVECTION DO d=1,ND DO j=JstrV,Jend DO i=Istr,Iend cy(i,j,d)=0.5_r8*(cy_rho(i,j-1,d)+cy_rho(i,j,d)) # ifdef MASKING cy(i,j,d)=cy(i,j,d)*vmask(i,j) # endif # ifdef WET_DRY cy(i,j,d)=cy(i,j,d)*vmask_wet(i,j) # endif END DO END DO ENDDO # endif # ifdef ACT_ADVECTION DO d=2,ND DO j=Jstr,Jend DO i=Istr,Iend ct(i,j,d)=0.5_r8*(ct_r(i,j,d-1)+ ct_r(i,j,d)) # ifdef MASKING ! Apply Land/Sea mask ct(i,j,d)=ct(i,j,d)*rmask(i,j) # endif # ifdef WET_DRY ct(i,j,d)=ct(i,j,d)*rmask_wet(i,j) # endif END DO END DO ENDDO ! DO j=Jstr,Jend DO i=Istr,Iend IF (Inwavecircle.eq.1) THEN ct(i,j,1)=0.5_r8*(ct_r(i,j,ND)+ ct_r(i,j,1)) ELSE ct(i,j,1)=0.0_r8 END IF # ifdef MASKING ! Apply Land/Sea mask ct(i,j,1)=ct(i,j,1)*rmask(i,j) # endif # ifdef WET_DRY ct(i,j,1)=ct(i,j,1)*rmask_wet(i,j) # endif ! ct(i,j,ND+1)=ct(i,j,1) # ifdef MASKING ! Apply Land/Sea mask ct(i,j,ND+1)=ct(i,j,ND+1)*rmask(i,j) # endif # ifdef WET_DRY ct(i,j,ND+1)=ct(i,j,ND+1)*rmask_wet(i,j) # endif END DO END DO # endif ! ! Apply nonperiodic boundary conditions in xi and etai space. ! # ifdef ACX_ADVECTION CALL cx3dbc_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & cx) # endif # ifdef ACY_ADVECTION CALL cy3dbc_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & cy) # endif # ifdef ACT_ADVECTION CALL ct3dbc_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & ct) CALL ct3dbc_dir_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & ct) # endif ! ! Apply periodic boundary conditions. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN # ifdef ACX_ADVECTION CALL exchange_u3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, ND, & & cx) # endif # ifdef ACY_ADVECTION CALL exchange_v3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, ND, & & cy) # endif # ifdef ACT_ADVECTION CALL exchange_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, ND+1, & & ct) # endif END IF # ifdef DISTRIBUTE ! ! Exchange boundary data. ! # ifdef ACX_ADVECTION CALL mp_exchange3d (ng, tile, iNLM, 1, & & LBi, UBi, LBj, UBj, 1, ND, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & cx) # endif # ifdef ACY_ADVECTION CALL mp_exchange3d (ng, tile, iNLM, 1, & & LBi, UBi, LBj, UBj, 1, ND, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & cy) # endif # ifdef ACT_ADVECTION CALL mp_exchange3d (ng, tile, iNLM, 1, & & LBi, UBi, LBj, UBj, 1, ND+1, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & ct) # endif # endif RETURN END SUBROUTINE celer_inw_tile #endif END MODULE celer_inw_mod