New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynspg_ts.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 1793

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

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

Download in other formats: