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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 5861

Last change on this file since 5861 was 5836, checked in by cetlod, 9 years ago
merge the simplification branch onto the trunk, see ticket #1612
Property svn:keywords set to `Id`
File size: 36.6 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
	8	!! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code
	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
	16	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
[4990]	17	!! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity
[5836]	18	!! - ! 2014-06 (G. Madec) suppression of velocity curl from in-core memory
[503]	19	!!----------------------------------------------------------------------
[3]	20
	21	!!----------------------------------------------------------------------
[2528]	22	!! dyn_vor : Update the momentum trend with the vorticity trend
	23	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
	24	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
	25	!! vor_een : energy and enstrophy conserving (ln_dynvor_een=T)
	26	!! dyn_vor_init : set and control of the different vorticity option
[3]	27	!!----------------------------------------------------------------------
[503]	28	USE oce ! ocean dynamics and tracers
	29	USE dom_oce ! ocean space and time domain
[3294]	30	USE dommsk ! ocean mask
[643]	31	USE dynadv ! momentum advection (use ln_dynadv_vec value)
[4990]	32	USE trd_oce ! trends: ocean variables
	33	USE trddyn ! trend manager: dynamics
[5836]	34	!
[503]	35	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	36	USE prtctl ! Print control
	37	USE in_out_manager ! I/O manager
[3294]	38	USE lib_mpp ! MPP library
	39	USE wrk_nemo ! Memory Allocation
	40	USE timing ! Timing
[3]	41
[3294]	42
[3]	43	IMPLICIT NONE
	44	PRIVATE
	45
[2528]	46	PUBLIC dyn_vor ! routine called by step.F90
[5836]	47	PUBLIC dyn_vor_init ! routine called by nemogcm.F90
[3]	48
[4147]	49	! !!* Namelist namdyn_vor: vorticity term
[5836]	50	LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme (ENE)
	51	LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS)
	52	LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX)
	53	LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme (EEN)
	54	INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1)
	55	LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes)
[3]	56
[5836]	57	INTEGER :: nvor_scheme ! choice of the type of advection scheme
	58	! ! associated indices:
	59	INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme
	60	INTEGER, PUBLIC, PARAMETER :: np_ENS = 2 ! ENS scheme
	61	INTEGER, PUBLIC, PARAMETER :: np_MIX = 3 ! MIX scheme
	62	INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme
[455]	63
[5836]	64	INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity
	65	! ! associated indices:
	66	INTEGER, PARAMETER :: np_COR = 1 ! Coriolis (planetary)
	67	INTEGER, PARAMETER :: np_RVO = 2 ! relative vorticity
	68	INTEGER, PARAMETER :: np_MET = 3 ! metric term
	69	INTEGER, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity)
	70	INTEGER, PARAMETER :: np_CME = 5 ! Coriolis + metric term
	71
	72	REAL(wp) :: r1_4 = 0.250_wp ! =1/4
	73	REAL(wp) :: r1_8 = 0.125_wp ! =1/8
	74	REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12
	75
[3]	76	!! * Substitutions
	77	# include "domzgr_substitute.h90"
	78	# include "vectopt_loop_substitute.h90"
	79	!!----------------------------------------------------------------------
[5836]	80	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
[1152]	81	!! $Id$
[2715]	82	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	83	!!----------------------------------------------------------------------
	84	CONTAINS
	85
[455]	86	SUBROUTINE dyn_vor( kt )
[3]	87	!!----------------------------------------------------------------------
	88	!!
[455]	89	!! ** Purpose : compute the lateral ocean tracer physics.
	90	!!
	91	!! ** Action : - Update (ua,va) with the now vorticity term trend
[503]	92	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[4990]	93	!! and planetary vorticity trends) and send them to trd_dyn
	94	!! for futher diagnostics (l_trddyn=T)
[503]	95	!!----------------------------------------------------------------------
[3294]	96	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	97	!
[3294]	98	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
[455]	99	!!----------------------------------------------------------------------
[2715]	100	!
[3294]	101	IF( nn_timing == 1 ) CALL timing_start('dyn_vor')
	102	!
	103	IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
	104	!
[5836]	105	SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==!
[643]	106	!
[5836]	107	CASE ( np_ENE ) !* energy conserving scheme
	108	IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two
[455]	109	ztrdu(:,:,:) = ua(:,:,:)
	110	ztrdv(:,:,:) = va(:,:,:)
[5836]	111	CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	112	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	113	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	114	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	115	ztrdu(:,:,:) = ua(:,:,:)
	116	ztrdv(:,:,:) = va(:,:,:)
[5836]	117	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	118	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	119	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	120	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	121	ELSE
[5836]	122	CALL vor_ene( kt, ntot, ua, va ) ! total vorticity trend
[455]	123	ENDIF
[643]	124	!
[5836]	125	CASE ( np_ENS ) !* enstrophy conserving scheme
	126	IF( l_trddyn ) THEN ! trend diagnostics: splitthe trend in two
[455]	127	ztrdu(:,:,:) = ua(:,:,:)
	128	ztrdv(:,:,:) = va(:,:,:)
[5836]	129	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	130	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	131	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	132	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	133	ztrdu(:,:,:) = ua(:,:,:)
	134	ztrdv(:,:,:) = va(:,:,:)
[5836]	135	CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	136	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	137	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	138	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	139	ELSE
[5836]	140	CALL vor_ens( kt, ntot, ua, va ) ! total vorticity trend
[455]	141	ENDIF
[643]	142	!
[5836]	143	CASE ( np_MIX ) !* mixed ene-ens scheme
	144	IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two
[455]	145	ztrdu(:,:,:) = ua(:,:,:)
	146	ztrdv(:,:,:) = va(:,:,:)
[5836]	147	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)
[455]	148	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	149	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	150	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	151	ztrdu(:,:,:) = ua(:,:,:)
	152	ztrdv(:,:,:) = va(:,:,:)
[5836]	153	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)
[455]	154	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	155	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	156	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	157	ELSE
[5836]	158	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)
	159	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)
	160	ENDIF
[643]	161	!
[5836]	162	CASE ( np_EEN ) !* energy and enstrophy conserving scheme
	163	IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two
[455]	164	ztrdu(:,:,:) = ua(:,:,:)
	165	ztrdv(:,:,:) = va(:,:,:)
[5836]	166	CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	167	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	168	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	169	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	170	ztrdu(:,:,:) = ua(:,:,:)
	171	ztrdv(:,:,:) = va(:,:,:)
[5836]	172	CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	173	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	174	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	175	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	176	ELSE
[5836]	177	CALL vor_een( kt, ntot, ua, va ) ! total vorticity trend
[455]	178	ENDIF
[643]	179	!
[455]	180	END SELECT
[2715]	181	!
[455]	182	! ! print sum trends (used for debugging)
[2715]	183	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, &
[455]	184	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]	185	!
[3294]	186	IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
	187	!
	188	IF( nn_timing == 1 ) CALL timing_stop('dyn_vor')
	189	!
[455]	190	END SUBROUTINE dyn_vor
	191
	192
[643]	193	SUBROUTINE vor_ene( kt, kvor, pua, pva )
[455]	194	!!----------------------------------------------------------------------
	195	!! * ROUTINE vor_ene *
	196	!!
[3]	197	!! ** Purpose : Compute the now total vorticity trend and add it to
	198	!! the general trend of the momentum equation.
	199	!!
	200	!! ** Method : Trend evaluated using now fields (centered in time)
[5836]	201	!! and the Sadourny (1975) flux form formulation : conserves the
	202	!! horizontal kinetic energy.
	203	!! The general trend of momentum is increased due to the vorticity
	204	!! term which is given by:
	205	!! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v vn) ]
	206	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u un) ]
	207	!! where rvor is the relative vorticity
[3]	208	!!
	209	!! ** Action : - Update (ua,va) with the now vorticity term trend
	210	!!
[503]	211	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	212	!!----------------------------------------------------------------------
[643]	213	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	214	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
[1438]	215	! ! =nrvm (relative vorticity or metric)
[643]	216	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	217	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[2715]	218	!
[5836]	219	INTEGER :: ji, jj, jk ! dummy loop indices
	220	REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
	221	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! 2D workspace
[3]	222	!!----------------------------------------------------------------------
[3294]	223	!
	224	IF( nn_timing == 1 ) CALL timing_start('vor_ene')
	225	!
	226	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
	227	!
[52]	228	IF( kt == nit000 ) THEN
	229	IF(lwp) WRITE(numout,*)
[455]	230	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
	231	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	232	ENDIF
[5836]	233	!
[3]	234	! ! ===============
	235	DO jk = 1, jpkm1 ! Horizontal slab
	236	! ! ===============
[1438]	237	!
[5836]	238	SELECT CASE( kvor ) !== vorticity considered ==!
	239	CASE ( np_COR ) !* Coriolis (planetary vorticity)
	240	zwz(:,:) = ff(:,:)
	241	CASE ( np_RVO ) !* relative vorticity
[643]	242	DO jj = 1, jpjm1
	243	DO ji = 1, fs_jpim1 ! vector opt.
[5836]	244	zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
	245	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	246	END DO
	247	END DO
	248	CASE ( np_MET ) !* metric term
	249	DO jj = 1, jpjm1
	250	DO ji = 1, fs_jpim1 ! vector opt.
[643]	251	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	252	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[5836]	253	& * 0.5 * r1_e1e2f(ji,jj)
[643]	254	END DO
	255	END DO
[5836]	256	CASE ( np_CRV ) !* Coriolis + relative vorticity
[643]	257	DO jj = 1, jpjm1
	258	DO ji = 1, fs_jpim1 ! vector opt.
[5836]	259	zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
	260	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) &
	261	& * r1_e1e2f(ji,jj)
[643]	262	END DO
	263	END DO
[5836]	264	CASE ( np_CME ) !* Coriolis + metric
	265	DO jj = 1, jpjm1
	266	DO ji = 1, fs_jpim1 ! vector opt.
	267	zwz(ji,jj) = ff(ji,jj) &
	268	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	269	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	270	& * 0.5 * r1_e1e2f(ji,jj)
	271	END DO
	272	END DO
	273	CASE DEFAULT ! error
	274	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	275	END SELECT
[5836]	276	!
	277	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
	278	DO jj = 1, jpjm1
	279	DO ji = 1, fs_jpim1 ! vector opt.
	280	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
	281	END DO
	282	END DO
	283	ENDIF
[455]	284
	285	IF( ln_sco ) THEN
	286	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
[3]	287	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	288	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	289	ELSE
	290	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	291	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
	292	ENDIF
[5836]	293	! !== compute and add the vorticity term trend =!
[3]	294	DO jj = 2, jpjm1
	295	DO ji = fs_2, fs_jpim1 ! vector opt.
	296	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	297	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	298	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	299	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
[5836]	300	pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	301	pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[3]	302	END DO
	303	END DO
	304	! ! ===============
	305	END DO ! End of slab
	306	! ! ===============
[3294]	307	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
[2715]	308	!
[3294]	309	IF( nn_timing == 1 ) CALL timing_stop('vor_ene')
	310	!
[455]	311	END SUBROUTINE vor_ene
[216]	312
	313
[643]	314	SUBROUTINE vor_ens( kt, kvor, pua, pva )
[3]	315	!!----------------------------------------------------------------------
[455]	316	!! * ROUTINE vor_ens *
[3]	317	!!
	318	!! ** Purpose : Compute the now total vorticity trend and add it to
	319	!! the general trend of the momentum equation.
	320	!!
	321	!! ** Method : Trend evaluated using now fields (centered in time)
	322	!! and the Sadourny (1975) flux FORM formulation : conserves the
	323	!! potential enstrophy of a horizontally non-divergent flow. the
	324	!! trend of the vorticity term is given by:
[5836]	325	!! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
	326	!! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
[3]	327	!! Add this trend to the general momentum trend (ua,va):
	328	!! (ua,va) = (ua,va) + ( voru , vorv )
	329	!!
	330	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
	331	!!
[503]	332	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	333	!!----------------------------------------------------------------------
[643]	334	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	335	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
	336	! ! =nrvm (relative vorticity or metric)
	337	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	338	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[2715]	339	!
[5836]	340	INTEGER :: ji, jj, jk ! dummy loop indices
	341	REAL(wp) :: zuav, zvau ! local scalars
	342	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww ! 2D workspace
[3]	343	!!----------------------------------------------------------------------
[3294]	344	!
	345	IF( nn_timing == 1 ) CALL timing_start('vor_ens')
	346	!
	347	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
	348	!
[52]	349	IF( kt == nit000 ) THEN
	350	IF(lwp) WRITE(numout,*)
[455]	351	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
	352	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	353	ENDIF
[3]	354	! ! ===============
	355	DO jk = 1, jpkm1 ! Horizontal slab
	356	! ! ===============
[1438]	357	!
[5836]	358	SELECT CASE( kvor ) !== vorticity considered ==!
	359	CASE ( np_COR ) !* Coriolis (planetary vorticity)
	360	zwz(:,:) = ff(:,:)
	361	CASE ( np_RVO ) !* relative vorticity
[643]	362	DO jj = 1, jpjm1
	363	DO ji = 1, fs_jpim1 ! vector opt.
[5836]	364	zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
	365	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj)
	366	END DO
	367	END DO
	368	CASE ( np_MET ) !* metric term
	369	DO jj = 1, jpjm1
	370	DO ji = 1, fs_jpim1 ! vector opt.
[643]	371	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	372	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[5836]	373	& * 0.5 * r1_e1e2f(ji,jj)
[643]	374	END DO
	375	END DO
[5836]	376	CASE ( np_CRV ) !* Coriolis + relative vorticity
[643]	377	DO jj = 1, jpjm1
	378	DO ji = 1, fs_jpim1 ! vector opt.
[5836]	379	zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
	380	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) &
	381	& * r1_e1e2f(ji,jj)
[643]	382	END DO
	383	END DO
[5836]	384	CASE ( np_CME ) !* Coriolis + metric
	385	DO jj = 1, jpjm1
	386	DO ji = 1, fs_jpim1 ! vector opt.
	387	zwz(ji,jj) = ff(ji,jj) &
	388	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	389	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	390	& * 0.5 * r1_e1e2f(ji,jj)
	391	END DO
	392	END DO
	393	CASE DEFAULT ! error
	394	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	395	END SELECT
[1438]	396	!
[5836]	397	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
	398	DO jj = 1, jpjm1
	399	DO ji = 1, fs_jpim1 ! vector opt.
	400	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
[3]	401	END DO
	402	END DO
[5836]	403	ENDIF
	404	!
	405	IF( ln_sco ) THEN !== horizontal fluxes ==!
	406	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
	407	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	408	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[3]	409	ELSE
[5836]	410	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	411	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
[3]	412	ENDIF
[5836]	413	! !== compute and add the vorticity term trend =!
[3]	414	DO jj = 2, jpjm1
	415	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	416	zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
	417	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
	418	zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
	419	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
[455]	420	pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	421	pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
[3]	422	END DO
	423	END DO
	424	! ! ===============
	425	END DO ! End of slab
	426	! ! ===============
[3294]	427	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
[2715]	428	!
[3294]	429	IF( nn_timing == 1 ) CALL timing_stop('vor_ens')
	430	!
[455]	431	END SUBROUTINE vor_ens
[216]	432
	433
[643]	434	SUBROUTINE vor_een( kt, kvor, pua, pva )
[108]	435	!!----------------------------------------------------------------------
[455]	436	!! * ROUTINE vor_een *
[108]	437	!!
	438	!! ** Purpose : Compute the now total vorticity trend and add it to
	439	!! the general trend of the momentum equation.
	440	!!
	441	!! ** Method : Trend evaluated using now fields (centered in time)
[1438]	442	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]	443	!! both the horizontal kinetic energy and the potential enstrophy
[1438]	444	!! when horizontal divergence is zero (see the NEMO documentation)
	445	!! Add this trend to the general momentum trend (ua,va).
[108]	446	!!
	447	!! ** Action : - Update (ua,va) with the now vorticity term trend
	448	!!
[503]	449	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	450	!!----------------------------------------------------------------------
[643]	451	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[5836]	452	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric)
[643]	453	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	454	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[5836]	455	!
	456	INTEGER :: ji, jj, jk ! dummy loop indices
	457	INTEGER :: ierr ! local integer
	458	REAL(wp) :: zua, zva ! local scalars
	459	REAL(wp) :: zmsk, ze3 ! local scalars
	460	!
	461	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, z1_e3f
	462	REAL(wp), POINTER, DIMENSION(:,:) :: ztnw, ztne, ztsw, ztse
[108]	463	!!----------------------------------------------------------------------
[3294]	464	!
	465	IF( nn_timing == 1 ) CALL timing_start('vor_een')
	466	!
[5836]	467	CALL wrk_alloc( jpi,jpj, zwx , zwy , zwz , z1_e3f )
	468	CALL wrk_alloc( jpi,jpj, ztnw, ztne, ztsw, ztse )
[3294]	469	!
[108]	470	IF( kt == nit000 ) THEN
	471	IF(lwp) WRITE(numout,*)
[455]	472	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	473	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[1438]	474	ENDIF
[5836]	475	!
	476	! ! ===============
	477	DO jk = 1, jpkm1 ! Horizontal slab
	478	! ! ===============
	479	!
	480	SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point
	481	CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
	482	DO jj = 1, jpjm1
	483	DO ji = 1, fs_jpim1 ! vector opt.
	484	ze3 = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
	485	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
	486	IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4.0_wp / ze3
	487	ELSE ; z1_e3f(ji,jj) = 0.0_wp
	488	ENDIF
[108]	489	END DO
	490	END DO
[5836]	491	CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
	492	DO jj = 1, jpjm1
	493	DO ji = 1, fs_jpim1 ! vector opt.
	494	ze3 = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
	495	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
	496	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	497	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
	498	IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3
	499	ELSE ; z1_e3f(ji,jj) = 0.0_wp
	500	ENDIF
[5029]	501	END DO
	502	END DO
[5836]	503	END SELECT
	504	!
	505	SELECT CASE( kvor ) !== vorticity considered ==!
	506	CASE ( np_COR ) !* Coriolis (planetary vorticity)
[643]	507	DO jj = 1, jpjm1
	508	DO ji = 1, fs_jpim1 ! vector opt.
[5836]	509	zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj)
	510	END DO
	511	END DO
	512	CASE ( np_RVO ) !* relative vorticity
	513	DO jj = 1, jpjm1
	514	DO ji = 1, fs_jpim1 ! vector opt.
	515	zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
	516	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) &
	517	& * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
	518	END DO
	519	END DO
	520	CASE ( np_MET ) !* metric term
	521	DO jj = 1, jpjm1
	522	DO ji = 1, fs_jpim1 ! vector opt.
[643]	523	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	524	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[5836]	525	& * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj)
[643]	526	END DO
	527	END DO
[5836]	528	CASE ( np_CRV ) !* Coriolis + relative vorticity
[643]	529	DO jj = 1, jpjm1
	530	DO ji = 1, fs_jpim1 ! vector opt.
[5836]	531	zwz(ji,jj) = ( ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) &
	532	& - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) &
	533	& * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
[643]	534	END DO
	535	END DO
[5836]	536	CASE ( np_CME ) !* Coriolis + metric
	537	DO jj = 1, jpjm1
	538	DO ji = 1, fs_jpim1 ! vector opt.
	539	zwz(ji,jj) = ( ff(ji,jj) &
	540	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	541	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	542	& * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
	543	END DO
	544	END DO
	545	CASE DEFAULT ! error
	546	CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
[455]	547	END SELECT
[5836]	548	!
	549	CALL lbc_lnk( zwz, 'F', 1. )
	550	!
	551	IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!
	552	DO jj = 1, jpjm1
	553	DO ji = 1, fs_jpim1 ! vector opt.
	554	zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
	555	END DO
	556	END DO
	557	ENDIF
	558	!
	559	! !== horizontal fluxes ==!
[108]	560	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	561	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	562
[5836]	563	! !== compute and add the vorticity term trend =!
[1438]	564	jj = 2
	565	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
[5836]	566	DO ji = 2, jpi ! split in 2 parts due to vector opt.
[108]	567	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	568	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	569	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	570	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	571	END DO
	572	DO jj = 3, jpj
[1694]	573	DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3
[108]	574	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	575	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	576	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	577	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	578	END DO
	579	END DO
	580	DO jj = 2, jpjm1
	581	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	582	zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	583	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	584	zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	585	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
[455]	586	pua(ji,jj,jk) = pua(ji,jj,jk) + zua
	587	pva(ji,jj,jk) = pva(ji,jj,jk) + zva
[108]	588	END DO
	589	END DO
	590	! ! ===============
	591	END DO ! End of slab
	592	! ! ===============
[2715]	593	!
[5836]	594	CALL wrk_dealloc( jpi,jpj, zwx , zwy , zwz , z1_e3f )
	595	CALL wrk_dealloc( jpi,jpj, ztnw, ztne, ztsw, ztse )
	596	!
[3294]	597	IF( nn_timing == 1 ) CALL timing_stop('vor_een')
	598	!
[455]	599	END SUBROUTINE vor_een
[216]	600
	601
[2528]	602	SUBROUTINE dyn_vor_init
[3]	603	!!---------------------------------------------------------------------
[2528]	604	!! * ROUTINE dyn_vor_init *
[3]	605	!!
	606	!! ** Purpose : Control the consistency between cpp options for
[1438]	607	!! tracer advection schemes
[3]	608	!!----------------------------------------------------------------------
[2715]	609	INTEGER :: ioptio ! local integer
[3294]	610	INTEGER :: ji, jj, jk ! dummy loop indices
[4147]	611	INTEGER :: ios ! Local integer output status for namelist read
[2715]	612	!!
[5836]	613	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk
[3]	614	!!----------------------------------------------------------------------
	615
[4147]	616	REWIND( numnam_ref ) ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
	617	READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
	618	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
[3]	619
[4147]	620	REWIND( numnam_cfg ) ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
	621	READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
	622	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
[4624]	623	IF(lwm) WRITE ( numond, namdyn_vor )
[4147]	624
[503]	625	IF(lwp) THEN ! Namelist print
[3]	626	WRITE(numout,*)
[2528]	627	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
	628	WRITE(numout,*) '~~~~~~~~~~~~'
[4147]	629	WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme'
[5836]	630	WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
	631	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
	632	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
	633	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
	634	WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f
	635	WRITE(numout,*) ' masked (=1) or unmasked(=0) vorticity ln_dynvor_msk = ', ln_dynvor_msk
[52]	636	ENDIF
	637
[5836]	638	!!gm this should be removed when choosing a unique strategy for fmask at the coast
[3294]	639	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
	640	! at angles with three ocean points and one land point
[5836]	641	IF(lwp) WRITE(numout,*)
	642	IF(lwp) WRITE(numout,*) ' namlbc: change fmask value in the angles (T) ln_vorlat = ', ln_vorlat
[3294]	643	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
	644	DO jk = 1, jpk
	645	DO jj = 2, jpjm1
	646	DO ji = 2, jpim1
	647	IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
	648	fmask(ji,jj,jk) = 1._wp
	649	END DO
	650	END DO
	651	END DO
	652	!
	653	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	654	!
	655	ENDIF
[5836]	656	!!gm end
[3294]	657
[5836]	658	ioptio = 0 ! type of scheme for vorticity (set nvor_scheme)
	659	IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF
	660	IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF
	661	IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF
	662	IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF
	663	!
[503]	664	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
[5836]	665	!
	666	IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot)
	667	ncor = np_COR
[643]	668	IF( ln_dynadv_vec ) THEN
	669	IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'
[5836]	670	nrvm = np_RVO ! relative vorticity
	671	ntot = np_CRV ! relative + planetary vorticity
[643]	672	ELSE
	673	IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'
[5836]	674	nrvm = np_MET ! metric term
	675	ntot = np_CME ! Coriolis + metric term
[643]	676	ENDIF
	677
[503]	678	IF(lwp) THEN ! Print the choice
	679	WRITE(numout,*)
[5836]	680	IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' vorticity scheme ==>> energy conserving scheme'
	681	IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' vorticity scheme ==>> enstrophy conserving scheme'
	682	IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' vorticity scheme ==>> mixed enstrophy/energy conserving scheme'
	683	IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' vorticity scheme ==>> energy and enstrophy conserving scheme'
[3]	684	ENDIF
[503]	685	!
[2528]	686	END SUBROUTINE dyn_vor_init
[3]	687
[503]	688	!!==============================================================================
[3]	689	END MODULE dynvor

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

Download in other formats: