New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

bdytides.F90 in NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/BDY/bdytides.F90 @ 11641

Last change on this file since 11641 was 11641, checked in by acc, 5 years ago
Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Further substantive changes to the BDY routines that read multiple copies of some namelists from namelist_cfg. Special treatment is required to replicate the original behaviour with internal files. These changes emulate previous behaviour but removal of the reliance on multiple, same-named namelist blocks would be preferable.
Property svn:keywords set to `Id`
File size: 27.8 KB

Rev	Line
[911]	1	MODULE bdytides
[1125]	2	!!======================================================================
[911]	3	!! * MODULE bdytides *
	4	!! Ocean dynamics: Tidal forcing at open boundaries
[1125]	5	!!======================================================================
	6	!! History : 2.0 ! 2007-01 (D.Storkey) Original code
	7	!! 2.3 ! 2008-01 (J.Holt) Add date correction. Origins POLCOMS v6.3 2007
	8	!! 3.0 ! 2008-04 (NEMO team) add in the reference version
[2528]	9	!! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes
[3651]	10	!! 3.4 ! 2012-09 (G. Reffray and J. Chanut) New inputs + mods
[4292]	11	!! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes
[1125]	12	!!----------------------------------------------------------------------
[6140]	13	!! bdytide_init : read of namelist and initialisation of tidal harmonics data
	14	!! tide_update : calculation of tidal forcing at each timestep
[1125]	15	!!----------------------------------------------------------------------
[6140]	16	USE oce ! ocean dynamics and tracers
	17	USE dom_oce ! ocean space and time domain
	18	USE phycst ! physical constants
	19	USE bdy_oce ! ocean open boundary conditions
	20	USE tideini !
	21	USE daymod ! calendar
	22	!
	23	USE in_out_manager ! I/O units
	24	USE iom ! xIO server
	25	USE fldread !
	26	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[1125]	27
[911]	28	IMPLICIT NONE
	29	PRIVATE
	30
[3651]	31	PUBLIC bdytide_init ! routine called in bdy_init
	32	PUBLIC bdytide_update ! routine called in bdy_dta
[4292]	33	PUBLIC bdy_dta_tides ! routine called in dyn_spg_ts
[911]	34
[3294]	35	TYPE, PUBLIC :: TIDES_DATA !: Storage for external tidal harmonics data
[6140]	36	REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh0 !: Tidal constituents : SSH0 (read in file)
	37	REAL(wp), POINTER, DIMENSION(:,:,:) :: u0, v0 !: Tidal constituents : U0, V0 (read in file)
	38	REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh !: Tidal constituents : SSH (after nodal cor.)
	39	REAL(wp), POINTER, DIMENSION(:,:,:) :: u , v !: Tidal constituents : U , V (after nodal cor.)
[3294]	40	END TYPE TIDES_DATA
[911]	41
[4354]	42	!$AGRIF_DO_NOT_TREAT
[3651]	43	TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data
[4354]	44	!$AGRIF_END_DO_NOT_TREAT
[5930]	45	TYPE(OBC_DATA) , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s !: bdy external data (slow component)
[911]	46
[1125]	47	!!----------------------------------------------------------------------
[9598]	48	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[6140]	49	!! $Id$
[10068]	50	!! Software governed by the CeCILL license (see ./LICENSE)
[1125]	51	!!----------------------------------------------------------------------
[911]	52	CONTAINS
	53
[3651]	54	SUBROUTINE bdytide_init
[1125]	55	!!----------------------------------------------------------------------
[3651]	56	!! * SUBROUTINE bdytide_init *
[1125]	57	!!
[3294]	58	!! ** Purpose : - Read in namelist for tides and initialise external
	59	!! tidal harmonics data
[911]	60	!!
	61	!!----------------------------------------------------------------------
[3294]	62	!! namelist variables
	63	!!-------------------
[3651]	64	CHARACTER(len=80) :: filtide !: Filename root for tidal input files
	65	LOGICAL :: ln_bdytide_2ddta !: If true, read 2d harmonic data
	66	LOGICAL :: ln_bdytide_conj !: If true, assume complex conjugate tidal data
[1125]	67	!!
[3651]	68	INTEGER :: ib_bdy, itide, ib !: dummy loop indices
	69	INTEGER :: ii, ij !: dummy loop indices
[3294]	70	INTEGER :: inum, igrd
[3651]	71	INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays)
[4147]	72	INTEGER :: ios ! Local integer output status for namelist read
[11641]	73	INTEGER :: nbdy_rdstart, nbdy_loc
	74	CHARACTER(LEN=50) :: cerrmsg !: error string
[3651]	75	CHARACTER(len=80) :: clfile !: full file name for tidal input file
	76	REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data
[9125]	77	REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: ztr, zti !: " " " " " " " "
[3294]	78	!!
[3651]	79	TYPE(TIDES_DATA), POINTER :: td !: local short cut
[3294]	80	!!
[3651]	81	NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
[1125]	82	!!----------------------------------------------------------------------
[6140]	83	!
[11536]	84	IF(lwp) WRITE(numout,*)
	85	IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
	86	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
[911]	87
[3651]	88
[11641]	89	nbdy_rdstart = 1
[3294]	90	DO ib_bdy = 1, nb_bdy
[7646]	91	IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
	92	!
[3294]	93	td => tides(ib_bdy)
[911]	94
[3294]	95	! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
	96	filtide(:) = ''
[911]	97
[11536]	98	READ ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
	99	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
[11641]	100	!
	101	! Need to support possibility of reading more than one
	102	! nambdy_tide from the namelist_cfg internal file.
	103	! Do this by finding the ib_bdy'th occurence of nambdy_tide in the
	104	! character buffer as the starting point.
	105	!
	106	nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_tide' )
	107	IF( nbdy_loc .GT. 0 ) THEN
	108	nbdy_rdstart = nbdy_rdstart + nbdy_loc
	109	ELSE
	110	WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',ib_bdy,' of nambdy_tide not found'
	111	ios = -1
	112	CALL ctl_nam ( ios , cerrmsg )
	113	ENDIF
	114	nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
	115	READ ( numnam_cfg( nbdy_rdstart: ), nambdy_tide, IOSTAT = ios, ERR = 902)
[11536]	116	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
[4624]	117	IF(lwm) WRITE ( numond, nambdy_tide )
[3651]	118	! ! Parameter control and print
	119	IF(lwp) WRITE(numout,*) ' '
	120	IF(lwp) WRITE(numout,*) ' Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
	121	IF(lwp) WRITE(numout,*) ' read tidal data in 2d files: ', ln_bdytide_2ddta
	122	IF(lwp) WRITE(numout,*) ' assume complex conjugate : ', ln_bdytide_conj
	123	IF(lwp) WRITE(numout,*) ' Number of tidal components to read: ', nb_harmo
	124	IF(lwp) THEN
[5084]	125	WRITE(numout,*) ' Tidal components: '
[3651]	126	DO itide = 1, nb_harmo
[5084]	127	WRITE(numout,*) ' ', Wave(ntide(itide))%cname_tide
[3651]	128	END DO
	129	ENDIF
	130	IF(lwp) WRITE(numout,*) ' '
[1125]	131
[3651]	132	! Allocate space for tidal harmonics data - get size from OBC data arrays
	133	! -----------------------------------------------------------------------
[911]	134
[3651]	135	! JC: If FRS scheme is used, we assume that tidal is needed over the whole
	136	! relaxation area
[11536]	137	IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = idx_bdy(ib_bdy)%nblen (:)
	138	ELSE ; ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
[3294]	139	ENDIF
[911]	140
[3651]	141	ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
	142	ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
[911]	143
[3651]	144	ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
	145	ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
[911]	146
[3651]	147	ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
	148	ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
[911]	149
[4292]	150	td%ssh0(:,:,:) = 0._wp
	151	td%ssh (:,:,:) = 0._wp
	152	td%u0 (:,:,:) = 0._wp
	153	td%u (:,:,:) = 0._wp
	154	td%v0 (:,:,:) = 0._wp
	155	td%v (:,:,:) = 0._wp
[3294]	156
[7646]	157	IF( ln_bdytide_2ddta ) THEN
[3651]	158	! It is assumed that each data file contains all complex harmonic amplitudes
[7646]	159	! given on the global domain (ie global, jpiglo x jpjglo)
[3651]	160	!
[9125]	161	ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
[3651]	162	!
	163	! SSH fields
	164	clfile = TRIM(filtide)//'_grid_T.nc'
[7646]	165	CALL iom_open( clfile , inum )
[3651]	166	igrd = 1 ! Everything is at T-points here
	167	DO itide = 1, nb_harmo
[7646]	168	CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
	169	CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )
[3651]	170	DO ib = 1, ilen0(igrd)
	171	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	172	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
[11536]	173	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?
[3651]	174	td%ssh0(ib,itide,1) = ztr(ii,ij)
	175	td%ssh0(ib,itide,2) = zti(ii,ij)
	176	END DO
	177	END DO
[3294]	178	CALL iom_close( inum )
[3651]	179	!
	180	! U fields
	181	clfile = TRIM(filtide)//'_grid_U.nc'
[7646]	182	CALL iom_open( clfile , inum )
[3651]	183	igrd = 2 ! Everything is at U-points here
	184	DO itide = 1, nb_harmo
[7646]	185	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
	186	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
[3651]	187	DO ib = 1, ilen0(igrd)
	188	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	189	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
[11536]	190	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?
[3651]	191	td%u0(ib,itide,1) = ztr(ii,ij)
	192	td%u0(ib,itide,2) = zti(ii,ij)
	193	END DO
	194	END DO
[3294]	195	CALL iom_close( inum )
[3651]	196	!
	197	! V fields
	198	clfile = TRIM(filtide)//'_grid_V.nc'
[7646]	199	CALL iom_open( clfile , inum )
[3651]	200	igrd = 3 ! Everything is at V-points here
	201	DO itide = 1, nb_harmo
[7646]	202	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
	203	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
[3651]	204	DO ib = 1, ilen0(igrd)
	205	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	206	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
[11536]	207	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?
[3651]	208	td%v0(ib,itide,1) = ztr(ii,ij)
	209	td%v0(ib,itide,2) = zti(ii,ij)
	210	END DO
	211	END DO
[3294]	212	CALL iom_close( inum )
	213	!
[9125]	214	DEALLOCATE( ztr, zti )
[3651]	215	!
	216	ELSE
	217	!
	218	! Read tidal data only on bdy segments
	219	!
	220	ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
[5132]	221	!
[3651]	222	! Open files and read in tidal forcing data
	223	! -----------------------------------------
[3294]	224
[3651]	225	DO itide = 1, nb_harmo
	226	! ! SSH fields
	227	clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
	228	CALL iom_open( clfile, inum )
[11536]	229	CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
[3651]	230	td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
[11536]	231	CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
[3651]	232	td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
	233	CALL iom_close( inum )
	234	! ! U fields
	235	clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
	236	CALL iom_open( clfile, inum )
[11536]	237	CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
[3651]	238	td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
[11536]	239	CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
[3651]	240	td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
	241	CALL iom_close( inum )
	242	! ! V fields
	243	clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
	244	CALL iom_open( clfile, inum )
[11536]	245	CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
[3651]	246	td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
[11536]	247	CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
[3651]	248	td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
	249	CALL iom_close( inum )
[3294]	250	!
[3651]	251	END DO ! end loop on tidal components
[3294]	252	!
[3651]	253	DEALLOCATE( dta_read )
[7646]	254	!
[3651]	255	ENDIF ! ln_bdytide_2ddta=.true.
[1125]	256	!
[6140]	257	IF( ln_bdytide_conj ) THEN ! assume complex conjugate in data files
[3651]	258	td%ssh0(:,:,2) = - td%ssh0(:,:,2)
	259	td%u0 (:,:,2) = - td%u0 (:,:,2)
	260	td%v0 (:,:,2) = - td%v0 (:,:,2)
	261	ENDIF
	262	!
[5930]	263	! Allocate slow varying data in the case of time splitting:
	264	! Do it anyway because at this stage knowledge of free surface scheme is unknown
	265	ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
	266	ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
	267	ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
[6140]	268	dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
	269	dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
	270	dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
[4292]	271	!
[7646]	272	ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
[1125]	273	!
[3294]	274	END DO ! loop on ib_bdy
[6140]	275	!
	276	END SUBROUTINE bdytide_init
[911]	277
[1125]	278
[11536]	279	SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset )
[1125]	280	!!----------------------------------------------------------------------
[3651]	281	!! * SUBROUTINE bdytide_update *
[1125]	282	!!
[3294]	283	!! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
[911]	284	!!
[1125]	285	!!----------------------------------------------------------------------
[6140]	286	INTEGER , INTENT(in ) :: kt ! Main timestep counter
	287	TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
	288	TYPE(OBC_DATA) , INTENT(inout) :: dta ! OBC external data
	289	TYPE(TIDES_DATA) , INTENT(inout) :: td ! tidal harmonics data
[11536]	290	INTEGER, OPTIONAL, INTENT(in ) :: kit ! Barotropic timestep counter (for timesplitting option)
	291	INTEGER, OPTIONAL, INTENT(in ) :: kt_offset ! time offset in units of timesteps. NB. if kit
[6140]	292	! ! is present then units = subcycle timesteps.
[11536]	293	! ! kt_offset = 0 => get data at "now" time level
	294	! ! kt_offset = -1 => get data at "before" time level
	295	! ! kt_offset = +1 => get data at "after" time level
[6140]	296	! ! etc.
	297	!
	298	INTEGER :: itide, igrd, ib ! dummy loop indices
	299	INTEGER :: time_add ! time offset in units of timesteps
	300	INTEGER, DIMENSION(3) :: ilen0 ! length of boundary data (from OBC arrays)
	301	REAL(wp) :: z_arg, z_sarg, zflag, zramp ! local scalars
[3651]	302	REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
[1125]	303	!!----------------------------------------------------------------------
[6140]	304	!
[3651]	305	ilen0(1) = SIZE(td%ssh(:,1,1))
	306	ilen0(2) = SIZE(td%u(:,1,1))
	307	ilen0(3) = SIZE(td%v(:,1,1))
	308
	309	zflag=1
[11536]	310	IF ( PRESENT(kit) ) THEN
	311	IF ( kit /= 1 ) zflag=0
[3651]	312	ENDIF
	313
[6140]	314	IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
[3651]	315	!
[6140]	316	kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
[3651]	317	!
	318	IF(lwp) THEN
	319	WRITE(numout,*)
	320	WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
	321	WRITE(numout,*) '~~~~~~~~~~~~~~ '
	322	ENDIF
	323	!
	324	CALL tide_init_elevation ( idx, td )
	325	CALL tide_init_velocities( idx, td )
	326	!
	327	ENDIF
	328
[3294]	329	time_add = 0
[11536]	330	IF( PRESENT(kt_offset) ) THEN
	331	time_add = kt_offset
[3294]	332	ENDIF
	333
[11536]	334	IF( PRESENT(kit) ) THEN
	335	z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp(time_add-1)) rdt / REAL(nn_baro,wp) )
[3294]	336	ELSE
[3651]	337	z_arg = ((kt-kt_tide)+time_add) * rdt
[1125]	338	ENDIF
[911]	339
[3651]	340	! Linear ramp on tidal component at open boundaries
[4292]	341	zramp = 1._wp
	342	IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)rdt)/(rdttideramprday),0._wp),1._wp)
[3651]	343
	344	DO itide = 1, nb_harmo
	345	z_sarg = z_arg * omega_tide(itide)
[1125]	346	z_cost(itide) = COS( z_sarg )
	347	z_sist(itide) = SIN( z_sarg )
	348	END DO
[911]	349
[3651]	350	DO itide = 1, nb_harmo
	351	igrd=1 ! SSH on tracer grid
	352	DO ib = 1, ilen0(igrd)
	353	dta%ssh(ib) = dta%ssh(ib) + zramp(td%ssh(ib,itide,1)z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
[1125]	354	END DO
[3294]	355	igrd=2 ! U grid
[3651]	356	DO ib = 1, ilen0(igrd)
	357	dta%u2d(ib) = dta%u2d(ib) + zramp(td%u (ib,itide,1)z_cost(itide) + td%u (ib,itide,2)*z_sist(itide))
[1125]	358	END DO
[3294]	359	igrd=3 ! V grid
[3651]	360	DO ib = 1, ilen0(igrd)
	361	dta%v2d(ib) = dta%v2d(ib) + zramp(td%v (ib,itide,1)z_cost(itide) + td%v (ib,itide,2)*z_sist(itide))
[1125]	362	END DO
	363	END DO
	364	!
[3651]	365	END SUBROUTINE bdytide_update
[911]	366
[7646]	367
[11536]	368	SUBROUTINE bdy_dta_tides( kt, kit, kt_offset )
[4292]	369	!!----------------------------------------------------------------------
	370	!! * SUBROUTINE bdy_dta_tides *
	371	!!
	372	!! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
	373	!!
	374	!!----------------------------------------------------------------------
[6140]	375	INTEGER, INTENT(in) :: kt ! Main timestep counter
	376	INTEGER, OPTIONAL, INTENT(in) :: kit ! Barotropic timestep counter (for timesplitting option)
[11536]	377	INTEGER, OPTIONAL, INTENT(in) :: kt_offset ! time offset in units of timesteps. NB. if kit
[6140]	378	! ! is present then units = subcycle timesteps.
[11536]	379	! ! kt_offset = 0 => get data at "now" time level
	380	! ! kt_offset = -1 => get data at "before" time level
	381	! ! kt_offset = +1 => get data at "after" time level
[6140]	382	! ! etc.
	383	!
	384	LOGICAL :: lk_first_btstp ! =.TRUE. if time splitting and first barotropic step
	385	INTEGER :: itide, ib_bdy, ib, igrd ! loop indices
	386	INTEGER :: time_add ! time offset in units of timesteps
	387	INTEGER, DIMENSION(jpbgrd) :: ilen0
	388	INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim ! short cuts
	389	REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist
[4292]	390	!!----------------------------------------------------------------------
[6140]	391	!
[4292]	392	lk_first_btstp=.TRUE.
	393	IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
	394
	395	time_add = 0
[11536]	396	IF( PRESENT(kt_offset) ) THEN
	397	time_add = kt_offset
[4292]	398	ENDIF
	399
	400	! Absolute time from model initialization:
	401	IF( PRESENT(kit) ) THEN
[5913]	402	z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
[4292]	403	ELSE
	404	z_arg = ( kt + time_add ) * rdt
	405	ENDIF
	406
	407	! Linear ramp on tidal component at open boundaries
	408	zramp = 1.
	409	IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000rdt)/(rdttideramprday),0.),1.)
	410
	411	DO ib_bdy = 1,nb_bdy
[7646]	412	!
	413	IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
	414	!
[4292]	415	nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
	416	nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
[7646]	417	!
	418	IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = nblen (:)
	419	ELSE ; ilen0(:) = nblenrim(:)
[4292]	420	ENDIF
[7646]	421	!
[4292]	422	! We refresh nodal factors every day below
	423	! This should be done somewhere else
[6140]	424	IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
[4292]	425	!
[6140]	426	kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
[4292]	427	!
	428	IF(lwp) THEN
	429	WRITE(numout,*)
	430	WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
	431	WRITE(numout,*) '~~~~~~~~~~~~~~ '
	432	ENDIF
	433	!
	434	CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
	435	CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
	436	!
	437	ENDIF
	438	zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
	439	!
[5930]	440	! If time splitting, initialize arrays from slow varying open boundary data:
	441	IF ( PRESENT(kit) ) THEN
[11536]	442	IF ( dta_bdy(ib_bdy)%lneed_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
	443	IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
	444	IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
[4292]	445	ENDIF
	446	!
	447	! Update open boundary data arrays:
	448	DO itide = 1, nb_harmo
	449	!
	450	z_sarg = (z_arg + zoff) * omega_tide(itide)
	451	z_cost = zramp * COS( z_sarg )
	452	z_sist = zramp * SIN( z_sarg )
	453	!
[11536]	454	IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
[4758]	455	igrd=1 ! SSH on tracer grid
	456	DO ib = 1, ilen0(igrd)
	457	dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
	458	& ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
	459	& tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
	460	END DO
	461	ENDIF
[4292]	462	!
[11536]	463	IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
[4758]	464	igrd=2 ! U grid
	465	DO ib = 1, ilen0(igrd)
	466	dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
	467	& ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
	468	& tides(ib_bdy)%u(ib,itide,2)*z_sist )
	469	END DO
	470	igrd=3 ! V grid
	471	DO ib = 1, ilen0(igrd)
	472	dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
	473	& ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
	474	& tides(ib_bdy)%v(ib,itide,2)*z_sist )
	475	END DO
	476	ENDIF
	477	END DO
[4292]	478	END IF
	479	END DO
	480	!
	481	END SUBROUTINE bdy_dta_tides
	482
[6140]	483
[3651]	484	SUBROUTINE tide_init_elevation( idx, td )
[1125]	485	!!----------------------------------------------------------------------
[3651]	486	!! * ROUTINE tide_init_elevation *
[1125]	487	!!----------------------------------------------------------------------
[6140]	488	TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
	489	TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data
	490	!
	491	INTEGER :: itide, igrd, ib ! dummy loop indices
	492	INTEGER, DIMENSION(1) :: ilen0 ! length of boundary data (from OBC arrays)
[3651]	493	REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide
[6140]	494	!!----------------------------------------------------------------------
	495	!
[3651]	496	igrd=1
	497	! SSH on tracer grid.
	498	ilen0(1) = SIZE(td%ssh0(:,1,1))
[6140]	499	!
	500	ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
	501	!
[3651]	502	DO itide = 1, nb_harmo
	503	DO ib = 1, ilen0(igrd)
	504	mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)2.+td%ssh0(ib,itide,2)2.)
	505	phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
[1125]	506	END DO
[3651]	507	DO ib = 1 , ilen0(igrd)
	508	mod_tide(ib)=mod_tide(ib)*ftide(itide)
	509	phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
	510	ENDDO
	511	DO ib = 1 , ilen0(igrd)
	512	td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
	513	td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
	514	ENDDO
	515	END DO
[6140]	516	!
	517	DEALLOCATE( mod_tide, phi_tide )
	518	!
	519	END SUBROUTINE tide_init_elevation
[911]	520
	521
[3651]	522	SUBROUTINE tide_init_velocities( idx, td )
[1125]	523	!!----------------------------------------------------------------------
[3651]	524	!! * ROUTINE tide_init_elevation *
[1125]	525	!!----------------------------------------------------------------------
[6140]	526	TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
	527	TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data
	528	!
	529	INTEGER :: itide, igrd, ib ! dummy loop indices
	530	INTEGER, DIMENSION(3) :: ilen0 ! length of boundary data (from OBC arrays)
[3651]	531	REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide
[6140]	532	!!----------------------------------------------------------------------
	533	!
[3651]	534	ilen0(2) = SIZE(td%u0(:,1,1))
	535	ilen0(3) = SIZE(td%v0(:,1,1))
[6140]	536	!
[3651]	537	igrd=2 ! U grid.
[6140]	538	!
	539	ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
	540	!
[3651]	541	DO itide = 1, nb_harmo
	542	DO ib = 1, ilen0(igrd)
	543	mod_tide(ib)=SQRT(td%u0(ib,itide,1)2.+td%u0(ib,itide,2)2.)
	544	phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
	545	END DO
	546	DO ib = 1, ilen0(igrd)
	547	mod_tide(ib)=mod_tide(ib)*ftide(itide)
	548	phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
	549	ENDDO
	550	DO ib = 1, ilen0(igrd)
	551	td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
	552	td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
	553	ENDDO
	554	END DO
[6140]	555	!
	556	DEALLOCATE( mod_tide , phi_tide )
	557	!
[3651]	558	igrd=3 ! V grid.
[6140]	559	!
	560	ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
[911]	561
[3651]	562	DO itide = 1, nb_harmo
	563	DO ib = 1, ilen0(igrd)
	564	mod_tide(ib)=SQRT(td%v0(ib,itide,1)2.+td%v0(ib,itide,2)2.)
	565	phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
	566	END DO
	567	DO ib = 1, ilen0(igrd)
	568	mod_tide(ib)=mod_tide(ib)*ftide(itide)
	569	phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
	570	ENDDO
	571	DO ib = 1, ilen0(igrd)
	572	td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
	573	td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
	574	ENDDO
	575	END DO
[6140]	576	!
	577	DEALLOCATE( mod_tide, phi_tide )
	578	!
	579	END SUBROUTINE tide_init_velocities
[1125]	580
	581	!!======================================================================
[911]	582	END MODULE bdytides
[4292]	583

Note: See TracBrowser for help on using the repository browser.

Download in other formats: