Context Navigation

dynvor.F90 @ 15795

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DYN/dynvor.F90 @ 15795

Download in other formats: