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_r12563_ASINTER-06_ABL_improvement/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/DYN/dynvor.F90 @ 13159

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

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

Download in other formats: