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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3649

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

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

Download in other formats: