Context Navigation

dynvor.F90 @ 2382

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

Download in other formats: