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 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge !! 3.4 ! 2013 (J. Harle) rewite to used tide_mod for phase and nodal !! corrections every day !!---------------------------------------------------------------------- #if defined key_bdy !!---------------------------------------------------------------------- !! 'key_bdy' Open Boundary Condition !!---------------------------------------------------------------------- !! PUBLIC !! tide_init : read of namelist and initialisation of tidal harmonics data !! tide_update : calculation of tidal forcing at each timestep !!---------------------------------------------------------------------- USE timing ! Timing USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE iom USE in_out_manager ! I/O units USE phycst ! physical constants USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE bdy_par ! Unstructured boundary parameters USE bdy_oce ! ocean open boundary conditions USE fldread, ONLY: fld_map USE daymod ! calendar USE tide_mod USE ioipsl, ONLY : ymds2ju ! for calendar IMPLICIT NONE PRIVATE PUBLIC tide_init ! routine called in nemo_init PUBLIC tide_update ! routine called in bdydyn TYPE, PUBLIC :: TIDES_DATA !: Storage for external tidal harmonics data INTEGER :: ncpt !: Actual number of tidal components REAL(wp), POINTER, DIMENSION(:) :: speed !: Phase speed of tidal constituent (deg/hr) REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh !: Tidal constituents : SSH REAL(wp), POINTER, DIMENSION(:,:,:) :: u !: Tidal constituents : U REAL(wp), POINTER, DIMENSION(:,:,:) :: v !: Tidal constituents : V REAL(wp), POINTER, DIMENSION(:,:,:) :: sshr !: Tidal constituents : SSH (reference) REAL(wp), POINTER, DIMENSION(:,:,:) :: ur !: Tidal constituents : U (reference) REAL(wp), POINTER, DIMENSION(:,:,:) :: vr !: Tidal constituents : V (reference) END TYPE TIDES_DATA TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data INTEGER, ALLOCATABLE, DIMENSION(:) :: bdy_ntide REAL(wp), ALLOCATABLE, DIMENSION(:) :: bdy_omega_tide REAL(wp), ALLOCATABLE, DIMENSION(:) :: bdy_v0tide, & bdy_blank, & bdy_utide, & bdy_ftide, & rbdy_ftide LOGICAL :: ln_tide_date !: =T correct tide phases and amplitude for model start date LOGICAL :: ln_tide_v0 !: =T correct tide phases and amplitude for model start date INTEGER :: nn_tide_date !: yyyymmdd reference date of tidal data INTEGER :: bdy_nn_tide INTEGER :: bdy_kt_tide ! Main tide timestep counter INTEGER :: bdy_tide_offset ! Main tide timestep counter !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tide_init !!---------------------------------------------------------------------- !! *** SUBROUTINE tide_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 CHARACTER(len= 4), DIMENSION(jpmax_harmo) :: tide_cpt !: Names of tidal components used. !! INTEGER :: ib_bdy, itide, ib, ji !: dummy loop indices INTEGER :: inum, igrd INTEGER :: lcl_ryear, lcl_rmonth, lcl_rday INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) CHARACTER(len=80) :: clfile !: full file name for tidal input file REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t, fdayn, fdayr REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data !! TYPE(TIDES_DATA), POINTER :: td !: local short cut !! NAMELIST/nambdy_tide/filtide, tide_cpt, ln_tide_date, nn_tide_date, ln_tide_v0 !!---------------------------------------------------------------------- IF( nn_timing == 1 ) CALL timing_start('tide_init') IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' IF(lwp) WRITE(numout,*) '~~~~~~~~~' REWIND(numnam) DO ib_bdy = 1, nb_bdy IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN td => tides(ib_bdy) ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries ln_tide_date = .false. ln_tide_v0 = .false. filtide(:) = '' tide_cpt(:) = '' ! Initialise bdy_ky_tide: updated in tide_update if using time correction otherwise defaults to 1 bdy_kt_tide=1 ! Don't REWIND here - may need to read more than one of these namelists. READ ( numnam, nambdy_tide ) ! ! Count number of components specified td%ncpt = 0 DO itide = 1, jpmax_harmo IF( tide_cpt(itide) /= '' ) THEN td%ncpt = td%ncpt + 1 ENDIF END DO CALL tide_init_Wave ! Find constituents in standard list ALLOCATE(bdy_ntide (td%ncpt)) DO itide=1,td%ncpt bdy_ntide(itide)=0 DO ji=1,jpmax_harmo IF ( TRIM( tide_cpt(itide) ) .eq. Wave(ji)%cname_tide) THEN bdy_ntide(itide) = ji EXIT END IF END DO IF (bdy_ntide(itide).eq.0) THEN CALL ctl_stop( 'BDYTIDE tidal components do not match up with tide.h90' ) ENDIF END DO ! Fill in phase speeds from tide_pulse ALLOCATE(bdy_omega_tide(td%ncpt)) CALL tide_pulse( bdy_omega_tide, bdy_ntide ,td%ncpt) ALLOCATE( td%speed(td%ncpt) ) td%speed = bdy_omega_tide(1:td%ncpt) ! ! Parameter control and print IF( td%ncpt < 1 ) THEN CALL ctl_stop( ' Did not find any tidal components in namelist nambdy_tide' ) ELSE IF(lwp) WRITE(numout,*) ' Namelist nambdy_tide : tidal harmonic forcing at open boundaries' IF(lwp) WRITE(numout,*) ' tidal components specified ', td%ncpt IF(lwp) WRITE(numout,*) ' ', tide_cpt(1:td%ncpt) IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : ' IF(lwp) WRITE(numout,*) ' ', td%speed(1:td%ncpt) ENDIF ! Allocate space for tidal harmonics data - ! get size from OBC data arrays ! --------------------------------------- ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) ALLOCATE( td%sshr( ilen0(1), td%ncpt, 2 ) ) ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) ALLOCATE( td%ur( ilen0(2), td%ncpt, 2 ) ) ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) ALLOCATE( td%vr( ilen0(3), td%ncpt, 2 ) ) ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) ! Set day length in timesteps for use if making phase and nodal corrections bdy_nn_tide=NINT(rday/rdt) ALLOCATE(bdy_v0tide (td%ncpt)) ALLOCATE(bdy_blank (td%ncpt)) ALLOCATE(bdy_utide (td%ncpt)) ALLOCATE(bdy_ftide (td%ncpt)) ALLOCATE(rbdy_ftide (td%ncpt)) ! Open files and read in tidal forcing data ! ----------------------------------------- DO itide = 1, td%ncpt ! ! SSH fields clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile CALL iom_open( clfile, inum ) CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1) CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1) CALL iom_close( inum ) ! ! U fields clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile CALL iom_open( clfile, inum ) CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) td%u(:,itide,1) = dta_read(1:ilen0(2),1,1) CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) td%u(:,itide,2) = dta_read(1:ilen0(2),1,1) CALL iom_close( inum ) ! ! V fields clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile CALL iom_open( clfile, inum ) CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) td%v(:,itide,1) = dta_read(1:ilen0(3),1,1) CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) td%v(:,itide,2) = dta_read(1:ilen0(3),1,1) CALL iom_close( inum ) ! END DO ! end loop on tidal components IF( ln_tide_date .and. ln_tide_v0 ) THEN ! correct for date factors: gather v0 CALL tide_harmo(bdy_omega_tide, bdy_v0tide, bdy_utide, bdy_ftide, bdy_ntide, td%ncpt, nn_tide_date) lcl_ryear = INT(nn_tide_date / 10000 ) lcl_rmonth = INT((nn_tide_date - lcl_ryear * 10000 ) / 100 ) lcl_rday = INT(nn_tide_date - lcl_ryear * 10000 - lcl_rmonth * 100) nyear = int(ndate0 / 10000 ) ! initial year nmonth = int((ndate0 - nyear * 10000 ) / 100 ) ! initial month nday = int(ndate0 - nyear * 10000 - nmonth * 100) CALL ymds2ju( nyear, nmonth, nday, 0._wp, fdayn ) CALL ymds2ju( lcl_ryear, lcl_rmonth, lcl_rday, 0._wp, fdayr ) bdy_tide_offset = NINT( fdayn - fdayr ) * 86400 IF(lwp) WRITE(numout,*) ' BDYTIDE offset ' IF(lwp) WRITE(numout,*) ' ', lcl_ryear, lcl_rmonth, lcl_rday IF(lwp) WRITE(numout,*) ' ', nyear, nmonth, nday IF(lwp) WRITE(numout,*) ' ', fdayn, fdayr, bdy_tide_offset ELSE bdy_v0tide(:)=0 bdy_utide(:)=0 bdy_ftide(:)=1 bdy_tide_offset = 0 IF(lwp) WRITE(numout,*) ' BDYTIDE offset ', bdy_tide_offset ENDIF ! Pass tidal forcing data to reference arrays for date correction to tidal harmonics DO itide = 1, td%ncpt ! loop on tidal components ! ! elevation igrd = 1 DO ib = 1, ilen0(igrd) td%sshr(ib,itide,1) = td%ssh(ib,itide,1) td%sshr(ib,itide,2) = td%ssh(ib,itide,2) END DO ! ! u igrd = 2 DO ib = 1, ilen0(igrd) td%ur(ib,itide,1) = td%u(ib,itide,1) td%ur(ib,itide,2) = td%u(ib,itide,2) END DO ! ! v igrd = 3 DO ib = 1, ilen0(igrd) td%vr(ib,itide,1) = td%v(ib,itide,1) td%vr(ib,itide,2) = td%v(ib,itide,2) ENDDO ENDDO ! loop on tidal components IF(lwp) WRITE(numout,*) 'BDYTIDE: summary of mappings' DO itide = 1, td%ncpt ! loop on tidal components IF(lwp) WRITE(numout,'(2i3,x,a)') itide, bdy_ntide(itide), tide_cpt(itide) ENDDO ! ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 ! END DO ! loop on ib_bdy IF( nn_timing == 1 ) CALL timing_stop('tide_init') END SUBROUTINE tide_init SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) !!---------------------------------------------------------------------- !! *** SUBROUTINE tide_update *** !! !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. !! !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! Main timestep counter !!gm doctor jit ==> kit TYPE(OBC_INDEX), INTENT( in ) :: idx ! OBC indices TYPE(OBC_DATA), INTENT(inout) :: dta ! OBC external data TYPE(TIDES_DATA),INTENT(inout) :: td ! tidal harmonics data INTEGER,INTENT(in),OPTIONAL :: jit ! Barotropic timestep counter (for timesplitting option) INTEGER,INTENT( in ), OPTIONAL :: time_offset ! time offset in units of timesteps. NB. if jit ! is present then units = subcycle timesteps. ! time_offset = 0 => get data at "now" time level ! time_offset = -1 => get data at "before" time level ! time_offset = +1 => get data at "after" time level ! etc. !! INTEGER :: itide, igrd, ib ! dummy loop indices INTEGER :: time_add ! time offset in units of timesteps INTEGER :: sub_step ! dummy for jit (probably not required as ! timesplitting always used?) REAL(wp) :: z_arg, z_sarg REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost REAL(wp) :: z_atde, z_btde REAL(wp) :: z1t, z2t !!---------------------------------------------------------------------- IF( nn_timing == 1 ) CALL timing_start('tide_update') time_add = 0 IF( PRESENT(time_offset) ) THEN time_add = time_offset ENDIF ! Phase corrections for the current day sub_step = 1 IF( PRESENT(jit) ) THEN sub_step = jit ENDIF IF( ln_tide_date ) THEN ! correct for date factors IF ( ( MOD( kt - 1, bdy_nn_tide ) == 0 ) .and. (sub_step==1) ) THEN IF ( ln_tide_v0 ) THEN bdy_kt_tide = 1 CALL tide_harmo(bdy_omega_tide, bdy_blank, bdy_utide, bdy_ftide, bdy_ntide, td%ncpt, ndastp) ELSE bdy_kt_tide = kt CALL tide_harmo(bdy_omega_tide, bdy_v0tide, bdy_utide, bdy_ftide, bdy_ntide, td%ncpt, ndastp) ENDIF DO itide = 1, td%ncpt ! loop on tidal components IF(lwp) WRITE(numout,*) 'BDYTIDE CORR:', itide, bdy_omega_tide(itide), bdy_v0tide(itide), & & bdy_utide(itide), bdy_ftide(itide) ENDDO ! Make adjustment for reference date in tidal harmonic data IF(lwp) WRITE(numout,*) 'BDYTIDE: nodal and phase correction at the start of day ', & & (kt-1)*rdt/rday + 1 DO itide = 1, td%ncpt ! loop on tidal components z_arg = bdy_utide(itide)+bdy_v0tide(itide) z_atde= bdy_ftide(itide)* cos(z_arg) z_btde= bdy_ftide(itide)* sin(z_arg) ! ! elevation igrd = 1 DO ib = 1, idx%nblenrim(igrd) z1t = z_atde * td%sshr(ib,itide,1) + z_btde * td%sshr(ib,itide,2) z2t = z_atde * td%sshr(ib,itide,2) - z_btde * td%sshr(ib,itide,1) td%ssh(ib,itide,1) = z1t td%ssh(ib,itide,2) = z2t END DO ! ! u igrd = 2 DO ib = 1, idx%nblenrim(igrd) z1t = z_atde * td%ur(ib,itide,1) + z_btde * td%ur(ib,itide,2) z2t = z_atde * td%ur(ib,itide,2) - z_btde * td%ur(ib,itide,1) td%u(ib,itide,1) = z1t td%u(ib,itide,2) = z2t END DO ! ! v igrd = 3 DO ib = 1, idx%nblenrim(igrd) z1t = z_atde * td%vr(ib,itide,1) + z_btde * td%vr(ib,itide,2) z2t = z_atde * td%vr(ib,itide,2) - z_btde * td%vr(ib,itide,1) td%v(ib,itide,1) = z1t td%v(ib,itide,2) = z2t ENDDO ENDDO ! loop on tidal components ENDIF ENDIF ! correct for date factors IF( PRESENT(jit) ) THEN IF( ln_tide_date ) THEN ! correct for date factors z_arg = ( (kt-bdy_kt_tide) * rdt + bdy_tide_offset + (jit+time_add) * rdt / REAL(nn_baro,wp) ) ELSE z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) ENDIF ELSE IF(lwp) WRITE(numout,*) 'BDYTIDE: should I be in here?' IF( ln_tide_date ) THEN ! correct for date factors z_arg = (kt+time_add-bdy_kt_tide+1) * rdt + bdy_tide_offset ELSE z_arg = (kt+time_add) * rdt ENDIF ENDIF DO itide = 1, td%ncpt z_sarg = z_arg * td%speed(itide) z_cost(itide) = COS( z_sarg ) z_sist(itide) = SIN( z_sarg ) END DO DO itide = 1, td%ncpt igrd=1 ! SSH on tracer grid. DO ib = 1, idx%nblenrim(igrd) dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) IF ( (idx%nbmap(ib,igrd) == 100 .and. (itide==10)) .and. (sub_step==1) ) THEN write(numout,*) 'z', ib, idx%nbmap(ib,igrd), idx%nbi(ib,igrd), idx%nbj(ib,igrd), & & itide, (td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) ENDIF ! IF ( (idx%nbmap(ib,igrd) == 100 .and. (itide==10)) ) THEN ! write(numout,*) 'z', ib, idx%nbmap(ib,igrd), idx%nbi(ib,igrd), idx%nbj(ib,igrd), & ! & itide, (td%ssh(ib,itide,1)*z_cost(itide) - td%ssh(ib,itide,2)*z_sist(itide)) ! ENDIF END DO igrd=2 ! U grid DO ib=1, idx%nblenrim(igrd) dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) ! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) END DO igrd=3 ! V grid DO ib=1, idx%nblenrim(igrd) dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) ! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) END DO END DO ! IF( nn_timing == 1 ) CALL timing_stop('tide_update') ! END SUBROUTINE tide_update #else !!---------------------------------------------------------------------- !! Dummy module NO Unstruct Open Boundary Conditions for tides !!---------------------------------------------------------------------- !!gm are you sure we need to define filtide and tide_cpt ? CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files CHARACTER(len=4 ), PUBLIC, DIMENSION(1) :: tide_cpt !: Names of tidal components used. CONTAINS SUBROUTINE tide_init ! Empty routine END SUBROUTINE tide_init SUBROUTINE tide_data ! Empty routine END SUBROUTINE tide_data SUBROUTINE tide_update( kt, kit ) ! Empty routine WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit END SUBROUTINE tide_update #endif !!====================================================================== END MODULE bdytides