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.

dynspg_ts.F90 in branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4380

Last change on this file since 4380 was 4380, checked in by hliu, 10 years ago
wetting and drying: some bugs removed from dynspg_ts.F90
Property svn:keywords set to `Id`
File size: 44.6 KB

Rev	Line
[358]	1	MODULE dynspg_ts
	2	!!======================================================================
[1502]	3	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
	4	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
	5	!! - ! 2006-08 (S. Masson) distributed restart using iom
	6	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
	7	!! - ! 2008-01 (R. Benshila) change averaging method
	8	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
[2528]	9	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
[2724]	10	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
	11	!!---------------------------------------------------------------------
[575]	12	#if defined key_dynspg_ts \|\| defined key_esopa
[358]	13	!!----------------------------------------------------------------------
[455]	14	!! 'key_dynspg_ts' free surface cst volume with time splitting
[358]	15	!!----------------------------------------------------------------------
	16	!! dyn_spg_ts : compute surface pressure gradient trend using a time-
	17	!! splitting scheme and add to the general trend
[508]	18	!! ts_rst : read/write the time-splitting restart fields in the ocean restart file
[358]	19	!!----------------------------------------------------------------------
	20	USE oce ! ocean dynamics and tracers
	21	USE dom_oce ! ocean space and time domain
[888]	22	USE sbc_oce ! surface boundary condition: ocean
	23	USE dynspg_oce ! surface pressure gradient variables
[358]	24	USE phycst ! physical constants
[888]	25	USE domvvl ! variable volume
[1662]	26	USE zdfbfr ! bottom friction
[358]	27	USE dynvor ! vorticity term
	28	USE obc_oce ! Lateral open boundary condition
[371]	29	USE obc_par ! open boundary condition parameters
[3294]	30	USE obcdta ! open boundary condition data
	31	USE obcfla ! Flather open boundary condition
	32	USE bdy_par ! for lk_bdy
	33	USE bdy_oce ! Lateral open boundary condition
	34	USE bdydta ! open boundary condition data
	35	USE bdydyn2d ! open boundary conditions on barotropic variables
	36	USE sbctide
	37	USE updtide
[358]	38	USE lib_mpp ! distributed memory computing library
	39	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	40	USE prtctl ! Print control
	41	USE in_out_manager ! I/O manager
[2715]	42	USE iom ! IOM library
[3294]	43	USE zdf_oce ! Vertical diffusion
	44	USE wrk_nemo ! Memory Allocation
	45	USE timing ! Timing
[358]	46
[3294]	47
[358]	48	IMPLICIT NONE
	49	PRIVATE
	50
[2715]	51	PUBLIC dyn_spg_ts ! routine called by step.F90
	52	PUBLIC ts_rst ! routine called by istate.F90
	53	PUBLIC dyn_spg_ts_alloc ! routine called by dynspg.F90
[358]	54
[1502]	55
[2715]	56	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
	57	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
[508]	58
[2715]	59	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_b, vn_b ! now averaged velocity
	60	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub_b, vb_b ! before averaged velocity
[1502]	61
[358]	62	!! * Substitutions
	63	# include "domzgr_substitute.h90"
	64	# include "vectopt_loop_substitute.h90"
[2715]	65	!!----------------------------------------------------------------------
	66	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[888]	67	!! $Id$
[2715]	68	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	69	!!----------------------------------------------------------------------
[358]	70	CONTAINS
	71
[2715]	72	INTEGER FUNCTION dyn_spg_ts_alloc()
	73	!!----------------------------------------------------------------------
	74	!! * routine dyn_spg_ts_alloc *
	75	!!----------------------------------------------------------------------
	76	ALLOCATE( ftnw (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) , &
	77	& ftsw (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc )
	78	!
	79	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
	80	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
	81	!
	82	END FUNCTION dyn_spg_ts_alloc
	83
	84
[358]	85	SUBROUTINE dyn_spg_ts( kt )
	86	!!----------------------------------------------------------------------
	87	!! * routine dyn_spg_ts *
	88	!!
	89	!! ** Purpose : Compute the now trend due to the surface pressure
	90	!! gradient in case of free surface formulation with time-splitting.
	91	!! Add it to the general trend of momentum equation.
	92	!!
	93	!! ** Method : Free surface formulation with time-splitting
	94	!! -1- Save the vertically integrated trend. This general trend is
	95	!! held constant over the barotropic integration.
	96	!! The Coriolis force is removed from the general trend as the
	97	!! surface gradient and the Coriolis force are updated within
	98	!! the barotropic integration.
[367]	99	!! -2- Barotropic loop : updates of sea surface height (ssha_e) and
[1502]	100	!! barotropic velocity (ua_e and va_e) through barotropic
[358]	101	!! momentum and continuity integration. Barotropic former
	102	!! variables are time averaging over the full barotropic cycle
[2528]	103	!! (= 2 * baroclinic time step) and saved in uX_b
	104	!! and vX_b (X specifying after, now or before).
[1438]	105	!! -3- The new general trend becomes :
[2528]	106	!! ua = ua - sum_k(ua)/H + ( un_b - ub_b )
[358]	107	!!
	108	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
	109	!!
[508]	110	!! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
[358]	111	!!---------------------------------------------------------------------
[2715]	112	!
[1502]	113	INTEGER, INTENT(in) :: kt ! ocean time-step index
[2715]	114	!
[1662]	115	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	116	INTEGER :: icycle ! local scalar
[3294]	117	INTEGER :: ikbu, ikbv ! local scalar
	118	REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars
	119	REAL(wp) :: z1_8, zx1, zy1 ! - -
	120	REAL(wp) :: z1_4, zx2, zy2 ! - -
	121	REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - -
	122	REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - -
	123	REAL(wp) :: ua_btm, va_btm ! - -
	124	!
	125	REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv
	126	REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e
	127	REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum
[4376]	128
[4380]	129	REAL(wp), POINTER, DIMENSION(:,:) :: zwadfltu, zwadfltv
[358]	130	!!----------------------------------------------------------------------
[3294]	131	!
	132	IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts')
	133	!
	134	CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv )
	135	CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e )
	136	CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
[4376]	137
[4380]	138	IF(ln_wad) CALL wrk_alloc( jpi, jpj, zwadfltu, zwadfltv)
[3294]	139	!
[1502]	140	IF( kt == nit000 ) THEN !* initialisation
[508]	141	!
[358]	142	IF(lwp) WRITE(numout,*)
	143	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
	144	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
[1241]	145	IF(lwp) WRITE(numout,) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2nn_baro
[1502]	146	!
[1708]	147	CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b
[1502]	148	!
	149	ua_e (:,:) = un_b (:,:)
	150	va_e (:,:) = vn_b (:,:)
	151	hu_e (:,:) = hu (:,:)
	152	hv_e (:,:) = hv (:,:)
	153	hur_e (:,:) = hur (:,:)
	154	hvr_e (:,:) = hvr (:,:)
[358]	155	IF( ln_dynvor_een ) THEN
[3294]	156	ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
[358]	157	DO jj = 2, jpj
[1708]	158	DO ji = fs_2, jpi ! vector opt.
[3294]	159	ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp
	160	ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp
	161	ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp
	162	ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp
[358]	163	END DO
	164	END DO
	165	ENDIF
[508]	166	!
	167	ENDIF
[358]	168
[3294]	169	! !* Local constant initialization
	170	z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! reciprocal of baroclinic time step
	171	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! reciprocal of baroclinic
	172	! time step (euler timestep)
	173	z1_8 = 0.125_wp ! coefficient for vorticity estimates
	174	z1_4 = 0.25_wp
	175	zraur = 1._wp / rau0 ! 1 / volumic mass
[1502]	176	!
[3294]	177	zhdiv(:,:) = 0._wp ! barotropic divergence
	178	zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F)
	179	zv_sld = 0._wp ; zv_asp = 0._wp
[1438]	180
[3294]	181	IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction
	182	z2dt_bf = rdt
	183	ELSE
	184	z2dt_bf = 2.0_wp * rdt
	185	ENDIF
	186
[4380]	187	IF(ln_wad) THEN
	188	zwadfltu(:,:) = 1._wp
	189	zwadfltv(:,:) = 1._wp
	190	END IF
[4376]	191
[358]	192	! -----------------------------------------------------------------------------
	193	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
	194	! -----------------------------------------------------------------------------
[1502]	195	!
	196	! !* e3d/dt(Ua), e3Ub, e3*Vn (Vertically integrated)
	197	! ! --------------------------
[3294]	198	zua(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp
	199	zva(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp
[1502]	200	!
	201	DO jk = 1, jpkm1
	202	#if defined key_vectopt_loop
	203	DO jj = 1, 1 !Vector opt. => forced unrolling
[358]	204	DO ji = 1, jpij
[1502]	205	#else
	206	DO jj = 1, jpj
	207	DO ji = 1, jpi
	208	#endif
	209	! ! now trend
	210	zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
	211	zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
	212	! ! now velocity
	213	zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk)
	214	zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk)
[2724]	215	!
	216	#if defined key_vvl
[3294]	217	ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk) *umask(ji,jj,jk)
	218	vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk) *vmask(ji,jj,jk)
[2724]	219	#else
	220	ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
	221	vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
	222	#endif
[358]	223	END DO
	224	END DO
[1502]	225	END DO
	226
	227	! !* baroclinic momentum trend (remove the vertical mean trend)
	228	DO jk = 1, jpkm1 ! --------------------------
	229	DO jj = 2, jpjm1
	230	DO ji = fs_2, fs_jpim1 ! vector opt.
	231	ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
	232	va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
	233	END DO
[358]	234	END DO
[1502]	235	END DO
[358]	236
[1502]	237	! !* barotropic Coriolis trends * H (vorticity scheme dependent)
	238	! ! ---------------------------====
	239	zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport
	240	zwy(:,:) = zvn(:,:) * e1v(:,:)
	241	!
[358]	242	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
	243	DO jj = 2, jpjm1
	244	DO ji = fs_2, fs_jpim1 ! vector opt.
	245	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	246	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	247	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	248	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	249	! energy conserving formulation for planetary vorticity term
[1502]	250	zcu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )
	251	zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )
[358]	252	END DO
	253	END DO
[508]	254	!
[358]	255	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
	256	DO jj = 2, jpjm1
	257	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	258	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	259	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[358]	260	zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) )
	261	zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) )
	262	END DO
	263	END DO
[508]	264	!
[358]	265	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
	266	DO jj = 2, jpjm1
	267	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	268	zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	269	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
	270	zcv(ji,jj) = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &
	271	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) )
[358]	272	END DO
	273	END DO
[508]	274	!
[358]	275	ENDIF
	276
[1502]	277	! !* Right-Hand-Side of the barotropic momentum equation
	278	! ! ----------------------------------------------------
	279	IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient
	280	DO jj = 2, jpjm1
[358]	281	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	282	zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) &
	283	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)
	284	zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) &
	285	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)
[358]	286	END DO
	287	END DO
[1502]	288	ENDIF
[358]	289
[1502]	290	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
[358]	291	DO ji = fs_2, fs_jpim1
[3294]	292	zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
	293	zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
	294	END DO
[358]	295	END DO
	296
[1708]	297
	298	! ! Remove barotropic contribution of bottom friction
	299	! ! from the barotropic transport trend
[3294]	300	zcoef = -1._wp * z1_2dt_b
	301
	302	IF(ln_bfrimp) THEN
	303	! ! Remove the bottom stress trend from 3-D sea surface level gradient
	304	! ! and Coriolis forcing in case of 3D semi-implicit bottom friction
	305	DO jj = 2, jpjm1
	306	DO ji = fs_2, fs_jpim1
	307	ikbu = mbku(ji,jj)
	308	ikbv = mbkv(ji,jj)
	309	ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu)
	310	va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv)
	311
	312	zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm
	313	zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm
	314	END DO
	315	END DO
	316
	317	ELSE
	318
[1708]	319	# if defined key_vectopt_loop
[3294]	320	DO jj = 1, 1
	321	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[1708]	322	# else
[3294]	323	DO jj = 2, jpjm1
	324	DO ji = 2, jpim1
[1708]	325	# endif
	326	! Apply stability criteria for bottom friction
[2528]	327	!RBbug for vvl and external mode we may need to use varying fse3
	328	!!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
[3294]	329	zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef )
	330	zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef )
	331	END DO
	332	END DO
[1708]	333
[3294]	334	IF( lk_vvl ) THEN
	335	DO jj = 2, jpjm1
	336	DO ji = fs_2, fs_jpim1 ! vector opt.
	337	zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) &
	338	& / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) )
	339	zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) &
	340	& / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) )
	341	END DO
	342	END DO
	343	ELSE
	344	DO jj = 2, jpjm1
	345	DO ji = fs_2, fs_jpim1 ! vector opt.
	346	zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
	347	zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
	348	END DO
	349	END DO
	350	ENDIF
	351	END IF ! end (ln_bfrimp)
[1662]	352
[3294]	353
[1502]	354	! !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
	355	! ! --------------------------
	356	zua(:,:) = zua(:,:) * hur(:,:)
	357	zva(:,:) = zva(:,:) * hvr(:,:)
	358	!
[2724]	359	IF( lk_vvl ) THEN
[3294]	360	ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
	361	vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
[2724]	362	ELSE
	363	ub_b(:,:) = ub_b(:,:) * hur(:,:)
	364	vb_b(:,:) = vb_b(:,:) * hvr(:,:)
	365	ENDIF
[1502]	366
[358]	367	! -----------------------------------------------------------------------
	368	! Phase 2 : Integration of the barotropic equations with time splitting
	369	! -----------------------------------------------------------------------
[1502]	370	!
	371	! ! ==================== !
	372	! ! Initialisations !
	373	! ! ==================== !
	374	icycle = 2 * nn_baro ! Number of barotropic sub time-step
	375
	376	! ! Start from NOW field
	377	hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points
	378	hv_e (:,:) = hv (:,:)
	379	hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points
	380	hvr_e (:,:) = hvr (:,:)
[1662]	381	!RBbug zsshb_e(:,:) = sshn (:,:)
	382	zsshb_e(:,:) = sshn_b(:,:) ! sea surface height (before and now)
[1502]	383	sshn_e (:,:) = sshn (:,:)
	384
	385	zun_e (:,:) = un_b (:,:) ! barotropic velocity (external)
	386	zvn_e (:,:) = vn_b (:,:)
	387	zub_e (:,:) = un_b (:,:)
	388	zvb_e (:,:) = vn_b (:,:)
[358]	389
[1502]	390	zu_sum (:,:) = un_b (:,:) ! summation
	391	zv_sum (:,:) = vn_b (:,:)
	392	zssh_sum(:,:) = sshn (:,:)
[358]	393
[1502]	394	#if defined key_obc
[367]	395	! set ssh corrections to 0
	396	! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
[3294]	397	IF( lp_obc_east ) sshfoe_b(:,:) = 0._wp
	398	IF( lp_obc_west ) sshfow_b(:,:) = 0._wp
	399	IF( lp_obc_south ) sshfos_b(:,:) = 0._wp
	400	IF( lp_obc_north ) sshfon_b(:,:) = 0._wp
[367]	401	#endif
	402
[1502]	403	! ! ==================== !
	404	DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1)
	405	! ! ==================== !
[1241]	406	z2dt_e = 2. * ( rdt / nn_baro )
[1502]	407	IF( jn == 1 ) z2dt_e = rdt / nn_baro
[358]	408
[3294]	409	! !* Update the forcing (BDY and tides)
[1502]	410	! ! ------------------
[2528]	411	IF( lk_obc ) CALL obc_dta_bt ( kt, jn )
[3294]	412	IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=+1 )
[3651]	413	IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn )
[367]	414
[1502]	415	! !* after ssh_e
	416	! ! -----------
[4380]	417	IF(ln_wad) THEN
	418	DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports
	419	DO ji = fs_2, fs_jpim1 ! vector opt.
	420	zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) * zwadfltu(ji, jj ) &
	421	& - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) * zwadfltu(ji-1, jj ) &
	422	& + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) * zwadfltv(ji, jj ) &
	423	& - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) * zwadfltv(ji, jj-1) ) &
	424	& / ( e1t(ji,jj) * e2t(ji,jj) )
	425	END DO
[358]	426	END DO
[4380]	427	ELSE
	428	DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports
	429	DO ji = fs_2, fs_jpim1 ! vector opt.
	430	zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) &
	431	& - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) &
	432	& + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) &
	433	& - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) &
	434	& / ( e1t(ji,jj) * e2t(ji,jj) )
	435	END DO
	436	END DO
	437	END IF
[1502]	438	!
[358]	439	#if defined key_obc
[1502]	440	! ! OBC : zhdiv must be zero behind the open boundary
	441	!! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
[3294]	442	IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0._wp ! east
	443	IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0._wp ! west
	444	IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0._wp ! north
	445	IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0._wp ! south
[358]	446	#endif
[1170]	447	#if defined key_bdy
[1502]	448	zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask
[1170]	449	#endif
[1502]	450	!
[4380]	451	DO jj = 2, jpjm1 ! leap-frog on ssh_e
	452	DO ji = fs_2, fs_jpim1 ! vector opt.
	453	ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * &
	454	& ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
	455	END DO
	456	END DO
	457
	458	!! generate W/D filter
[4376]	459	IF(ln_wad) THEN
[4380]	460	DO jj = 2, jpjm1
[4376]	461	DO ji = fs_2, fs_jpim1 ! vector opt.
[4380]	462	IF(ssha_e(ji,jj) + bathy(ji,jj) <= rn_wadmin) THEN
	463	zwadfltu(ji, jj ) = 0._wp
	464	zwadfltu(ji-1, jj ) = 0._wp
	465	zwadfltv(ji, jj ) = 0._wp
	466	zwadfltv(ji, jj-1) = 0._wp
[4376]	467	END IF
	468	END DO
[358]	469	END DO
[4380]	470	CALL lbc_lnk(zwadfltu, 'U', 1._wp)
	471	CALL lbc_lnk(zwadfltv, 'V', 1._wp)
[4376]	472	END IF
[358]	473
[4380]	474
[1502]	475	! !* after barotropic velocities (vorticity scheme dependent)
	476	! ! ---------------------------
[3294]	477	zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport
[1502]	478	zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
	479	!
	480	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
[358]	481	DO jj = 2, jpjm1
	482	DO ji = fs_2, fs_jpim1 ! vector opt.
	483	! surface pressure gradient
[592]	484	IF( lk_vvl) THEN
[1662]	485	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	486	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	487	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	488	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	489	ELSE
[1662]	490	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	491	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	492	ENDIF
[3294]	493	! add tidal astronomical forcing
[3651]	494	IF ( ln_tide_pot .AND. lk_tide ) THEN
[3294]	495	zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
	496	zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
	497	ENDIF
[358]	498	! energy conserving formulation for planetary vorticity term
	499	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	500	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	501	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	502	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[1662]	503	zu_cor = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
	504	zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
	505	! after velocities with implicit bottom friction
[3294]	506
	507	IF( ln_bfrimp ) THEN ! implicit bottom friction
	508	! A new method to implement the implicit bottom friction.
	509	! H. Liu
	510	! Sept 2011
	511	ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + &
	512	& z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) &
	513	& / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
	514	ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1)
	515	!
	516	va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + &
	517	& z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) &
	518	& / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
	519	va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1)
	520	!
	521	ELSE
	522	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) &
	523	& / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
	524	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) &
	525	& / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
	526	ENDIF
[358]	527	END DO
	528	END DO
[508]	529	!
[1502]	530	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
[358]	531	DO jj = 2, jpjm1
	532	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	533	! surface pressure gradient
[592]	534	IF( lk_vvl) THEN
[1662]	535	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	536	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	537	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	538	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	539	ELSE
[1662]	540	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	541	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	542	ENDIF
[3294]	543	! add tidal astronomical forcing
[3651]	544	IF ( ln_tide_pot .AND. lk_tide ) THEN
[3294]	545	zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
	546	zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
	547	ENDIF
[358]	548	! enstrophy conserving formulation for planetary vorticity term
[1502]	549	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	550	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[1662]	551	zu_cor = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
	552	zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj)
	553	! after velocities with implicit bottom friction
[3294]	554	IF( ln_bfrimp ) THEN
	555	! A new method to implement the implicit bottom friction.
	556	! H. Liu
	557	! Sept 2011
	558	ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + &
	559	& z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) &
	560	& / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
	561	ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1)
	562	!
	563	va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + &
	564	& z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) &
	565	& / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
	566	va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1)
	567	!
	568	ELSE
	569	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) &
	570	& / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
	571	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) &
	572	& / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
	573	ENDIF
[358]	574	END DO
	575	END DO
[508]	576	!
[1502]	577	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
[358]	578	DO jj = 2, jpjm1
	579	DO ji = fs_2, fs_jpim1 ! vector opt.
	580	! surface pressure gradient
[592]	581	IF( lk_vvl) THEN
[1662]	582	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	583	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	584	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	585	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	586	ELSE
[1662]	587	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	588	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	589	ENDIF
[3294]	590	! add tidal astronomical forcing
[3651]	591	IF ( ln_tide_pot .AND. lk_tide ) THEN
[3294]	592	zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
	593	zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
	594	ENDIF
[358]	595	! energy/enstrophy conserving formulation for planetary vorticity term
[1662]	596	zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
[1502]	597	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
[1662]	598	zv_cor = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &
[1502]	599	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj)
[1662]	600	! after velocities with implicit bottom friction
[3294]	601	IF( ln_bfrimp ) THEN
	602	! A new method to implement the implicit bottom friction.
	603	! H. Liu
	604	! Sept 2011
	605	ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + &
	606	& z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) &
	607	& / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) )
	608	ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1)
	609	!
	610	va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + &
	611	& z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) &
	612	& / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) )
	613	va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1)
	614	!
	615	ELSE
	616	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) &
	617	& / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
	618	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) &
	619	& / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
	620	ENDIF
[358]	621	END DO
	622	END DO
[508]	623	!
[358]	624	ENDIF
[3294]	625	! !* domain lateral boundary
	626	! ! -----------------------
[358]	627
[3294]	628	! OBC open boundaries
[2715]	629	IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
[3294]	630
	631	! BDY open boundaries
	632	#if defined key_bdy
	633	pssh => sshn_e
	634	phur => hur_e
	635	phvr => hvr_e
	636	pu2d => ua_e
	637	pv2d => va_e
	638
	639	IF( lk_bdy ) CALL bdy_dyn2d( kt )
	640	#endif
	641
[4380]	642	IF(ln_wad) THEN
	643	DO jj = 2, jpjm1
	644	DO ji = fs_2, fs_jpim1 ! vector opt.
	645	ua_e(ji,jj) = ua_e(ji,jj) * zwadfltu(ji,jj)
	646	va_e(ji,jj) = va_e(ji,jj) * zwadfltv(ji,jj)
	647	END DO
	648	END DO
	649	END IF
	650
[1502]	651	!
	652	CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries
	653	CALL lbc_lnk( va_e , 'V', -1. )
	654	CALL lbc_lnk( ssha_e, 'T', 1. )
[358]	655
[1502]	656	zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps
	657	zv_sum (:,:) = zv_sum (:,:) + va_e (:,:)
	658	zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)
[367]	659
[1502]	660	! !* Time filter and swap
	661	! ! --------------------
	662	IF( jn == 1 ) THEN ! Swap only (1st Euler time step)
	663	zsshb_e(:,:) = sshn_e(:,:)
	664	zub_e (:,:) = zun_e (:,:)
	665	zvb_e (:,:) = zvn_e (:,:)
	666	sshn_e (:,:) = ssha_e(:,:)
	667	zun_e (:,:) = ua_e (:,:)
	668	zvn_e (:,:) = va_e (:,:)
	669	ELSE ! Swap + Filter
	670	zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
	671	zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)
	672	zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)
	673	sshn_e (:,:) = ssha_e(:,:)
	674	zun_e (:,:) = ua_e (:,:)
	675	zvn_e (:,:) = va_e (:,:)
[358]	676	ENDIF
	677
[1502]	678	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
	679	! ! ------------------
	680	DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points
	681	DO ji = 1, fs_jpim1 ! Vector opt.
[3294]	682	zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	683	& * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &
	684	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
	685	zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	686	& * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &
	687	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
[592]	688	END DO
	689	END DO
[1502]	690	CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions
	691	CALL lbc_lnk( zsshvn_e, 'V', 1. )
[1438]	692	!
[1502]	693	hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points
	694	hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
[3294]	695	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
	696	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
[1502]	697	!
[1438]	698	ENDIF
[358]	699	! ! ==================== !
	700	END DO ! end loop !
	701	! ! ==================== !
	702
[367]	703	#if defined key_obc
[1502]	704	IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ?????
[1241]	705	IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)
	706	IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)
	707	IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)
[367]	708	#endif
[358]	709
[1438]	710	! -----------------------------------------------------------------------------
[1502]	711	! Phase 3. update the general trend with the barotropic trend
[1438]	712	! -----------------------------------------------------------------------------
[1502]	713	!
	714	! !* Time average ==> after barotropic u, v, ssh
[3294]	715	zcoef = 1._wp / ( 2 * nn_baro + 1 )
[2528]	716	zu_sum(:,:) = zcoef * zu_sum (:,:)
	717	zv_sum(:,:) = zcoef * zv_sum (:,:)
[1502]	718	!
	719	! !* update the general momentum trend
[4376]	720	IF(ln_wad) THEN
	721	DO jk=1,jpkm1
[4380]	722	ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b * zwadfltu(:,:)
	723	va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b * zwadfltv(:,:)
[4376]	724	END DO
	725	ELSE
	726	DO jk=1,jpkm1
	727	ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b
	728	va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b
	729	END DO
	730	END IF
	731
[2528]	732	un_b (:,:) = zu_sum(:,:)
	733	vn_b (:,:) = zv_sum(:,:)
	734	sshn_b(:,:) = zcoef * zssh_sum(:,:)
[1502]	735	!
	736	! !* write time-spliting arrays in the restart
[508]	737	IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' )
	738	!
[3294]	739	CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv )
	740	CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e )
	741	CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum )
[4376]	742
[4380]	743	IF(ln_wad) CALL wrk_dealloc( jpi, jpj, zwadfltu, zwadfltv)
[1662]	744	!
[3294]	745	IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
[2715]	746	!
[508]	747	END SUBROUTINE dyn_spg_ts
	748
	749
	750	SUBROUTINE ts_rst( kt, cdrw )
	751	!!---------------------------------------------------------------------
	752	!! * ROUTINE ts_rst *
	753	!!
	754	!! ** Purpose : Read or write time-splitting arrays in restart file
	755	!!----------------------------------------------------------------------
	756	INTEGER , INTENT(in) :: kt ! ocean time-step
	757	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	758	!
	759	INTEGER :: ji, jk ! dummy loop indices
	760	!!----------------------------------------------------------------------
	761	!
	762	IF( TRIM(cdrw) == 'READ' ) THEN
[1502]	763	IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
	764	CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued
	765	CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop
[508]	766	ELSE
[3294]	767	un_b (:,:) = 0._wp
	768	vn_b (:,:) = 0._wp
[508]	769	! vertical sum
	770	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	771	DO jk = 1, jpkm1
	772	DO ji = 1, jpij
[1502]	773	un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
	774	vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
[508]	775	END DO
	776	END DO
	777	ELSE ! No vector opt.
	778	DO jk = 1, jpkm1
[1502]	779	un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
	780	vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
[508]	781	END DO
	782	ENDIF
[1502]	783	un_b (:,:) = un_b(:,:) * hur(:,:)
	784	vn_b (:,:) = vn_b(:,:) * hvr(:,:)
[508]	785	ENDIF
[2528]	786
	787	! Vertically integrated velocity (before)
	788	IF (neuler/=0) THEN
[3294]	789	ub_b (:,:) = 0._wp
	790	vb_b (:,:) = 0._wp
[2528]	791
	792	! vertical sum
	793	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	794	DO jk = 1, jpkm1
	795	DO ji = 1, jpij
	796	ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
	797	vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
	798	END DO
	799	END DO
	800	ELSE ! No vector opt.
	801	DO jk = 1, jpkm1
	802	ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
	803	vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
	804	END DO
	805	ENDIF
	806
	807	IF( lk_vvl ) THEN
[3294]	808	ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
	809	vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
[2528]	810	ELSE
	811	ub_b(:,:) = ub_b(:,:) * hur(:,:)
	812	vb_b(:,:) = vb_b(:,:) * hvr(:,:)
	813	ENDIF
	814	ELSE ! neuler==0
	815	ub_b (:,:) = un_b (:,:)
	816	vb_b (:,:) = vn_b (:,:)
	817	ENDIF
	818
[2145]	819	IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
[2715]	820	CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) ) ! filtered ssh
[2145]	821	ELSE
[2715]	822	sshn_b(:,:) = sshb(:,:) ! if not in restart set previous time mean to current baroclinic before value
[2145]	823	ENDIF
[508]	824	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
[2145]	825	CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity and ssh
	826	CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! issued from barotropic loop
	827	CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) ) !
[358]	828	ENDIF
[508]	829	!
	830	END SUBROUTINE ts_rst
	831
[358]	832	#else
	833	!!----------------------------------------------------------------------
	834	!! Default case : Empty module No standart free surface cst volume
	835	!!----------------------------------------------------------------------
	836	CONTAINS
[2715]	837	INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function
	838	dyn_spg_ts_alloc = 0
	839	END FUNCTION dyn_spg_ts_alloc
	840	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
	841	INTEGER, INTENT(in) :: kt
[358]	842	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
	843	END SUBROUTINE dyn_spg_ts
[2715]	844	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
[657]	845	INTEGER , INTENT(in) :: kt ! ocean time-step
	846	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	847	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
	848	END SUBROUTINE ts_rst
[358]	849	#endif
	850
	851	!!======================================================================
	852	END MODULE dynspg_ts

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

Download in other formats: