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 @ 2305

Last change on this file since 2305 was 2305, checked in by davestorkey, 14 years ago
Updates to BDY - DOCTOR standards
Property svn:keywords set to `Id`
File size: 33.6 KB

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

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

Download in other formats: