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.

dynzad.F90 in branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90 @ 5189

Last change on this file since 5189 was 5189, checked in by mathiot, 9 years ago
ISF cleaning branch: simplification and bug correction in hpg_isf, zps_hde_isf, mixed layer definition, slope, diffusion, vertical advection and top friction
Property svn:keywords set to `Id`
File size: 12.7 KB

Rev	Line
[3]	1	MODULE dynzad
	2	!!======================================================================
	3	!! * MODULE dynzad *
	4	!! Ocean dynamics : vertical advection trend
	5	!!======================================================================
[2715]	6	!! History : OPA ! 1991-01 (G. Madec) Original code
	7	!! 7.0 ! 1991-11 (G. Madec)
	8	!! 7.5 ! 1996-01 (G. Madec) statement function for e3
	9	!! NEMO 0.5 ! 2002-07 (G. Madec) Free form, F90
[503]	10	!!----------------------------------------------------------------------
[3]	11
	12	!!----------------------------------------------------------------------
[503]	13	!! dyn_zad : vertical advection momentum trend
[3]	14	!!----------------------------------------------------------------------
[503]	15	USE oce ! ocean dynamics and tracers
	16	USE dom_oce ! ocean space and time domain
[888]	17	USE sbc_oce ! surface boundary condition: ocean
[4990]	18	USE trd_oce ! trends: ocean variables
	19	USE trddyn ! trend manager: dynamics
	20	!
[719]	21	USE in_out_manager ! I/O manager
[4990]	22	USE lib_mpp ! MPP library
[503]	23	USE prtctl ! Print control
[4990]	24	USE wrk_nemo ! Memory Allocation
	25	USE timing ! Timing
[3]	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
[4990]	30	PUBLIC dyn_zad ! routine called by dynadv.F90
	31	PUBLIC dyn_zad_zts ! routine called by dynadv.F90
[3]	32
	33	!! * Substitutions
	34	# include "domzgr_substitute.h90"
	35	# include "vectopt_loop_substitute.h90"
	36	!!----------------------------------------------------------------------
[2528]	37	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	38	!! $Id$
[2715]	39	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	40	!!----------------------------------------------------------------------
	41	CONTAINS
	42
	43	SUBROUTINE dyn_zad ( kt )
	44	!!----------------------------------------------------------------------
	45	!! * ROUTINE dynzad *
	46	!!
	47	!! ** Purpose : Compute the now vertical momentum advection trend and
	48	!! add it to the general trend of momentum equation.
	49	!!
	50	!! ** Method : The now vertical advection of momentum is given by:
	51	!! w dz(u) = ua + 1/(e1ue2ue3u) mk+1[ mi(e1te2twn) dk(un) ]
	52	!! w dz(v) = va + 1/(e1ve2ve3v) mk+1[ mj(e1te2twn) dk(vn) ]
	53	!! Add this trend to the general trend (ua,va):
	54	!! (ua,va) = (ua,va) + w dz(u,v)
	55	!!
	56	!! ** Action : - Update (ua,va) with the vert. momentum adv. trends
[4990]	57	!! - Send the trends to trddyn for diagnostics (l_trddyn=T)
[3294]	58	!!----------------------------------------------------------------------
[503]	59	INTEGER, INTENT(in) :: kt ! ocean time-step inedx
[2715]	60	!
[503]	61	INTEGER :: ji, jj, jk ! dummy loop indices
	62	REAL(wp) :: zua, zva ! temporary scalars
[3294]	63	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw
	64	REAL(wp), POINTER, DIMENSION(:,: ) :: zww
	65	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
[3]	66	!!----------------------------------------------------------------------
[3294]	67	!
	68	IF( nn_timing == 1 ) CALL timing_start('dyn_zad')
	69	!
	70	CALL wrk_alloc( jpi,jpj, zww )
	71	CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw )
	72	!
[3]	73	IF( kt == nit000 ) THEN
	74	IF(lwp)WRITE(numout,*)
	75	IF(lwp)WRITE(numout,*) 'dyn_zad : arakawa advection scheme'
	76	ENDIF
[216]	77
[503]	78	IF( l_trddyn ) THEN ! Save ua and va trends
[3294]	79	CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
[503]	80	ztrdu(:,:,:) = ua(:,:,:)
	81	ztrdv(:,:,:) = va(:,:,:)
[216]	82	ENDIF
[3]	83
[5189]	84	! ! Vertical momentum advection at uw and vw-points
	85	DO jk = 2, jpk
[503]	86	DO jj = 2, jpj ! vertical fluxes
[5189]	87	DO ji = fs_2, jpi
	88	zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
[3]	89	END DO
	90	END DO
[5189]	91	DO jj = 2, jpjm1 ! interior vertical momentum advection at w-point
	92	DO ji = fs_2, fs_jpim1
	93	zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )
	94	zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )
[3]	95	END DO
	96	END DO
	97	END DO
[5189]	98	DO jj = 2, jpjm1 ! First/last level advective fluxes
	99	DO ji = fs_2, fs_jpim1 ! w is zero at k=1 and k=jpk
	100	zwuw(ji,jj, 1 ) = 0._wp
	101	zwvw(ji,jj, 1 ) = 0._wp
	102	zwuw(ji,jj,jpk) = 0._wp
	103	zwvw(ji,jj,jpk) = 0._wp
[5120]	104	END DO
[5189]	105	END DO
[3]	106
[503]	107	DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points
[3]	108	DO jj = 2, jpjm1
[503]	109	DO ji = fs_2, fs_jpim1 ! vector opt.
	110	! ! vertical momentum advective trends
[3]	111	zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
	112	zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
[503]	113	! ! add the trends to the general momentum trends
[3]	114	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
	115	va(ji,jj,jk) = va(ji,jj,jk) + zva
	116	END DO
	117	END DO
	118	END DO
	119
[503]	120	IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic
	121	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	122	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[4990]	123	CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
[3294]	124	CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )
[216]	125	ENDIF
[503]	126	! ! Control print
	127	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, &
	128	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
	129	!
[3294]	130	CALL wrk_dealloc( jpi,jpj, zww )
	131	CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw )
[2715]	132	!
[3294]	133	IF( nn_timing == 1 ) CALL timing_stop('dyn_zad')
	134	!
[3]	135	END SUBROUTINE dyn_zad
	136
[4990]	137	SUBROUTINE dyn_zad_zts ( kt )
	138	!!----------------------------------------------------------------------
	139	!! * ROUTINE dynzad_zts *
	140	!!
	141	!! ** Purpose : Compute the now vertical momentum advection trend and
	142	!! add it to the general trend of momentum equation. This version
	143	!! uses sub-timesteps for improved numerical stability with small
	144	!! vertical grid sizes. This is especially relevant when using
	145	!! embedded ice with thin surface boxes.
	146	!!
	147	!! ** Method : The now vertical advection of momentum is given by:
	148	!! w dz(u) = ua + 1/(e1ue2ue3u) mk+1[ mi(e1te2twn) dk(un) ]
	149	!! w dz(v) = va + 1/(e1ve2ve3v) mk+1[ mj(e1te2twn) dk(vn) ]
	150	!! Add this trend to the general trend (ua,va):
	151	!! (ua,va) = (ua,va) + w dz(u,v)
	152	!!
	153	!! ** Action : - Update (ua,va) with the vert. momentum adv. trends
	154	!! - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
	155	!!----------------------------------------------------------------------
	156	INTEGER, INTENT(in) :: kt ! ocean time-step inedx
	157	!
	158	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	159	INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection
	160	INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps
	161	REAL(wp) :: zua, zva ! temporary scalars
	162	REAL(wp) :: zr_rdt ! temporary scalar
	163	REAL(wp) :: z2dtzts ! length of Euler forward sub-timestep for vertical advection
	164	REAL(wp) :: zts ! length of sub-timestep for vertical advection
	165	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw, zww
	166	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
	167	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zus , zvs
	168	!!----------------------------------------------------------------------
	169	!
	170	IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts')
	171	!
	172	CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww )
	173	CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs )
	174	!
	175	IF( kt == nit000 ) THEN
	176	IF(lwp)WRITE(numout,*)
	177	IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
	178	ENDIF
	179
	180	IF( l_trddyn ) THEN ! Save ua and va trends
	181	CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
	182	ztrdu(:,:,:) = ua(:,:,:)
	183	ztrdv(:,:,:) = va(:,:,:)
	184	ENDIF
	185
	186	IF( neuler == 0 .AND. kt == nit000 ) THEN
	187	z2dtzts = rdt / REAL( jnzts, wp ) ! = rdt (restart with Euler time stepping)
	188	ELSE
	189	z2dtzts = 2._wp * rdt / REAL( jnzts, wp ) ! = 2 rdt (leapfrog)
	190	ENDIF
	191
	192	DO jk = 2, jpkm1 ! Calculate and store vertical fluxes
	193	DO jj = 2, jpj
	194	DO ji = fs_2, jpi ! vector opt.
[5189]	195	zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
[4990]	196	END DO
	197	END DO
	198	END DO
[5120]	199	!
	200	! Surface and bottom advective fluxes set to zero
	201	DO jj = 2, jpjm1
[4990]	202	DO ji = fs_2, fs_jpim1 ! vector opt.
[5120]	203	zwuw(ji,jj, 1 ) = 0._wp
	204	zwvw(ji,jj, 1 ) = 0._wp
[4990]	205	zwuw(ji,jj,jpk) = 0._wp
	206	zwvw(ji,jj,jpk) = 0._wp
	207	END DO
	208	END DO
	209
	210	! Start with before values and use sub timestepping to reach after values
	211
	212	zus(:,:,:,1) = ub(:,:,:)
	213	zvs(:,:,:,1) = vb(:,:,:)
	214
	215	DO jl = 1, jnzts ! Start of sub timestepping loop
	216
	217	IF( jl == 1 ) THEN ! Euler forward to kick things off
	218	jtb = 1 ; jtn = 1 ; jta = 2
	219	zts = z2dtzts
	220	ELSEIF( jl == 2 ) THEN ! First leapfrog step
	221	jtb = 1 ; jtn = 2 ; jta = 3
	222	zts = 2._wp * z2dtzts
	223	ELSE ! Shuffle pointers for subsequent leapfrog steps
	224	jtb = MOD(jtb,3) + 1
	225	jtn = MOD(jtn,3) + 1
	226	jta = MOD(jta,3) + 1
	227	ENDIF
	228
	229	DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical
	230	DO jj = 2, jpjm1 ! vertical momentum advection at w-point
	231	DO ji = fs_2, fs_jpim1 ! vector opt.
[5120]	232	zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk)
	233	zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk)
[4990]	234	END DO
	235	END DO
	236	END DO
	237	DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points
	238	DO jj = 2, jpjm1
	239	DO ji = fs_2, fs_jpim1 ! vector opt.
	240	! ! vertical momentum advective trends
	241	zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
	242	zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
	243	zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
	244	zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
	245	END DO
	246	END DO
	247	END DO
	248
	249	END DO ! End of sub timestepping loop
	250
	251	zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
	252	DO jk = 1, jpkm1 ! Recover trends over the outer timestep
	253	DO jj = 2, jpjm1
	254	DO ji = fs_2, fs_jpim1 ! vector opt.
	255	! ! vertical momentum advective trends
	256	! ! add the trends to the general momentum trends
	257	ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
	258	va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
	259	END DO
	260	END DO
	261	END DO
	262
	263	IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic
	264	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	265	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
	266	CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
	267	CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )
	268	ENDIF
	269	! ! Control print
	270	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, &
	271	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
	272	!
	273	CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww )
	274	CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs )
	275	!
	276	IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts')
	277	!
	278	END SUBROUTINE dyn_zad_zts
	279
[503]	280	!!======================================================================
[3]	281	END MODULE dynzad

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

Download in other formats: