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

Context Navigation

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3317

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

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

Download in other formats: