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

Last change on this file since 800 was 800, checked in by rblod, 16 years ago
Correct a mistake in jki routine suppression, see ticket #49
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 33.7 KB

Rev	Line
[358]	1	MODULE dynspg_ts
	2	!!======================================================================
	3	!! * MODULE dynspg_ts *
	4	!! Ocean dynamics: surface pressure gradient trend
	5	!!======================================================================
[508]	6	!! History : 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code
	7	!! " " ! 05-11 (V. Garnier, G. Madec) optimization
	8	!! 9.0 ! 06-08 (S. Masson) distributed restart using iom
[788]	9	!! 9.0 ! 08-01 (R. Benshila) change averaging method
[508]	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	!!----------------------------------------------------------------------
[508]	15	!!----------------------------------------------------------------------
[358]	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	!! * Modules used
	21	USE oce ! ocean dynamics and tracers
	22	USE dom_oce ! ocean space and time domain
	23	USE phycst ! physical constants
[719]	24	USE ocesbc ! ocean surface boundary condition
[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
[358]	30	USE lib_mpp ! distributed memory computing library
	31	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	32	USE prtctl ! Print control
[719]	33	USE dynspg_oce ! surface pressure gradient variables
[358]	34	USE in_out_manager ! I/O manager
[508]	35	USE iom
	36	USE restart ! only for lrst_oce
[719]	37	USE domvvl ! variable volume
[358]	38
	39	IMPLICIT NONE
	40	PRIVATE
	41
	42	PUBLIC dyn_spg_ts ! routine called by step.F90
[800]	43	PUBLIC ts_rst ! routine called by istate.F90
[358]	44
[508]	45	REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne, & ! triad of coriolis parameter
	46	& ftsw, ftse ! (only used with een vorticity scheme)
	47
	48
[358]	49	!! * Substitutions
	50	# include "domzgr_substitute.h90"
	51	# include "vectopt_loop_substitute.h90"
	52	!!----------------------------------------------------------------------
[359]	53	!! OPA 9.0 , LOCEAN-IPSL (2005)
[719]	54	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/DYN/dynspg_ts.F90,v 1.16 2007/06/05 10:38:27 opalod Exp $
[359]	55	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
[358]	56	!!----------------------------------------------------------------------
	57
	58	CONTAINS
	59
	60	SUBROUTINE dyn_spg_ts( kt )
	61	!!----------------------------------------------------------------------
	62	!! * routine dyn_spg_ts *
	63	!!
	64	!! ** Purpose : Compute the now trend due to the surface pressure
	65	!! gradient in case of free surface formulation with time-splitting.
	66	!! Add it to the general trend of momentum equation.
	67	!! Compute the free surface.
	68	!!
	69	!! ** Method : Free surface formulation with time-splitting
	70	!! -1- Save the vertically integrated trend. This general trend is
	71	!! held constant over the barotropic integration.
	72	!! The Coriolis force is removed from the general trend as the
	73	!! surface gradient and the Coriolis force are updated within
	74	!! the barotropic integration.
[367]	75	!! -2- Barotropic loop : updates of sea surface height (ssha_e) and
	76	!! barotropic transports (ua_e and va_e) through barotropic
[358]	77	!! momentum and continuity integration. Barotropic former
	78	!! variables are time averaging over the full barotropic cycle
	79	!! (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b
	80	!! and zvX_b (X specifying after, now or before).
	81	!! -3- Update of sea surface height from time averaged barotropic
	82	!! variables.
	83	!! - apply lateral boundary conditions on sshn.
	84	!! -4- The new general trend becomes :
	85	!! ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H
	86	!!
	87	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
	88	!!
[508]	89	!! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
[358]	90	!!---------------------------------------------------------------------
	91	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	92
	93	!! * Local declarations
	94	INTEGER :: ji, jj, jk, jit ! dummy loop indices
[788]	95	INTEGER :: icycle, ibaro ! temporary scalar
[358]	96	REAL(wp) :: &
	97	zraur, zcoef, z2dt_e, z2dt_b, zfac25, & ! temporary scalars
	98	zfact1, zspgu, zcubt, zx1, zy1, & ! " "
	99	zfact2, zspgv, zcvbt, zx2, zy2 ! " "
	100	REAL(wp), DIMENSION(jpi,jpj) :: &
	101	zcu, zcv, zwx, zwy, zhdiv, & ! temporary arrays
	102	zua, zva, zub, zvb, & ! " "
	103	zssha_b, zua_b, zva_b, & ! " "
	104	zsshb_e, zub_e, zvb_e, & ! " "
[367]	105	zun_e, zvn_e ! " "
[592]	106	!! Variable volume
	107	REAL(wp), DIMENSION(jpi,jpj) :: &
	108	zspgu_1, zspgv_1, zsshun_e, zsshvn_e, hu_e, hv_e ! 2D workspace
	109	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3un_e, zfse3vn_e ! 3D workspace
[358]	110	!!----------------------------------------------------------------------
	111
	112	! Arrays initialization
	113	! ---------------------
[367]	114	zua_b(:,:) = 0.e0 ; zub_e(:,:) = 0.e0 ; zun_e(:,:) = 0.e0
	115	zva_b(:,:) = 0.e0 ; zvb_e(:,:) = 0.e0 ; zvn_e(:,:) = 0.e0
[363]	116	zhdiv(:,:) = 0.e0
[358]	117
	118
	119	IF( kt == nit000 ) THEN
[508]	120	!
[358]	121	IF(lwp) WRITE(numout,*)
	122	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
	123	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
	124	IF(lwp) WRITE(numout,) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', FLOOR( 2rdt/rdtbt )
	125
[508]	126	CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields:
	127	! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b
	128
[367]	129	ssha_e(:,:) = sshn(:,:)
	130	ua_e(:,:) = un_b(:,:)
	131	va_e(:,:) = vn_b(:,:)
[358]	132
	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
	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
[592]	147	! Substract the surface pressure gradient from the momentum
	148	! ---------------------------------------------------------
	149	IF( lk_vvl ) THEN ! Variable volume
	150
	151	! 0. Local initialization
	152	! -----------------------
	153	zspgu_1(:,:) = 0.e0
	154	zspgv_1(:,:) = 0.e0
	155
	156	! 1. Surface pressure gradient (now)
	157	! ----------------------------
	158	DO jj = 2, jpjm1
	159	DO ji = fs_2, fs_jpim1 ! vector opt.
	160	zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) &
	161	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e1u(ji,jj)
	162	zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) &
	163	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj)
	164	END DO
	165	END DO
	166
	167	! 2. Substract the surface pressure trend from the general trend
	168	! ------------------------------------------------------
	169	DO jk = 1, jpkm1
	170	DO jj = 2, jpjm1
	171	DO ji = fs_2, fs_jpim1 ! vector opt.
	172	ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj)
	173	va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj)
	174	END DO
	175	END DO
	176	END DO
	177
	178	ENDIF
	179
[358]	180	! Local constant initialization
	181	! --------------------------------
	182	z2dt_b = 2.0 * rdt ! baroclinic time step
	183	IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt
[374]	184	zfact1 = 0.5 * 0.25 ! coefficient for vorticity estimates
[358]	185	zfact2 = 0.5 * 0.5
[374]	186	zraur = 1. / rauw ! 1 / volumic mass of pure water
[358]	187
	188	! -----------------------------------------------------------------------------
	189	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
	190	! -----------------------------------------------------------------------------
	191
	192	! Vertically integrated quantities
	193	! --------------------------------
	194	zua(:,:) = 0.e0
	195	zva(:,:) = 0.e0
	196	zub(:,:) = 0.e0
	197	zvb(:,:) = 0.e0
	198	zwx(:,:) = 0.e0
	199	zwy(:,:) = 0.e0
	200
	201	! vertical sum
	202	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	203	DO jk = 1, jpkm1
	204	DO ji = 1, jpij
	205	! ! Vertically integrated momentum trends
	206	zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk)
	207	zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk)
	208	! ! Vertically integrated transports (before)
	209	zub(ji,1) = zub(ji,1) + fse3u(ji,1,jk) * ub(ji,1,jk)
	210	zvb(ji,1) = zvb(ji,1) + fse3v(ji,1,jk) * vb(ji,1,jk)
	211	! ! Planetary vorticity transport fluxes (now)
	212	zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk)
	213	zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk)
	214	END DO
	215	END DO
	216	ELSE ! No vector opt.
	217	DO jk = 1, jpkm1
	218	! ! Vertically integrated momentum trends
	219	zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
	220	zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
	221	! ! Vertically integrated transports (before)
	222	zub(:,:) = zub(:,:) + fse3u(:,:,jk) * ub(:,:,jk)
	223	zvb(:,:) = zvb(:,:) + fse3v(:,:,jk) * vb(:,:,jk)
	224	! ! Planetary vorticity (now)
	225	zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	226	zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	227	END DO
	228	ENDIF
	229
	230	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
	231	DO jj = 2, jpjm1
	232	DO ji = fs_2, fs_jpim1 ! vector opt.
	233	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	234	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	235	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	236	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	237	! energy conserving formulation for planetary vorticity term
	238	zcu(ji,jj) = zfact2 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )
	239	zcv(ji,jj) =-zfact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )
	240	END DO
	241	END DO
[508]	242	!
[358]	243	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
	244	DO jj = 2, jpjm1
	245	DO ji = fs_2, fs_jpim1 ! vector opt.
	246	zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
	247	+ zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	248	zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
	249	+ zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	250	zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) )
	251	zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) )
	252	END DO
	253	END DO
[508]	254	!
[358]	255	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
	256	zfac25 = 0.25
	257	DO jj = 2, jpjm1
	258	DO ji = fs_2, fs_jpim1 ! vector opt.
	259	zcu(ji,jj) = + zfac25 / e1u(ji,jj) &
[553]	260	& * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	261	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
[358]	262	zcv(ji,jj) = - zfac25 / e2v(ji,jj) &
[553]	263	& * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &
	264	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) )
[358]	265	END DO
	266	END DO
[508]	267	!
[358]	268	ENDIF
	269
	270
	271	! Remove barotropic trend from general momentum trend
	272	! ---------------------------------------------------
	273	DO jk = 1 , jpkm1
	274	DO jj = 2, jpjm1
	275	DO ji = fs_2, fs_jpim1 ! vector opt.
	276	ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
	277	va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
	278	END DO
	279	END DO
	280	END DO
	281
	282	! Remove coriolis term from barotropic trend
	283	! ------------------------------------------
	284	DO jj = 2, jpjm1
	285	DO ji = fs_2, fs_jpim1
	286	zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
	287	zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
	288	END DO
	289	END DO
	290
	291	! -----------------------------------------------------------------------
	292	! Phase 2 : Integration of the barotropic equations with time splitting
	293	! -----------------------------------------------------------------------
	294
	295	! Initialisations
	296	!----------------
	297	! Number of iteration of the barotropic loop
[788]	298	ibaro = FLOOR( rdt / rdtbt )
	299	icycle = 3 /2 * ibaro
[358]	300
	301	! variables for the barotropic equations
	302	zsshb_e(:,:) = sshn_b(:,:) ! (barotropic) sea surface height (before and now)
[367]	303	sshn_e (:,:) = sshn_b(:,:)
	304	zub_e (:,:) = un_b (:,:) ! barotropic transports issued from the barotropic equations (before and now)
	305	zvb_e (:,:) = vn_b (:,:)
	306	zun_e (:,:) = un_b (:,:)
	307	zvn_e (:,:) = vn_b (:,:)
[788]	308	zssha_b(:,:) = 0.e0
	309	zua_b (:,:) = 0.e0
	310	zva_b (:,:) = 0.e0
[592]	311	IF( lk_vvl ) THEN
	312	zsshun_e(:,:) = sshu (:,:) ! (barotropic) sea surface height (now) at u-point
	313	zsshvn_e(:,:) = sshv (:,:) ! (barotropic) sea surface height (now) at v-point
	314	hu_e(:,:) = hu (:,:) ! (barotropic) ocean depth at u-point
	315	hv_e(:,:) = hv (:,:) ! (barotropic) ocean depth at v-point
	316	zfse3un_e(:,:,:) = fse3u(:,:,:) ! (barotropic) scale factors at u-point
	317	zfse3un_e(:,:,:) = fse3v(:,:,:) ! (barotropic) scale factors at v-point
	318	ENDIF
[358]	319
[367]	320	! set ssh corrections to 0
	321	! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
	322	#if defined key_obc
	323	IF( lp_obc_east ) sshfoe_b(:,:) = 0.e0
	324	IF( lp_obc_west ) sshfow_b(:,:) = 0.e0
	325	IF( lp_obc_south ) sshfos_b(:,:) = 0.e0
	326	IF( lp_obc_north ) sshfon_b(:,:) = 0.e0
	327	#endif
	328
[358]	329	! Barotropic integration over 2 baroclinic time steps
	330	! ---------------------------------------------------
	331
	332	! ! ==================== !
	333	DO jit = 1, icycle ! sub-time-step loop !
	334	! ! ==================== !
	335	z2dt_e = 2. * rdtbt
	336	IF ( jit == 1 ) z2dt_e = rdtbt
	337
[367]	338	! Time interpolation of open boundary condition data
	339	IF( lk_obc ) CALL obc_dta_bt( kt, jit )
	340
[358]	341	! Horizontal divergence of barotropic transports
	342	!--------------------------------------------------
	343	DO jj = 2, jpjm1
	344	DO ji = fs_2, fs_jpim1 ! vector opt.
	345	zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun_e(ji ,jj) &
	346	& -e2u(ji-1,jj ) * zun_e(ji-1,jj) &
	347	& +e1v(ji ,jj ) * zvn_e(ji ,jj) &
	348	& -e1v(ji ,jj-1) * zvn_e(ji ,jj-1) ) &
	349	& / (e1t(ji,jj)*e2t(ji,jj))
	350	END DO
	351	END DO
	352
	353	#if defined key_obc
	354	! open boundaries (div must be zero behind the open boundary)
	355	! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
[367]	356	IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east
	357	IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west
	358	IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north
	359	IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south
[358]	360	#endif
	361
	362	! Sea surface height from the barotropic system
	363	!----------------------------------------------
	364	DO jj = 2, jpjm1
	365	DO ji = fs_2, fs_jpim1 ! vector opt.
[367]	366	ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) &
[358]	367	& + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
	368	END DO
	369	END DO
	370
	371	! evolution of the barotropic transport ( following the vorticity scheme used)
	372	! ----------------------------------------------------------------------------
	373	zwx(:,:) = e2u(:,:) * zun_e(:,:)
	374	zwy(:,:) = e1v(:,:) * zvn_e(:,:)
	375
	376	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
	377	DO jj = 2, jpjm1
	378	DO ji = fs_2, fs_jpim1 ! vector opt.
	379	! surface pressure gradient
[592]	380	IF( lk_vvl) THEN
	381	zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) &
	382	& - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)
	383	zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &
	384	& - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)
	385	ELSE
	386	zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
	387	zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
	388	ENDIF
[358]	389	! energy conserving formulation for planetary vorticity term
	390	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	391	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	392	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	393	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	394	zcubt = zfact2 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )
	395	zcvbt =-zfact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )
	396	! after transports
[367]	397	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
	398	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
[358]	399	END DO
	400	END DO
[508]	401	!
[358]	402	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
	403	DO jj = 2, jpjm1
	404	DO ji = fs_2, fs_jpim1 ! vector opt.
	405	! surface pressure gradient
[592]	406	IF( lk_vvl) THEN
	407	zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) &
	408	& - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)
	409	zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &
	410	& - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)
	411	ELSE
	412	zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
	413	zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
	414	ENDIF
[358]	415	! enstrophy conserving formulation for planetary vorticity term
	416	zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
	417	+ zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	418	zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
	419	+ zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	420	zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) )
	421	zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) )
	422	! after transports
[367]	423	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
	424	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
[358]	425	END DO
	426	END DO
[508]	427	!
[358]	428	ELSEIF ( ln_dynvor_een ) THEN ! energy and enstrophy conserving scheme
	429	zfac25 = 0.25
	430	DO jj = 2, jpjm1
	431	DO ji = fs_2, fs_jpim1 ! vector opt.
	432	! surface pressure gradient
[592]	433	IF( lk_vvl) THEN
	434	zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) &
	435	& - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)
	436	zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &
	437	& - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)
	438	ELSE
	439	zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
	440	zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
	441	ENDIF
[358]	442	! energy/enstrophy conserving formulation for planetary vorticity term
[553]	443	zcubt = + zfac25 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	444	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
	445	zcvbt = - zfac25 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &
	446	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) )
[358]	447	! after transports
[367]	448	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
	449	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
[358]	450	END DO
	451	END DO
[508]	452	!
[358]	453	ENDIF
	454
[367]	455	! Flather's boundary condition for the barotropic loop :
	456	! - Update sea surface height on each open boundary
	457	! - Correct the barotropic transports
[371]	458	IF( lk_obc ) CALL obc_fla_ts
[358]	459
[367]	460
	461	! ... Boundary conditions on ua_e, va_e, ssha_e
	462	CALL lbc_lnk( ua_e , 'U', -1. )
	463	CALL lbc_lnk( va_e , 'V', -1. )
	464	CALL lbc_lnk( ssha_e, 'T', 1. )
	465
[358]	466	! temporal sum
	467	!-------------
[788]	468	IF( jit >= ibaro/2 ) THEN
	469	zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:)
	470	zua_b (:,:) = zua_b (:,:) + ua_e (:,:)
	471	zva_b (:,:) = zva_b (:,:) + va_e (:,:)
	472	ENDIF
[358]	473
	474	! Time filter and swap of dynamics arrays
	475	! ---------------------------------------
[592]	476	IF( neuler == 0 .AND. jit == 1 ) THEN ! Euler (forward) time stepping
[367]	477	zsshb_e(:,:) = sshn_e(:,:)
	478	zub_e (:,:) = zun_e (:,:)
	479	zvb_e (:,:) = zvn_e (:,:)
	480	sshn_e (:,:) = ssha_e(:,:)
	481	zun_e (:,:) = ua_e (:,:)
	482	zvn_e (:,:) = va_e (:,:)
[358]	483	ELSE ! Asselin filtering
[367]	484	zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
	485	zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)
	486	zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)
	487	sshn_e (:,:) = ssha_e(:,:)
	488	zun_e (:,:) = ua_e (:,:)
	489	zvn_e (:,:) = va_e (:,:)
[358]	490	ENDIF
	491
[592]	492	IF( lk_vvl ) THEN ! Variable volume
	493
	494	! Sea surface elevation
	495	! ---------------------
	496	DO jj = 1, jpjm1
	497	DO ji = 1,jpim1
	498
	499	! Sea Surface Height at u-point before
	500	zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	501	& * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &
	502	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
	503
	504	! Sea Surface Height at v-point before
	505	zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	506	& * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &
	507	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
	508
	509	END DO
	510	END DO
	511
	512	! Boundaries conditions
	513	CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. )
	514
	515	! Scale factors at before and after time step
	516	! -------------------------------------------
[661]	517	CALL dom_vvl_sf( zsshun_e, 'U', zfse3un_e ) ; CALL dom_vvl_sf( zsshvn_e, 'V', zfse3vn_e )
[592]	518
	519	! Ocean depth at U- and V-points
	520	hu_e(:,:) = 0.e0
	521	hv_e(:,:) = 0.e0
	522
	523	DO jk = 1, jpk
	524	hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk)
	525	hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk)
	526	END DO
	527
	528	ENDIF ! End variable volume case
	529
[358]	530	! ! ==================== !
	531	END DO ! end loop !
	532	! ! ==================== !
	533
	534
	535	! Time average of after barotropic variables
[788]	536	zcoef = 1.e0 / ( ibaro + 1 )
[358]	537	zssha_b(:,:) = zcoef * zssha_b(:,:)
[367]	538	zua_b (:,:) = zcoef * zua_b (:,:)
	539	zva_b (:,:) = zcoef * zva_b (:,:)
	540	#if defined key_obc
	541	IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
	542	IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)
	543	IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)
	544	IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)
	545	#endif
[358]	546
	547
	548	! ---------------------------------------------------------------------------
	549	! Phase 3 : Update sea surface height from time averaged barotropic variables
	550	! ---------------------------------------------------------------------------
	551
	552
	553	! Horizontal divergence of time averaged barotropic transports
	554	!-------------------------------------------------------------
[592]	555	IF( .NOT. lk_vvl ) THEN
	556	DO jj = 2, jpjm1
	557	DO ji = fs_2, fs_jpim1 ! vector opt.
	558	zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) &
	559	& +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) &
	560	& / ( e1t(ji,jj) * e2t(ji,jj) )
	561	END DO
[358]	562	END DO
[592]	563	ENDIF
[358]	564
[592]	565	#if defined key_obc && ! defined key_vvl
[358]	566	! open boundaries (div must be zero behind the open boundary)
	567	! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
[367]	568	IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east
	569	IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west
[358]	570	IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north
[367]	571	IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south
[358]	572	#endif
	573
	574	! sea surface height
	575	!-------------------
[592]	576	IF( lk_vvl ) THEN
	577	sshbb(:,:) = sshb(:,:)
	578	sshb (:,:) = sshn(:,:)
	579	sshn (:,:) = ssha(:,:)
	580	ELSE
	581	sshb (:,:) = sshn(:,:)
	582	sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)
	583	ENDIF
[358]	584
	585	! ... Boundary conditions on sshn
[367]	586	IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
[358]	587
	588
	589	! -----------------------------------------------------------------------------
	590	! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)
	591	! -----------------------------------------------------------------------------
	592
	593	! Swap on time averaged barotropic variables
	594	! ------------------------------------------
	595	sshb_b(:,:) = sshn_b (:,:)
[592]	596	IF ( kt == nit000 ) zssha_b(:,:) = sshn(:,:)
[358]	597	sshn_b(:,:) = zssha_b(:,:)
	598	un_b (:,:) = zua_b (:,:)
	599	vn_b (:,:) = zva_b (:,:)
	600
	601	! add time averaged barotropic coriolis and surface pressure gradient
	602	! terms to the general momentum trend
	603	! --------------------------------------------------------------------
	604	DO jk=1,jpkm1
	605	ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b
	606	va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b
	607	END DO
	608
[508]	609	! write filtered free surface arrays in restart file
	610	! --------------------------------------------------
	611	IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' )
	612
	613	! print sum trends (used for debugging)
	614	IF( ln_ctl ) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask )
	615	!
	616	END SUBROUTINE dyn_spg_ts
	617
	618
	619	SUBROUTINE ts_rst( kt, cdrw )
	620	!!---------------------------------------------------------------------
	621	!! * ROUTINE ts_rst *
	622	!!
	623	!! ** Purpose : Read or write time-splitting arrays in restart file
	624	!!----------------------------------------------------------------------
	625	INTEGER , INTENT(in) :: kt ! ocean time-step
	626	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	627	!
	628	INTEGER :: ji, jk ! dummy loop indices
	629	!!----------------------------------------------------------------------
	630	!
	631	IF( TRIM(cdrw) == 'READ' ) THEN
[746]	632	IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
[683]	633	CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) )
	634	CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) )
[508]	635	IF( neuler == 0 ) sshb(:,:) = sshn(:,:)
	636	ELSE
[558]	637	IF( nn_rstssh == 1 ) THEN
	638	sshb(:,:) = 0.e0
	639	sshn(:,:) = 0.e0
	640	ENDIF
[508]	641	ENDIF
[746]	642	IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
[683]	643	CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) ) ! free surface issued
	644	CALL iom_get( numror, jpdom_autoglo, 'sshn_b', sshn_b(:,:) ) ! from time-splitting loop
	645	CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! horizontal transports issued
	646	CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop
[508]	647	IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:)
	648	ELSE
	649	sshb_b(:,:) = sshb(:,:)
	650	sshn_b(:,:) = sshn(:,:)
	651	un_b (:,:) = 0.e0
	652	vn_b (:,:) = 0.e0
	653	! vertical sum
	654	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	655	DO jk = 1, jpkm1
	656	DO ji = 1, jpij
	657	un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
	658	vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
	659	END DO
	660	END DO
	661	ELSE ! No vector opt.
	662	DO jk = 1, jpkm1
	663	un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
	664	vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
	665	END DO
	666	ENDIF
	667	ENDIF
	668	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
	669	CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) )
	670	CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) )
	671	CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) ) ! free surface issued
	672	CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) ) ! from barotropic loop
	673	CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! horizontal transports issued
	674	CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop
[358]	675	ENDIF
[508]	676	!
	677	END SUBROUTINE ts_rst
	678
[358]	679	#else
	680	!!----------------------------------------------------------------------
	681	!! Default case : Empty module No standart free surface cst volume
	682	!!----------------------------------------------------------------------
	683	CONTAINS
	684	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
	685	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
	686	END SUBROUTINE dyn_spg_ts
[657]	687	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
	688	INTEGER , INTENT(in) :: kt ! ocean time-step
	689	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	690	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
	691	END SUBROUTINE ts_rst
[358]	692	#endif
	693
	694	!!======================================================================
	695	END MODULE dynspg_ts

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

Download in other formats: