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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 4428

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

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

Download in other formats: