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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2364

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

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

Download in other formats: