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_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3970

Last change on this file since 3970 was 3970, checked in by cbricaud, 11 years ago
Time splitting update, see ticket #1079
Property svn:keywords set to `Id`
File size: 18.2 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
[3970]	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
[367]	32	USE obc_oce ! ocean open boundary conditions
[3]	33	USE obcdyn ! open boundary condition for momentum (obc_dyn routine)
[367]	34	USE obcdyn_bt ! 2D open boundary condition for momentum (obc_dyn_bt routine)
	35	USE obcvol ! ocean open boundary condition (obc_vol routines)
[3294]	36	USE bdy_oce ! ocean open boundary conditions
	37	USE bdydta ! ocean open boundary conditions
	38	USE bdydyn ! ocean open boundary conditions
	39	USE bdyvol ! ocean open boundary condition (bdy_vol routines)
[1502]	40	USE in_out_manager ! I/O manager
[3]	41	USE lbclnk ! lateral boundary condition (or mpp link)
[2715]	42	USE lib_mpp ! MPP library
[3294]	43	USE wrk_nemo ! Memory Allocation
[258]	44	USE prtctl ! Print control
[3970]	45	USE dynspg_ts ! Barotropic velocities
	46
[2528]	47	#if defined key_agrif
	48	USE agrif_opa_interp
	49	#endif
[3294]	50	USE timing ! Timing
[3]	51
	52	IMPLICIT NONE
	53	PRIVATE
	54
[1438]	55	PUBLIC dyn_nxt ! routine called by step.F90
	56
[592]	57	!! * Substitutions
	58	# include "domzgr_substitute.h90"
[2715]	59	!!----------------------------------------------------------------------
[2528]	60	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1438]	61	!! $Id$
[2715]	62	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	63	!!----------------------------------------------------------------------
[3]	64	CONTAINS
	65
	66	SUBROUTINE dyn_nxt ( kt )
	67	!!----------------------------------------------------------------------
	68	!! * ROUTINE dyn_nxt *
	69	!!
[1502]	70	!! ** Purpose : Compute the after horizontal velocity. Apply the boundary
	71	!! condition on the after velocity, achieved the time stepping
	72	!! by applying the Asselin filter on now fields and swapping
	73	!! the fields.
[3]	74	!!
[1502]	75	!! ** Method : * After velocity is compute using a leap-frog scheme:
	76	!! (ua,va) = (ub,vb) + 2 rdt (ua,va)
	77	!! Note that with flux form advection and variable volume layer
	78	!! (lk_vvl=T), the leap-frog is applied on thickness weighted
	79	!! velocity.
	80	!! Note also that in filtered free surface (lk_dynspg_flt=T),
	81	!! the time stepping has already been done in dynspg module
[3]	82	!!
[1502]	83	!! * Apply lateral boundary conditions on after velocity
	84	!! at the local domain boundaries through lbc_lnk call,
[3294]	85	!! at the one-way open boundaries (lk_obc=T),
[1502]	86	!! at the AGRIF zoom boundaries (lk_agrif=T)
[3]	87	!!
[1502]	88	!! * Apply the time filter applied and swap of the dynamics
	89	!! arrays to start the next time step:
	90	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
	91	!! (un,vn) = (ua,va).
	92	!! Note that with flux form advection and variable volume layer
	93	!! (lk_vvl=T), the time filter is applied on thickness weighted
	94	!! velocity.
	95	!!
	96	!! ** Action : ub,vb filtered before horizontal velocity of next time-step
	97	!! un,vn now horizontal velocity of next time-step
[3]	98	!!----------------------------------------------------------------------
	99	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	100	!
[3]	101	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	102	INTEGER :: iku, ikv ! local integers
[1566]	103	#if ! defined key_dynspg_flt
[3]	104	REAL(wp) :: z2dt ! temporary scalar
[1566]	105	#endif
[3294]	106	REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars
	107	REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - -
[3970]	108	REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva
[3294]	109	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f
[1502]	110	!!----------------------------------------------------------------------
[3294]	111	!
	112	IF( nn_timing == 1 ) CALL timing_start('dyn_nxt')
	113	!
	114	CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
[3970]	115	IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
[3294]	116	!
[3]	117	IF( kt == nit000 ) THEN
	118	IF(lwp) WRITE(numout,*)
	119	IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
	120	IF(lwp) WRITE(numout,*) '~~~~~~~'
	121	ENDIF
	122
[1502]	123	#if defined key_dynspg_flt
	124	!
	125	! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine
	126	! -------------
[3]	127
[1502]	128	! Update after velocity on domain lateral boundaries (only local domain required)
	129	! --------------------------------------------------
	130	CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries
	131	CALL lbc_lnk( va, 'V', -1. )
	132	!
	133	#else
[3970]	134
	135	#if defined key_dynspg_exp
[1502]	136	! Next velocity : Leap-frog time stepping
[1438]	137	! -------------
[1502]	138	z2dt = 2. * rdt ! Euler or leap-frog time step
	139	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
	140	!
	141	IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity
[1438]	142	DO jk = 1, jpkm1
[1502]	143	ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
	144	va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
	145	END DO
	146	ELSE ! applied on thickness weighted velocity
	147	DO jk = 1, jpkm1
	148	ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &
	149	& + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &
[1438]	150	& / fse3u_a(:,:,jk) * umask(:,:,jk)
[1502]	151	va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &
	152	& + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &
[1438]	153	& / fse3v_a(:,:,jk) * vmask(:,:,jk)
[592]	154	END DO
	155	ENDIF
[3970]	156	#endif
[592]	157
[3970]	158	#if defined key_dynspg_ts
	159	! Ensure below that barotropic velocities match time splitting estimate
	160	! Compute actual transport and replace it with ts estimate at "after" time step
	161	zua(:,:) = 0._wp
	162	zva(:,:) = 0._wp
	163	IF (lk_vvl) THEN
	164	DO jk = 1, jpkm1
	165	zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
	166	zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
	167	END DO
	168	DO jk = 1, jpkm1
	169	ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) &
	170	& / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) + ua_b(:,:) ) * umask(:,:,jk)
	171	va(:,:,jk) = ( va(:,:,jk) - zva(:,:) &
	172	& / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) + va_b(:,:) ) * vmask(:,:,jk)
	173	END DO
	174	ELSE
	175	DO jk = 1, jpkm1
	176	zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
	177	zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
	178	END DO
	179	DO jk = 1, jpkm1
	180	ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk)
	181	va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk)
	182	END DO
	183	ENDIF
[1502]	184
[3970]	185	IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
	186	! Remove advective velocity from "now velocities"
	187	! prior to asselin filtering
	188	! In the forward case, this is done below after asselin filtering
	189	DO jk = 1, jpkm1
	190	un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
	191	vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
	192	END DO
	193	ENDIF
	194	#endif
	195
[1502]	196	! Update after velocity on domain lateral boundaries
	197	! --------------------------------------------------
	198	CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries
	199	CALL lbc_lnk( va, 'V', -1. )
	200	!
[3]	201	# if defined key_obc
[1502]	202	! !* OBC open boundaries
[3764]	203	IF( lk_obc ) CALL obc_dyn( kt )
[1502]	204	!
[2715]	205	IF( .NOT. lk_dynspg_flt ) THEN
[1502]	206	! Flather boundary condition : - Update sea surface height on each open boundary
[2715]	207	! sshn (= after ssh ) for explicit case (lk_dynspg_exp=T)
	208	! sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
[1502]	209	! - Correct the barotropic velocities
[3764]	210	IF( lk_obc ) CALL obc_dyn_bt( kt )
[1438]	211	!
[1502]	212	!!gm ERROR - potential BUG: sshn should not be modified at this stage !! ssh_nxt not alrady called
	213	CALL lbc_lnk( sshn, 'T', 1. ) ! Boundary conditions on sshn
[1438]	214	!
[3764]	215	IF( lk_obc .AND. ln_vol_cst ) CALL obc_vol( kt )
[1438]	216	!
	217	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask )
[367]	218	ENDIF
[1502]	219	!
[3294]	220	# elif defined key_bdy
[1502]	221	! !* BDY open boundaries
[3764]	222	IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
	223	IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )
[3294]	224
	225	!!$ Do we need a call to bdy_vol here??
	226	!
[1438]	227	# endif
[1502]	228	!
[392]	229	# if defined key_agrif
[1502]	230	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
[389]	231	# endif
[3]	232	#endif
[592]	233
[1438]	234	! Time filter and swap of dynamics arrays
	235	! ------------------------------------------
[1502]	236	IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap
	237	DO jk = 1, jpkm1
	238	un(:,:,jk) = ua(:,:,jk) ! un <-- ua
[1438]	239	vn(:,:,jk) = va(:,:,jk)
	240	END DO
[3970]	241	IF (lk_vvl) THEN
	242	DO jk = 1, jpkm1
	243	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
	244	fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
	245	fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
	246	ENDDO
	247	ENDIF
[1502]	248	ELSE !* Leap-Frog : Asselin filter and swap
[2528]	249	! ! =============!
	250	IF( .NOT. lk_vvl ) THEN ! Fixed volume !
	251	! ! =============!
[1502]	252	DO jk = 1, jpkm1
[592]	253	DO jj = 1, jpj
[1502]	254	DO ji = 1, jpi
[2528]	255	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
	256	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
[1502]	257	!
	258	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	259	vb(ji,jj,jk) = zvf
	260	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	261	vn(ji,jj,jk) = va(ji,jj,jk)
	262	END DO
	263	END DO
	264	END DO
[2528]	265	! ! ================!
	266	ELSE ! Variable volume !
	267	! ! ================!
[3294]	268	!
[3970]	269	IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
	270	! Remove asselin filtering on thicknesses if forward time splitting
	271	DO jk = 1, jpkm1
	272	fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
	273	ENDDO
	274	ELSE
	275
	276	DO jk = 1, jpkm1 ! Before scale factor at t-points
	277	fse3t_b(:,:,jk) = fse3t_n(:,:,jk) &
	278	& + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) &
	279	& - 2._wp * fse3t_n(:,:,jk) )
	280	END DO
	281	zec = atfp * rdt / rau0 ! Add filter correction only at the 1st level of t-point scale factors
	282	fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
	283	ENDIF
[2528]	284	!
[3294]	285	IF( ln_dynadv_vec ) THEN ! vector invariant form (no thickness weighted calulation)
	286	!
	287	! ! before scale factors at u- & v-pts (computed from fse3t_b)
	288	CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
	289	!
	290	DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: applied on velocity
	291	DO jj = 1, jpj ! --------
[2528]	292	DO ji = 1, jpi
	293	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
	294	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
	295	!
	296	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	297	vb(ji,jj,jk) = zvf
	298	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	299	vn(ji,jj,jk) = va(ji,jj,jk)
	300	END DO
	301	END DO
	302	END DO
	303	!
[3294]	304	ELSE ! flux form (thickness weighted calulation)
	305	!
	306	CALL dom_vvl_2( kt, ze3u_f, ze3v_f ) ! before scale factors at u- & v-pts (computed from fse3t_b)
	307	!
	308	DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap:
	309	DO jj = 1, jpj ! applied on thickness weighted velocity
[3764]	310	DO ji = 1, jpi ! ---------------------------
[2528]	311	zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
	312	zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
	313	zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
	314	zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
	315	zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
	316	zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
	317	!
[3294]	318	zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)
	319	zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)
[2528]	320	!
[3294]	321	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
[2528]	322	vb(ji,jj,jk) = zvf
[3294]	323	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
[2528]	324	vn(ji,jj,jk) = va(ji,jj,jk)
	325	END DO
	326	END DO
	327	END DO
[3294]	328	fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor
	329	fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
[2528]	330	ENDIF
	331	!
[3]	332	ENDIF
[2528]	333	!
[3970]	334	#if defined key_dynspg_ts
	335	IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
	336	! Remove asselin filtering of barotropic velocities if forward time splitting
	337	! note that we replace barotropic velocities by advective velocities
	338	zua(:,:) = 0._wp
	339	zva(:,:) = 0._wp
	340	IF (lk_vvl) THEN
	341	DO jk = 1, jpkm1
	342	zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
	343	zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
	344	END DO
	345	DO jk = 1, jpkm1
	346	ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) &
	347	& / ( hu_0(:,:) + sshu_n(:,:) + 1._wp - umask(:,:,1) ) - un_b(:,:)) * umask(:,:,jk)
	348	vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) &
	349	& / ( hv_0(:,:) + sshv_n(:,:) + 1._wp - vmask(:,:,1) ) - vn_b(:,:)) * vmask(:,:,jk)
	350	END DO
	351	ELSE
	352	DO jk = 1, jpkm1
	353	zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
	354	zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
	355	END DO
	356	DO jk = 1, jpkm1
	357	ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
	358	vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
	359	END DO
	360	ENDIF
	361	ENDIF
	362	#endif
	363	!
	364	ENDIF ! neuler =/0
[3]	365
[1438]	366	IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
	367	& tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask )
	368	!
[3294]	369	CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
[3970]	370	IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
[2715]	371	!
[3294]	372	IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
	373	!
[3]	374	END SUBROUTINE dyn_nxt
	375
[1502]	376	!!=========================================================================
[3]	377	END MODULE dynnxt

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

Download in other formats: