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

Last change on this file since 1502 was 1502, checked in by rblod, 15 years ago
Update dynnxt and dynspg_ts for variable volume, see ticket #474
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 30.3 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
[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
[911]	29	USE bdy_oce ! unstructured open boundaries
	30	USE bdy_par ! unstructured open boundaries
	31	USE bdydta ! unstructured open boundaries
	32	USE bdydyn ! unstructured open boundaries
	33	USE bdytides ! tidal forcing at unstructured open boundaries.
[358]	34	USE lib_mpp ! distributed memory computing library
	35	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	36	USE prtctl ! Print control
	37	USE in_out_manager ! I/O manager
[508]	38	USE iom
	39	USE restart ! only for lrst_oce
[358]	40
	41	IMPLICIT NONE
	42	PRIVATE
	43
	44	PUBLIC dyn_spg_ts ! routine called by step.F90
[800]	45	PUBLIC ts_rst ! routine called by istate.F90
[358]	46
[1502]	47
[1438]	48	REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne ! triad of coriolis parameter
	49	REAL(wp), DIMENSION(jpi,jpj) :: ftsw, ftse ! (only used with een vorticity scheme)
[508]	50
[1502]	51	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: un_b, vn_b ! averaged velocity
	52
[358]	53	!! * Substitutions
	54	# include "domzgr_substitute.h90"
	55	# include "vectopt_loop_substitute.h90"
[1438]	56	!!-------------------------------------------------------------------------
	57	!! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
[888]	58	!! $Id$
[1438]	59	!! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	60	!!-------------------------------------------------------------------------
[358]	61
	62	CONTAINS
	63
	64	SUBROUTINE dyn_spg_ts( kt )
	65	!!----------------------------------------------------------------------
	66	!! * routine dyn_spg_ts *
	67	!!
	68	!! ** Purpose : Compute the now trend due to the surface pressure
	69	!! gradient in case of free surface formulation with time-splitting.
	70	!! Add it to the general trend of momentum equation.
	71	!!
	72	!! ** Method : Free surface formulation with time-splitting
	73	!! -1- Save the vertically integrated trend. This general trend is
	74	!! held constant over the barotropic integration.
	75	!! The Coriolis force is removed from the general trend as the
	76	!! surface gradient and the Coriolis force are updated within
	77	!! the barotropic integration.
[367]	78	!! -2- Barotropic loop : updates of sea surface height (ssha_e) and
[1502]	79	!! barotropic velocity (ua_e and va_e) through barotropic
[358]	80	!! momentum and continuity integration. Barotropic former
	81	!! variables are time averaging over the full barotropic cycle
[1502]	82	!! (= 2 * baroclinic time step) and saved in zuX_b
[358]	83	!! and zvX_b (X specifying after, now or before).
[1438]	84	!! -3- The new general trend becomes :
[1502]	85	!! ua = ua - sum_k(ua)/H + ( ua_e - sum_k(ub) )
[358]	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	!!---------------------------------------------------------------------
[1502]	91	INTEGER, INTENT(in) :: kt ! ocean time-step index
[1438]	92	!!
[1502]	93	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	94	INTEGER :: icycle ! temporary scalar
	95	REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b, & ! temporary scalars
	96	z1_8, zspgu, zcubt, zx1, zy1, & ! " "
	97	z1_4, zspgv, zcvbt, zx2, zy2 ! " "
	98	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv, zsshb_e
	99	!!
	100	REAL(wp), DIMENSION(jpi,jpj) :: zsshun_e, zsshvn_e ! 2D workspace
	101	!!
	102	REAL(wp), DIMENSION(jpi,jpj) :: zcu, zwx, zua, zun, zub ! 2D workspace
	103	REAL(wp), DIMENSION(jpi,jpj) :: zcv, zwy, zva, zvn, zvb ! - -
	104	REAL(wp), DIMENSION(jpi,jpj) :: zun_e, zub_e, zu_sum ! 2D workspace
	105	REAL(wp), DIMENSION(jpi,jpj) :: zvn_e, zvb_e, zv_sum ! - -
	106	REAL(wp), DIMENSION(jpi,jpj) :: zssh_sum ! - -
[358]	107	!!----------------------------------------------------------------------
	108
[1502]	109	IF( kt == nit000 ) THEN !* initialisation
[508]	110	!
[358]	111	IF(lwp) WRITE(numout,*)
	112	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
	113	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
[1241]	114	IF(lwp) WRITE(numout,) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2nn_baro
[1502]	115	!
[508]	116	CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields:
[1502]	117	! ! un_b, vn_b
	118	!
	119	ua_e (:,:) = un_b (:,:)
	120	va_e (:,:) = vn_b (:,:)
	121	hu_e (:,:) = hu (:,:)
	122	hv_e (:,:) = hv (:,:)
	123	hur_e (:,:) = hur (:,:)
	124	hvr_e (:,:) = hvr (:,:)
[358]	125	IF( ln_dynvor_een ) THEN
[508]	126	ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0
[358]	127	DO jj = 2, jpj
	128	DO ji = fs_2, jpi ! vector opt.
[508]	129	ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3.
	130	ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3.
	131	ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3.
	132	ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3.
[358]	133	END DO
	134	END DO
	135	ENDIF
[508]	136	!
	137	ENDIF
[358]	138
[1502]	139	! !* Local constant initialization
[358]	140	z2dt_b = 2.0 * rdt ! baroclinic time step
[1502]	141	z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates
	142	z1_4 = 0.5 * 0.5
[374]	143	zraur = 1. / rauw ! 1 / volumic mass of pure water
[1502]	144	!
	145	zhdiv(:,:) = 0.e0 ! barotropic divergence
[1438]	146
[358]	147	! -----------------------------------------------------------------------------
	148	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
	149	! -----------------------------------------------------------------------------
[1502]	150	!
	151	! !* e3d/dt(Ua), e3Ub, e3*Vn (Vertically integrated)
	152	! ! --------------------------
	153	zua(:,:) = 0.e0 ; zun(:,:) = 0.e0 ; zub(:,:) = 0.e0
	154	zva(:,:) = 0.e0 ; zvn(:,:) = 0.e0 ; zvb(:,:) = 0.e0
	155	!
	156	DO jk = 1, jpkm1
	157	#if defined key_vectopt_loop
	158	DO jj = 1, 1 !Vector opt. => forced unrolling
[358]	159	DO ji = 1, jpij
[1502]	160	#else
	161	DO jj = 1, jpj
	162	DO ji = 1, jpi
	163	#endif
	164	! ! now trend
	165	zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
	166	zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
	167	! ! now velocity
	168	zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk)
	169	zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk)
	170	! ! before velocity
	171	zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)
	172	zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
[358]	173	END DO
	174	END DO
[1502]	175	END DO
	176
	177	! !* baroclinic momentum trend (remove the vertical mean trend)
	178	DO jk = 1, jpkm1 ! --------------------------
	179	DO jj = 2, jpjm1
	180	DO ji = fs_2, fs_jpim1 ! vector opt.
	181	ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
	182	va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
	183	END DO
[358]	184	END DO
[1502]	185	END DO
[358]	186
[1502]	187	! !* barotropic Coriolis trends * H (vorticity scheme dependent)
	188	! ! ---------------------------====
	189	zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport
	190	zwy(:,:) = zvn(:,:) * e1v(:,:)
	191	!
[358]	192	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
	193	DO jj = 2, jpjm1
	194	DO ji = fs_2, fs_jpim1 ! vector opt.
	195	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	196	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	197	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	198	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	199	! energy conserving formulation for planetary vorticity term
[1502]	200	zcu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )
	201	zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )
[358]	202	END DO
	203	END DO
[508]	204	!
[358]	205	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
	206	DO jj = 2, jpjm1
	207	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	208	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	209	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[358]	210	zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) )
	211	zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) )
	212	END DO
	213	END DO
[508]	214	!
[358]	215	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
	216	DO jj = 2, jpjm1
	217	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	218	zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	219	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
	220	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) &
	221	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) )
[358]	222	END DO
	223	END DO
[508]	224	!
[358]	225	ENDIF
	226
[1502]	227	! !* Right-Hand-Side of the barotropic momentum equation
	228	! ! ----------------------------------------------------
	229	IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient
	230	DO jj = 2, jpjm1
[358]	231	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	232	zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) &
	233	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)
	234	zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) &
	235	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)
[358]	236	END DO
	237	END DO
[1502]	238	ENDIF
[358]	239
[1502]	240	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
[358]	241	DO ji = fs_2, fs_jpim1
	242	zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
	243	zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
	244	END DO
	245	END DO
	246
[1502]	247	! !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
	248	! ! --------------------------
	249	zua(:,:) = zua(:,:) * hur(:,:)
	250	zva(:,:) = zva(:,:) * hvr(:,:)
	251	!
	252	IF( lk_vvl ) THEN
	253	zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
	254	zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
	255	ELSE
	256	zub(:,:) = zub(:,:) * hur(:,:)
	257	zvb(:,:) = zvb(:,:) * hvr(:,:)
	258	ENDIF
	259
[358]	260	! -----------------------------------------------------------------------
	261	! Phase 2 : Integration of the barotropic equations with time splitting
	262	! -----------------------------------------------------------------------
[1502]	263	!
	264	! ! ==================== !
	265	! ! Initialisations !
	266	! ! ==================== !
	267	icycle = 2 * nn_baro ! Number of barotropic sub time-step
	268
	269	! ! Start from NOW field
	270	hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points
	271	hv_e (:,:) = hv (:,:)
	272	hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points
	273	hvr_e (:,:) = hvr (:,:)
	274	zsshb_e(:,:) = sshn (:,:) ! sea surface height (before and now)
	275	sshn_e (:,:) = sshn (:,:)
	276
	277	zun_e (:,:) = un_b (:,:) ! barotropic velocity (external)
	278	zvn_e (:,:) = vn_b (:,:)
	279	zub_e (:,:) = un_b (:,:)
	280	zvb_e (:,:) = vn_b (:,:)
[358]	281
[1502]	282	zu_sum (:,:) = un_b (:,:) ! summation
	283	zv_sum (:,:) = vn_b (:,:)
	284	zssh_sum(:,:) = sshn (:,:)
[358]	285
[1502]	286	#if defined key_obc
[367]	287	! set ssh corrections to 0
	288	! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
	289	IF( lp_obc_east ) sshfoe_b(:,:) = 0.e0
	290	IF( lp_obc_west ) sshfow_b(:,:) = 0.e0
	291	IF( lp_obc_south ) sshfos_b(:,:) = 0.e0
	292	IF( lp_obc_north ) sshfon_b(:,:) = 0.e0
	293	#endif
	294
[1502]	295	! ! ==================== !
	296	DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1)
	297	! ! ==================== !
[1241]	298	z2dt_e = 2. * ( rdt / nn_baro )
[1502]	299	IF( jn == 1 ) z2dt_e = rdt / nn_baro
[358]	300
[1502]	301	! !* Update the forcing (OBC, BDY and tides)
	302	! ! ------------------
	303	IF( lk_obc ) CALL obc_dta_bt( kt, jn )
	304	IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dta_bt( kt, jn+1 )
[367]	305
[1502]	306	! !* after ssh_e
	307	! ! -----------
	308	DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports
[358]	309	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	310	zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) &
	311	& - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) &
	312	& + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) &
	313	& - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
[358]	314	END DO
	315	END DO
[1502]	316	!
[358]	317	#if defined key_obc
[1502]	318	! ! OBC : zhdiv must be zero behind the open boundary
	319	!! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
	320	IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east
	321	IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west
[367]	322	IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north
[1502]	323	IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south
[358]	324	#endif
[1170]	325	#if defined key_bdy
[1502]	326	zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask
[1170]	327	#endif
[1502]	328	!
	329	DO jj = 2, jpjm1 ! leap-frog on ssh_e
[358]	330	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	331	ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
[358]	332	END DO
	333	END DO
	334
[1502]	335	! !* after barotropic velocities (vorticity scheme dependent)
	336	! ! ---------------------------
	337	zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport
	338	zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
	339	!
	340	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
[358]	341	DO jj = 2, jpjm1
	342	DO ji = fs_2, fs_jpim1 ! vector opt.
	343	! surface pressure gradient
[592]	344	IF( lk_vvl) THEN
[1502]	345	zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	346	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	347	zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	348	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	349	ELSE
[1502]	350	zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	351	zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	352	ENDIF
[358]	353	! energy conserving formulation for planetary vorticity term
	354	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	355	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	356	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	357	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
[1502]	358	zcubt = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
	359	zcvbt =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
	360	! after barotropic velocity
	361	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
	362	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
[358]	363	END DO
	364	END DO
[508]	365	!
[1502]	366	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
[358]	367	DO jj = 2, jpjm1
	368	DO ji = fs_2, fs_jpim1 ! vector opt.
[1502]	369	! surface pressure gradient
[592]	370	IF( lk_vvl) THEN
[1502]	371	zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
	372	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
	373	zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
	374	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
[592]	375	ELSE
[1502]	376	zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	377	zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	378	ENDIF
[358]	379	! enstrophy conserving formulation for planetary vorticity term
[1502]	380	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	381	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	382	zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
	383	zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj)
	384	! after barotropic velocity
	385	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
	386	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
[358]	387	END DO
	388	END DO
[508]	389	!
[1502]	390	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving 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
[1502]	395	zspgu = -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	zspgv = -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
[1502]	400	zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
	401	zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
[592]	402	ENDIF
[358]	403	! energy/enstrophy conserving formulation for planetary vorticity term
[1502]	404	zcubt = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
	405	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
	406	zcvbt = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &
	407	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj)
	408	! after barotropic velocity
	409	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
	410	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
[358]	411	END DO
	412	END DO
[508]	413	!
[358]	414	ENDIF
[1502]	415	! !* domain lateral boundary
	416	! ! -----------------------
	417	! ! Flather's boundary condition for the barotropic loop :
	418	! ! - Update sea surface height on each open boundary
	419	! ! - Correct the velocity
[358]	420
[1502]	421	IF( lk_obc ) CALL obc_fla_ts
	422	IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dyn_fla( sshn_e )
	423	!
	424	CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries
	425	CALL lbc_lnk( va_e , 'V', -1. )
	426	CALL lbc_lnk( ssha_e, 'T', 1. )
[358]	427
[1502]	428	zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps
	429	zv_sum (:,:) = zv_sum (:,:) + va_e (:,:)
	430	zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)
[367]	431
[1502]	432	! !* Time filter and swap
	433	! ! --------------------
	434	IF( jn == 1 ) THEN ! Swap only (1st Euler time step)
	435	zsshb_e(:,:) = sshn_e(:,:)
	436	zub_e (:,:) = zun_e (:,:)
	437	zvb_e (:,:) = zvn_e (:,:)
	438	sshn_e (:,:) = ssha_e(:,:)
	439	zun_e (:,:) = ua_e (:,:)
	440	zvn_e (:,:) = va_e (:,:)
	441	ELSE ! Swap + Filter
	442	zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
	443	zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)
	444	zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)
	445	sshn_e (:,:) = ssha_e(:,:)
	446	zun_e (:,:) = ua_e (:,:)
	447	zvn_e (:,:) = va_e (:,:)
[358]	448	ENDIF
	449
[1502]	450	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
	451	! ! ------------------
	452	DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points
	453	DO ji = 1, fs_jpim1 ! Vector opt.
	454	zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	455	& * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &
	456	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
	457	zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	458	& * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &
	459	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
[592]	460	END DO
	461	END DO
[1502]	462	CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions
	463	CALL lbc_lnk( zsshvn_e, 'V', 1. )
[1438]	464	!
[1502]	465	hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points
	466	hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
	467	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
	468	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
	469	!
[1438]	470	ENDIF
[358]	471	! ! ==================== !
	472	END DO ! end loop !
	473	! ! ==================== !
	474
[367]	475	#if defined key_obc
[1502]	476	IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ?????
[1241]	477	IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)
	478	IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)
	479	IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)
[367]	480	#endif
[358]	481
[1438]	482	! -----------------------------------------------------------------------------
[1502]	483	! Phase 3. update the general trend with the barotropic trend
[1438]	484	! -----------------------------------------------------------------------------
[1502]	485	!
	486	! !* Time average ==> after barotropic u, v, ssh
	487	zcoef = 1.e0 / ( 2 * nn_baro + 1 )
	488	un_b (:,:) = zcoef * zu_sum (:,:)
	489	vn_b (:,:) = zcoef * zv_sum (:,:)
	490	sshn_b(:,:) = zcoef * zssh_sum(:,:)
	491	!
	492	! !* update the general momentum trend
[358]	493	DO jk=1,jpkm1
[1502]	494	ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b
	495	va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b
[358]	496	END DO
[1502]	497	!
	498	! !* write time-spliting arrays in the restart
[508]	499	IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' )
	500	!
	501	END SUBROUTINE dyn_spg_ts
	502
	503
	504	SUBROUTINE ts_rst( kt, cdrw )
	505	!!---------------------------------------------------------------------
	506	!! * ROUTINE ts_rst *
	507	!!
	508	!! ** Purpose : Read or write time-splitting arrays in restart file
	509	!!----------------------------------------------------------------------
	510	INTEGER , INTENT(in) :: kt ! ocean time-step
	511	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	512	!
	513	INTEGER :: ji, jk ! dummy loop indices
	514	!!----------------------------------------------------------------------
	515	!
	516	IF( TRIM(cdrw) == 'READ' ) THEN
[1502]	517	IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
	518	CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued
	519	CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop
[508]	520	ELSE
[1502]	521	un_b (:,:) = 0.e0
	522	vn_b (:,:) = 0.e0
[508]	523	! vertical sum
	524	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
	525	DO jk = 1, jpkm1
	526	DO ji = 1, jpij
[1502]	527	un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
	528	vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
[508]	529	END DO
	530	END DO
	531	ELSE ! No vector opt.
	532	DO jk = 1, jpkm1
[1502]	533	un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
	534	vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
[508]	535	END DO
	536	ENDIF
[1502]	537	un_b (:,:) = un_b(:,:) * hur(:,:)
	538	vn_b (:,:) = vn_b(:,:) * hvr(:,:)
[508]	539	ENDIF
	540	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
[1502]	541	CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity issued
	542	CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop
[358]	543	ENDIF
[508]	544	!
	545	END SUBROUTINE ts_rst
	546
[358]	547	#else
	548	!!----------------------------------------------------------------------
	549	!! Default case : Empty module No standart free surface cst volume
	550	!!----------------------------------------------------------------------
	551	CONTAINS
	552	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
	553	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
	554	END SUBROUTINE dyn_spg_ts
[657]	555	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
	556	INTEGER , INTENT(in) :: kt ! ocean time-step
	557	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	558	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
	559	END SUBROUTINE ts_rst
[358]	560	#endif
	561
	562	!!======================================================================
	563	END MODULE dynspg_ts

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

Download in other formats: