MODULE bdytides
!!======================================================================
!! *** MODULE bdytides ***
!! Ocean dynamics: Tidal forcing at open boundaries
!!======================================================================
!! History : 2.0 ! 2007-01 (D.Storkey) Original code
!! 2.3 ! 2008-01 (J.Holt) Add date correction. Origins POLCOMS v6.3 2007
!! 3.0 ! 2008-04 (NEMO team) add in the reference version
!! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes
!! 3.4 ! 2012-09 (G. Reffray and J. Chanut) New inputs + mods
!! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes
!!----------------------------------------------------------------------
!! bdytide_init : read of namelist and initialisation of tidal harmonics data
!! tide_update : calculation of tidal forcing at each timestep
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
USE bdy_oce ! ocean open boundary conditions
USE tide_mod !
USE daymod ! calendar
!
USE in_out_manager ! I/O units
USE iom ! xIO server
USE fldread !
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
IMPLICIT NONE
PRIVATE
PUBLIC bdytide_init ! routine called in bdy_init
PUBLIC bdy_dta_tides ! routine called in dyn_spg_ts
TYPE, PUBLIC :: TIDES_DATA !: Storage for external tidal harmonics data
REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh0 !: Tidal constituents : SSH0 (read in file)
REAL(wp), POINTER, DIMENSION(:,:,:) :: u0, v0 !: Tidal constituents : U0, V0 (read in file)
REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh !: Tidal constituents : SSH (after nodal cor.)
REAL(wp), POINTER, DIMENSION(:,:,:) :: u , v !: Tidal constituents : U , V (after nodal cor.)
END TYPE TIDES_DATA
!$AGRIF_DO_NOT_TREAT
TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data
!$AGRIF_END_DO_NOT_TREAT
TYPE(OBC_DATA) , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s !: bdy external data (slow component)
INTEGER :: kt_tide
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: bdytides.F90 14169 2020-12-14 18:32:36Z jchanut $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE bdytide_init
!!----------------------------------------------------------------------
!! *** SUBROUTINE bdytide_init ***
!!
!! ** Purpose : - Read in namelist for tides and initialise external
!! tidal harmonics data
!!
!!----------------------------------------------------------------------
!! namelist variables
!!-------------------
CHARACTER(len=80) :: filtide ! Filename root for tidal input files
LOGICAL :: ln_bdytide_2ddta ! If true, read 2d harmonic data
!!
INTEGER :: ib_bdy, itide, ib ! dummy loop indices
INTEGER :: ii, ij ! dummy loop indices
INTEGER :: inum, igrd
INTEGER :: isz ! bdy data size
INTEGER :: ios ! Local integer output status for namelist read
INTEGER :: nbdy_rdstart, nbdy_loc
CHARACTER(LEN=50) :: cerrmsg ! error string
CHARACTER(len=80) :: clfile ! full file name for tidal input file
REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read ! work space to read in tidal harmonics data
REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: ztr, zti ! " " " " " " " "
!!
TYPE(TIDES_DATA), POINTER :: td ! local short cut
TYPE( OBC_DATA), POINTER :: dta ! local short cut
!!
NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
!!----------------------------------------------------------------------
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
nbdy_rdstart = 1
DO ib_bdy = 1, nb_bdy
IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
!
td => tides(ib_bdy)
dta => dta_bdy(ib_bdy)
! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
filtide(:) = ''
READ ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
!
! Need to support possibility of reading more than one
! nambdy_tide from the namelist_cfg internal file.
! Do this by finding the ib_bdy'th occurence of nambdy_tide in the
! character buffer as the starting point.
!
nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_tide' )
IF( nbdy_loc .GT. 0 ) THEN
nbdy_rdstart = nbdy_rdstart + nbdy_loc
ELSE
WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',ib_bdy,' of nambdy_tide not found'
ios = -1
CALL ctl_nam ( ios , cerrmsg )
ENDIF
READ ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_tide, IOSTAT = ios, ERR = 902)
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
IF(lwm) WRITE ( numond, nambdy_tide )
! ! Parameter control and print
IF(lwp) WRITE(numout,*) ' '
IF(lwp) WRITE(numout,*) ' Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
IF(lwp) WRITE(numout,*) ' read tidal data in 2d files: ', ln_bdytide_2ddta
IF(lwp) WRITE(numout,*) ' Number of tidal components to read: ', nb_harmo
IF(lwp) THEN
WRITE(numout,*) ' Tidal components: '
DO itide = 1, nb_harmo
WRITE(numout,*) ' ', tide_harmonics(itide)%cname_tide
END DO
ENDIF
IF(lwp) WRITE(numout,*) ' '
! Allocate space for tidal harmonics data - get size from BDY data arrays
! Allocate also slow varying data in the case of time splitting:
! Do it anyway because at this stage knowledge of free surface scheme is unknown
! -----------------------------------------------------------------------
IF( ASSOCIATED(dta%ssh) ) THEN ! we use bdy ssh on this mpi subdomain
isz = SIZE(dta%ssh)
ALLOCATE( td%ssh0( isz, nb_harmo, 2 ), td%ssh( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%ssh( isz ) )
dta_bdy_s(ib_bdy)%ssh(:) = 0._wp ! needed?
ENDIF
IF( ASSOCIATED(dta%u2d) ) THEN ! we use bdy u2d on this mpi subdomain
isz = SIZE(dta%u2d)
ALLOCATE( td%u0 ( isz, nb_harmo, 2 ), td%u ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%u2d( isz ) )
dta_bdy_s(ib_bdy)%u2d(:) = 0._wp ! needed?
ENDIF
IF( ASSOCIATED(dta%v2d) ) THEN ! we use bdy v2d on this mpi subdomain
isz = SIZE(dta%v2d)
ALLOCATE( td%v0 ( isz, nb_harmo, 2 ), td%v ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%v2d( isz ) )
dta_bdy_s(ib_bdy)%v2d(:) = 0._wp ! needed?
ENDIF
! fill td%ssh0, td%u0, td%v0
! -----------------------------------------------------------------------
IF( ln_bdytide_2ddta ) THEN
!
! It is assumed that each data file contains all complex harmonic amplitudes
! given on the global domain (ie global, jpiglo x jpjglo)
!
ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
!
! SSH fields
clfile = TRIM(filtide)//'_grid_T.nc'
CALL iom_open( clfile , inum )
igrd = 1 ! Everything is at T-points here
DO itide = 1, nb_harmo
CALL iom_get( inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
CALL iom_get( inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) )
IF( ASSOCIATED(dta%ssh) ) THEN ! we use bdy ssh on this mpi subdomain
DO ib = 1, SIZE(dta%ssh)
ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
td%ssh0(ib,itide,1) = ztr(ii,ij)
td%ssh0(ib,itide,2) = zti(ii,ij)
END DO
ENDIF
END DO
CALL iom_close( inum )
!
! U fields
clfile = TRIM(filtide)//'_grid_U.nc'
CALL iom_open( clfile , inum )
igrd = 2 ! Everything is at U-points here
DO itide = 1, nb_harmo
CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:),cd_type='U',psgn=-1._wp)
CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:),cd_type='U',psgn=-1._wp)
IF( ASSOCIATED(dta%u2d) ) THEN ! we use bdy u2d on this mpi subdomain
DO ib = 1, SIZE(dta%u2d)
ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
td%u0(ib,itide,1) = ztr(ii,ij)
td%u0(ib,itide,2) = zti(ii,ij)
END DO
ENDIF
END DO
CALL iom_close( inum )
!
! V fields
clfile = TRIM(filtide)//'_grid_V.nc'
CALL iom_open( clfile , inum )
igrd = 3 ! Everything is at V-points here
DO itide = 1, nb_harmo
CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:),cd_type='V',psgn=-1._wp)
CALL iom_get(inum, jpdom_auto, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:),cd_type='V',psgn=-1._wp)
IF( ASSOCIATED(dta%v2d) ) THEN ! we use bdy v2d on this mpi subdomain
DO ib = 1, SIZE(dta%v2d)
ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
td%v0(ib,itide,1) = ztr(ii,ij)
td%v0(ib,itide,2) = zti(ii,ij)
END DO
ENDIF
END DO
CALL iom_close( inum )
!
DEALLOCATE( ztr, zti )
!
ELSE
!
! Read tidal data only on bdy segments
!
ALLOCATE( dta_read( MAXVAL( idx_bdy(ib_bdy)%nblen(:) ), 1, 1 ) )
!
! Open files and read in tidal forcing data
! -----------------------------------------
DO itide = 1, nb_harmo
! ! SSH fields
IF( ASSOCIATED(dta%ssh) ) THEN ! we use bdy ssh on this mpi subdomain
isz = SIZE(dta%ssh)
clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
CALL iom_open( clfile, inum )
CALL fld_map( inum, 'z1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
td%ssh0(:,itide,1) = dta_read(1:isz,1,1)
CALL fld_map( inum, 'z2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
td%ssh0(:,itide,2) = dta_read(1:isz,1,1)
CALL iom_close( inum )
ENDIF
! ! U fields
IF( ASSOCIATED(dta%u2d) ) THEN ! we use bdy u2d on this mpi subdomain
isz = SIZE(dta%u2d)
clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
CALL iom_open( clfile, inum )
CALL fld_map( inum, 'u1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
td%u0(:,itide,1) = dta_read(1:isz,1,1)
CALL fld_map( inum, 'u2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
td%u0(:,itide,2) = dta_read(1:isz,1,1)
CALL iom_close( inum )
ENDIF
! ! V fields
IF( ASSOCIATED(dta%v2d) ) THEN ! we use bdy v2d on this mpi subdomain
isz = SIZE(dta%v2d)
clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
CALL iom_open( clfile, inum )
CALL fld_map( inum, 'v1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
td%v0(:,itide,1) = dta_read(1:isz,1,1)
CALL fld_map( inum, 'v2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
td%v0(:,itide,2) = dta_read(1:isz,1,1)
CALL iom_close( inum )
ENDIF
!
END DO ! end loop on tidal components
!
DEALLOCATE( dta_read )
!
ENDIF ! ln_bdytide_2ddta=.true.
!
ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
!
END DO ! loop on ib_bdy
!
END SUBROUTINE bdytide_init
SUBROUTINE bdy_dta_tides( kt, kit, pt_offset )
!!----------------------------------------------------------------------
!! *** SUBROUTINE bdy_dta_tides ***
!!
!! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
!!
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! Main timestep counter
INTEGER, OPTIONAL, INTENT(in) :: kit ! Barotropic timestep counter (for timesplitting option)
REAL(wp),OPTIONAL, INTENT(in) :: pt_offset ! time offset in units of timesteps
!
LOGICAL :: lk_first_btstp ! =.TRUE. if time splitting and first barotropic step
INTEGER :: itide, ib_bdy, ib ! loop indices
REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset
!!----------------------------------------------------------------------
!
lk_first_btstp=.TRUE.
IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
zt_offset = 0._wp
IF( PRESENT(pt_offset) ) zt_offset = pt_offset
! Absolute time from model initialization:
IF( PRESENT(kit) ) THEN
z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_e, wp) ) * rn_Dt
ELSE
z_arg = ( REAL(kt, wp) + zt_offset ) * rn_Dt
ENDIF
! Linear ramp on tidal component at open boundaries
zramp = 1.
IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - REAL(nit000,wp)*rn_Dt)/(rn_tide_ramp_dt*rday),0.),1.)
DO ib_bdy = 1,nb_bdy
!
IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
!
! We refresh nodal factors every day below
! This should be done somewhere else
IF ( ( nsec_day == NINT(0.5_wp * rn_Dt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
!
kt_tide = kt - NINT((REAL(nsec_day,wp) - 0.5_wp * rn_Dt)/rn_Dt)
!
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
WRITE(numout,*) '~~~~~~~~~~~~~~ '
ENDIF
!
CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
!
ENDIF
zoff = REAL(-kt_tide,wp) * rn_Dt ! time offset relative to nodal factor computation time
!
! If time splitting, initialize arrays from slow varying open boundary data:
IF ( PRESENT(kit) ) THEN
IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) dta_bdy(ib_bdy)%ssh(:) = dta_bdy_s(ib_bdy)%ssh(:)
IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) dta_bdy(ib_bdy)%u2d(:) = dta_bdy_s(ib_bdy)%u2d(:)
IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) dta_bdy(ib_bdy)%v2d(:) = dta_bdy_s(ib_bdy)%v2d(:)
ENDIF
!
! Update open boundary data arrays:
DO itide = 1, nb_harmo
!
z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega
z_cost = zramp * COS( z_sarg )
z_sist = zramp * SIN( z_sarg )
!
IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) THEN ! SSH on tracer grid
DO ib = 1, SIZE(dta_bdy(ib_bdy)%ssh)
dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
& ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
& tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
END DO
ENDIF
!
IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) THEN ! U grid
DO ib = 1, SIZE(dta_bdy(ib_bdy)%u2d)
dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
& ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
& tides(ib_bdy)%u(ib,itide,2)*z_sist )
END DO
ENDIF
!
IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) THEN ! V grid
DO ib = 1, SIZE(dta_bdy(ib_bdy)%v2d)
dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
& ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
& tides(ib_bdy)%v(ib,itide,2)*z_sist )
END DO
ENDIF
!
END DO
ENDIF
END DO
!
END SUBROUTINE bdy_dta_tides
SUBROUTINE tide_init_elevation( idx, td )
!!----------------------------------------------------------------------
!! *** ROUTINE tide_init_elevation ***
!!----------------------------------------------------------------------
TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data
!
INTEGER :: itide, isz, ib ! dummy loop indices
REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide
!!----------------------------------------------------------------------
!
IF( ASSOCIATED(td%ssh0) ) THEN ! SSH on tracer grid.
!
isz = SIZE( td%ssh0, dim = 1 )
ALLOCATE( mod_tide(isz), phi_tide(isz) )
!
DO itide = 1, nb_harmo
DO ib = 1, isz
mod_tide(ib)=SQRT( td%ssh0(ib,itide,1)*td%ssh0(ib,itide,1) + td%ssh0(ib,itide,2)*td%ssh0(ib,itide,2) )
phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
END DO
DO ib = 1, isz
mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
END DO
DO ib = 1, isz
td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
END DO
END DO
!
DEALLOCATE( mod_tide, phi_tide )
!
ENDIF
!
END SUBROUTINE tide_init_elevation
SUBROUTINE tide_init_velocities( idx, td )
!!----------------------------------------------------------------------
!! *** ROUTINE tide_init_elevation ***
!!----------------------------------------------------------------------
TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data
!
INTEGER :: itide, isz, ib ! dummy loop indices
REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide
!!----------------------------------------------------------------------
!
IF( ASSOCIATED(td%u0) ) THEN ! U grid. we use bdy u2d on this mpi subdomain
!
isz = SIZE( td%u0, dim = 1 )
ALLOCATE( mod_tide(isz), phi_tide(isz) )
!
DO itide = 1, nb_harmo
DO ib = 1, isz
mod_tide(ib)=SQRT( td%u0(ib,itide,1)*td%u0(ib,itide,1) + td%u0(ib,itide,2)*td%u0(ib,itide,2) )
phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
END DO
DO ib = 1, isz
mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
END DO
DO ib = 1, isz
td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
END DO
END DO
!
DEALLOCATE( mod_tide, phi_tide )
!
ENDIF
!
IF( ASSOCIATED(td%v0) ) THEN ! V grid. we use bdy u2d on this mpi subdomain
!
isz = SIZE( td%v0, dim = 1 )
ALLOCATE( mod_tide(isz), phi_tide(isz) )
!
DO itide = 1, nb_harmo
DO ib = 1, isz
mod_tide(ib)=SQRT( td%v0(ib,itide,1)*td%v0(ib,itide,1) + td%v0(ib,itide,2)*td%v0(ib,itide,2) )
phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
END DO
DO ib = 1, isz
mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
END DO
DO ib = 1, isz
td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
END DO
END DO
!
DEALLOCATE( mod_tide, phi_tide )
!
ENDIF
!
END SUBROUTINE tide_init_velocities
!!======================================================================
END MODULE bdytides