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

Context Navigation

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3723

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

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

Download in other formats: