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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 10759

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

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

Download in other formats: