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

Context Navigation

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3802

Last change on this file since 3802 was 3764, checked in by smasson, 12 years ago
dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012
Property svn:keywords set to `Id`
File size: 14.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
[1502]	19	!!-------------------------------------------------------------------------
[1438]	20
[1502]	21	!!-------------------------------------------------------------------------
	22	!! dyn_nxt : obtain the next (after) horizontal velocity
	23	!!-------------------------------------------------------------------------
[3]	24	USE oce ! ocean dynamics and tracers
	25	USE dom_oce ! ocean space and time domain
[2528]	26	USE sbc_oce ! Surface boundary condition: ocean fields
	27	USE phycst ! physical constants
[1502]	28	USE dynspg_oce ! type of surface pressure gradient
	29	USE dynadv ! dynamics: vector invariant versus flux form
	30	USE domvvl ! variable volume
[367]	31	USE obc_oce ! ocean open boundary conditions
[3]	32	USE obcdyn ! open boundary condition for momentum (obc_dyn routine)
[367]	33	USE obcdyn_bt ! 2D open boundary condition for momentum (obc_dyn_bt routine)
	34	USE obcvol ! ocean open boundary condition (obc_vol routines)
[3294]	35	USE bdy_oce ! ocean open boundary conditions
	36	USE bdydta ! ocean open boundary conditions
	37	USE bdydyn ! ocean open boundary conditions
	38	USE bdyvol ! ocean open boundary condition (bdy_vol routines)
[1502]	39	USE in_out_manager ! I/O manager
[3]	40	USE lbclnk ! lateral boundary condition (or mpp link)
[2715]	41	USE lib_mpp ! MPP library
[3294]	42	USE wrk_nemo ! Memory Allocation
[258]	43	USE prtctl ! Print control
[2528]	44	#if defined key_agrif
	45	USE agrif_opa_interp
	46	#endif
[3294]	47	USE timing ! Timing
[3]	48
	49	IMPLICIT NONE
	50	PRIVATE
	51
[1438]	52	PUBLIC dyn_nxt ! routine called by step.F90
	53
[592]	54	!! * Substitutions
	55	# include "domzgr_substitute.h90"
[2715]	56	!!----------------------------------------------------------------------
[2528]	57	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1438]	58	!! $Id$
[2715]	59	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	60	!!----------------------------------------------------------------------
[3]	61	CONTAINS
	62
	63	SUBROUTINE dyn_nxt ( kt )
	64	!!----------------------------------------------------------------------
	65	!! * ROUTINE dyn_nxt *
	66	!!
[1502]	67	!! ** Purpose : Compute the after horizontal velocity. Apply the boundary
	68	!! condition on the after velocity, achieved the time stepping
	69	!! by applying the Asselin filter on now fields and swapping
	70	!! the fields.
[3]	71	!!
[1502]	72	!! ** Method : * After velocity is compute using a leap-frog scheme:
	73	!! (ua,va) = (ub,vb) + 2 rdt (ua,va)
	74	!! Note that with flux form advection and variable volume layer
	75	!! (lk_vvl=T), the leap-frog is applied on thickness weighted
	76	!! velocity.
	77	!! Note also that in filtered free surface (lk_dynspg_flt=T),
	78	!! the time stepping has already been done in dynspg module
[3]	79	!!
[1502]	80	!! * Apply lateral boundary conditions on after velocity
	81	!! at the local domain boundaries through lbc_lnk call,
[3294]	82	!! at the one-way open boundaries (lk_obc=T),
[1502]	83	!! at the AGRIF zoom boundaries (lk_agrif=T)
[3]	84	!!
[1502]	85	!! * Apply the time filter applied and swap of the dynamics
	86	!! arrays to start the next time step:
	87	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
	88	!! (un,vn) = (ua,va).
	89	!! Note that with flux form advection and variable volume layer
	90	!! (lk_vvl=T), the time filter is applied on thickness weighted
	91	!! velocity.
	92	!!
	93	!! ** Action : ub,vb filtered before horizontal velocity of next time-step
	94	!! un,vn now horizontal velocity of next time-step
[3]	95	!!----------------------------------------------------------------------
	96	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	97	!
[3]	98	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	99	INTEGER :: iku, ikv ! local integers
[1566]	100	#if ! defined key_dynspg_flt
[3]	101	REAL(wp) :: z2dt ! temporary scalar
[1566]	102	#endif
[3294]	103	REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars
	104	REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - -
	105	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f
[1502]	106	!!----------------------------------------------------------------------
[3294]	107	!
	108	IF( nn_timing == 1 ) CALL timing_start('dyn_nxt')
	109	!
	110	CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
	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
	129	! Next velocity : Leap-frog time stepping
[1438]	130	! -------------
[1502]	131	z2dt = 2. * rdt ! Euler or leap-frog time step
	132	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
	133	!
	134	IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity
[1438]	135	DO jk = 1, jpkm1
[1502]	136	ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
	137	va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
	138	END DO
	139	ELSE ! applied on thickness weighted velocity
	140	DO jk = 1, jpkm1
	141	ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &
	142	& + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &
[1438]	143	& / fse3u_a(:,:,jk) * umask(:,:,jk)
[1502]	144	va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &
	145	& + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &
[1438]	146	& / fse3v_a(:,:,jk) * vmask(:,:,jk)
[592]	147	END DO
	148	ENDIF
	149
[1502]	150
	151	! Update after velocity on domain lateral boundaries
	152	! --------------------------------------------------
	153	CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries
	154	CALL lbc_lnk( va, 'V', -1. )
	155	!
[3]	156	# if defined key_obc
[1502]	157	! !* OBC open boundaries
[3764]	158	IF( lk_obc ) CALL obc_dyn( kt )
[1502]	159	!
[2715]	160	IF( .NOT. lk_dynspg_flt ) THEN
[1502]	161	! Flather boundary condition : - Update sea surface height on each open boundary
[2715]	162	! sshn (= after ssh ) for explicit case (lk_dynspg_exp=T)
	163	! sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
[1502]	164	! - Correct the barotropic velocities
[3764]	165	IF( lk_obc ) CALL obc_dyn_bt( kt )
[1438]	166	!
[1502]	167	!!gm ERROR - potential BUG: sshn should not be modified at this stage !! ssh_nxt not alrady called
	168	CALL lbc_lnk( sshn, 'T', 1. ) ! Boundary conditions on sshn
[1438]	169	!
[3764]	170	IF( lk_obc .AND. ln_vol_cst ) CALL obc_vol( kt )
[1438]	171	!
	172	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask )
[367]	173	ENDIF
[1502]	174	!
[3294]	175	# elif defined key_bdy
[1502]	176	! !* BDY open boundaries
[3764]	177	IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
	178	IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )
[3294]	179
	180	!!$ Do we need a call to bdy_vol here??
	181	!
[1438]	182	# endif
[1502]	183	!
[392]	184	# if defined key_agrif
[1502]	185	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
[389]	186	# endif
[3]	187	#endif
[592]	188
[1438]	189	! Time filter and swap of dynamics arrays
	190	! ------------------------------------------
[1502]	191	IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap
	192	DO jk = 1, jpkm1
	193	un(:,:,jk) = ua(:,:,jk) ! un <-- ua
[1438]	194	vn(:,:,jk) = va(:,:,jk)
	195	END DO
[1502]	196	ELSE !* Leap-Frog : Asselin filter and swap
[2528]	197	! ! =============!
	198	IF( .NOT. lk_vvl ) THEN ! Fixed volume !
	199	! ! =============!
[1502]	200	DO jk = 1, jpkm1
[592]	201	DO jj = 1, jpj
[1502]	202	DO ji = 1, jpi
[2528]	203	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
	204	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
[1502]	205	!
	206	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	207	vb(ji,jj,jk) = zvf
	208	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	209	vn(ji,jj,jk) = va(ji,jj,jk)
	210	END DO
	211	END DO
	212	END DO
[2528]	213	! ! ================!
	214	ELSE ! Variable volume !
	215	! ! ================!
[3294]	216	!
	217	DO jk = 1, jpkm1 ! Before scale factor at t-points
[2528]	218	fse3t_b(:,:,jk) = fse3t_n(:,:,jk) &
	219	& + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) &
[3294]	220	& - 2._wp * fse3t_n(:,:,jk) )
	221	END DO
	222	zec = atfp * rdt / rau0 ! Add filter correction only at the 1st level of t-point scale factors
[2528]	223	fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
	224	!
[3294]	225	IF( ln_dynadv_vec ) THEN ! vector invariant form (no thickness weighted calulation)
	226	!
	227	! ! before scale factors at u- & v-pts (computed from fse3t_b)
	228	CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
	229	!
	230	DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: applied on velocity
	231	DO jj = 1, jpj ! --------
[2528]	232	DO ji = 1, jpi
	233	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
	234	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
	235	!
	236	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	237	vb(ji,jj,jk) = zvf
	238	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	239	vn(ji,jj,jk) = va(ji,jj,jk)
	240	END DO
	241	END DO
	242	END DO
	243	!
[3294]	244	ELSE ! flux form (thickness weighted calulation)
	245	!
	246	CALL dom_vvl_2( kt, ze3u_f, ze3v_f ) ! before scale factors at u- & v-pts (computed from fse3t_b)
	247	!
	248	DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap:
	249	DO jj = 1, jpj ! applied on thickness weighted velocity
[3764]	250	DO ji = 1, jpi ! ---------------------------
[2528]	251	zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
	252	zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
	253	zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
	254	zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
	255	zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
	256	zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
	257	!
[3294]	258	zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)
	259	zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)
[2528]	260	!
[3294]	261	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
[2528]	262	vb(ji,jj,jk) = zvf
[3294]	263	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
[2528]	264	vn(ji,jj,jk) = va(ji,jj,jk)
	265	END DO
	266	END DO
	267	END DO
[3294]	268	fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor
	269	fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
[2528]	270	ENDIF
	271	!
[3]	272	ENDIF
[2528]	273	!
[258]	274	ENDIF
[3]	275
[1438]	276	IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
	277	& tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask )
	278	!
[3294]	279	CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
[2715]	280	!
[3294]	281	IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
	282	!
[3]	283	END SUBROUTINE dyn_nxt
	284
[1502]	285	!!=========================================================================
[3]	286	END MODULE dynnxt

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

Download in other formats: