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_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdytides.F90 @ 11223

Last change on this file since 11223 was 11223, checked in by smasson, 5 years ago
dev_r10984_HPC-13 : cleaning of rewriting of bdydta, see #2285
Property svn:keywords set to `Id`
File size: 27.0 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
[3651]	73	CHARACTER(len=80) :: clfile !: full file name for tidal input file
	74	REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data
[9125]	75	REAL(wp),ALLOCATABLE, DIMENSION(:,:) :: ztr, zti !: " " " " " " " "
[3294]	76	!!
[3651]	77	TYPE(TIDES_DATA), POINTER :: td !: local short cut
[3294]	78	!!
[3651]	79	NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
[1125]	80	!!----------------------------------------------------------------------
[6140]	81	!
[11223]	82	IF(lwp) WRITE(numout,*)
	83	IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
	84	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
[911]	85
[4147]	86	REWIND(numnam_cfg)
[3651]	87
[3294]	88	DO ib_bdy = 1, nb_bdy
[7646]	89	IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
	90	!
[3294]	91	td => tides(ib_bdy)
[911]	92
[3294]	93	! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
	94	filtide(:) = ''
[911]	95
[11223]	96	REWIND( numnam_ref )
[4147]	97	READ ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
[9168]	98	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in reference namelist', lwp )
[11223]	99	! Don't REWIND here - may need to read more than one of these namelists.
[4147]	100	READ ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
[9168]	101	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist', lwp )
[4624]	102	IF(lwm) WRITE ( numond, nambdy_tide )
[3651]	103	! ! Parameter control and print
	104	IF(lwp) WRITE(numout,*) ' '
	105	IF(lwp) WRITE(numout,*) ' Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
	106	IF(lwp) WRITE(numout,*) ' read tidal data in 2d files: ', ln_bdytide_2ddta
	107	IF(lwp) WRITE(numout,*) ' assume complex conjugate : ', ln_bdytide_conj
	108	IF(lwp) WRITE(numout,*) ' Number of tidal components to read: ', nb_harmo
	109	IF(lwp) THEN
[5084]	110	WRITE(numout,*) ' Tidal components: '
[3651]	111	DO itide = 1, nb_harmo
[5084]	112	WRITE(numout,*) ' ', Wave(ntide(itide))%cname_tide
[3651]	113	END DO
	114	ENDIF
	115	IF(lwp) WRITE(numout,*) ' '
[1125]	116
[3651]	117	! Allocate space for tidal harmonics data - get size from OBC data arrays
	118	! -----------------------------------------------------------------------
[911]	119
[3651]	120	! JC: If FRS scheme is used, we assume that tidal is needed over the whole
	121	! relaxation area
[11223]	122	IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = idx_bdy(ib_bdy)%nblen (:)
	123	ELSE ; ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
[3294]	124	ENDIF
[911]	125
[3651]	126	ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
	127	ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
[911]	128
[3651]	129	ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
	130	ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
[911]	131
[3651]	132	ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
	133	ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
[911]	134
[4292]	135	td%ssh0(:,:,:) = 0._wp
	136	td%ssh (:,:,:) = 0._wp
	137	td%u0 (:,:,:) = 0._wp
	138	td%u (:,:,:) = 0._wp
	139	td%v0 (:,:,:) = 0._wp
	140	td%v (:,:,:) = 0._wp
[3294]	141
[7646]	142	IF( ln_bdytide_2ddta ) THEN
[3651]	143	! It is assumed that each data file contains all complex harmonic amplitudes
[7646]	144	! given on the global domain (ie global, jpiglo x jpjglo)
[3651]	145	!
[9125]	146	ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
[3651]	147	!
	148	! SSH fields
	149	clfile = TRIM(filtide)//'_grid_T.nc'
[7646]	150	CALL iom_open( clfile , inum )
[3651]	151	igrd = 1 ! Everything is at T-points here
	152	DO itide = 1, nb_harmo
[7646]	153	CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
	154	CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )
[3651]	155	DO ib = 1, ilen0(igrd)
	156	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	157	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
[11048]	158	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?
[3651]	159	td%ssh0(ib,itide,1) = ztr(ii,ij)
	160	td%ssh0(ib,itide,2) = zti(ii,ij)
	161	END DO
	162	END DO
[3294]	163	CALL iom_close( inum )
[3651]	164	!
	165	! U fields
	166	clfile = TRIM(filtide)//'_grid_U.nc'
[7646]	167	CALL iom_open( clfile , inum )
[3651]	168	igrd = 2 ! Everything is at U-points here
	169	DO itide = 1, nb_harmo
[7646]	170	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
	171	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
[3651]	172	DO ib = 1, ilen0(igrd)
	173	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	174	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
[11048]	175	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?
[3651]	176	td%u0(ib,itide,1) = ztr(ii,ij)
	177	td%u0(ib,itide,2) = zti(ii,ij)
	178	END DO
	179	END DO
[3294]	180	CALL iom_close( inum )
[3651]	181	!
	182	! V fields
	183	clfile = TRIM(filtide)//'_grid_V.nc'
[7646]	184	CALL iom_open( clfile , inum )
[3651]	185	igrd = 3 ! Everything is at V-points here
	186	DO itide = 1, nb_harmo
[7646]	187	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
	188	CALL iom_get ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
[3651]	189	DO ib = 1, ilen0(igrd)
	190	ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
	191	ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
[11048]	192	IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?
[3651]	193	td%v0(ib,itide,1) = ztr(ii,ij)
	194	td%v0(ib,itide,2) = zti(ii,ij)
	195	END DO
	196	END DO
[3294]	197	CALL iom_close( inum )
	198	!
[9125]	199	DEALLOCATE( ztr, zti )
[3651]	200	!
	201	ELSE
	202	!
	203	! Read tidal data only on bdy segments
	204	!
	205	ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
[5132]	206	!
[3651]	207	! Open files and read in tidal forcing data
	208	! -----------------------------------------
[3294]	209
[3651]	210	DO itide = 1, nb_harmo
	211	! ! SSH fields
	212	clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
	213	CALL iom_open( clfile, inum )
[11223]	214	CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
[3651]	215	td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
[11223]	216	CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
[3651]	217	td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
	218	CALL iom_close( inum )
	219	! ! U fields
	220	clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
	221	CALL iom_open( clfile, inum )
[11223]	222	CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
[3651]	223	td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
[11223]	224	CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
[3651]	225	td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
	226	CALL iom_close( inum )
	227	! ! V fields
	228	clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
	229	CALL iom_open( clfile, inum )
[11223]	230	CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
[3651]	231	td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
[11223]	232	CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
[3651]	233	td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
	234	CALL iom_close( inum )
[3294]	235	!
[3651]	236	END DO ! end loop on tidal components
[3294]	237	!
[3651]	238	DEALLOCATE( dta_read )
[7646]	239	!
[3651]	240	ENDIF ! ln_bdytide_2ddta=.true.
[1125]	241	!
[6140]	242	IF( ln_bdytide_conj ) THEN ! assume complex conjugate in data files
[3651]	243	td%ssh0(:,:,2) = - td%ssh0(:,:,2)
	244	td%u0 (:,:,2) = - td%u0 (:,:,2)
	245	td%v0 (:,:,2) = - td%v0 (:,:,2)
	246	ENDIF
	247	!
[5930]	248	! Allocate slow varying data in the case of time splitting:
	249	! Do it anyway because at this stage knowledge of free surface scheme is unknown
	250	ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
	251	ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
	252	ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
[6140]	253	dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
	254	dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
	255	dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
[4292]	256	!
[7646]	257	ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
[1125]	258	!
[3294]	259	END DO ! loop on ib_bdy
[6140]	260	!
	261	END SUBROUTINE bdytide_init
[911]	262
[1125]	263
[11223]	264	SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset )
[1125]	265	!!----------------------------------------------------------------------
[3651]	266	!! * SUBROUTINE bdytide_update *
[1125]	267	!!
[3294]	268	!! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
[911]	269	!!
[1125]	270	!!----------------------------------------------------------------------
[6140]	271	INTEGER , INTENT(in ) :: kt ! Main timestep counter
	272	TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
	273	TYPE(OBC_DATA) , INTENT(inout) :: dta ! OBC external data
	274	TYPE(TIDES_DATA) , INTENT(inout) :: td ! tidal harmonics data
[11223]	275	INTEGER, OPTIONAL, INTENT(in ) :: kit ! Barotropic timestep counter (for timesplitting option)
	276	INTEGER, OPTIONAL, INTENT(in ) :: kt_offset ! time offset in units of timesteps. NB. if kit
[6140]	277	! ! is present then units = subcycle timesteps.
[11223]	278	! ! kt_offset = 0 => get data at "now" time level
	279	! ! kt_offset = -1 => get data at "before" time level
	280	! ! kt_offset = +1 => get data at "after" time level
[6140]	281	! ! etc.
	282	!
	283	INTEGER :: itide, igrd, ib ! dummy loop indices
	284	INTEGER :: time_add ! time offset in units of timesteps
	285	INTEGER, DIMENSION(3) :: ilen0 ! length of boundary data (from OBC arrays)
	286	REAL(wp) :: z_arg, z_sarg, zflag, zramp ! local scalars
[3651]	287	REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
[1125]	288	!!----------------------------------------------------------------------
[6140]	289	!
[3651]	290	ilen0(1) = SIZE(td%ssh(:,1,1))
	291	ilen0(2) = SIZE(td%u(:,1,1))
	292	ilen0(3) = SIZE(td%v(:,1,1))
	293
	294	zflag=1
[11223]	295	IF ( PRESENT(kit) ) THEN
	296	IF ( kit /= 1 ) zflag=0
[3651]	297	ENDIF
	298
[6140]	299	IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
[3651]	300	!
[6140]	301	kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
[3651]	302	!
	303	IF(lwp) THEN
	304	WRITE(numout,*)
	305	WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
	306	WRITE(numout,*) '~~~~~~~~~~~~~~ '
	307	ENDIF
	308	!
	309	CALL tide_init_elevation ( idx, td )
	310	CALL tide_init_velocities( idx, td )
	311	!
	312	ENDIF
	313
[3294]	314	time_add = 0
[11223]	315	IF( PRESENT(kt_offset) ) THEN
	316	time_add = kt_offset
[3294]	317	ENDIF
	318
[11223]	319	IF( PRESENT(kit) ) THEN
	320	z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp(time_add-1)) rdt / REAL(nn_baro,wp) )
[3294]	321	ELSE
[3651]	322	z_arg = ((kt-kt_tide)+time_add) * rdt
[1125]	323	ENDIF
[911]	324
[3651]	325	! Linear ramp on tidal component at open boundaries
[4292]	326	zramp = 1._wp
	327	IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)rdt)/(rdttideramprday),0._wp),1._wp)
[3651]	328
	329	DO itide = 1, nb_harmo
	330	z_sarg = z_arg * omega_tide(itide)
[1125]	331	z_cost(itide) = COS( z_sarg )
	332	z_sist(itide) = SIN( z_sarg )
	333	END DO
[911]	334
[3651]	335	DO itide = 1, nb_harmo
	336	igrd=1 ! SSH on tracer grid
	337	DO ib = 1, ilen0(igrd)
	338	dta%ssh(ib) = dta%ssh(ib) + zramp(td%ssh(ib,itide,1)z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
[1125]	339	END DO
[3294]	340	igrd=2 ! U grid
[3651]	341	DO ib = 1, ilen0(igrd)
	342	dta%u2d(ib) = dta%u2d(ib) + zramp(td%u (ib,itide,1)z_cost(itide) + td%u (ib,itide,2)*z_sist(itide))
[1125]	343	END DO
[3294]	344	igrd=3 ! V grid
[3651]	345	DO ib = 1, ilen0(igrd)
	346	dta%v2d(ib) = dta%v2d(ib) + zramp(td%v (ib,itide,1)z_cost(itide) + td%v (ib,itide,2)*z_sist(itide))
[1125]	347	END DO
	348	END DO
	349	!
[3651]	350	END SUBROUTINE bdytide_update
[911]	351
[7646]	352
[11223]	353	SUBROUTINE bdy_dta_tides( kt, kit, kt_offset )
[4292]	354	!!----------------------------------------------------------------------
	355	!! * SUBROUTINE bdy_dta_tides *
	356	!!
	357	!! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
	358	!!
	359	!!----------------------------------------------------------------------
[6140]	360	INTEGER, INTENT(in) :: kt ! Main timestep counter
	361	INTEGER, OPTIONAL, INTENT(in) :: kit ! Barotropic timestep counter (for timesplitting option)
[11223]	362	INTEGER, OPTIONAL, INTENT(in) :: kt_offset ! time offset in units of timesteps. NB. if kit
[6140]	363	! ! is present then units = subcycle timesteps.
[11223]	364	! ! kt_offset = 0 => get data at "now" time level
	365	! ! kt_offset = -1 => get data at "before" time level
	366	! ! kt_offset = +1 => get data at "after" time level
[6140]	367	! ! etc.
	368	!
	369	LOGICAL :: lk_first_btstp ! =.TRUE. if time splitting and first barotropic step
	370	INTEGER :: itide, ib_bdy, ib, igrd ! loop indices
	371	INTEGER :: time_add ! time offset in units of timesteps
	372	INTEGER, DIMENSION(jpbgrd) :: ilen0
	373	INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim ! short cuts
	374	REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist
[4292]	375	!!----------------------------------------------------------------------
[6140]	376	!
[4292]	377	lk_first_btstp=.TRUE.
	378	IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
	379
	380	time_add = 0
[11223]	381	IF( PRESENT(kt_offset) ) THEN
	382	time_add = kt_offset
[4292]	383	ENDIF
	384
	385	! Absolute time from model initialization:
	386	IF( PRESENT(kit) ) THEN
[5913]	387	z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
[4292]	388	ELSE
	389	z_arg = ( kt + time_add ) * rdt
	390	ENDIF
	391
	392	! Linear ramp on tidal component at open boundaries
	393	zramp = 1.
	394	IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000rdt)/(rdttideramprday),0.),1.)
	395
	396	DO ib_bdy = 1,nb_bdy
[7646]	397	!
	398	IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
	399	!
[4292]	400	nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
	401	nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
[7646]	402	!
	403	IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN ; ilen0(:) = nblen (:)
	404	ELSE ; ilen0(:) = nblenrim(:)
[4292]	405	ENDIF
[7646]	406	!
[4292]	407	! We refresh nodal factors every day below
	408	! This should be done somewhere else
[6140]	409	IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
[4292]	410	!
[6140]	411	kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
[4292]	412	!
	413	IF(lwp) THEN
	414	WRITE(numout,*)
	415	WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
	416	WRITE(numout,*) '~~~~~~~~~~~~~~ '
	417	ENDIF
	418	!
	419	CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
	420	CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
	421	!
	422	ENDIF
	423	zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
	424	!
[5930]	425	! If time splitting, initialize arrays from slow varying open boundary data:
	426	IF ( PRESENT(kit) ) THEN
[11223]	427	IF ( dta_bdy(ib_bdy)%lneed_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
	428	IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
	429	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]	430	ENDIF
	431	!
	432	! Update open boundary data arrays:
	433	DO itide = 1, nb_harmo
	434	!
	435	z_sarg = (z_arg + zoff) * omega_tide(itide)
	436	z_cost = zramp * COS( z_sarg )
	437	z_sist = zramp * SIN( z_sarg )
	438	!
[11223]	439	IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
[4758]	440	igrd=1 ! SSH on tracer grid
	441	DO ib = 1, ilen0(igrd)
	442	dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
	443	& ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
	444	& tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
	445	END DO
	446	ENDIF
[4292]	447	!
[11223]	448	IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
[4758]	449	igrd=2 ! U grid
	450	DO ib = 1, ilen0(igrd)
	451	dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
	452	& ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
	453	& tides(ib_bdy)%u(ib,itide,2)*z_sist )
	454	END DO
	455	igrd=3 ! V grid
	456	DO ib = 1, ilen0(igrd)
	457	dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
	458	& ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
	459	& tides(ib_bdy)%v(ib,itide,2)*z_sist )
	460	END DO
	461	ENDIF
	462	END DO
[4292]	463	END IF
	464	END DO
	465	!
	466	END SUBROUTINE bdy_dta_tides
	467
[6140]	468
[3651]	469	SUBROUTINE tide_init_elevation( idx, td )
[1125]	470	!!----------------------------------------------------------------------
[3651]	471	!! * ROUTINE tide_init_elevation *
[1125]	472	!!----------------------------------------------------------------------
[6140]	473	TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
	474	TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data
	475	!
	476	INTEGER :: itide, igrd, ib ! dummy loop indices
	477	INTEGER, DIMENSION(1) :: ilen0 ! length of boundary data (from OBC arrays)
[3651]	478	REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide
[6140]	479	!!----------------------------------------------------------------------
	480	!
[3651]	481	igrd=1
	482	! SSH on tracer grid.
	483	ilen0(1) = SIZE(td%ssh0(:,1,1))
[6140]	484	!
	485	ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
	486	!
[3651]	487	DO itide = 1, nb_harmo
	488	DO ib = 1, ilen0(igrd)
	489	mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)2.+td%ssh0(ib,itide,2)2.)
	490	phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
[1125]	491	END DO
[3651]	492	DO ib = 1 , ilen0(igrd)
	493	mod_tide(ib)=mod_tide(ib)*ftide(itide)
	494	phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
	495	ENDDO
	496	DO ib = 1 , ilen0(igrd)
	497	td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
	498	td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
	499	ENDDO
	500	END DO
[6140]	501	!
	502	DEALLOCATE( mod_tide, phi_tide )
	503	!
	504	END SUBROUTINE tide_init_elevation
[911]	505
	506
[3651]	507	SUBROUTINE tide_init_velocities( idx, td )
[1125]	508	!!----------------------------------------------------------------------
[3651]	509	!! * ROUTINE tide_init_elevation *
[1125]	510	!!----------------------------------------------------------------------
[6140]	511	TYPE(OBC_INDEX) , INTENT(in ) :: idx ! OBC indices
	512	TYPE(TIDES_DATA), INTENT(inout) :: td ! tidal harmonics data
	513	!
	514	INTEGER :: itide, igrd, ib ! dummy loop indices
	515	INTEGER, DIMENSION(3) :: ilen0 ! length of boundary data (from OBC arrays)
[3651]	516	REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide
[6140]	517	!!----------------------------------------------------------------------
	518	!
[3651]	519	ilen0(2) = SIZE(td%u0(:,1,1))
	520	ilen0(3) = SIZE(td%v0(:,1,1))
[6140]	521	!
[3651]	522	igrd=2 ! U grid.
[6140]	523	!
	524	ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
	525	!
[3651]	526	DO itide = 1, nb_harmo
	527	DO ib = 1, ilen0(igrd)
	528	mod_tide(ib)=SQRT(td%u0(ib,itide,1)2.+td%u0(ib,itide,2)2.)
	529	phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
	530	END DO
	531	DO ib = 1, ilen0(igrd)
	532	mod_tide(ib)=mod_tide(ib)*ftide(itide)
	533	phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
	534	ENDDO
	535	DO ib = 1, ilen0(igrd)
	536	td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
	537	td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
	538	ENDDO
	539	END DO
[6140]	540	!
	541	DEALLOCATE( mod_tide , phi_tide )
	542	!
[3651]	543	igrd=3 ! V grid.
[6140]	544	!
	545	ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
[911]	546
[3651]	547	DO itide = 1, nb_harmo
	548	DO ib = 1, ilen0(igrd)
	549	mod_tide(ib)=SQRT(td%v0(ib,itide,1)2.+td%v0(ib,itide,2)2.)
	550	phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
	551	END DO
	552	DO ib = 1, ilen0(igrd)
	553	mod_tide(ib)=mod_tide(ib)*ftide(itide)
	554	phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
	555	ENDDO
	556	DO ib = 1, ilen0(igrd)
	557	td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
	558	td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
	559	ENDDO
	560	END DO
[6140]	561	!
	562	DEALLOCATE( mod_tide, phi_tide )
	563	!
	564	END SUBROUTINE tide_init_velocities
[1125]	565
	566	!!======================================================================
[911]	567	END MODULE bdytides
[4292]	568

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

Download in other formats: