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.

dynvor.F90 in NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/DYN/dynvor.F90 @ 13553

Last change on this file since 13553 was 13553, checked in by hadcv, 4 years ago
Merge in trunk up to [13550]
Property svn:keywords set to `Id`
File size: 53.1 KB

Rev	Line
[3]	1	MODULE dynvor
	2	!!======================================================================
	3	!! * MODULE dynvor *
	4	!! Ocean dynamics: Update the momentum trend with the relative and
	5	!! planetary vorticity trends
	6	!!======================================================================
[2715]	7	!! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code
[9528]	8	!! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code
[2715]	9	!! 6.0 ! 1996-01 (G. Madec) s-coord, suppress work arrays
	10	!! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module
	11	!! 1.0 ! 2004-02 (G. Madec) vor_een: Original code
	12	!! - ! 2003-08 (G. Madec) add vor_ctl
	13	!! - ! 2005-11 (G. Madec) add dyn_vor (new step architecture)
	14	!! 2.0 ! 2006-11 (G. Madec) flux form advection: add metric term
	15	!! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme
[9019]	16	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
	17	!! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity
	18	!! - ! 2014-06 (G. Madec) suppression of velocity curl from in-core memory
[7646]	19	!! - ! 2016-12 (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T)
[9019]	20	!! 4.0 ! 2017-07 (G. Madec) linear dynamics + trends diag. with Stokes-Coriolis
[9528]	21	!! - ! 2018-03 (G. Madec) add two new schemes (ln_dynvor_enT and ln_dynvor_eet)
	22	!! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation
[503]	23	!!----------------------------------------------------------------------
[3]	24
	25	!!----------------------------------------------------------------------
[9019]	26	!! dyn_vor : Update the momentum trend with the vorticity trend
	27	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
	28	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
	29	!! vor_een : energy and enstrophy conserving (ln_dynvor_een=T)
	30	!! dyn_vor_init : set and control of the different vorticity option
[3]	31	!!----------------------------------------------------------------------
[503]	32	USE oce ! ocean dynamics and tracers
	33	USE dom_oce ! ocean space and time domain
[3294]	34	USE dommsk ! ocean mask
[9019]	35	USE dynadv ! momentum advection
[4990]	36	USE trd_oce ! trends: ocean variables
	37	USE trddyn ! trend manager: dynamics
[7646]	38	USE sbcwave ! Surface Waves (add Stokes-Coriolis force)
	39	USE sbc_oce , ONLY : ln_stcor ! use Stoke-Coriolis force
[5836]	40	!
[503]	41	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	42	USE prtctl ! Print control
	43	USE in_out_manager ! I/O manager
[3294]	44	USE lib_mpp ! MPP library
	45	USE timing ! Timing
[3]	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
[2528]	50	PUBLIC dyn_vor ! routine called by step.F90
[5836]	51	PUBLIC dyn_vor_init ! routine called by nemogcm.F90
[3]	52
[4147]	53	! !!* Namelist namdyn_vor: vorticity term
[9528]	54	LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS)
	55	LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE)
	56	LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT)
	57	LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET)
	58	LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN)
[5836]	59	INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
[9528]	60	LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX)
[5836]	61	LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
[3]	62
[9528]	63	INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme
	64	! ! associated indices:
	65	INTEGER, PUBLIC, PARAMETER :: np_ENS = 0 ! ENS scheme
[5836]	66	INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme
[9528]	67	INTEGER, PUBLIC, PARAMETER :: np_ENT = 2 ! ENT scheme (t-point vorticity)
	68	INTEGER, PUBLIC, PARAMETER :: np_EET = 3 ! EET scheme (EEN using e3t)
[5836]	69	INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme
[9528]	70	INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme
[455]	71
[5836]	72	INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity
	73	! ! associated indices:
[9528]	74	INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary)
	75	INTEGER, PUBLIC, PARAMETER :: np_RVO = 2 ! relative vorticity
	76	INTEGER, PUBLIC, PARAMETER :: np_MET = 3 ! metric term
	77	INTEGER, PUBLIC, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity)
	78	INTEGER, PUBLIC, PARAMETER :: np_CME = 5 ! Coriolis + metric term
	79
	80	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2u_2 ! = di(e2u)/2 used in T-point metric term calculation
	81	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - -
[13286]	82	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2v)/(2*e1e2f) used in F-point metric term calculation
	83	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1u)/(2*e1e2f) - - - -
[5836]	84
	85	REAL(wp) :: r1_4 = 0.250_wp ! =1/4
	86	REAL(wp) :: r1_8 = 0.125_wp ! =1/8
	87	REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12
	88
[3]	89	!! * Substitutions
[12377]	90	# include "do_loop_substitute.h90"
[13237]	91	# include "domzgr_substitute.h90"
	92
[3]	93	!!----------------------------------------------------------------------
[9598]	94	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	95	!! $Id$
[10068]	96	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	97	!!----------------------------------------------------------------------
	98	CONTAINS
	99
[12377]	100	SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
[3]	101	!!----------------------------------------------------------------------
	102	!!
[455]	103	!! ** Purpose : compute the lateral ocean tracer physics.
	104	!!
[12377]	105	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
[503]	106	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[4990]	107	!! and planetary vorticity trends) and send them to trd_dyn
	108	!! for futher diagnostics (l_trddyn=T)
[503]	109	!!----------------------------------------------------------------------
[12377]	110	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	111	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
	112	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation
[2715]	113	!
[9019]	114	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
[455]	115	!!----------------------------------------------------------------------
[2715]	116	!
[9019]	117	IF( ln_timing ) CALL timing_start('dyn_vor')
[3294]	118	!
[9019]	119	IF( l_trddyn ) THEN !== trend diagnostics case : split the added trend in two parts ==!
	120	!
	121	ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
	122	!
[12377]	123	ztrdu(:,:,:) = puu(:,:,:,Krhs) !* planetary vorticity trend (including Stokes-Coriolis force)
	124	ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[9019]	125	SELECT CASE( nvor_scheme )
[12377]	126	CASE( np_ENS ) ; CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! enstrophy conserving scheme
	127	IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	128	CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme
	129	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	130	CASE( np_ENT ) ; CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (T-pts)
	131	IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	132	CASE( np_EET ) ; CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (een with e3t)
	133	IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
	134	CASE( np_EEN ) ; CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy & enstrophy scheme
	135	IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	136	END SELECT
[12377]	137	ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
	138	ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
	139	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt, Kmm )
[9019]	140	!
	141	IF( n_dynadv /= np_LIN_dyn ) THEN !* relative vorticity or metric trend (only in non-linear case)
[12377]	142	ztrdu(:,:,:) = puu(:,:,:,Krhs)
	143	ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[9019]	144	SELECT CASE( nvor_scheme )
[12377]	145	CASE( np_ENT ) ; CALL vor_enT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (T-pts)
	146	CASE( np_EET ) ; CALL vor_eeT( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (een with e3t)
	147	CASE( np_ENE ) ; CALL vor_ene( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme
	148	CASE( np_ENS, np_MIX ) ; CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! enstrophy conserving scheme
	149	CASE( np_EEN ) ; CALL vor_een( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy & enstrophy scheme
[9019]	150	END SELECT
[12377]	151	ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
	152	ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
	153	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt, Kmm )
[9019]	154	ENDIF
	155	!
	156	DEALLOCATE( ztrdu, ztrdv )
	157	!
	158	ELSE !== total vorticity trend added to the general trend ==!
	159	!
	160	SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==!
[9528]	161	CASE( np_ENT ) !* energy conserving scheme (T-pts)
[12377]	162	CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	163	IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9528]	164	CASE( np_EET ) !* energy conserving scheme (een scheme using e3t)
[12377]	165	CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	166	IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	167	CASE( np_ENE ) !* energy conserving scheme
[12377]	168	CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	169	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	170	CASE( np_ENS ) !* enstrophy conserving scheme
[12377]	171	CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	172	IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	173	CASE( np_MIX ) !* mixed ene-ens scheme
[12377]	174	CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! relative vorticity or metric trend (ens)
	175	CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! planetary vorticity trend (ene)
	176	IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	177	CASE( np_EEN ) !* energy and enstrophy conserving scheme
[12377]	178	CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
	179	IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
[9019]	180	END SELECT
[643]	181	!
[9019]	182	ENDIF
[2715]	183	!
[455]	184	! ! print sum trends (used for debugging)
[12377]	185	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, &
	186	& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]	187	!
[9019]	188	IF( ln_timing ) CALL timing_stop('dyn_vor')
[3294]	189	!
[455]	190	END SUBROUTINE dyn_vor
	191
	192
[12377]	193	SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[9528]	194	!!----------------------------------------------------------------------
	195	!! * ROUTINE vor_enT *
	196	!!
	197	!! ** Purpose : Compute the now total vorticity trend and add it to
	198	!! the general trend of the momentum equation.
	199	!!
	200	!! ** Method : Trend evaluated using now fields (centered in time)
	201	!! and t-point evaluation of vorticity (planetary and relative).
	202	!! conserves the horizontal kinetic energy.
	203	!! The general trend of momentum is increased due to the vorticity
	204	!! term which is given by:
	205	!! voru = 1/bu mj[ ( mi(mj(bfrvor))+btf_t)/e3t mj[vn] ]
	206	!! vorv = 1/bv mi[ ( mi(mj(bfrvor))+btf_t)/e3f mj[un] ]
	207	!! where rvor is the relative vorticity at f-point
	208	!!
[12377]	209	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[9528]	210	!!----------------------------------------------------------------------
	211	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	212	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9528]	213	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
	214	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	215	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
	216	!
	217	INTEGER :: ji, jj, jk ! dummy loop indices
	218	REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
[13553]	219	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace
	220	REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
[9528]	221	!!----------------------------------------------------------------------
	222	!
	223	IF( kt == nit000 ) THEN
	224	IF(lwp) WRITE(numout,*)
	225	IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
	226	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	227	ENDIF
	228	!
[10425]	229	!
	230	SELECT CASE( kvor ) !== volume weighted vorticity considered ==!
	231	CASE ( np_RVO ) !* relative vorticity
	232	DO jk = 1, jpkm1 ! Horizontal slab
[13295]	233	DO_2D( 1, 0, 1, 0 )
[12377]	234	zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	235	& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	236	END_2D
[9528]	237	IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity
[13295]	238	DO_2D( 1, 0, 1, 0 )
[12377]	239	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	240	END_2D
[9528]	241	ENDIF
[10425]	242	END DO
	243
[13226]	244	CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
[10425]	245
	246	CASE ( np_CRV ) !* Coriolis + relative vorticity
	247	DO jk = 1, jpkm1 ! Horizontal slab
[13553]	248	DO_2D( 1, 0, 1, 0 ) ! relative vorticity
[12377]	249	zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	250	& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	251	END_2D
[10425]	252	IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity
[13295]	253	DO_2D( 1, 0, 1, 0 )
[12377]	254	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	255	END_2D
[10425]	256	ENDIF
	257	END DO
	258
[13226]	259	CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
[10425]	260
	261	END SELECT
	262
	263	! ! ===============
	264	DO jk = 1, jpkm1 ! Horizontal slab
	265	! ! ===============
	266
	267	SELECT CASE( kvor ) !== volume weighted vorticity considered ==!
	268	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[12377]	269	zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm)
[10425]	270	CASE ( np_RVO ) !* relative vorticity
[13295]	271	DO_2D( 0, 1, 0, 1 )
[12377]	272	zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
[13237]	273	& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) &
	274	& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
[12377]	275	END_2D
[9528]	276	CASE ( np_MET ) !* metric term
[13295]	277	DO_2D( 0, 1, 0, 1 )
[12377]	278	zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
[13237]	279	& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
	280	& * e3t(ji,jj,jk,Kmm)
[12377]	281	END_2D
[9528]	282	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	283	DO_2D( 0, 1, 0, 1 )
[12377]	284	zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
[13237]	285	& + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) &
	286	& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
[12377]	287	END_2D
[9528]	288	CASE ( np_CME ) !* Coriolis + metric
[13295]	289	DO_2D( 0, 1, 0, 1 )
[12377]	290	zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) &
	291	& + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
[13237]	292	& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
	293	& * e3t(ji,jj,jk,Kmm)
[12377]	294	END_2D
[9528]	295	CASE DEFAULT ! error
	296	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
	297	END SELECT
	298	!
	299	! !== compute and add the vorticity term trend =!
[13295]	300	DO_2D( 0, 0, 0, 0 )
[12377]	301	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
	302	& * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) &
	303	& + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) )
	304	!
	305	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
	306	& * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) &
	307	& + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) )
	308	END_2D
[9528]	309	! ! ===============
	310	END DO ! End of slab
	311	! ! ===============
	312	END SUBROUTINE vor_enT
	313
	314
[12377]	315	SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[455]	316	!!----------------------------------------------------------------------
	317	!! * ROUTINE vor_ene *
	318	!!
[3]	319	!! ** Purpose : Compute the now total vorticity trend and add it to
	320	!! the general trend of the momentum equation.
	321	!!
	322	!! ** Method : Trend evaluated using now fields (centered in time)
[5836]	323	!! and the Sadourny (1975) flux form formulation : conserves the
	324	!! horizontal kinetic energy.
	325	!! The general trend of momentum is increased due to the vorticity
	326	!! term which is given by:
[12377]	327	!! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v pvv(:,:,:,Kmm)) ]
	328	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u puu(:,:,:,Kmm)) ]
[5836]	329	!! where rvor is the relative vorticity
[3]	330	!!
[12377]	331	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[3]	332	!!
[503]	333	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	334	!!----------------------------------------------------------------------
[9019]	335	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	336	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9019]	337	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	338	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	339	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[2715]	340	!
[5836]	341	INTEGER :: ji, jj, jk ! dummy loop indices
	342	REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
[9019]	343	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! 2D workspace
[3]	344	!!----------------------------------------------------------------------
[3294]	345	!
[52]	346	IF( kt == nit000 ) THEN
	347	IF(lwp) WRITE(numout,*)
[455]	348	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
	349	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	350	ENDIF
[5836]	351	!
[3]	352	! ! ===============
	353	DO jk = 1, jpkm1 ! Horizontal slab
	354	! ! ===============
[1438]	355	!
[5836]	356	SELECT CASE( kvor ) !== vorticity considered ==!
	357	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[7753]	358	zwz(:,:) = ff_f(:,:)
[5836]	359	CASE ( np_RVO ) !* relative vorticity
[13295]	360	DO_2D( 1, 0, 1, 0 )
[12377]	361	zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	362	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	363	END_2D
[5836]	364	CASE ( np_MET ) !* metric term
[13295]	365	DO_2D( 1, 0, 1, 0 )
[12377]	366	zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	367	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	368	END_2D
[5836]	369	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	370	DO_2D( 1, 0, 1, 0 )
[12377]	371	zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	372	& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	373	END_2D
[5836]	374	CASE ( np_CME ) !* Coriolis + metric
[13295]	375	DO_2D( 1, 0, 1, 0 )
[12377]	376	zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	377	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	378	END_2D
[5836]	379	CASE DEFAULT ! error
	380	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	381	END SELECT
[5836]	382	!
	383	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
[13295]	384	DO_2D( 1, 0, 1, 0 )
[12377]	385	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
	386	END_2D
[5836]	387	ENDIF
[455]	388
	389	IF( ln_sco ) THEN
[12377]	390	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	391	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	392	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[3]	393	ELSE
[12377]	394	zwx(:,:) = e2u(:,:) * pu(:,:,jk)
	395	zwy(:,:) = e1v(:,:) * pv(:,:,jk)
[3]	396	ENDIF
[5836]	397	! !== compute and add the vorticity term trend =!
[13295]	398	DO_2D( 0, 0, 0, 0 )
[12377]	399	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	400	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	401	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	402	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
	403	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	404	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
	405	END_2D
[3]	406	! ! ===============
	407	END DO ! End of slab
	408	! ! ===============
[455]	409	END SUBROUTINE vor_ene
[216]	410
	411
[12377]	412	SUBROUTINE vor_ens( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[3]	413	!!----------------------------------------------------------------------
[455]	414	!! * ROUTINE vor_ens *
[3]	415	!!
	416	!! ** Purpose : Compute the now total vorticity trend and add it to
	417	!! the general trend of the momentum equation.
	418	!!
	419	!! ** Method : Trend evaluated using now fields (centered in time)
	420	!! and the Sadourny (1975) flux FORM formulation : conserves the
	421	!! potential enstrophy of a horizontally non-divergent flow. the
	422	!! trend of the vorticity term is given by:
[12377]	423	!! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v pvv(:,:,:,Kmm)) ]
	424	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u puu(:,:,:,Kmm)) ]
	425	!! Add this trend to the general momentum trend:
	426	!! (u(rhs),v(Krhs)) = (u(rhs),v(Krhs)) + ( voru , vorv )
[3]	427	!!
[12377]	428	!! ** Action : - Update (pu_rhs,pv_rhs)) arrays with the now vorticity term trend
[3]	429	!!
[503]	430	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	431	!!----------------------------------------------------------------------
[9019]	432	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	433	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9019]	434	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	435	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	436	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[2715]	437	!
[5836]	438	INTEGER :: ji, jj, jk ! dummy loop indices
	439	REAL(wp) :: zuav, zvau ! local scalars
[9019]	440	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zww ! 2D workspace
[3]	441	!!----------------------------------------------------------------------
[3294]	442	!
[52]	443	IF( kt == nit000 ) THEN
	444	IF(lwp) WRITE(numout,*)
[455]	445	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
	446	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	447	ENDIF
[3]	448	! ! ===============
	449	DO jk = 1, jpkm1 ! Horizontal slab
	450	! ! ===============
[1438]	451	!
[5836]	452	SELECT CASE( kvor ) !== vorticity considered ==!
	453	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[7646]	454	zwz(:,:) = ff_f(:,:)
[5836]	455	CASE ( np_RVO ) !* relative vorticity
[13295]	456	DO_2D( 1, 0, 1, 0 )
[12377]	457	zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	458	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	459	END_2D
[5836]	460	CASE ( np_MET ) !* metric term
[13295]	461	DO_2D( 1, 0, 1, 0 )
[12377]	462	zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	463	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	464	END_2D
[5836]	465	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	466	DO_2D( 1, 0, 1, 0 )
[12377]	467	zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	468	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	469	END_2D
[5836]	470	CASE ( np_CME ) !* Coriolis + metric
[13295]	471	DO_2D( 1, 0, 1, 0 )
[12377]	472	zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	473	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	474	END_2D
[5836]	475	CASE DEFAULT ! error
	476	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	477	END SELECT
[1438]	478	!
[5836]	479	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
[13295]	480	DO_2D( 1, 0, 1, 0 )
[12377]	481	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
	482	END_2D
[5836]	483	ENDIF
	484	!
	485	IF( ln_sco ) THEN !== horizontal fluxes ==!
[12377]	486	zwz(:,:) = zwz(:,:) / e3f(:,:,jk)
	487	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	488	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[3]	489	ELSE
[12377]	490	zwx(:,:) = e2u(:,:) * pu(:,:,jk)
	491	zwy(:,:) = e1v(:,:) * pv(:,:,jk)
[3]	492	ENDIF
[5836]	493	! !== compute and add the vorticity term trend =!
[13295]	494	DO_2D( 0, 0, 0, 0 )
[12377]	495	zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
	496	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
	497	zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
	498	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
	499	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	500	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
	501	END_2D
[3]	502	! ! ===============
	503	END DO ! End of slab
	504	! ! ===============
[455]	505	END SUBROUTINE vor_ens
[216]	506
	507
[12377]	508	SUBROUTINE vor_een( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[108]	509	!!----------------------------------------------------------------------
[455]	510	!! * ROUTINE vor_een *
[108]	511	!!
	512	!! ** Purpose : Compute the now total vorticity trend and add it to
	513	!! the general trend of the momentum equation.
	514	!!
	515	!! ** Method : Trend evaluated using now fields (centered in time)
[1438]	516	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]	517	!! both the horizontal kinetic energy and the potential enstrophy
[1438]	518	!! when horizontal divergence is zero (see the NEMO documentation)
[12377]	519	!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[108]	520	!!
[12377]	521	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[108]	522	!!
[503]	523	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	524	!!----------------------------------------------------------------------
[9019]	525	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	526	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9019]	527	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	528	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	529	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[5836]	530	!
	531	INTEGER :: ji, jj, jk ! dummy loop indices
	532	INTEGER :: ierr ! local integer
	533	REAL(wp) :: zua, zva ! local scalars
[9528]	534	REAL(wp) :: zmsk, ze3f ! local scalars
[13553]	535	REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f
	536	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse
	537	REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
[108]	538	!!----------------------------------------------------------------------
[3294]	539	!
[108]	540	IF( kt == nit000 ) THEN
	541	IF(lwp) WRITE(numout,*)
[455]	542	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	543	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[1438]	544	ENDIF
[5836]	545	!
	546	! ! ===============
	547	DO jk = 1, jpkm1 ! Horizontal slab
	548	! ! ===============
	549	!
	550	SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point
	551	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
[13295]	552	DO_2D( 1, 0, 1, 0 )
[13237]	553	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	554	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	555	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	556	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
[12377]	557	IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f
	558	ELSE ; z1_e3f(ji,jj) = 0._wp
	559	ENDIF
	560	END_2D
[5836]	561	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
[13295]	562	DO_2D( 1, 0, 1, 0 )
[13237]	563	ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
	564	& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) &
	565	& + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
	566	& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) )
[12377]	567	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	568	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
	569	IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f
	570	ELSE ; z1_e3f(ji,jj) = 0._wp
	571	ENDIF
	572	END_2D
[5836]	573	END SELECT
	574	!
	575	SELECT CASE( kvor ) !== vorticity considered ==!
	576	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[13295]	577	DO_2D( 1, 0, 1, 0 )
[12377]	578	zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
	579	END_2D
[5836]	580	CASE ( np_RVO ) !* relative vorticity
[13295]	581	DO_2D( 1, 0, 1, 0 )
[12377]	582	zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	583	& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
	584	END_2D
[5836]	585	CASE ( np_MET ) !* metric term
[13295]	586	DO_2D( 1, 0, 1, 0 )
[12377]	587	zwz(ji,jj,jk) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	588	& - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	589	END_2D
[5836]	590	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	591	DO_2D( 1, 0, 1, 0 )
[12377]	592	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	593	& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	594	& * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	595	END_2D
[5836]	596	CASE ( np_CME ) !* Coriolis + metric
[13295]	597	DO_2D( 1, 0, 1, 0 )
[12377]	598	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	599	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	600	END_2D
[5836]	601	CASE DEFAULT ! error
	602	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	603	END SELECT
[5836]	604	!
	605	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
[13295]	606	DO_2D( 1, 0, 1, 0 )
[12377]	607	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	608	END_2D
[5836]	609	ENDIF
[10425]	610	END DO ! End of slab
[5836]	611	!
[13226]	612	CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
[10425]	613
	614	DO jk = 1, jpkm1 ! Horizontal slab
[5907]	615	!
[5836]	616	! !== horizontal fluxes ==!
[12377]	617	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	618	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[108]	619
[5836]	620	! !== compute and add the vorticity term trend =!
[1438]	621	jj = 2
	622	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
[5836]	623	DO ji = 2, jpi ! split in 2 parts due to vector opt.
[10425]	624	ztne(ji,jj) = zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk)
	625	ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk)
	626	ztse(ji,jj) = zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk)
	627	ztsw(ji,jj) = zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk)
[108]	628	END DO
	629	DO jj = 3, jpj
[12377]	630	DO ji = 2, jpi ! vector opt. ok because we start at jj = 3
[10425]	631	ztne(ji,jj) = zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk)
	632	ztnw(ji,jj) = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk)
	633	ztse(ji,jj) = zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk)
	634	ztsw(ji,jj) = zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk)
[108]	635	END DO
	636	END DO
[13295]	637	DO_2D( 0, 0, 0, 0 )
[12377]	638	zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	639	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	640	zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	641	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
	642	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
	643	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
	644	END_2D
[108]	645	! ! ===============
	646	END DO ! End of slab
	647	! ! ===============
[455]	648	END SUBROUTINE vor_een
[216]	649
	650
[9528]	651
[12377]	652	SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
[9528]	653	!!----------------------------------------------------------------------
	654	!! * ROUTINE vor_eeT *
	655	!!
	656	!! ** Purpose : Compute the now total vorticity trend and add it to
	657	!! the general trend of the momentum equation.
	658	!!
	659	!! ** Method : Trend evaluated using now fields (centered in time)
	660	!! and the Arakawa and Lamb (1980) vector form formulation using
	661	!! a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
	662	!! The change consists in
[12377]	663	!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
[9528]	664	!!
[12377]	665	!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
[9528]	666	!!
	667	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	668	!!----------------------------------------------------------------------
	669	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[12377]	670	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[9528]	671	INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
[12377]	672	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
	673	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
[9528]	674	!
	675	INTEGER :: ji, jj, jk ! dummy loop indices
	676	INTEGER :: ierr ! local integer
	677	REAL(wp) :: zua, zva ! local scalars
	678	REAL(wp) :: zmsk, z1_e3t ! local scalars
[13553]	679	REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy
	680	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse
	681	REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
[9528]	682	!!----------------------------------------------------------------------
	683	!
	684	IF( kt == nit000 ) THEN
	685	IF(lwp) WRITE(numout,*)
	686	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	687	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	688	ENDIF
	689	!
	690	! ! ===============
	691	DO jk = 1, jpkm1 ! Horizontal slab
	692	! ! ===============
	693	!
	694	!
	695	SELECT CASE( kvor ) !== vorticity considered ==!
	696	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[13295]	697	DO_2D( 1, 0, 1, 0 )
[12377]	698	zwz(ji,jj,jk) = ff_f(ji,jj)
	699	END_2D
[9528]	700	CASE ( np_RVO ) !* relative vorticity
[13295]	701	DO_2D( 1, 0, 1, 0 )
[12377]	702	zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	703	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	704	& * r1_e1e2f(ji,jj)
	705	END_2D
[9528]	706	CASE ( np_MET ) !* metric term
[13295]	707	DO_2D( 1, 0, 1, 0 )
[12377]	708	zwz(ji,jj,jk) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	709	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	710	END_2D
[9528]	711	CASE ( np_CRV ) !* Coriolis + relative vorticity
[13295]	712	DO_2D( 1, 0, 1, 0 )
[12377]	713	zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
	714	& - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &
	715	& * r1_e1e2f(ji,jj) )
	716	END_2D
[9528]	717	CASE ( np_CME ) !* Coriolis + metric
[13295]	718	DO_2D( 1, 0, 1, 0 )
[12377]	719	zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
	720	& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
	721	END_2D
[9528]	722	CASE DEFAULT ! error
	723	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
	724	END SELECT
	725	!
	726	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
[13295]	727	DO_2D( 1, 0, 1, 0 )
[12377]	728	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	729	END_2D
[9528]	730	ENDIF
[10425]	731	END DO
	732	!
[13226]	733	CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )
[10425]	734	!
	735	DO jk = 1, jpkm1 ! Horizontal slab
	736
	737	! !== horizontal fluxes ==!
[12377]	738	zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk)
	739	zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk)
[9528]	740
	741	! !== compute and add the vorticity term trend =!
	742	jj = 2
	743	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
	744	DO ji = 2, jpi ! split in 2 parts due to vector opt.
[12377]	745	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[10425]	746	ztne(ji,jj) = ( zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) ) * z1_e3t
	747	ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) ) * z1_e3t
	748	ztse(ji,jj) = ( zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
	749	ztsw(ji,jj) = ( zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) ) * z1_e3t
[9528]	750	END DO
	751	DO jj = 3, jpj
[12377]	752	DO ji = 2, jpi ! vector opt. ok because we start at jj = 3
	753	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
[10425]	754	ztne(ji,jj) = ( zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) ) * z1_e3t
	755	ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) ) * z1_e3t
	756	ztse(ji,jj) = ( zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) ) * z1_e3t
	757	ztsw(ji,jj) = ( zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) ) * z1_e3t
[9528]	758	END DO
	759	END DO
[13295]	760	DO_2D( 0, 0, 0, 0 )
[12377]	761	zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	762	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	763	zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	764	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
	765	pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
	766	pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
	767	END_2D
[9528]	768	! ! ===============
	769	END DO ! End of slab
	770	! ! ===============
	771	END SUBROUTINE vor_eeT
	772
	773
[2528]	774	SUBROUTINE dyn_vor_init
[3]	775	!!---------------------------------------------------------------------
[2528]	776	!! * ROUTINE dyn_vor_init *
[3]	777	!!
	778	!! ** Purpose : Control the consistency between cpp options for
[1438]	779	!! tracer advection schemes
[3]	780	!!----------------------------------------------------------------------
[9528]	781	INTEGER :: ji, jj, jk ! dummy loop indices
	782	INTEGER :: ioptio, ios ! local integer
[2715]	783	!!
[9528]	784	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, &
	785	& ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk
[3]	786	!!----------------------------------------------------------------------
[9528]	787	!
	788	IF(lwp) THEN
	789	WRITE(numout,*)
	790	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
	791	WRITE(numout,*) '~~~~~~~~~~~~'
	792	ENDIF
	793	!
[4147]	794	READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
[11536]	795	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist' )
[4147]	796	READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
[11536]	797	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist' )
[4624]	798	IF(lwm) WRITE ( numond, namdyn_vor )
[9528]	799	!
[503]	800	IF(lwp) THEN ! Namelist print
[7646]	801	WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme'
	802	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
[9528]	803	WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
	804	WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT
	805	WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT
[7646]	806	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
	807	WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f
[9528]	808	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
[7646]	809	WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk
[52]	810	ENDIF
	811
[9528]	812	IF( ln_dynvor_msk ) CALL ctl_stop( 'dyn_vor_init: masked vorticity is not currently not available')
	813
[5836]	814	!!gm this should be removed when choosing a unique strategy for fmask at the coast
[3294]	815	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
	816	! at angles with three ocean points and one land point
[5836]	817	IF(lwp) WRITE(numout,*)
[7646]	818	IF(lwp) WRITE(numout,*) ' change fmask value in the angles (T) ln_vorlat = ', ln_vorlat
[3294]	819	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
[13295]	820	DO_3D( 1, 0, 1, 0, 1, jpk )
[12377]	821	IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
[12793]	822	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp
[12377]	823	END_3D
[9528]	824	!
[10425]	825	CALL lbc_lnk( 'dynvor', fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
[9528]	826	!
[3294]	827	ENDIF
[5836]	828	!!gm end
[3294]	829
[5836]	830	ioptio = 0 ! type of scheme for vorticity (set nvor_scheme)
[9528]	831	IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF
	832	IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF
	833	IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF
	834	IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF
	835	IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF
	836	IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF
[5836]	837	!
[6140]	838	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
[5836]	839	!
	840	IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot)
[9019]	841	ncor = np_COR ! planetary vorticity
	842	SELECT CASE( n_dynadv )
	843	CASE( np_LIN_dyn )
[9190]	844	IF(lwp) WRITE(numout,*) ' ==>>> linear dynamics : total vorticity = Coriolis'
[9019]	845	nrvm = np_COR ! planetary vorticity
	846	ntot = np_COR ! - -
	847	CASE( np_VEC_c2 )
[9190]	848	IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity'
[5836]	849	nrvm = np_RVO ! relative vorticity
[9019]	850	ntot = np_CRV ! relative + planetary vorticity
	851	CASE( np_FLX_c2 , np_FLX_ubs )
[9190]	852	IF(lwp) WRITE(numout,*) ' ==>>> flux form dynamics : total vorticity = Coriolis + metric term'
[5836]	853	nrvm = np_MET ! metric term
	854	ntot = np_CME ! Coriolis + metric term
[9528]	855	!
	856	SELECT CASE( nvor_scheme ) ! pre-computed gradients for the metric term:
	857	CASE( np_ENT ) !* T-point metric term : pre-compute di(e2u)/2 and dj(e1v)/2
	858	ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) )
[13295]	859	DO_2D( 0, 0, 0, 0 )
[12377]	860	di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj ) ) * 0.5_wp
	861	dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp
	862	END_2D
[13226]	863	CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp ) ! Lateral boundary conditions
[9528]	864	!
	865	CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2e1e2f) and dj(e1v)/(2e1e2f)
	866	ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) )
[13295]	867	DO_2D( 1, 0, 1, 0 )
[12377]	868	di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj ) - e2v(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj)
	869	dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj)
	870	END_2D
[13226]	871	CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp ) ! Lateral boundary conditions
[9528]	872	END SELECT
	873	!
[9019]	874	END SELECT
[643]	875
[503]	876	IF(lwp) THEN ! Print the choice
	877	WRITE(numout,*)
[9019]	878	SELECT CASE( nvor_scheme )
[9528]	879	CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)'
	880	CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)'
	881	CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)'
	882	CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)'
	883	CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)'
	884	CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)'
[9019]	885	END SELECT
[3]	886	ENDIF
[503]	887	!
[2528]	888	END SUBROUTINE dyn_vor_init
[3]	889
[503]	890	!!==============================================================================
[3]	891	END MODULE dynvor

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

Download in other formats: