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 branches/UKMO/r5936_restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/r5936_restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90 @ 7114

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

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

Download in other formats: