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 @ 684

Last change on this file since 684 was 661, checked in by opalod, 17 years ago
nemo_v2_bugfix_039:RB: change dom_vvl functions to corresponding subroutines for compatibility with AGRIF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 13.9 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)
	18	USE dynspg_oce ! type of surface pressure gradient
[3]	19	USE lbclnk ! lateral boundary condition (or mpp link)
[258]	20	USE prtctl ! Print control
[389]	21	USE agrif_opa_update
	22	USE agrif_opa_interp
[592]	23	USE domvvl ! variable volume
[3]	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
	28	!! * Accessibility
	29	PUBLIC dyn_nxt ! routine called by step.F90
[592]	30	!! * Substitutions
	31	# include "domzgr_substitute.h90"
[3]	32	!!----------------------------------------------------------------------
	33
	34	CONTAINS
	35
	36	SUBROUTINE dyn_nxt ( kt )
	37	!!----------------------------------------------------------------------
	38	!! * ROUTINE dyn_nxt *
	39	!!
	40	!! ** Purpose : Compute the after horizontal velocity from the
	41	!! momentum trend.
	42	!!
	43	!! ** Method : Apply lateral boundary conditions on the trends (ua,va)
	44	!! through calls to routine lbc_lnk.
	45	!! After velocity is compute using a leap-frog scheme environment:
	46	!! (ua,va) = (ub,vb) + 2 rdt (ua,va)
[359]	47	!! Note that if lk_dynspg_flt=T, the time stepping has already been
[32]	48	!! performed in dynspg module
[3]	49	!! Time filter applied on now horizontal velocity to avoid the
	50	!! divergence of two consecutive time-steps and swap of dynamics
	51	!! arrays to start the next time step:
	52	!! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
	53	!! (un,vn) = (ua,va)
	54	!!
	55	!! ** Action : - Update ub,vb arrays, the before horizontal velocity
	56	!! - Update un,vn arrays, the now horizontal velocity
	57	!!
	58	!! History :
	59	!! ! 87-02 (P. Andrich, D. L Hostis) Original code
	60	!! ! 90-10 (C. Levy, G. Madec)
	61	!! ! 93-03 (M. Guyon) symetrical conditions
	62	!! ! 97-02 (G. Madec & M. Imbard) opa, release 8.0
	63	!! ! 97-04 (A. Weaver) Euler forward step
	64	!! ! 97-06 (G. Madec) lateral boudary cond., lbc routine
	65	!! 8.5 ! 02-08 (G. Madec) F90: Free form and module
	66	!! ! 02-10 (C. Talandier, A-M. Treguier) Open boundary cond.
[359]	67	!! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization
[3]	68	!!----------------------------------------------------------------------
	69	!! * Arguments
	70	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	71
	72	!! * Local declarations
	73	INTEGER :: ji, jj, jk ! dummy loop indices
	74	REAL(wp) :: z2dt ! temporary scalar
[642]	75	REAL(wp) :: zsshun1, zsshvn1 ! temporary scalar
[592]	76	!! Variable volume
[642]	77	REAL(wp), DIMENSION(jpi,jpj) :: & ! 2D workspace
	78	zsshub, zsshun, zsshua, &
	79	zsshvb, zsshvn, zsshva
[592]	80	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
	81	zfse3ub, zfse3un, zfse3ua, & ! 3D workspace
	82	zfse3vb, zfse3vn, zfse3va
[3]	83	!!----------------------------------------------------------------------
[247]	84	!! OPA 9.0 , LOCEAN-IPSL (2005)
	85	!! $Header$
	86	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
[3]	87	!!----------------------------------------------------------------------
	88
	89	IF( kt == nit000 ) THEN
	90	IF(lwp) WRITE(numout,*)
	91	IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
	92	IF(lwp) WRITE(numout,*) '~~~~~~~'
	93	ENDIF
	94
	95	! Local constant initialization
	96	z2dt = 2. * rdt
	97	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
	98
[592]	99	!! Explicit physics with thickness weighted updates
	100	IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN
	101
	102	! Sea surface elevation time stepping
	103	! -----------------------------------
	104	!
	105	DO jj = 1, jpjm1
	106	DO ji = 1,jpim1
	107
	108	! Sea Surface Height at u-point before
	109	zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	110	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) &
	111	& + e1t(ji+1,jj ) * e2t(ji+1,jj ) * sshbb(ji+1,jj ) )
	112
	113	! Sea Surface Height at v-point before
	114	zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	115	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) &
	116	& + e1t(ji ,jj+1) * e2t(ji ,jj+1) * sshbb(ji ,jj+1) )
	117
	118	! Sea Surface Height at u-point after
	119	zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
	120	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) &
	121	& + e1t(ji+1,jj ) * e2t(ji+1,jj ) * ssha(ji+1,jj ) )
	122
	123	! Sea Surface Height at v-point after
	124	zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
	125	& * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) &
	126	& + e1t(ji ,jj+1) * e2t(ji ,jj+1) * ssha(ji ,jj+1) )
	127
	128	END DO
	129	END DO
	130
	131	! Boundaries conditions
	132	CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. )
	133	CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. )
	134
	135	! Scale factors at before and after time step
	136	! -------------------------------------------
[661]	137	CALL dom_vvl_sf( zsshub, 'U', zfse3ub ) ; CALL dom_vvl_sf( zsshua, 'U', zfse3ua )
	138	CALL dom_vvl_sf( zsshvb, 'V', zfse3vb ) ; CALL dom_vvl_sf( zsshva, 'V', zfse3va )
[592]	139
	140	! Asselin filtered scale factor at now time step
	141	! ----------------------------------------------
	142	IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
[661]	143	CALL dom_vvl_sf_ini( 'U', zfse3un ) ; CALL dom_vvl_sf_ini( 'V', zfse3vn )
[592]	144	ELSE
[642]	145	zsshun(:,:) = atfp * ( zsshub(:,:) + zsshua(:,:) ) + atfp1 * sshu(:,:)
	146	zsshvn(:,:) = atfp * ( zsshvb(:,:) + zsshva(:,:) ) + atfp1 * sshv(:,:)
[661]	147	CALL dom_vvl_sf( zsshun, 'U', zfse3un ) ; CALL dom_vvl_sf( zsshvn, 'V', zfse3vn )
[592]	148	ENDIF
	149
	150	! Thickness weighting
	151	! -------------------
[642]	152	DO jk = 1, jpkm1
	153	DO jj = 1, jpj
	154	DO ji = 1, jpi
	155	ua(ji,jj,jk) = ua(ji,jj,jk) * fse3u(ji,jj,jk)
	156	va(ji,jj,jk) = va(ji,jj,jk) * fse3v(ji,jj,jk)
[592]	157
[642]	158	un(ji,jj,jk) = un(ji,jj,jk) * fse3u(ji,jj,jk)
	159	vn(ji,jj,jk) = vn(ji,jj,jk) * fse3v(ji,jj,jk)
[592]	160
[642]	161	ub(ji,jj,jk) = ub(ji,jj,jk) * zfse3ub(ji,jj,jk)
	162	vb(ji,jj,jk) = vb(ji,jj,jk) * zfse3vb(ji,jj,jk)
	163	END DO
	164	END DO
	165	END DO
[592]	166
	167	ENDIF
	168
[3]	169	! Lateral boundary conditions on ( ua, va )
	170	CALL lbc_lnk( ua, 'U', -1. )
	171	CALL lbc_lnk( va, 'V', -1. )
	172
	173	! ! ===============
	174	DO jk = 1, jpkm1 ! Horizontal slab
	175	! ! ===============
	176	! Next velocity
	177	! -------------
[359]	178	#if defined key_dynspg_flt
[3]	179	! Leap-frog time stepping already done in dynspg.F routine
	180	#else
	181	DO jj = 1, jpj ! caution: don't use (:,:) for this loop
	182	DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking
	183	! Leap-frog time stepping
	184	ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk)
	185	va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk)
	186	END DO
	187	END DO
	188	# if defined key_obc
	189	! ! ===============
	190	END DO ! End of slab
	191	! ! ===============
	192	! Update (ua,va) along open boundaries (only in the rigid-lid case)
	193	CALL obc_dyn( kt )
[367]	194
	195	IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
	196	!Flather boundary condition :
	197	! - Update sea surface height on each open boundary
	198	! sshn (= after ssh) for explicit case
	199	! sshn_b (= after ssha_b) for time-splitting case
	200	! - Correct the barotropic velocities
	201	CALL obc_dyn_bt( kt )
	202
	203	!Boundary conditions on sshn ( after ssh)
	204	CALL lbc_lnk( sshn, 'T', 1. )
	205
	206	IF(ln_ctl) THEN ! print sum trends (used for debugging)
	207	CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask)
	208	ENDIF
	209
	210	IF ( ln_vol_cst ) CALL obc_vol( kt )
	211
	212	ENDIF
	213
[3]	214	! ! ===============
	215	DO jk = 1, jpkm1 ! Horizontal slab
	216	! ! ===============
	217	# endif
[392]	218	# if defined key_agrif
[389]	219	! ! ===============
	220	END DO ! End of slab
	221	! ! ===============
	222	! Update (ua,va) along open boundaries (only in the rigid-lid case)
	223	CALL Agrif_dyn( kt )
	224	! ! ===============
	225	DO jk = 1, jpkm1 ! Horizontal slab
	226	! ! ===============
	227	# endif
[3]	228	#endif
[592]	229
[3]	230	! Time filter and swap of dynamics arrays
	231	! ------------------------------------------
	232	IF( neuler == 0 .AND. kt == nit000 ) THEN
[592]	233	IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels
	234	! caution: don't use (:,:) for this loop
	235	! it causes optimization problems on NEC in auto-tasking
	236	DO jj = 1, jpj
	237	DO ji = 1, jpi
[642]	238	zsshun1 = umask(ji,jj,jk) / fse3u(ji,jj,jk)
	239	zsshvn1 = vmask(ji,jj,jk) / fse3v(ji,jj,jk)
	240	ub(ji,jj,jk) = un(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
	241	vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
	242	zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
	243	zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
	244	un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1 * umask(ji,jj,jk)
	245	vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1 * vmask(ji,jj,jk)
[592]	246	END DO
[3]	247	END DO
[592]	248	ELSE ! Fixed levels
	249	DO jj = 1, jpj
	250	DO ji = 1, jpi
	251	! Euler (forward) time stepping
	252	ub(ji,jj,jk) = un(ji,jj,jk)
	253	vb(ji,jj,jk) = vn(ji,jj,jk)
	254	un(ji,jj,jk) = ua(ji,jj,jk)
	255	vn(ji,jj,jk) = va(ji,jj,jk)
	256	END DO
	257	END DO
	258	ENDIF
[3]	259	ELSE
[592]	260	IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels
	261	! caution: don't use (:,:) for this loop
	262	! it causes optimization problems on NEC in auto-tasking
	263	DO jj = 1, jpj
	264	DO ji = 1, jpi
[642]	265	zsshun1 = umask(ji,jj,jk) / zfse3un(ji,jj,jk)
	266	zsshvn1 = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk)
[592]	267	ub(ji,jj,jk) = ( atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) &
[642]	268	& + atfp1 * un(ji,jj,jk) ) * zsshun1
[592]	269	vb(ji,jj,jk) = ( atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) &
[642]	270	& + atfp1 * vn(ji,jj,jk) ) * zsshvn1
	271	zsshun1 = umask(ji,jj,jk) / zfse3ua(ji,jj,jk)
	272	zsshvn1 = vmask(ji,jj,jk) / zfse3va(ji,jj,jk)
	273	un(ji,jj,jk) = ua(ji,jj,jk) * zsshun1
	274	vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn1
[592]	275	END DO
[3]	276	END DO
[592]	277	ELSE ! Fixed levels
	278	DO jj = 1, jpj
	279	DO ji = 1, jpi
	280	! Leap-frog time stepping
	281	ub(ji,jj,jk) = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)
	282	vb(ji,jj,jk) = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)
	283	un(ji,jj,jk) = ua(ji,jj,jk)
	284	vn(ji,jj,jk) = va(ji,jj,jk)
	285	END DO
	286	END DO
	287	ENDIF
[3]	288	ENDIF
	289	! ! ===============
	290	END DO ! End of slab
	291	! ! ===============
	292
[258]	293	IF(ln_ctl) THEN
	294	CALL prt_ctl(tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, &
	295	& tab3d_2=vn, clinfo2=' Vn: ', mask2=vmask)
	296	ENDIF
[3]	297
[392]	298	#if defined key_agrif
[389]	299	IF (.NOT.Agrif_Root()) CALL Agrif_Update_Dyn( kt )
	300	#endif
	301
[3]	302	END SUBROUTINE dyn_nxt
	303
	304	!!======================================================================
	305	END MODULE dynnxt

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

Download in other formats: