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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynvor.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

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

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

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

Download in other formats: