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_r12527_Gurvan_ShallowWater/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynvor.F90 @ 13151

Last change on this file since 13151 was 13151, checked in by gm, 4 years ago
result from merge with qco r12983
Property svn:keywords set to `Id`
File size: 52.5 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
[13151]	17	!! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity
[9019]	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
[13151]	72	INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity
[5836]	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
[13151]	81	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - -
[9528]	82	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2u)/(2*e1e2f) used in F-point metric term calculation
[13151]	83	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1v)/(2*e1e2f) - - - -
	84
[5836]	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
[13151]	88
[3]	89	!! * Substitutions
[12377]	90	# include "do_loop_substitute.h90"
[13151]	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
[13151]	107	!! and planetary vorticity trends) and send them to trd_dyn
[4990]	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	!!
[13151]	197	!! ** Purpose : Compute the now total vorticity trend and add it to
[9528]	198	!! the general trend of the momentum equation.
	199	!!
[13151]	200	!! ** Method : Trend evaluated using now fields (centered in time)
[9528]	201	!! and t-point evaluation of vorticity (planetary and relative).
	202	!! conserves the horizontal kinetic energy.
[13151]	203	!! The general trend of momentum is increased due to the vorticity
[9528]	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
[10425]	219	REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace
	220	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz ! 3D workspace
[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
[12377]	233	DO_2D_10_10
	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
[13151]	237	IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity
[12377]	238	DO_2D_10_10
	239	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	240	END_2D
[9528]	241	ENDIF
[10425]	242	END DO
	243
	244	CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
	245
	246	CASE ( np_CRV ) !* Coriolis + relative vorticity
	247	DO jk = 1, jpkm1 ! Horizontal slab
[12377]	248	DO_2D_10_10
	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
[13151]	252	IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity
[12377]	253	DO_2D_10_10
	254	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	255	END_2D
[10425]	256	ENDIF
	257	END DO
	258
	259	CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
	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
[12377]	271	DO_2D_01_01
	272	zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
[13151]	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
[12377]	277	DO_2D_01_01
[13151]	278	zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
	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
[12377]	283	DO_2D_01_01
[13151]	284	zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) &
	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
[12377]	289	DO_2D_01_01
[13151]	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) &
	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 =!
[12377]	300	DO_2D_00_00
	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) &
[13151]	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) ) )
[12377]	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	!!
[13151]	319	!! ** Purpose : Compute the now total vorticity trend and add it to
[3]	320	!! the general trend of the momentum equation.
	321	!!
[13151]	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.
[13151]	325	!! The general trend of momentum is increased due to the vorticity
[5836]	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)
[13151]	358	zwz(:,:) = ff_f(:,:)
[5836]	359	CASE ( np_RVO ) !* relative vorticity
[12377]	360	DO_2D_10_10
	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
[12377]	365	DO_2D_10_10
	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
[12377]	370	DO_2D_10_10
	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
[12377]	375	DO_2D_10_10
	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 ==!
[12377]	384	DO_2D_10_10
	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 =!
[12377]	398	DO_2D_00_00
	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 )
[13151]	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 )
[12377]	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)
[13151]	454	zwz(:,:) = ff_f(:,:)
[5836]	455	CASE ( np_RVO ) !* relative vorticity
[12377]	456	DO_2D_10_10
	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
[12377]	461	DO_2D_10_10
	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
[12377]	466	DO_2D_10_10
	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
[12377]	471	DO_2D_10_10
	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 ==!
[12377]	480	DO_2D_10_10
	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 =!
[12377]	494	DO_2D_00_00
	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	!!
[13151]	512	!! ** Purpose : Compute the now total vorticity trend and add it to
[108]	513	!! the general trend of the momentum equation.
	514	!!
[13151]	515	!! ** Method : Trend evaluated using now fields (centered in time)
	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
[10425]	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,jpk) :: zwz
[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)
[12377]	552	DO_2D_10_10
[13151]	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)
[12377]	562	DO_2D_10_10
[13151]	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)
[12377]	577	DO_2D_10_10
	578	zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj)
	579	END_2D
[5836]	580	CASE ( np_RVO ) !* relative vorticity
[12377]	581	DO_2D_10_10
	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
[12377]	586	DO_2D_10_10
	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
[12377]	591	DO_2D_10_10
	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
[12377]	597	DO_2D_10_10
	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 ==!
[12377]	606	DO_2D_10_10
	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	!
[10425]	612	CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
	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
[12377]	637	DO_2D_00_00
	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	!!
[13151]	656	!! ** Purpose : Compute the now total vorticity trend and add it to
[9528]	657	!! the general trend of the momentum equation.
	658	!!
[13151]	659	!! ** Method : Trend evaluated using now fields (centered in time)
	660	!! and the Arakawa and Lamb (1980) vector form formulation using
[9528]	661	!! a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
[13151]	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
[13151]	679	REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy
[10425]	680	REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse
	681	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz
[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)
[12377]	697	DO_2D_10_10
	698	zwz(ji,jj,jk) = ff_f(ji,jj)
	699	END_2D
[9528]	700	CASE ( np_RVO ) !* relative vorticity
[12377]	701	DO_2D_10_10
	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
[12377]	707	DO_2D_10_10
	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
[12377]	712	DO_2D_10_10
	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
[12377]	718	DO_2D_10_10
	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 ==!
[12377]	727	DO_2D_10_10
	728	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk)
	729	END_2D
[9528]	730	ENDIF
[10425]	731	END DO
	732	!
	733	CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )
	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
[12377]	760	DO_2D_00_00
	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
[12377]	820	DO_3D_10_10( 1, jpk )
	821	IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	822	& + tmask(ji,jj ,jk) + tmask(ji+1,jj+1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp
	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' )
[13151]	839	!
[5836]	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 )
[13151]	848	IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity'
[5836]	849	nrvm = np_RVO ! relative vorticity
[13151]	850	ntot = np_CRV ! relative + planetary vorticity
[9019]	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) )
[12377]	859	DO_2D_00_00
	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
[10425]	863	CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. ) ! 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) )
[12377]	867	DO_2D_10_10
	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
[10425]	871	CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. ) ! Lateral boundary conditions
[9528]	872	END SELECT
	873	!
[9019]	874	END SELECT
[13151]	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)'
[13151]	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: