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

Context Navigation

source: trunk/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 973

Last change on this file since 973 was 911, checked in by ctlod, 16 years ago
Implementation of the BDY package, see ticket: #126
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.4 KB

Rev	Line
[3]	1	MODULE dynnxt
	2	!!======================================================================
	3	!! * MODULE dynnxt *
	4	!! Ocean dynamics: time stepping
	5	!!======================================================================
	6
	7	!!----------------------------------------------------------------------
	8	!! dyn_nxt : update the horizontal velocity from the momentum trend
	9	!!----------------------------------------------------------------------
	10	!! * Modules used
	11	USE oce ! ocean dynamics and tracers
	12	USE dom_oce ! ocean space and time domain
	13	USE in_out_manager ! I/O manager
[367]	14	USE obc_oce ! ocean open boundary conditions
[3]	15	USE obcdyn ! open boundary condition for momentum (obc_dyn routine)
[367]	16	USE obcdyn_bt ! 2D open boundary condition for momentum (obc_dyn_bt routine)
	17	USE obcvol ! ocean open boundary condition (obc_vol routines)
[911]	18	USE bdy_oce ! unstructured open boundary conditions
	19	USE bdydta ! unstructured open boundary conditions
	20	USE bdydyn ! unstructured open boundary conditions
[367]	21	USE dynspg_oce ! type of surface pressure gradient
[3]	22	USE lbclnk ! lateral boundary condition (or mpp link)
[258]	23	USE prtctl ! Print control
[389]	24	USE agrif_opa_update
	25	USE agrif_opa_interp
[592]	26	USE domvvl ! variable volume
[3]	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
	31	!! * Accessibility
	32	PUBLIC dyn_nxt ! routine called by step.F90
[592]	33	!! * Substitutions
	34	# include "domzgr_substitute.h90"
[3]	35	!!----------------------------------------------------------------------
	36
	37	CONTAINS
	38
	39	SUBROUTINE dyn_nxt ( kt )
	40	!!----------------------------------------------------------------------
	41	!! * ROUTINE dyn_nxt *
	42	!!
	43	!! ** Purpose : Compute the after horizontal velocity from the
	44	!! momentum trend.
	45	!!
	46	!! ** Method : Apply lateral boundary conditions on the trends (ua,va)
	47	!! through calls to routine lbc_lnk.
	48	!! After velocity is compute using a leap-frog scheme environment:
	49	!! (ua,va) = (ub,vb) + 2 rdt (ua,va)
[359]	50	!! Note that if lk_dynspg_flt=T, the time stepping has already been
[32]	51	!! performed in dynspg module
[3]	52	!! Time filter applied on now horizontal velocity to avoid the
	53	!! divergence of two consecutive time-steps and swap of dynamics
	54	!! arrays to start the next time step:
	55	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
	56	!! (un,vn) = (ua,va)
	57	!!
	58	!! ** Action : - Update ub,vb arrays, the before horizontal velocity
	59	!! - Update un,vn arrays, the now horizontal velocity
	60	!!
	61	!! History :
	62	!! ! 87-02 (P. Andrich, D. L Hostis) Original code
	63	!! ! 90-10 (C. Levy, G. Madec)
	64	!! ! 93-03 (M. Guyon) symetrical conditions
	65	!! ! 97-02 (G. Madec & M. Imbard) opa, release 8.0
	66	!! ! 97-04 (A. Weaver) Euler forward step
	67	!! ! 97-06 (G. Madec) lateral boudary cond., lbc routine
	68	!! 8.5 ! 02-08 (G. Madec) F90: Free form and module
	69	!! ! 02-10 (C. Talandier, A-M. Treguier) Open boundary cond.
[359]	70	!! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization
[911]	71	!! " ! 07-07 (D. Storkey) Calls to BDY routines.
[3]	72	!!----------------------------------------------------------------------
	73	!! * Arguments
	74	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	75
	76	!! * Local declarations
	77	INTEGER :: ji, jj, jk ! dummy loop indices
	78	REAL(wp) :: z2dt ! temporary scalar
[642]	79	REAL(wp) :: zsshun1, zsshvn1 ! temporary scalar
[592]	80	!! Variable volume
[642]	81	REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace
	82	zsshub, zsshun, zsshua, &
	83	zsshvb, zsshvn, zsshva
[592]	84	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
	85	zfse3ub, zfse3un, zfse3ua, & ! 3D workspace
	86	zfse3vb, zfse3vn, zfse3va
[3]	87	!!----------------------------------------------------------------------
[247]	88	!! OPA 9.0 , LOCEAN-IPSL (2005)
[782]	89	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/DYN/dynnxt.F90,v 1.13 2007/05/25 15:51:50 opalod Exp $
[247]	90	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
[3]	91	!!----------------------------------------------------------------------
	92
	93	IF( kt == nit000 ) THEN
	94	IF(lwp) WRITE(numout,*)
	95	IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
	96	IF(lwp) WRITE(numout,*) '~~~~~~~'
	97	ENDIF
	98
	99	! Local constant initialization
	100	z2dt = 2. * rdt
	101	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
	102
[592]	103	!! Explicit physics with thickness weighted updates
	104	IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN
	105
	106	! Sea surface elevation time stepping
	107	! -----------------------------------
	108	!
	109	DO jj = 1, jpjm1
	110	DO ji = 1,jpim1
	111
	112	! Sea Surface Height at u-point before
	113	zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	114	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) &
	115	& + e1t(ji+1,jj ) * e2t(ji+1,jj ) * sshbb(ji+1,jj ) )
	116
	117	! Sea Surface Height at v-point before
	118	zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	119	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) &
	120	& + e1t(ji ,jj+1) * e2t(ji ,jj+1) * sshbb(ji ,jj+1) )
	121
	122	! Sea Surface Height at u-point after
	123	zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	124	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) &
	125	& + e1t(ji+1,jj ) * e2t(ji+1,jj ) * ssha(ji+1,jj ) )
	126
	127	! Sea Surface Height at v-point after
	128	zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	129	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) &
	130	& + e1t(ji ,jj+1) * e2t(ji ,jj+1) * ssha(ji ,jj+1) )
	131
	132	END DO
	133	END DO
	134
	135	! Boundaries conditions
	136	CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. )
	137	CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. )
	138
	139	! Scale factors at before and after time step
	140	! -------------------------------------------
[661]	141	CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ; CALL dom_vvl_sf( zsshua, 'U', zfse3ua )
	142	CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ; CALL dom_vvl_sf( zsshva, 'V', zfse3va )
[592]	143
	144	! Asselin filtered scale factor at now time step
	145	! ----------------------------------------------
	146	IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
[661]	147	CALL dom_vvl_sf_ini( 'U', zfse3un ) ; CALL dom_vvl_sf_ini( 'V', zfse3vn )
[592]	148	ELSE
[642]	149	zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:)
	150	zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:)
[661]	151	CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ; CALL dom_vvl_sf( zsshvn, 'V', zfse3vn )
[592]	152	ENDIF
	153
	154	! Thickness weighting
	155	! -------------------
[642]	156	DO jk = 1, jpkm1
	157	DO jj = 1, jpj
	158	DO ji = 1, jpi
	159	ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)
	160	va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)
[592]	161
[642]	162	un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk)
	163	vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk)
[592]	164
[642]	165	ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)
	166	vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)
	167	END DO
	168	END DO
	169	END DO
[592]	170
	171	ENDIF
	172
[3]	173	! Lateral boundary conditions on ( ua, va )
	174	CALL lbc_lnk( ua, 'U', -1. )
	175	CALL lbc_lnk( va, 'V', -1. )
	176
	177	! ! ===============
	178	DO jk = 1, jpkm1 ! Horizontal slab
	179	! ! ===============
	180	! Next velocity
	181	! -------------
[359]	182	#if defined key_dynspg_flt
[3]	183	! Leap-frog time stepping already done in dynspg.F routine
	184	#else
	185	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	186	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
	187	! Leap-frog time stepping
	188	ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
	189	va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
	190	END DO
	191	END DO
[911]	192
	193	IF( lk_vvl ) THEN
	194	! Unweight velocities prior to updating open boundaries.
	195
	196	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	197	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
	198	ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk)
	199	va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk)
	200
	201	un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk)
	202	vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk)
	203
	204	ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk)
	205	vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk)
	206	END DO
	207	END DO
	208
	209	ENDIF
	210
[3]	211	# if defined key_obc
	212	! ! ===============
	213	END DO ! End of slab
	214	! ! ===============
	215	! Update (ua,va) along open boundaries (only in the rigid-lid case)
	216	CALL obc_dyn( kt )
[367]	217
	218	IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
	219	!Flather boundary condition :
	220	! - Update sea surface height on each open boundary
	221	! sshn (= after ssh) for explicit case
	222	! sshn_b (= after ssha_b) for time-splitting case
	223	! - Correct the barotropic velocities
	224	CALL obc_dyn_bt( kt )
	225
	226	!Boundary conditions on sshn ( after ssh)
	227	CALL lbc_lnk( sshn, 'T', 1. )
	228
	229	IF(ln_ctl) THEN ! print sum trends (used for debugging)
	230	CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask)
	231	ENDIF
	232
	233	IF ( ln_vol_cst ) CALL obc_vol( kt )
	234
	235	ENDIF
	236
[3]	237	! ! ===============
	238	DO jk = 1, jpkm1 ! Horizontal slab
	239	! ! ===============
[911]	240	# elif defined key_bdy \|\| key_bdy_tides
	241	! ! ===============
	242	END DO ! End of slab
	243	! ! ===============
	244	! Update (ua,va) along open boundaries (for exp or ts options).
	245	IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN
	246
	247	CALL bdy_dyn_frs( kt )
	248
	249	IF ( ln_bdy_fla ) THEN
	250
	251	ua_e(:,:)=0.0
	252	va_e(:,:)=0.0
	253
	254	! Set these variables for use in bdy_dyn_fla
	255	hu_e(:,:) = hu(:,:)
	256	hv_e(:,:) = hv(:,:)
	257
	258	DO jk = 1, jpkm1
	259	!! Vertically integrated momentum trends
	260	ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
	261	va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
	262	END DO
	263
	264	DO jk = 1 , jpkm1
	265	ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)
	266	va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)
	267	END DO
	268
	269	CALL bdy_dta_bt( kt+1, 0)
	270	CALL bdy_dyn_fla
	271
	272	ENDIF
	273
	274	DO jk = 1 , jpkm1
	275	ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)
	276	va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)
	277	END DO
	278
	279	ENDIF
	280
	281	! ! ===============
	282	DO jk = 1, jpkm1 ! Horizontal slab
	283	! ! ===============
[3]	284	# endif
[392]	285	# if defined key_agrif
[389]	286	! ! ===============
	287	END DO ! End of slab
	288	! ! ===============
	289	! Update (ua,va) along open boundaries (only in the rigid-lid case)
	290	CALL Agrif_dyn( kt )
	291	! ! ===============
	292	DO jk = 1, jpkm1 ! Horizontal slab
	293	! ! ===============
	294	# endif
[3]	295	#endif
[592]	296
[3]	297	! Time filter and swap of dynamics arrays
	298	! ------------------------------------------
	299	IF( neuler == 0 .AND. kt == nit000 ) THEN
[592]	300	IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels
	301	! caution: don't use (:,:) for this loop
	302	! it causes optimization problems on NEC in auto-tasking
	303	DO jj = 1, jpj
	304	DO ji = 1, jpi
[642]	305	zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk)
	306	zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk)
	307	ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
	308	vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
	309	zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
	310	zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
	311	un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
	312	vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
[592]	313	END DO
[3]	314	END DO
[592]	315	ELSE ! Fixed levels
	316	DO jj = 1, jpj
	317	DO ji = 1, jpi
	318	! Euler (forward) time stepping
	319	ub(ji,jj,jk) = un(ji,jj,jk)
	320	vb(ji,jj,jk) = vn(ji,jj,jk)
	321	un(ji,jj,jk) = ua(ji,jj,jk)
	322	vn(ji,jj,jk) = va(ji,jj,jk)
	323	END DO
	324	END DO
	325	ENDIF
[3]	326	ELSE
[592]	327	IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels
	328	! caution: don't use (:,:) for this loop
	329	! it causes optimization problems on NEC in auto-tasking
	330	DO jj = 1, jpj
	331	DO ji = 1, jpi
[642]	332	zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk)
	333	zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk)
[592]	334	ub(ji,jj,jk) = ( atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) &
[642]	335	& + atfp1 * un(ji,jj,jk) ) * zsshun1
[592]	336	vb(ji,jj,jk) = ( atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) &
[642]	337	& + atfp1 * vn(ji,jj,jk) ) * zsshvn1
	338	zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
	339	zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
	340	un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1
	341	vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1
[592]	342	END DO
[3]	343	END DO
[592]	344	ELSE ! Fixed levels
	345	DO jj = 1, jpj
	346	DO ji = 1, jpi
	347	! Leap-frog time stepping
	348	ub(ji,jj,jk) = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
	349	vb(ji,jj,jk) = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
	350	un(ji,jj,jk) = ua(ji,jj,jk)
	351	vn(ji,jj,jk) = va(ji,jj,jk)
	352	END DO
	353	END DO
	354	ENDIF
[3]	355	ENDIF
	356	! ! ===============
	357	END DO ! End of slab
	358	! ! ===============
	359
[258]	360	IF(ln_ctl) THEN
	361	CALL prt_ctl(tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
	362	& tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)
	363	ENDIF
[3]	364
[392]	365	#if defined key_agrif
[389]	366	IF (.NOT.Agrif_Root()) CALL Agrif_Update_Dyn( kt )
	367	#endif
	368
[3]	369	END SUBROUTINE dyn_nxt
	370
	371	!!======================================================================
	372	END MODULE dynnxt

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

Download in other formats: