Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynvor.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 43.4 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
[503]	18	!!----------------------------------------------------------------------
[3]	19
	20	!!----------------------------------------------------------------------
[2528]	21	!! dyn_vor : Update the momentum trend with the vorticity trend
	22	!! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T)
	23	!! vor_ene : energy conserving scheme (ln_dynvor_ene=T)
	24	!! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=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
[503]	34	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	35	USE prtctl ! Print control
	36	USE in_out_manager ! I/O manager
[3294]	37	USE lib_mpp ! MPP library
	38	USE wrk_nemo ! Memory Allocation
	39	USE timing ! Timing
[8280]	40	USE lib_fortran
[3]	41
[3294]	42
[11738]	43	USE yomhook, ONLY: lhook, dr_hook
	44	USE parkind1, ONLY: jprb, jpim
	45
[3]	46	IMPLICIT NONE
	47	PRIVATE
	48
[2528]	49	PUBLIC dyn_vor ! routine called by step.F90
	50	PUBLIC dyn_vor_init ! routine called by opa.F90
[3]	51
[4147]	52	! !!* Namelist namdyn_vor: vorticity term
	53	LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme
	54	LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme
	55	LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme
	56	LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme
[5029]	57	LOGICAL, PUBLIC :: ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation)
[3]	58
[503]	59	INTEGER :: nvor = 0 ! type of vorticity trend used
[643]	60	INTEGER :: ncor = 1 ! coriolis
	61	INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term
	62	INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term
[455]	63
[3]	64	!! * Substitutions
	65	# include "domzgr_substitute.h90"
	66	# include "vectopt_loop_substitute.h90"
	67	!!----------------------------------------------------------------------
[2528]	68	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	69	!! $Id$
[2715]	70	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	71	!!----------------------------------------------------------------------
	72	CONTAINS
	73
[455]	74	SUBROUTINE dyn_vor( kt )
[3]	75	!!----------------------------------------------------------------------
	76	!!
[455]	77	!! ** Purpose : compute the lateral ocean tracer physics.
	78	!!
	79	!! ** Action : - Update (ua,va) with the now vorticity term trend
[503]	80	!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
[4990]	81	!! and planetary vorticity trends) and send them to trd_dyn
	82	!! for futher diagnostics (l_trddyn=T)
[503]	83	!!----------------------------------------------------------------------
[3294]	84	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	85	!
[3294]	86	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
[11738]	87	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	88	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	89	REAL(KIND=jprb) :: zhook_handle
	90
	91	CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_VOR'
	92
	93	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	94
[455]	95	!!----------------------------------------------------------------------
[2715]	96	!
[3294]	97	IF( nn_timing == 1 ) CALL timing_start('dyn_vor')
	98	!
	99	IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
	100	!
[643]	101	! ! vorticity term
[455]	102	SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend
[643]	103	!
[455]	104	CASE ( -1 ) ! esopa: test all possibility with control print
[643]	105	CALL vor_ene( kt, ntot, ua, va )
[503]	106	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, &
	107	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]	108	CALL vor_ens( kt, ntot, ua, va )
[503]	109	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, &
	110	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[455]	111	CALL vor_mix( kt )
[503]	112	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, &
	113	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]	114	CALL vor_een( kt, ntot, ua, va )
[503]	115	CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, &
	116	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[643]	117	!
[455]	118	CASE ( 0 ) ! energy conserving scheme
	119	IF( l_trddyn ) THEN
	120	ztrdu(:,:,:) = ua(:,:,:)
	121	ztrdv(:,:,:) = va(:,:,:)
[643]	122	CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	123	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	124	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	125	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	126	ztrdu(:,:,:) = ua(:,:,:)
	127	ztrdv(:,:,:) = va(:,:,:)
[643]	128	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	129	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	130	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	131	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	132	ELSE
[643]	133	CALL vor_ene( kt, ntot, ua, va ) ! total vorticity
[455]	134	ENDIF
[643]	135	!
[455]	136	CASE ( 1 ) ! enstrophy conserving scheme
	137	IF( l_trddyn ) THEN
	138	ztrdu(:,:,:) = ua(:,:,:)
	139	ztrdv(:,:,:) = va(:,:,:)
[643]	140	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	141	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	142	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	143	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	144	ztrdu(:,:,:) = ua(:,:,:)
	145	ztrdv(:,:,:) = va(:,:,:)
[643]	146	CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	147	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	148	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	149	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	150	ELSE
[643]	151	CALL vor_ens( kt, ntot, ua, va ) ! total vorticity
[455]	152	ENDIF
[643]	153	!
[455]	154	CASE ( 2 ) ! mixed ene-ens scheme
	155	IF( l_trddyn ) THEN
	156	ztrdu(:,:,:) = ua(:,:,:)
	157	ztrdv(:,:,:) = va(:,:,:)
[643]	158	CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)
[455]	159	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	160	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	161	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	162	ztrdu(:,:,:) = ua(:,:,:)
	163	ztrdv(:,:,:) = va(:,:,:)
[643]	164	CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)
[455]	165	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	166	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	167	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	168	ELSE
	169	CALL vor_mix( kt ) ! total vorticity (mix=ens-ene)
	170	ENDIF
[643]	171	!
[455]	172	CASE ( 3 ) ! energy and enstrophy conserving scheme
	173	IF( l_trddyn ) THEN
	174	ztrdu(:,:,:) = ua(:,:,:)
	175	ztrdv(:,:,:) = va(:,:,:)
[643]	176	CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend
[455]	177	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	178	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	179	CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
[455]	180	ztrdu(:,:,:) = ua(:,:,:)
	181	ztrdv(:,:,:) = va(:,:,:)
[643]	182	CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend
[455]	183	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	184	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	185	CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
[455]	186	ELSE
[643]	187	CALL vor_een( kt, ntot, ua, va ) ! total vorticity
[455]	188	ENDIF
[643]	189	!
[455]	190	END SELECT
[2715]	191	!
[455]	192	! ! print sum trends (used for debugging)
[2715]	193	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, &
[455]	194	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[1438]	195	!
[3294]	196	IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
	197	!
	198	IF( nn_timing == 1 ) CALL timing_stop('dyn_vor')
	199	!
[11738]	200	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[455]	201	END SUBROUTINE dyn_vor
	202
	203
[643]	204	SUBROUTINE vor_ene( kt, kvor, pua, pva )
[455]	205	!!----------------------------------------------------------------------
	206	!! * ROUTINE vor_ene *
	207	!!
[3]	208	!! ** Purpose : Compute the now total vorticity trend and add it to
	209	!! the general trend of the momentum equation.
	210	!!
	211	!! ** Method : Trend evaluated using now fields (centered in time)
	212	!! and the Sadourny (1975) flux form formulation : conserves the
	213	!! horizontal kinetic energy.
	214	!! The trend of the vorticity term is given by:
[455]	215	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
[3]	216	!! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ]
	217	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ]
	218	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	219	!! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ]
	220	!! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ]
	221	!! Add this trend to the general momentum trend (ua,va):
	222	!! (ua,va) = (ua,va) + ( voru , vorv )
	223	!!
	224	!! ** Action : - Update (ua,va) with the now vorticity term trend
	225	!!
[503]	226	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	227	!!----------------------------------------------------------------------
[2715]	228	!
[643]	229	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	230	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
[1438]	231	! ! =nrvm (relative vorticity or metric)
[643]	232	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	233	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[2715]	234	!
	235	INTEGER :: ji, jj, jk ! dummy loop indices
	236	REAL(wp) :: zx1, zy1, zfact2, zx2, zy2 ! local scalars
[3294]	237	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz
[11738]	238	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	239	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	240	REAL(KIND=jprb) :: zhook_handle
	241
	242	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_ENE'
	243
	244	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	245
[3]	246	!!----------------------------------------------------------------------
[3294]	247	!
	248	IF( nn_timing == 1 ) CALL timing_start('vor_ene')
	249	!
	250	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
	251	!
[52]	252	IF( kt == nit000 ) THEN
	253	IF(lwp) WRITE(numout,*)
[455]	254	IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme'
	255	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	256	ENDIF
[3]	257
[1438]	258	zfact2 = 0.5 * 0.5 ! Local constant initialization
[216]	259
[455]	260	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]	261	! ! ===============
	262	DO jk = 1, jpkm1 ! Horizontal slab
	263	! ! ===============
[1438]	264	!
[3]	265	! Potential vorticity and horizontal fluxes
	266	! -----------------------------------------
[643]	267	SELECT CASE( kvor ) ! vorticity considered
	268	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
	269	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
	270	CASE ( 3 ) ! metric term
	271	DO jj = 1, jpjm1
	272	DO ji = 1, fs_jpim1 ! vector opt.
	273	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	274	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	275	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
	276	END DO
	277	END DO
	278	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
	279	CASE ( 5 ) ! total (coriolis + metric)
	280	DO jj = 1, jpjm1
	281	DO ji = 1, fs_jpim1 ! vector opt.
	282	zwz(ji,jj) = ( ff (ji,jj) &
	283	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	284	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	285	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
	286	& )
	287	END DO
	288	END DO
[455]	289	END SELECT
	290
	291	IF( ln_sco ) THEN
	292	zwz(:,:) = zwz(:,:) / fse3f(:,:,jk)
[3]	293	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	294	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	295	ELSE
	296	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	297	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
	298	ENDIF
	299
	300	! Compute and add the vorticity term trend
	301	! ----------------------------------------
	302	DO jj = 2, jpjm1
	303	DO ji = fs_2, fs_jpim1 ! vector opt.
	304	zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1)
	305	zy2 = zwy(ji,jj ) + zwy(ji+1,jj )
	306	zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1)
	307	zx2 = zwx(ji ,jj) + zwx(ji ,jj+1)
[455]	308	pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	309	pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[3]	310	END DO
	311	END DO
	312	! ! ===============
	313	END DO ! End of slab
	314	! ! ===============
[3294]	315	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
[2715]	316	!
[3294]	317	IF( nn_timing == 1 ) CALL timing_stop('vor_ene')
	318	!
[11738]	319	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[455]	320	END SUBROUTINE vor_ene
[216]	321
	322
[455]	323	SUBROUTINE vor_mix( kt )
[3]	324	!!----------------------------------------------------------------------
[455]	325	!! * ROUTINE vor_mix *
[3]	326	!!
	327	!! ** Purpose : Compute the now total vorticity trend and add it to
	328	!! the general trend of the momentum equation.
	329	!!
	330	!! ** Method : Trend evaluated using now fields (centered in time)
	331	!! Mixte formulation : conserves the potential enstrophy of a hori-
	332	!! zontally non-divergent flow for (rotzu x uh), the relative vor-
	333	!! ticity term and the horizontal kinetic energy for (f x uh), the
	334	!! coriolis term. the now trend of the vorticity term is given by:
[455]	335	!! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:
[3]	336	!! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]
	337	!! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]
	338	!! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]
	339	!! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]
	340	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	341	!! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]
	342	!! +1/e1u mj-1[ f mi(e1v vn) ]
	343	!! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]
	344	!! +1/e2v mi-1[ f mj(e2u un) ]
	345	!! Add this now trend to the general momentum trend (ua,va):
	346	!! (ua,va) = (ua,va) + ( voru , vorv )
	347	!!
	348	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
	349	!!
[503]	350	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	351	!!----------------------------------------------------------------------
[2715]	352	!
[503]	353	INTEGER, INTENT(in) :: kt ! ocean timestep index
[2715]	354	!
[1438]	355	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	356	REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars
	357	REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - -
[3294]	358	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
[11738]	359	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	360	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	361	REAL(KIND=jprb) :: zhook_handle
	362
	363	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_MIX'
	364
	365	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	366
[3]	367	!!----------------------------------------------------------------------
[3294]	368	!
	369	IF( nn_timing == 1 ) CALL timing_start('vor_mix')
	370	!
	371	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )
	372	!
[52]	373	IF( kt == nit000 ) THEN
	374	IF(lwp) WRITE(numout,*)
[455]	375	IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'
	376	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	377	ENDIF
[3]	378
[1438]	379	zfact1 = 0.5 * 0.25 ! Local constant initialization
[3]	380	zfact2 = 0.5 * 0.5
	381
[455]	382	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )
[3]	383	! ! ===============
	384	DO jk = 1, jpkm1 ! Horizontal slab
	385	! ! ===============
[1438]	386	!
[3]	387	! Relative and planetary potential vorticity and horizontal fluxes
	388	! ----------------------------------------------------------------
[455]	389	IF( ln_sco ) THEN
[643]	390	IF( ln_dynadv_vec ) THEN
	391	zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)
	392	ELSE
	393	DO jj = 1, jpjm1
	394	DO ji = 1, fs_jpim1 ! vector opt.
	395	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	396	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	397	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )
	398	END DO
	399	END DO
	400	ENDIF
[3]	401	zwz(:,:) = ff (:,:) / fse3f(:,:,jk)
	402	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	403	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	404	ELSE
[643]	405	IF( ln_dynadv_vec ) THEN
	406	zww(:,:) = rotn(:,:,jk)
	407	ELSE
	408	DO jj = 1, jpjm1
	409	DO ji = 1, fs_jpim1 ! vector opt.
	410	zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	411	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	412	& * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )
	413	END DO
	414	END DO
	415	ENDIF
	416	zwz(:,:) = ff (:,:)
[3]	417	zwx(:,:) = e2u(:,:) * un(:,:,jk)
	418	zwy(:,:) = e1v(:,:) * vn(:,:,jk)
	419	ENDIF
	420
	421	! Compute and add the vorticity term trend
	422	! ----------------------------------------
	423	DO jj = 2, jpjm1
	424	DO ji = fs_2, fs_jpim1 ! vector opt.
	425	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
	426	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
	427	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
	428	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
	429	! enstrophy conserving formulation for relative vorticity term
	430	zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )
	431	zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )
	432	! energy conserving formulation for planetary vorticity term
	433	zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
	434	zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )
[503]	435	! mixed vorticity trend added to the momentum trends
[3]	436	ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua
	437	va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva
	438	END DO
	439	END DO
	440	! ! ===============
	441	END DO ! End of slab
	442	! ! ===============
[3294]	443	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )
[2715]	444	!
[3294]	445	IF( nn_timing == 1 ) CALL timing_stop('vor_mix')
	446	!
[11738]	447	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[455]	448	END SUBROUTINE vor_mix
[216]	449
	450
[643]	451	SUBROUTINE vor_ens( kt, kvor, pua, pva )
[3]	452	!!----------------------------------------------------------------------
[455]	453	!! * ROUTINE vor_ens *
[3]	454	!!
	455	!! ** Purpose : Compute the now total vorticity trend and add it to
	456	!! the general trend of the momentum equation.
	457	!!
	458	!! ** Method : Trend evaluated using now fields (centered in time)
	459	!! and the Sadourny (1975) flux FORM formulation : conserves the
	460	!! potential enstrophy of a horizontally non-divergent flow. the
	461	!! trend of the vorticity term is given by:
[455]	462	!! * s-coordinate (ln_sco=T), the e3. are inside the derivative:
[3]	463	!! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ]
	464	!! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ]
	465	!! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
	466	!! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ]
	467	!! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ]
	468	!! Add this trend to the general momentum trend (ua,va):
	469	!! (ua,va) = (ua,va) + ( voru , vorv )
	470	!!
	471	!! ** Action : - Update (ua,va) arrays with the now vorticity term trend
	472	!!
[503]	473	!! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
[3]	474	!!----------------------------------------------------------------------
[2715]	475	!
[643]	476	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	477	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
	478	! ! =nrvm (relative vorticity or metric)
	479	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	480	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[2715]	481	!
[503]	482	INTEGER :: ji, jj, jk ! dummy loop indices
	483	REAL(wp) :: zfact1, zuav, zvau ! temporary scalars
[3294]	484	REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww
[11738]	485	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	486	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	487	REAL(KIND=jprb) :: zhook_handle
	488
	489	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_ENS'
	490
	491	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	492
[3]	493	!!----------------------------------------------------------------------
[3294]	494	!
	495	IF( nn_timing == 1 ) CALL timing_start('vor_ens')
	496	!
	497	CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )
	498	!
[52]	499	IF( kt == nit000 ) THEN
	500	IF(lwp) WRITE(numout,*)
[455]	501	IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme'
	502	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[52]	503	ENDIF
[3]	504
[1438]	505	zfact1 = 0.5 * 0.25 ! Local constant initialization
[3]	506
[455]	507	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )
[3]	508	! ! ===============
	509	DO jk = 1, jpkm1 ! Horizontal slab
	510	! ! ===============
[1438]	511	!
[3]	512	! Potential vorticity and horizontal fluxes
	513	! -----------------------------------------
[643]	514	SELECT CASE( kvor ) ! vorticity considered
	515	CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis)
	516	CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity
	517	CASE ( 3 ) ! metric term
	518	DO jj = 1, jpjm1
	519	DO ji = 1, fs_jpim1 ! vector opt.
	520	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	521	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	522	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) )
	523	END DO
	524	END DO
	525	CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity)
	526	CASE ( 5 ) ! total (coriolis + metric)
	527	DO jj = 1, jpjm1
	528	DO ji = 1, fs_jpim1 ! vector opt.
	529	zwz(ji,jj) = ( ff (ji,jj) &
	530	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	531	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[1438]	532	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
[643]	533	& )
	534	END DO
	535	END DO
[455]	536	END SELECT
[1438]	537	!
[455]	538	IF( ln_sco ) THEN
[3]	539	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	540	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
[455]	541	zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk)
	542	zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk)
	543	zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk)
[3]	544	END DO
	545	END DO
	546	ELSE
	547	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	548	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
[455]	549	zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk)
	550	zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk)
[3]	551	END DO
	552	END DO
	553	ENDIF
[1438]	554	!
[3]	555	! Compute and add the vorticity term trend
	556	! ----------------------------------------
	557	DO jj = 2, jpjm1
	558	DO ji = fs_2, fs_jpim1 ! vector opt.
[455]	559	zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &
[503]	560	& + zwy(ji ,jj ) + zwy(ji+1,jj ) )
[455]	561	zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &
[503]	562	& + zwx(ji ,jj ) + zwx(ji ,jj+1) )
[455]	563	pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )
	564	pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )
[3]	565	END DO
	566	END DO
	567	! ! ===============
	568	END DO ! End of slab
	569	! ! ===============
[3294]	570	CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )
[2715]	571	!
[3294]	572	IF( nn_timing == 1 ) CALL timing_stop('vor_ens')
	573	!
[11738]	574	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[455]	575	END SUBROUTINE vor_ens
[216]	576
	577
[643]	578	SUBROUTINE vor_een( kt, kvor, pua, pva )
[108]	579	!!----------------------------------------------------------------------
[455]	580	!! * ROUTINE vor_een *
[108]	581	!!
	582	!! ** Purpose : Compute the now total vorticity trend and add it to
	583	!! the general trend of the momentum equation.
	584	!!
	585	!! ** Method : Trend evaluated using now fields (centered in time)
[1438]	586	!! and the Arakawa and Lamb (1980) flux form formulation : conserves
[108]	587	!! both the horizontal kinetic energy and the potential enstrophy
[1438]	588	!! when horizontal divergence is zero (see the NEMO documentation)
	589	!! Add this trend to the general momentum trend (ua,va).
[108]	590	!!
	591	!! ** Action : - Update (ua,va) with the now vorticity term trend
	592	!!
[503]	593	!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
	594	!!----------------------------------------------------------------------
[2715]	595	!
[643]	596	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	597	INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ;
[1438]	598	! ! =nrvm (relative vorticity or metric)
[643]	599	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend
	600	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend
[218]	601	!!
[3294]	602	INTEGER :: ji, jj, jk ! dummy loop indices
	603	INTEGER :: ierr ! local integer
	604	REAL(wp) :: zfac12, zua, zva ! local scalars
[4292]	605	REAL(wp) :: zmsk, ze3 ! local scalars
[3294]	606	! ! 3D workspace
	607	REAL(wp), POINTER , DIMENSION(:,: ) :: zwx, zwy, zwz
	608	REAL(wp), POINTER , DIMENSION(:,: ) :: ztnw, ztne, ztsw, ztse
	609	#if defined key_vvl
	610	REAL(wp), POINTER , DIMENSION(:,:,:) :: ze3f ! 3D workspace (lk_vvl=T)
[11738]	611	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	612	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	613	REAL(KIND=jprb) :: zhook_handle
	614
	615	CHARACTER(LEN=*), PARAMETER :: RoutineName='VOR_EEN'
	616
	617	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	618
[4292]	619	#else
[3294]	620	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all
[1438]	621	#endif
[108]	622	!!----------------------------------------------------------------------
[3294]	623	!
	624	IF( nn_timing == 1 ) CALL timing_start('vor_een')
	625	!
	626	CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz )
	627	CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
	628	#if defined key_vvl
	629	CALL wrk_alloc( jpi, jpj, jpk, ze3f )
	630	#endif
	631	!
[108]	632	IF( kt == nit000 ) THEN
	633	IF(lwp) WRITE(numout,*)
[455]	634	IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
	635	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[3802]	636	#if ! defined key_vvl
	637	IF( .NOT.ALLOCATED(ze3f) ) THEN
[2715]	638	ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )
	639	IF( lk_mpp ) CALL mpp_sum ( ierr )
	640	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
	641	ENDIF
[4990]	642	ze3f(:,:,:) = 0._wp
[3802]	643	#endif
[1438]	644	ENDIF
[108]	645
[4292]	646	IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
[5029]	647
	648	IF( ln_dynvor_een_old ) THEN ! original formulation
	649	DO jk = 1, jpk
	650	DO jj = 1, jpjm1
	651	DO ji = 1, jpim1
	652	ze3 = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
	653	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
	654	IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = 4.0_wp / ze3
	655	END DO
[108]	656	END DO
	657	END DO
[5029]	658	ELSE ! new formulation from NEMO 3.6
	659	DO jk = 1, jpk
	660	DO jj = 1, jpjm1
	661	DO ji = 1, jpim1
	662	ze3 = ( fse3t(ji,jj+1,jk)tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)tmask(ji+1,jj+1,jk) &
	663	& + fse3t(ji,jj ,jk)tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)tmask(ji+1,jj ,jk) )
	664	zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
	665	& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
	666	IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = zmsk / ze3
	667	END DO
	668	END DO
	669	END DO
	670	ENDIF
	671
[108]	672	CALL lbc_lnk( ze3f, 'F', 1. )
	673	ENDIF
	674
[2715]	675	zfac12 = 1._wp / 12._wp ! Local constant initialization
[216]	676
[108]	677
[455]	678	!CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse )
[108]	679	! ! ===============
	680	DO jk = 1, jpkm1 ! Horizontal slab
	681	! ! ===============
	682
	683	! Potential vorticity and horizontal fluxes
	684	! -----------------------------------------
[643]	685	SELECT CASE( kvor ) ! vorticity considered
[1438]	686	CASE ( 1 ) ! planetary vorticity (Coriolis)
	687	zwz(:,:) = ff(:,:) * ze3f(:,:,jk)
	688	CASE ( 2 ) ! relative vorticity
	689	zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk)
[643]	690	CASE ( 3 ) ! metric term
	691	DO jj = 1, jpjm1
	692	DO ji = 1, fs_jpim1 ! vector opt.
	693	zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	694	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
	695	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk)
	696	END DO
	697	END DO
[1516]	698	CALL lbc_lnk( zwz, 'F', 1. )
	699	CASE ( 4 ) ! total (relative + planetary vorticity)
[1438]	700	zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk)
[643]	701	CASE ( 5 ) ! total (coriolis + metric)
	702	DO jj = 1, jpjm1
	703	DO ji = 1, fs_jpim1 ! vector opt.
	704	zwz(ji,jj) = ( ff (ji,jj) &
	705	& + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &
	706	& - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &
[1438]	707	& * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) &
[643]	708	& ) * ze3f(ji,jj,jk)
	709	END DO
	710	END DO
[1516]	711	CALL lbc_lnk( zwz, 'F', 1. )
[455]	712	END SELECT
	713
[108]	714	zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
	715	zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
	716
	717	! Compute and add the vorticity term trend
	718	! ----------------------------------------
[1438]	719	jj = 2
	720	ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0
[108]	721	DO ji = 2, jpi
	722	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	723	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	724	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	725	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	726	END DO
	727	DO jj = 3, jpj
[1694]	728	DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3
[108]	729	ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
	730	ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
	731	ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
	732	ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
	733	END DO
	734	END DO
	735	DO jj = 2, jpjm1
	736	DO ji = fs_2, fs_jpim1 ! vector opt.
	737	zua = + zfac12 / e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
	738	& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
	739	zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
	740	& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
[455]	741	pua(ji,jj,jk) = pua(ji,jj,jk) + zua
	742	pva(ji,jj,jk) = pva(ji,jj,jk) + zva
[108]	743	END DO
	744	END DO
	745	! ! ===============
	746	END DO ! End of slab
	747	! ! ===============
[3294]	748	CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz )
	749	CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse )
	750	#if defined key_vvl
	751	CALL wrk_dealloc( jpi, jpj, jpk, ze3f )
	752	#endif
[2715]	753	!
[3294]	754	IF( nn_timing == 1 ) CALL timing_stop('vor_een')
	755	!
[11738]	756	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[455]	757	END SUBROUTINE vor_een
[216]	758
	759
[2528]	760	SUBROUTINE dyn_vor_init
[3]	761	!!---------------------------------------------------------------------
[2528]	762	!! * ROUTINE dyn_vor_init *
[3]	763	!!
	764	!! ** Purpose : Control the consistency between cpp options for
[1438]	765	!! tracer advection schemes
[3]	766	!!----------------------------------------------------------------------
[2715]	767	INTEGER :: ioptio ! local integer
[3294]	768	INTEGER :: ji, jj, jk ! dummy loop indices
[4147]	769	INTEGER :: ios ! Local integer output status for namelist read
[2715]	770	!!
[5029]	771	NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
[11738]	772	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
	773	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
	774	REAL(KIND=jprb) :: zhook_handle
	775
	776	CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_VOR_INIT'
	777
	778	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
	779
[3]	780	!!----------------------------------------------------------------------
	781
[4147]	782	REWIND( numnam_ref ) ! Namelist namdyn_vor in reference namelist : Vorticity scheme options
	783	READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901)
	784	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp )
[3]	785
[4147]	786	REWIND( numnam_cfg ) ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options
	787	READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 )
	788	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp )
[4624]	789	IF(lwm) WRITE ( numond, namdyn_vor )
[4147]	790
[503]	791	IF(lwp) THEN ! Namelist print
[3]	792	WRITE(numout,*)
[2528]	793	WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'
	794	WRITE(numout,*) '~~~~~~~~~~~~'
[4147]	795	WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme'
[503]	796	WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene
	797	WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens
	798	WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix
	799	WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een
[5029]	800	WRITE(numout,*) ' enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
[52]	801	ENDIF
	802
[3294]	803	! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks
	804	! at angles with three ocean points and one land point
	805	IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN
	806	DO jk = 1, jpk
	807	DO jj = 2, jpjm1
	808	DO ji = 2, jpim1
	809	IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) &
	810	fmask(ji,jj,jk) = 1._wp
	811	END DO
	812	END DO
	813	END DO
	814	!
	815	CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask
	816	!
	817	ENDIF
	818
[503]	819	ioptio = 0 ! Control of vorticity scheme options
	820	IF( ln_dynvor_ene ) ioptio = ioptio + 1
	821	IF( ln_dynvor_ens ) ioptio = ioptio + 1
	822	IF( ln_dynvor_mix ) ioptio = ioptio + 1
	823	IF( ln_dynvor_een ) ioptio = ioptio + 1
[5029]	824	IF( ln_dynvor_een_old ) ioptio = ioptio + 1
[503]	825	IF( lk_esopa ) ioptio = 1
	826
	827	IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )
	828
[643]	829	! ! Set nvor (type of scheme for vorticity)
[503]	830	IF( ln_dynvor_ene ) nvor = 0
	831	IF( ln_dynvor_ens ) nvor = 1
	832	IF( ln_dynvor_mix ) nvor = 2
[5029]	833	IF( ln_dynvor_een .or. ln_dynvor_een_old ) nvor = 3
[503]	834	IF( lk_esopa ) nvor = -1
	835
[643]	836	! ! Set ncor, nrvm, ntot (type of vorticity)
	837	IF(lwp) WRITE(numout,*)
	838	ncor = 1
	839	IF( ln_dynadv_vec ) THEN
	840	IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity'
	841	nrvm = 2
	842	ntot = 4
	843	ELSE
	844	IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term'
	845	nrvm = 3
	846	ntot = 5
	847	ENDIF
	848
[503]	849	IF(lwp) THEN ! Print the choice
	850	WRITE(numout,*)
[643]	851	IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme'
	852	IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme'
	853	IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme'
	854	IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme'
[503]	855	IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options'
[3]	856	ENDIF
[503]	857	!
[11738]	858	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[2528]	859	END SUBROUTINE dyn_vor_init
[3]	860
[503]	861	!!==============================================================================
[3]	862	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: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 11738

Download in other formats: