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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2382

Last change on this file since 2382 was 2338, checked in by rblod, 14 years ago
Fix pb with vvl and time splitting, see ticket #748
Property svn:keywords set to `Id`
File size: 34.4 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
[2236]	9	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
[1438]	10	!!---------------------------------------------------------------------
[575]	11	#if defined key_dynspg_ts \|\| defined key_esopa
[358]	12	!!----------------------------------------------------------------------
[455]	13	!! 'key_dynspg_ts' free surface cst volume with time splitting
[358]	14	!!----------------------------------------------------------------------
	15	!! dyn_spg_ts : compute surface pressure gradient trend using a time-
	16	!! splitting scheme and add to the general trend
[508]	17	!! ts_rst : read/write the time-splitting restart fields in the ocean restart file
[358]	18	!!----------------------------------------------------------------------
	19	USE oce ! ocean dynamics and tracers
	20	USE dom_oce ! ocean space and time domain
[888]	21	USE sbc_oce ! surface boundary condition: ocean
	22	USE dynspg_oce ! surface pressure gradient variables
[358]	23	USE phycst ! physical constants
[888]	24	USE domvvl ! variable volume
[1662]	25	USE zdfbfr ! bottom friction
[367]	26	USE obcdta ! open boundary condition data
	27	USE obcfla ! Flather open boundary condition
[358]	28	USE dynvor ! vorticity term
	29	USE obc_oce ! Lateral open boundary condition
[371]	30	USE obc_par ! open boundary condition parameters
[911]	31	USE bdy_oce ! unstructured open boundaries
	32	USE bdy_par ! unstructured open boundaries
	33	USE bdydta ! unstructured open boundaries
	34	USE bdydyn ! unstructured open boundaries
	35	USE bdytides ! tidal forcing at unstructured open boundaries.
[358]	36	USE lib_mpp ! distributed memory computing library
	37	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	38	USE prtctl ! Print control
	39	USE in_out_manager ! I/O manager
[508]	40	USE iom
	41	USE restart ! only for lrst_oce
[358]	42
	43	IMPLICIT NONE
	44	PRIVATE
	45
	46	PUBLIC dyn_spg_ts ! routine called by step.F90
[800]	47	PUBLIC ts_rst ! routine called by istate.F90
[358]	48
[1502]	49
[1438]	50	REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne ! triad of coriolis parameter
	51	REAL(wp), DIMENSION(jpi,jpj) :: ftsw, ftse ! (only used with een vorticity scheme)
[508]	52
[2338]	53	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: un_b, vn_b ! now averaged velocity
	54	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ub_b, vb_b ! before averaged velocity
[1502]	55
[2338]	56
[358]	57	!! * Substitutions
	58	# include "domzgr_substitute.h90"
	59	# include "vectopt_loop_substitute.h90"
[1438]	60	!!-------------------------------------------------------------------------
[2287]	61	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	62	!! $Id$
[2287]	63	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1438]	64	!!-------------------------------------------------------------------------
[358]	65
	66	CONTAINS
	67
	68	SUBROUTINE dyn_spg_ts( kt )
	69	!!----------------------------------------------------------------------
	70	!! * routine dyn_spg_ts *
	71	!!
	72	!! ** Purpose : Compute the now trend due to the surface pressure
	73	!! gradient in case of free surface formulation with time-splitting.
	74	!! Add it to the general trend of momentum equation.
	75	!!
	76	!! ** Method : Free surface formulation with time-splitting
	77	!! -1- Save the vertically integrated trend. This general trend is
	78	!! held constant over the barotropic integration.
	79	!! The Coriolis force is removed from the general trend as the
	80	!! surface gradient and the Coriolis force are updated within
	81	!! the barotropic integration.
[367]	82	!! -2- Barotropic loop : updates of sea surface height (ssha_e) and
[1502]	83	!! barotropic velocity (ua_e and va_e) through barotropic
[358]	84	!! momentum and continuity integration. Barotropic former
	85	!! variables are time averaging over the full barotropic cycle
[2338]	86	!! (= 2 * baroclinic time step) and saved in uX_b
	87	!! and vX_b (X specifying after, now or before).
[1438]	88	!! -3- The new general trend becomes :
[2338]	89	!! ua = ua - sum_k(ua)/H + ( un_b - ub_b )
[358]	90	!!
	91	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
	92	!!
[508]	93	!! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
[358]	94	!!---------------------------------------------------------------------
[1502]	95	INTEGER, INTENT(in) :: kt ! ocean time-step index
[1438]	96	!!
[1662]	97	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	98	INTEGER :: icycle ! temporary scalar
[1708]	99	INTEGER :: ikbu, ikbv ! temporary scalar
[1662]	100
	101	REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b ! temporary scalars
	102	REAL(wp) :: z1_8, zx1, zy1 ! - -
	103	REAL(wp) :: z1_4, zx2, zy2 ! - -
	104	REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - -
	105	REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - -
	106	!!
[1502]	107	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv, zsshb_e
[1662]	108	!!
[1708]	109	REAL(wp), DIMENSION(jpi,jpj) :: zbfru , zbfrv ! 2D workspace
	110	!!
[1502]	111	REAL(wp), DIMENSION(jpi,jpj) :: zsshun_e, zsshvn_e ! 2D workspace
	112	!!
[2338]	113	REAL(wp), DIMENSION(jpi,jpj) :: zcu, zwx, zua, zun ! 2D workspace
	114	REAL(wp), DIMENSION(jpi,jpj) :: zcv, zwy, zva, zvn ! - -
[1502]	115	REAL(wp), DIMENSION(jpi,jpj) :: zun_e, zub_e, zu_sum ! 2D workspace
	116	REAL(wp), DIMENSION(jpi,jpj) :: zvn_e, zvb_e, zv_sum ! - -
	117	REAL(wp), DIMENSION(jpi,jpj) :: zssh_sum ! - -
[358]	118	!!----------------------------------------------------------------------
	119
[1502]	120	IF( kt == nit000 ) THEN !* initialisation
[508]	121	!
[358]	122	IF(lwp) WRITE(numout,*)
	123	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
	124	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
[1241]	125	IF(lwp) WRITE(numout,) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2nn_baro
[1502]	126	!
[1708]	127	CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b
[1502]	128	!
	129	ua_e (:,:) = un_b (:,:)
	130	va_e (:,:) = vn_b (:,:)
	131	hu_e (:,:) = hu (:,:)
	132	hv_e (:,:) = hv (:,:)
	133	hur_e (:,:) = hur (:,:)
	134	hvr_e (:,:) = hvr (:,:)
[358]	135	IF( ln_dynvor_een ) THEN
[508]	136	ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0
[358]	137	DO jj = 2, jpj
[1708]	138	DO ji = fs_2, jpi ! vector opt.
[508]	139	ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3.
	140	ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3.
	141	ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3.
	142	ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3.
[358]	143	END DO
	144	END DO
	145	ENDIF
[508]	146	!
	147	ENDIF
[358]	148
[1502]	149	! !* Local constant initialization
[358]	150	z2dt_b = 2.0 * rdt ! baroclinic time step
[1502]	151	z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates
	152	z1_4 = 0.5 * 0.5
[1739]	153	zraur = 1. / rau0 ! 1 / volumic mass
[1502]	154	!
	155	zhdiv(:,:) = 0.e0 ! barotropic divergence
[1662]	156	zu_sld = 0.e0 ; zu_asp = 0.e0 ! tides trends (lk_tide=F)
	157	zv_sld = 0.e0 ; zv_asp = 0.e0
[1438]	158
[358]	159	! -----------------------------------------------------------------------------
	160	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
	161	! -----------------------------------------------------------------------------
[1502]	162	!
	163	! !* e3d/dt(Ua), e3Ub, e3*Vn (Vertically integrated)
	164	! ! --------------------------
[2338]	165	zua(:,:) = 0.e0 ; zun(:,:) = 0.e0
	166	zva(:,:) = 0.e0 ; zvn(:,:) = 0.e0
[1502]	167	!
	168	DO jk = 1, jpkm1
	169	#if defined key_vectopt_loop
	170	DO jj = 1, 1 !Vector opt. => forced unrolling
[358]	171	DO ji = 1, jpij
[1502]	172	#else
	173	DO jj = 1, jpj
	174	DO ji = 1, jpi
	175	#endif
	176	! ! now trend
	177	zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
	178	zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
	179	! ! now velocity
	180	zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk)
	181	zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk)
[358]	182	END DO
	183	END DO
[1502]	184	END DO
	185
	186	! !* baroclinic momentum trend (remove the vertical mean trend)
	187	DO jk = 1, jpkm1 ! --------------------------
	188	DO jj = 2, jpjm1
	189	DO ji = fs_2, fs_jpim1 ! vector opt.
	190	ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
	191	va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
	192	END DO
[358]	193	END DO
[1502]	194	END DO
[358]	195
[1502]	196	! !* barotropic Coriolis trends * H (vorticity scheme dependent)
	197	! ! ---------------------------====
	198	zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport
	199	zwy(:,:) = zvn(:,:) * e1v(:,:)
	200	!
[358]	201	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
	202	DO jj = 2, jpjm1
	203	DO ji = fs_2, fs_jpim1 ! vector opt.
	204	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	205	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	206	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	207	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	208	! energy conserving formulation for planetary vorticity term
[1502]	209	zcu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )
	210	zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )
[358]	211	END DO
	212	END DO
[508]	213	!
[358]	214	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
	215	DO jj = 2, jpjm1
	216	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	217	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	218	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[358]	219	zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) )
	220	zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) )
	221	END DO
	222	END DO
[508]	223	!
[358]	224	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
	225	DO jj = 2, jpjm1
	226	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	227	zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	228	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
	229	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) &
	230	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) )
[358]	231	END DO
	232	END DO
[508]	233	!
[358]	234	ENDIF
	235
[1502]	236	! !* Right-Hand-Side of the barotropic momentum equation
	237	! ! ----------------------------------------------------
	238	IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient
	239	DO jj = 2, jpjm1
[358]	240	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	241	zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) &
	242	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)
	243	zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) &
	244	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)
[358]	245	END DO
	246	END DO
[1502]	247	ENDIF
[358]	248
[1502]	249	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
[358]	250	DO ji = fs_2, fs_jpim1
	251	zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
	252	zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
	253	END DO
	254	END DO
	255
[1708]	256
	257	! ! Remove barotropic contribution of bottom friction
	258	! ! from the barotropic transport trend
	259	zcoef = -1. / z2dt_b
	260	# if defined key_vectopt_loop
	261	DO jj = 1, 1
[1779]	262	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
[1708]	263	# else
	264	DO jj = 2, jpjm1
	265	DO ji = 2, jpim1
	266	# endif
	267	ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )
	268	ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
	269	!
	270	! Apply stability criteria for bottom friction
	271	!RBbug for vvl and external mode we may need to
	272	! use varying fse3
	273	zbfru (ji,jj) = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbu)*zcoef )
	274	zbfrv (ji,jj) = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbv)*zcoef )
	275	END DO
	276	END DO
	277
[1662]	278	IF( lk_vvl ) THEN
[1708]	279	DO jj = 2, jpjm1
	280	DO ji = fs_2, fs_jpim1 ! vector opt.
[2338]	281	zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) &
[1708]	282	& / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
[2338]	283	zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) &
[1708]	284	& / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
	285	END DO
	286	END DO
[1662]	287	ELSE
[1708]	288	DO jj = 2, jpjm1
	289	DO ji = fs_2, fs_jpim1 ! vector opt.
[2338]	290	zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
	291	zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
[1708]	292	END DO
	293	END DO
[1662]	294	ENDIF
	295
[1502]	296	! !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
	297	! ! --------------------------
	298	zua(:,:) = zua(:,:) * hur(:,:)
	299	zva(:,:) = zva(:,:) * hvr(:,:)
	300	!
	301
[358]	302	! -----------------------------------------------------------------------
	303	! Phase 2 : Integration of the barotropic equations with time splitting
	304	! -----------------------------------------------------------------------
[1502]	305	!
	306	! ! ==================== !
	307	! ! Initialisations !
	308	! ! ==================== !
	309	icycle = 2 * nn_baro ! Number of barotropic sub time-step
	310
	311	! ! Start from NOW field
	312	hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points
	313	hv_e (:,:) = hv (:,:)
	314	hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points
	315	hvr_e (:,:) = hvr (:,:)
[1662]	316	!RBbug zsshb_e(:,:) = sshn (:,:)
	317	zsshb_e(:,:) = sshn_b(:,:) ! sea surface height (before and now)
[1502]	318	sshn_e (:,:) = sshn (:,:)
	319
	320	zun_e (:,:) = un_b (:,:) ! barotropic velocity (external)
	321	zvn_e (:,:) = vn_b (:,:)
	322	zub_e (:,:) = un_b (:,:)
	323	zvb_e (:,:) = vn_b (:,:)
[358]	324
[1502]	325	zu_sum (:,:) = un_b (:,:) ! summation
	326	zv_sum (:,:) = vn_b (:,:)
	327	zssh_sum(:,:) = sshn (:,:)
[358]	328
[1502]	329	#if defined key_obc
[367]	330	! set ssh corrections to 0
	331	! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
	332	IF( lp_obc_east ) sshfoe_b(:,:) = 0.e0
	333	IF( lp_obc_west ) sshfow_b(:,:) = 0.e0
	334	IF( lp_obc_south ) sshfos_b(:,:) = 0.e0
	335	IF( lp_obc_north ) sshfon_b(:,:) = 0.e0
	336	#endif
	337
[1502]	338	! ! ==================== !
	339	DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1)
	340	! ! ==================== !
[1241]	341	z2dt_e = 2. * ( rdt / nn_baro )
[1502]	342	IF( jn == 1 ) z2dt_e = rdt / nn_baro
[358]	343
[1502]	344	! !* Update the forcing (OBC, BDY and tides)
	345	! ! ------------------
[2338]	346	IF( lk_obc ) CALL obc_dta_bt ( kt, jn )
[2305]	347	IF( lk_bdy ) CALL bdy_dta_fla( kt, jn+1, icycle )
[367]	348
[1502]	349	! !* after ssh_e
	350	! ! -----------
	351	DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports
[358]	352	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	353	zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) &
	354	& - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) &
	355	& + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) &
	356	& - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
[358]	357	END DO
	358	END DO
[1502]	359	!
[358]	360	#if defined key_obc
[1502]	361	! ! OBC : zhdiv must be zero behind the open boundary
	362	!! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
	363	IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east
	364	IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west
[367]	365	IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north
[1502]	366	IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south
[358]	367	#endif
[1170]	368	#if defined key_bdy
[1502]	369	zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask
[1170]	370	#endif
[1502]	371	!
	372	DO jj = 2, jpjm1 ! leap-frog on ssh_e
[358]	373	DO ji = fs_2, fs_jpim1 ! vector opt.
[2000]	374	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]	375	END DO
	376	END DO
	377
[1502]	378	! !* after barotropic velocities (vorticity scheme dependent)
	379	! ! ---------------------------
	380	zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport
	381	zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
	382	!
	383	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
[358]	384	DO jj = 2, jpjm1
	385	DO ji = fs_2, fs_jpim1 ! vector opt.
	386	! surface pressure gradient
[592]	387	IF( lk_vvl) THEN
[1662]	388	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	389	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	390	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	391	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	392	ELSE
[1662]	393	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	394	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	395	ENDIF
[358]	396	! energy conserving formulation for planetary vorticity term
	397	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	398	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	399	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	400	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[1662]	401	zu_cor = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
	402	zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
	403	! after velocities with implicit bottom friction
	404	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) &
	405	& / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
	406	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) &
	407	& / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
[358]	408	END DO
	409	END DO
[508]	410	!
[1502]	411	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
[358]	412	DO jj = 2, jpjm1
	413	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	414	! surface pressure gradient
[592]	415	IF( lk_vvl) THEN
[1662]	416	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	417	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	418	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	419	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	420	ELSE
[1662]	421	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	422	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	423	ENDIF
[358]	424	! enstrophy conserving formulation for planetary vorticity term
[1502]	425	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	426	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[1662]	427	zu_cor = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
	428	zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj)
	429	! after velocities with implicit bottom friction
	430	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) &
	431	& / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
	432	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) &
	433	& / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
[358]	434	END DO
	435	END DO
[508]	436	!
[1502]	437	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
[358]	438	DO jj = 2, jpjm1
	439	DO ji = fs_2, fs_jpim1 ! vector opt.
	440	! surface pressure gradient
[592]	441	IF( lk_vvl) THEN
[1662]	442	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	443	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	444	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	445	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	446	ELSE
[1662]	447	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	448	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	449	ENDIF
[358]	450	! energy/enstrophy conserving formulation for planetary vorticity term
[1662]	451	zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
[1502]	452	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
[1662]	453	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]	454	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj)
[1662]	455	! after velocities with implicit bottom friction
	456	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) &
	457	& / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
	458	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) &
	459	& / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
[358]	460	END DO
	461	END DO
[508]	462	!
[358]	463	ENDIF
[1502]	464	! !* domain lateral boundary
	465	! ! -----------------------
	466	! ! Flather's boundary condition for the barotropic loop :
	467	! ! - Update sea surface height on each open boundary
	468	! ! - Correct the velocity
[358]	469
[2338]	470	IF( lk_obc ) CALL obc_fla_ts
	471	IF( lk_bdy .OR. ln_tides ) CALL bdy_dyn_fla( sshn_e )
[1502]	472	!
	473	CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries
	474	CALL lbc_lnk( va_e , 'V', -1. )
	475	CALL lbc_lnk( ssha_e, 'T', 1. )
[358]	476
[1502]	477	zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps
	478	zv_sum (:,:) = zv_sum (:,:) + va_e (:,:)
	479	zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)
[367]	480
[1502]	481	! !* Time filter and swap
	482	! ! --------------------
	483	IF( jn == 1 ) THEN ! Swap only (1st Euler time step)
	484	zsshb_e(:,:) = sshn_e(:,:)
	485	zub_e (:,:) = zun_e (:,:)
	486	zvb_e (:,:) = zvn_e (:,:)
	487	sshn_e (:,:) = ssha_e(:,:)
	488	zun_e (:,:) = ua_e (:,:)
	489	zvn_e (:,:) = va_e (:,:)
	490	ELSE ! Swap + Filter
	491	zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
	492	zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)
	493	zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)
	494	sshn_e (:,:) = ssha_e(:,:)
	495	zun_e (:,:) = ua_e (:,:)
	496	zvn_e (:,:) = va_e (:,:)
[358]	497	ENDIF
	498
[1502]	499	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
	500	! ! ------------------
	501	DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points
	502	DO ji = 1, fs_jpim1 ! Vector opt.
	503	zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	504	& * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &
	505	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
	506	zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	507	& * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &
	508	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
[592]	509	END DO
	510	END DO
[1502]	511	CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions
	512	CALL lbc_lnk( zsshvn_e, 'V', 1. )
[1438]	513	!
[1502]	514	hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points
	515	hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
	516	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
	517	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
	518	!
[1438]	519	ENDIF
[358]	520	! ! ==================== !
	521	END DO ! end loop !
	522	! ! ==================== !
	523
[367]	524	#if defined key_obc
[1502]	525	IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ?????
[1241]	526	IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)
	527	IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)
	528	IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)
[367]	529	#endif
[358]	530
[1438]	531	! -----------------------------------------------------------------------------
[1502]	532	! Phase 3. update the general trend with the barotropic trend
[1438]	533	! -----------------------------------------------------------------------------
[1502]	534	!
	535	! !* Time average ==> after barotropic u, v, ssh
	536	zcoef = 1.e0 / ( 2 * nn_baro + 1 )
[2338]	537	zu_sum(:,:) = zcoef * zu_sum (:,:)
	538	zv_sum(:,:) = zcoef * zv_sum (:,:)
[1502]	539	!
	540	! !* update the general momentum trend
[358]	541	DO jk=1,jpkm1
[2338]	542	ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
	543	va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
[358]	544	END DO
[2338]	545	ub_b (:,:) = un_b(:,:)
	546	vb_b (:,:) = vn_b(:,:)
	547	un_b (:,:) = zu_sum(:,:)
	548	vn_b (:,:) = zv_sum(:,:)
	549	sshn_b(:,:) = zcoef * zssh_sum(:,:)
[1502]	550	!
	551	! !* write time-spliting arrays in the restart
[508]	552	IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' )
	553	!
[1662]	554	!
[508]	555	END SUBROUTINE dyn_spg_ts
	556
	557
	558	SUBROUTINE ts_rst( kt, cdrw )
	559	!!---------------------------------------------------------------------
	560	!! * ROUTINE ts_rst *
	561	!!
	562	!! ** Purpose : Read or write time-splitting arrays in restart file
	563	!!----------------------------------------------------------------------
	564	INTEGER , INTENT(in) :: kt ! ocean time-step
	565	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	566	!
	567	INTEGER :: ji, jk ! dummy loop indices
	568	!!----------------------------------------------------------------------
	569	!
	570	IF( TRIM(cdrw) == 'READ' ) THEN
[1502]	571	IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
	572	CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued
	573	CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop
[508]	574	ELSE
[1502]	575	un_b (:,:) = 0.e0
	576	vn_b (:,:) = 0.e0
[508]	577	! vertical sum
	578	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	579	DO jk = 1, jpkm1
	580	DO ji = 1, jpij
[1502]	581	un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
	582	vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
[508]	583	END DO
	584	END DO
	585	ELSE ! No vector opt.
	586	DO jk = 1, jpkm1
[1502]	587	un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
	588	vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
[508]	589	END DO
	590	ENDIF
[1502]	591	un_b (:,:) = un_b(:,:) * hur(:,:)
	592	vn_b (:,:) = vn_b(:,:) * hvr(:,:)
[508]	593	ENDIF
[2338]	594
	595	! Vertically integrated velocity (before)
	596	IF (neuler/=0) THEN
	597	ub_b (:,:) = 0.e0
	598	vb_b (:,:) = 0.e0
	599
	600	! vertical sum
	601	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	602	DO jk = 1, jpkm1
	603	DO ji = 1, jpij
	604	ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
	605	vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
	606	END DO
	607	END DO
	608	ELSE ! No vector opt.
	609	DO jk = 1, jpkm1
	610	ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
	611	vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
	612	END DO
	613	ENDIF
	614
	615	IF( lk_vvl ) THEN
	616	ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
	617	vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
	618	ELSE
	619	ub_b(:,:) = ub_b(:,:) * hur(:,:)
	620	vb_b(:,:) = vb_b(:,:) * hvr(:,:)
	621	ENDIF
	622	ELSE ! neuler==0
	623	ub_b (:,:) = un_b (:,:)
	624	vb_b (:,:) = vn_b (:,:)
	625	ENDIF
	626
[2200]	627	IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
	628	CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) ) ! filtered extrenal ssh
	629	ELSE
	630	sshn_b(:,:)=sshb(:,:) ! if not in restart set previous time mean to current baroclinic before value
	631	ENDIF
[508]	632	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
[2200]	633	CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity and ssh
	634	CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! issued from barotropic loop
	635	CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) ) !
[358]	636	ENDIF
[508]	637	!
	638	END SUBROUTINE ts_rst
	639
[358]	640	#else
	641	!!----------------------------------------------------------------------
	642	!! Default case : Empty module No standart free surface cst volume
	643	!!----------------------------------------------------------------------
	644	CONTAINS
	645	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
	646	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
	647	END SUBROUTINE dyn_spg_ts
[657]	648	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
	649	INTEGER , INTENT(in) :: kt ! ocean time-step
	650	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	651	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
	652	END SUBROUTINE ts_rst
[358]	653	#endif
	654
	655	!!======================================================================
	656	END MODULE dynspg_ts

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

Download in other formats: