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 branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 7246

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

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

Download in other formats: