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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago
Update NEMOGCM from branch nemo_v3_3_beta
Property svn:keywords set to `Id`
File size: 16.3 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
[1502]	18	!!-------------------------------------------------------------------------
[1438]	19
[1502]	20	!!-------------------------------------------------------------------------
	21	!! dyn_nxt : obtain the next (after) horizontal velocity
	22	!!-------------------------------------------------------------------------
[3]	23	USE oce ! ocean dynamics and tracers
	24	USE dom_oce ! ocean space and time domain
[2528]	25	USE sbc_oce ! Surface boundary condition: ocean fields
	26	USE phycst ! physical constants
[1502]	27	USE dynspg_oce ! type of surface pressure gradient
	28	USE dynadv ! dynamics: vector invariant versus flux form
	29	USE domvvl ! variable volume
[367]	30	USE obc_oce ! ocean open boundary conditions
[3]	31	USE obcdyn ! open boundary condition for momentum (obc_dyn routine)
[367]	32	USE obcdyn_bt ! 2D open boundary condition for momentum (obc_dyn_bt routine)
	33	USE obcvol ! ocean open boundary condition (obc_vol routines)
[911]	34	USE bdy_oce ! unstructured open boundary conditions
	35	USE bdydta ! unstructured open boundary conditions
	36	USE bdydyn ! unstructured open boundary conditions
[1502]	37	USE in_out_manager ! I/O manager
[3]	38	USE lbclnk ! lateral boundary condition (or mpp link)
[258]	39	USE prtctl ! Print control
[2528]	40	#if defined key_agrif
	41	USE agrif_opa_interp
	42	#endif
[3]	43
	44	IMPLICIT NONE
	45	PRIVATE
	46
[1438]	47	PUBLIC dyn_nxt ! routine called by step.F90
	48
[592]	49	!! * Substitutions
	50	# include "domzgr_substitute.h90"
[1438]	51	!!-------------------------------------------------------------------------
[2528]	52	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1438]	53	!! $Id$
[2528]	54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1438]	55	!!-------------------------------------------------------------------------
[3]	56
	57	CONTAINS
	58
	59	SUBROUTINE dyn_nxt ( kt )
	60	!!----------------------------------------------------------------------
	61	!! * ROUTINE dyn_nxt *
	62	!!
[1502]	63	!! ** Purpose : Compute the after horizontal velocity. Apply the boundary
	64	!! condition on the after velocity, achieved the time stepping
	65	!! by applying the Asselin filter on now fields and swapping
	66	!! the fields.
[3]	67	!!
[1502]	68	!! ** Method : * After velocity is compute using a leap-frog scheme:
	69	!! (ua,va) = (ub,vb) + 2 rdt (ua,va)
	70	!! Note that with flux form advection and variable volume layer
	71	!! (lk_vvl=T), the leap-frog is applied on thickness weighted
	72	!! velocity.
	73	!! Note also that in filtered free surface (lk_dynspg_flt=T),
	74	!! the time stepping has already been done in dynspg module
[3]	75	!!
[1502]	76	!! * Apply lateral boundary conditions on after velocity
	77	!! at the local domain boundaries through lbc_lnk call,
	78	!! at the radiative open boundaries (lk_obc=T),
	79	!! at the relaxed open boundaries (lk_bdy=T), and
	80	!! at the AGRIF zoom boundaries (lk_agrif=T)
[3]	81	!!
[1502]	82	!! * Apply the time filter applied and swap of the dynamics
	83	!! arrays to start the next time step:
	84	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
	85	!! (un,vn) = (ua,va).
	86	!! Note that with flux form advection and variable volume layer
	87	!! (lk_vvl=T), the time filter is applied on thickness weighted
	88	!! velocity.
	89	!!
	90	!! ** Action : ub,vb filtered before horizontal velocity of next time-step
	91	!! un,vn now horizontal velocity of next time-step
[3]	92	!!----------------------------------------------------------------------
[2528]	93	USE oce, ONLY : ze3u_f => ta ! use ta as 3D workspace
	94	USE oce, ONLY : ze3v_f => sa ! use sa as 3D workspace
[3]	95	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1438]	96	!!
[3]	97	INTEGER :: ji, jj, jk ! dummy loop indices
[1566]	98	#if ! defined key_dynspg_flt
[3]	99	REAL(wp) :: z2dt ! temporary scalar
[1566]	100	#endif
[1438]	101	REAL(wp) :: zue3a , zue3n , zue3b ! temporary scalar
	102	REAL(wp) :: zve3a , zve3n , zve3b ! - -
	103	REAL(wp) :: zuf , zvf ! - -
[2528]	104	REAL(wp) :: zec ! - -
	105	REAL(wp) :: zv_t_ij , zv_t_ip1j ! - -
	106	REAL(wp) :: zv_t_ijp1 ! - -
	107	REAL(wp), DIMENSION(jpi,jpj) :: zs_t, zs_u_1, zs_v_1 ! temporary 2D workspace
[1502]	108	!!----------------------------------------------------------------------
[3]	109
	110	IF( kt == nit000 ) THEN
	111	IF(lwp) WRITE(numout,*)
	112	IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
	113	IF(lwp) WRITE(numout,*) '~~~~~~~'
	114	ENDIF
	115
[1502]	116	#if defined key_dynspg_flt
	117	!
	118	! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine
	119	! -------------
[3]	120
[1502]	121	! Update after velocity on domain lateral boundaries (only local domain required)
	122	! --------------------------------------------------
	123	CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries
	124	CALL lbc_lnk( va, 'V', -1. )
	125	!
	126	#else
	127	! Next velocity : Leap-frog time stepping
[1438]	128	! -------------
[1502]	129	z2dt = 2. * rdt ! Euler or leap-frog time step
	130	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
	131	!
	132	IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity
[1438]	133	DO jk = 1, jpkm1
[1502]	134	ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
	135	va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
	136	END DO
	137	ELSE ! applied on thickness weighted velocity
	138	DO jk = 1, jpkm1
	139	ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &
	140	& + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &
[1438]	141	& / fse3u_a(:,:,jk) * umask(:,:,jk)
[1502]	142	va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &
	143	& + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &
[1438]	144	& / fse3v_a(:,:,jk) * vmask(:,:,jk)
[592]	145	END DO
	146	ENDIF
	147
[1502]	148
	149	! Update after velocity on domain lateral boundaries
	150	! --------------------------------------------------
	151	CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries
	152	CALL lbc_lnk( va, 'V', -1. )
	153	!
[3]	154	# if defined key_obc
[1502]	155	! !* OBC open boundaries
[2528]	156	CALL obc_dyn( kt )
[1502]	157	!
[367]	158	IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
[1502]	159	! Flather boundary condition : - Update sea surface height on each open boundary
	160	! sshn (= after ssh ) for explicit case
	161	! sshn_b (= after ssha_b) for time-splitting case
	162	! - Correct the barotropic velocities
[367]	163	CALL obc_dyn_bt( kt )
[1438]	164	!
[1502]	165	!!gm ERROR - potential BUG: sshn should not be modified at this stage !! ssh_nxt not alrady called
	166	CALL lbc_lnk( sshn, 'T', 1. ) ! Boundary conditions on sshn
[1438]	167	!
	168	IF( ln_vol_cst ) CALL obc_vol( kt )
	169	!
	170	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask )
[367]	171	ENDIF
[1502]	172	!
[1125]	173	# elif defined key_bdy
[1502]	174	! !* BDY open boundaries
	175	IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN ! except for filtered option
[911]	176	CALL bdy_dyn_frs( kt )
	177	ENDIF
[1438]	178	# endif
[1502]	179	!
[392]	180	# if defined key_agrif
[1502]	181	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
[389]	182	# endif
[3]	183	#endif
[592]	184
[1438]	185	! Time filter and swap of dynamics arrays
	186	! ------------------------------------------
[1502]	187	IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap
	188	DO jk = 1, jpkm1
	189	un(:,:,jk) = ua(:,:,jk) ! un <-- ua
[1438]	190	vn(:,:,jk) = va(:,:,jk)
	191	END DO
[1502]	192	ELSE !* Leap-Frog : Asselin filter and swap
[2528]	193	! ! =============!
	194	IF( .NOT. lk_vvl ) THEN ! Fixed volume !
	195	! ! =============!
[1502]	196	DO jk = 1, jpkm1
[592]	197	DO jj = 1, jpj
[1502]	198	DO ji = 1, jpi
[2528]	199	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
	200	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
[1502]	201	!
	202	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	203	vb(ji,jj,jk) = zvf
	204	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	205	vn(ji,jj,jk) = va(ji,jj,jk)
	206	END DO
	207	END DO
	208	END DO
[2528]	209	! ! ================!
	210	ELSE ! Variable volume !
	211	! ! ================!
	212	! Before scale factor at t-points
	213	! -------------------------------
[1502]	214	DO jk = 1, jpkm1
[2528]	215	fse3t_b(:,:,jk) = fse3t_n(:,:,jk) &
	216	& + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) &
	217	& - 2.e0 * fse3t_n(:,:,jk) )
	218	ENDDO
	219	! Add volume filter correction only at the first level of t-point scale factors
	220	zec = atfp * rdt / rau0
	221	fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
	222	! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
	223	zs_t (:,:) = e1t(:,:) * e2t(:,:)
	224	zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
	225	zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
	226	!
	227	IF( ln_dynadv_vec ) THEN
	228	! Before scale factor at (u/v)-points
	229	! -----------------------------------
	230	! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
	231	DO jk = 1, jpkm1
	232	DO jj = 1, jpjm1
	233	DO ji = 1, jpim1
	234	zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk)
	235	zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk)
	236	zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk)
	237	fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
	238	fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
	239	END DO
[592]	240	END DO
	241	END DO
[2528]	242	! lateral boundary conditions
	243	CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
	244	CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
	245	! Add initial scale factor to scale factor anomaly
	246	fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
	247	fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
	248	! Leap-Frog - Asselin filter and swap: applied on velocity
	249	! -----------------------------------
	250	DO jk = 1, jpkm1
	251	DO jj = 1, jpj
	252	DO ji = 1, jpi
	253	zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
	254	zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
	255	!
	256	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	257	vb(ji,jj,jk) = zvf
	258	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	259	vn(ji,jj,jk) = va(ji,jj,jk)
	260	END DO
	261	END DO
	262	END DO
	263	!
	264	ELSE
	265	! Temporary filered scale factor at (u/v)-points (will become before scale factor)
	266	!-----------------------------------------------
	267	! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
	268	DO jk = 1, jpkm1
	269	DO jj = 1, jpjm1
	270	DO ji = 1, jpim1
	271	zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk)
	272	zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk)
	273	zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk)
	274	ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
	275	ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
	276	END DO
	277	END DO
	278	END DO
	279	! lateral boundary conditions
	280	CALL lbc_lnk( ze3u_f, 'U', 1. )
	281	CALL lbc_lnk( ze3v_f, 'V', 1. )
	282	! Add initial scale factor to scale factor anomaly
	283	ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
	284	ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
	285	! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
	286	! ----------------------------------- ===========================
	287	DO jk = 1, jpkm1
	288	DO jj = 1, jpj
	289	DO ji = 1, jpim1
	290	zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
	291	zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
	292	zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
	293	zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
	294	zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
	295	zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
	296	!
	297	zuf = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)
	298	zvf = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)
	299	!
	300	ub(ji,jj,jk) = zuf ! ub <-- filtered velocity
	301	vb(ji,jj,jk) = zvf
	302	un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua
	303	vn(ji,jj,jk) = va(ji,jj,jk)
	304	END DO
	305	END DO
	306	END DO
	307	fse3u_b(:,:,:) = ze3u_f(:,:,:) ! e3u_b <-- filtered scale factor
	308	fse3v_b(:,:,:) = ze3v_f(:,:,:)
	309	CALL lbc_lnk( ub, 'U', -1. ) ! lateral boundary conditions
	310	CALL lbc_lnk( vb, 'V', -1. )
	311	ENDIF
	312	!
[3]	313	ENDIF
[2528]	314	!
[258]	315	ENDIF
[3]	316
[1438]	317	IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
	318	& tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask )
	319	!
[3]	320	END SUBROUTINE dyn_nxt
	321
[1502]	322	!!=========================================================================
[3]	323	END MODULE dynnxt

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

Download in other formats: