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 @ 11049

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

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

Download in other formats: