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/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/UKMO/dev_5518_tide_analysis_restart/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 6116

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

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

Download in other formats: