Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynvor.F90 @ 7910

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

Download in other formats: