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.

dynnxt.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 4374

Last change on this file since 4374 was 4370, checked in by jchanut, 10 years ago
Time split cleaning / AMM12 restart / BDY and Agrif. See tickets #1228, #1227, #1225 and #1133
Property svn:keywords set to `Id`
File size: 17.1 KB

Rev	Line
[3]	1	MODULE dynnxt
[1502]	2	!!=========================================================================
[3]	3	!! * MODULE dynnxt *
	4	!! Ocean dynamics: time stepping
[1502]	5	!!=========================================================================
[1438]	6	!! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code
	7	!! ! 1990-10 (C. Levy, G. Madec)
	8	!! 7.0 ! 1993-03 (M. Guyon) symetrical conditions
	9	!! 8.0 ! 1997-02 (G. Madec & M. Imbard) opa, release 8.0
	10	!! 8.2 ! 1997-04 (A. Weaver) Euler forward step
	11	!! - ! 1997-06 (G. Madec) lateral boudary cond., lbc routine
	12	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	13	!! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond.
	14	!! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization
	15	!! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines.
[1502]	16	!! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option
[2528]	17	!! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module
[2723]	18	!! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
[4292]	19	!! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes
[1502]	20	!!-------------------------------------------------------------------------
[1438]	21
[1502]	22	!!-------------------------------------------------------------------------
	23	!! dyn_nxt : obtain the next (after) horizontal velocity
	24	!!-------------------------------------------------------------------------
[3]	25	USE oce ! ocean dynamics and tracers
	26	USE dom_oce ! ocean space and time domain
[2528]	27	USE sbc_oce ! Surface boundary condition: ocean fields
	28	USE phycst ! physical constants
[1502]	29	USE dynspg_oce ! type of surface pressure gradient
	30	USE dynadv ! dynamics: vector invariant versus flux form
	31	USE domvvl ! variable volume
[3294]	32	USE bdy_oce ! ocean open boundary conditions
	33	USE bdydta ! ocean open boundary conditions
	34	USE bdydyn ! ocean open boundary conditions
	35	USE bdyvol ! ocean open boundary condition (bdy_vol routines)
[1502]	36	USE in_out_manager ! I/O manager
[3]	37	USE lbclnk ! lateral boundary condition (or mpp link)
[2715]	38	USE lib_mpp ! MPP library
[3294]	39	USE wrk_nemo ! Memory Allocation
[258]	40	USE prtctl ! Print control
[4292]	41
[2528]	42	#if defined key_agrif
	43	USE agrif_opa_interp
	44	#endif
[3294]	45	USE timing ! Timing
[3]	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
[1438]	50	PUBLIC dyn_nxt ! routine called by step.F90
	51
[592]	52	!! * Substitutions
	53	# include "domzgr_substitute.h90"
[2715]	54	!!----------------------------------------------------------------------
[2528]	55	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1438]	56	!! $Id$
[2715]	57	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	58	!!----------------------------------------------------------------------
[3]	59	CONTAINS
	60
	61	SUBROUTINE dyn_nxt ( kt )
	62	!!----------------------------------------------------------------------
	63	!! * ROUTINE dyn_nxt *
	64	!!
[1502]	65	!! ** Purpose : Compute the after horizontal velocity. Apply the boundary
	66	!! condition on the after velocity, achieved the time stepping
	67	!! by applying the Asselin filter on now fields and swapping
	68	!! the fields.
[3]	69	!!
[1502]	70	!! ** Method : * After velocity is compute using a leap-frog scheme:
	71	!! (ua,va) = (ub,vb) + 2 rdt (ua,va)
	72	!! Note that with flux form advection and variable volume layer
	73	!! (lk_vvl=T), the leap-frog is applied on thickness weighted
	74	!! velocity.
	75	!! Note also that in filtered free surface (lk_dynspg_flt=T),
	76	!! the time stepping has already been done in dynspg module
[3]	77	!!
[1502]	78	!! * Apply lateral boundary conditions on after velocity
	79	!! at the local domain boundaries through lbc_lnk call,
[4328]	80	!! at the one-way open boundaries (lk_bdy=T),
[1502]	81	!! at the AGRIF zoom boundaries (lk_agrif=T)
[3]	82	!!
[1502]	83	!! * Apply the time filter applied and swap of the dynamics
	84	!! arrays to start the next time step:
	85	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
	86	!! (un,vn) = (ua,va).
	87	!! Note that with flux form advection and variable volume layer
	88	!! (lk_vvl=T), the time filter is applied on thickness weighted
	89	!! velocity.
	90	!!
	91	!! ** Action : ub,vb filtered before horizontal velocity of next time-step
	92	!! un,vn now horizontal velocity of next time-step
[3]	93	!!----------------------------------------------------------------------
	94	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	95	!
[3]	96	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	97	INTEGER :: iku, ikv ! local integers
[1566]	98	#if ! defined key_dynspg_flt
[3]	99	REAL(wp) :: z2dt ! temporary scalar
[1566]	100	#endif
[3294]	101	REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars
	102	REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - -
[4370]	103	REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva
[3294]	104	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f
[1502]	105	!!----------------------------------------------------------------------
[3294]	106	!
	107	IF( nn_timing == 1 ) CALL timing_start('dyn_nxt')
	108	!
	109	CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
[4370]	110	IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
[3294]	111	!
[3]	112	IF( kt == nit000 ) THEN
	113	IF(lwp) WRITE(numout,*)
	114	IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
	115	IF(lwp) WRITE(numout,*) '~~~~~~~'
	116	ENDIF
	117
[1502]	118	#if defined key_dynspg_flt
	119	!
	120	! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine
	121	! -------------
[3]	122
[1502]	123	! Update after velocity on domain lateral boundaries (only local domain required)
	124	! --------------------------------------------------
	125	CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries
	126	CALL lbc_lnk( va, 'V', -1. )
	127	!
	128	#else
[4292]	129
	130	# if defined key_dynspg_exp
[1502]	131	! Next velocity : Leap-frog time stepping
[1438]	132	! -------------
[1502]	133	z2dt = 2. * rdt ! Euler or leap-frog time step
	134	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
	135	!
	136	IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity
[1438]	137	DO jk = 1, jpkm1
[1502]	138	ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
	139	va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
	140	END DO
	141	ELSE ! applied on thickness weighted velocity
	142	DO jk = 1, jpkm1
	143	ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &
	144	& + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &
[1438]	145	& / fse3u_a(:,:,jk) * umask(:,:,jk)
[1502]	146	va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &
	147	& + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &
[1438]	148	& / fse3v_a(:,:,jk) * vmask(:,:,jk)
[592]	149	END DO
	150	ENDIF
[4292]	151	# endif
[592]	152
[4292]	153	# if defined key_dynspg_ts
	154	! Ensure below that barotropic velocities match time splitting estimate
	155	! Compute actual transport and replace it with ts estimate at "after" time step
	156	zua(:,:) = 0._wp
	157	zva(:,:) = 0._wp
[4370]	158	DO jk = 1, jpkm1
	159	zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
	160	zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
	161	END DO
	162	DO jk = 1, jpkm1
	163	ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
	164	va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
	165	END DO
[1502]	166
[4292]	167	IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
	168	! Remove advective velocity from "now velocities"
	169	! prior to asselin filtering
[4312]	170	! In the forward case, this is done below after asselin filtering
	171	! so that asselin contribution is removed at the same time
[4292]	172	DO jk = 1, jpkm1
	173	un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
	174	vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
	175	END DO
	176	ENDIF
	177	# endif
	178
[1502]	179	! Update after velocity on domain lateral boundaries
	180	! --------------------------------------------------
	181	CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries
	182	CALL lbc_lnk( va, 'V', -1. )
	183	!
[4328]	184	# if defined key_bdy
[1502]	185	! !* BDY open boundaries
[3764]	186	IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
	187	IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )
[3294]	188
	189	!!$ Do we need a call to bdy_vol here??
	190	!
[1438]	191	# endif
[1502]	192	!
[392]	193	# if defined key_agrif
[1502]	194	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
[389]	195	# endif
[3]	196	#endif
[592]	197
[1438]	198	! Time filter and swap of dynamics arrays
	199	! ------------------------------------------
[1502]	200	IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap
	201	DO jk = 1, jpkm1
	202	un(:,:,jk) = ua(:,:,jk) ! un <-- ua
[1438]	203	vn(:,:,jk) = va(:,:,jk)
	204	END DO
[4292]	205	IF (lk_vvl) THEN
	206	DO jk = 1, jpkm1
	207	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
	208	fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
	209	fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
	210	ENDDO
	211	ENDIF
[1502]	212	ELSE !* Leap-Frog : Asselin filter and swap
[2528]	213	! ! =============!
	214	IF( .NOT. lk_vvl ) THEN ! Fixed volume !
	215	! ! =============!
[1502]	216	DO jk = 1, jpkm1
[592]	217	DO jj = 1, jpj
[1502]	218	DO ji = 1, jpi
[4292]	219	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
	220	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
[1502]	221	!
	222	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	223	vb(ji,jj,jk) = zvf
	224	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	225	vn(ji,jj,jk) = va(ji,jj,jk)
	226	END DO
	227	END DO
	228	END DO
[2528]	229	! ! ================!
	230	ELSE ! Variable volume !
	231	! ! ================!
[4292]	232	! Before scale factor at t-points
	233	! (used as a now filtered scale factor until the swap)
	234	! ----------------------------------------------------
	235	IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
[4338]	236	! No asselin filtering on thicknesses if forward time splitting
[4292]	237	fse3t_b(:,:,:) = fse3t_n(:,:,:)
	238	ELSE
	239	fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
	240	! Add volume filter correction: compatibility with tracer advection scheme
	241	! => time filter + conservation correction (only at the first level)
	242	fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
[3294]	243	!
[4292]	244	ENDIF
[2528]	245	!
[4292]	246	IF( ln_dynadv_vec ) THEN
	247	! Before scale factor at (u/v)-points
	248	! -----------------------------------
	249	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
	250	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
	251	! Leap-Frog - Asselin filter and swap: applied on velocity
	252	! -----------------------------------
	253	DO jk = 1, jpkm1
	254	DO jj = 1, jpj
[2528]	255	DO ji = 1, jpi
[4292]	256	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
	257	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
[2528]	258	!
	259	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	260	vb(ji,jj,jk) = zvf
	261	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	262	vn(ji,jj,jk) = va(ji,jj,jk)
	263	END DO
	264	END DO
	265	END DO
	266	!
[4292]	267	ELSE
	268	! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
	269	!------------------------------------------------
	270	CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
	271	CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
	272	! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
	273	! ----------------------------------- ===========================
	274	DO jk = 1, jpkm1
	275	DO jj = 1, jpj
[4312]	276	DO ji = 1, jpi
[2528]	277	zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
	278	zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
	279	zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
	280	zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
	281	zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
	282	zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
	283	!
[3294]	284	zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)
	285	zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)
[2528]	286	!
[3294]	287	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
[2528]	288	vb(ji,jj,jk) = zvf
[3294]	289	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
[2528]	290	vn(ji,jj,jk) = va(ji,jj,jk)
	291	END DO
	292	END DO
	293	END DO
[3294]	294	fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor
	295	fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
[2528]	296	ENDIF
	297	!
[3]	298	ENDIF
[2528]	299	!
[4292]	300	IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
[4312]	301	! Revert "before" velocities to time split estimate
	302	! Doing it here also means that asselin filter contribution is removed
[4292]	303	zua(:,:) = 0._wp
	304	zva(:,:) = 0._wp
	305	DO jk = 1, jpkm1
[4370]	306	zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
	307	zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
	308	END DO
	309	DO jk = 1, jpkm1
[4292]	310	ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
	311	vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
	312	END DO
	313	ENDIF
	314	!
	315	ENDIF ! neuler =/0
[4354]	316	!
	317	! Set "now" and "before" barotropic velocities for next time step:
	318	! JC: Would be more clever to swap variables than to make a full vertical
	319	! integration
	320	!
[4370]	321	!
	322	IF (lk_vvl) THEN
	323	hu_b(:,:) = 0.
	324	hv_b(:,:) = 0.
	325	DO jk = 1, jpkm1
	326	hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
	327	hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
[4354]	328	END DO
[4370]	329	hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
	330	hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
[4354]	331	ENDIF
	332	!
	333	un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
	334	ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
	335	!
	336	DO jk = 1, jpkm1
	337	#if defined key_vectopt_loop
	338	DO jj = 1, 1 !Vector opt. => forced unrolling
	339	DO ji = 1, jpij
	340	#else
	341	DO jj = 1, jpj
	342	DO ji = 1, jpi
	343	#endif
[4370]	344	un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
	345	vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
[4354]	346	!
[4370]	347	ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
	348	vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
[4354]	349	END DO
	350	END DO
	351	END DO
	352	!
	353	!
[4370]	354	un_b(:,:) = un_b(:,:) * hur_a(:,:)
	355	vn_b(:,:) = vn_b(:,:) * hvr_a(:,:)
	356	ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
	357	vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
[4354]	358	!
	359	!
[1438]	360	IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
	361	& tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask )
	362	!
[3294]	363	CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
[4370]	364	IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
[2715]	365	!
[3294]	366	IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
	367	!
[3]	368	END SUBROUTINE dyn_nxt
	369
[1502]	370	!!=========================================================================
[3]	371	END MODULE dynnxt
[4354]	372

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

Download in other formats: