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

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

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

Download in other formats: