Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynzdf.F90 @ 11613

Last change on this file since 11613 was 11613, checked in by davestorkey, 5 years ago
UKMO/NEMO_4.0_momentum_trends branch : first set of code changes.
File size: 27.5 KB

Rev	Line
[456]	1	MODULE dynzdf
	2	!!==============================================================================
	3	!! * MODULE dynzdf *
	4	!! Ocean dynamics : vertical component of the momentum mixing trend
	5	!!==============================================================================
[2528]	6	!! History : 1.0 ! 2005-11 (G. Madec) Original code
	7	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
[9019]	8	!! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point
[456]	9	!!----------------------------------------------------------------------
[503]	10
	11	!!----------------------------------------------------------------------
[9019]	12	!! dyn_zdf : compute the after velocity through implicit calculation of vertical mixing
[456]	13	!!----------------------------------------------------------------------
[5836]	14	USE oce ! ocean dynamics and tracers variables
[9019]	15	USE phycst ! physical constants
[5836]	16	USE dom_oce ! ocean space and time domain variables
[9019]	17	USE sbc_oce ! surface boundary condition: ocean
[5836]	18	USE zdf_oce ! ocean vertical physics variables
[9019]	19	USE zdfdrg ! vertical physics: top/bottom drag coef.
	20	USE dynadv ,ONLY: ln_dynadv_vec ! dynamics: advection form
	21	USE dynldf_iso,ONLY: akzu, akzv ! dynamics: vertical component of rotated lateral mixing
[9490]	22	USE ldfdyn ! lateral diffusion: eddy viscosity coef. and type of operator
[5836]	23	USE trd_oce ! trends: ocean variables
	24	USE trddyn ! trend manager: dynamics
	25	!
	26	USE in_out_manager ! I/O manager
[11613]	27	USE iom ! IOM library
[5836]	28	USE lib_mpp ! MPP library
	29	USE prtctl ! Print control
	30	USE timing ! Timing
[456]	31
	32	IMPLICIT NONE
	33	PRIVATE
	34
[9019]	35	PUBLIC dyn_zdf ! routine called by step.F90
[456]	36
[9019]	37	REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
[456]	38
	39	!! * Substitutions
	40	# include "vectopt_loop_substitute.h90"
	41	!!----------------------------------------------------------------------
[9598]	42	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	43	!! $Id$
[10068]	44	!! Software governed by the CeCILL license (see ./LICENSE)
[456]	45	!!----------------------------------------------------------------------
	46	CONTAINS
	47
	48	SUBROUTINE dyn_zdf( kt )
	49	!!----------------------------------------------------------------------
	50	!! * ROUTINE dyn_zdf *
	51	!!
[9019]	52	!! ** Purpose : compute the trend due to the vert. momentum diffusion
	53	!! together with the Leap-Frog time stepping using an
	54	!! implicit scheme.
	55	!!
	56	!! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing
	57	!! ua = ub + 2dt ua vector form or linear free surf.
	58	!! ua = ( e3u_bub + 2dt * e3u_n*ua ) / e3u_a otherwise
	59	!! - update the after velocity with the implicit vertical mixing.
	60	!! This requires to solver the following system:
	61	!! ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
	62	!! with the following surface/top/bottom boundary condition:
	63	!! surface: wind stress input (averaged over kt-1/2 & kt+1/2)
	64	!! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
	65	!!
	66	!! ** Action : (ua,va) after velocity
[456]	67	!!---------------------------------------------------------------------
[9019]	68	INTEGER, INTENT(in) :: kt ! ocean time-step index
[3294]	69	!
[9019]	70	INTEGER :: ji, jj, jk ! dummy loop indices
	71	INTEGER :: iku, ikv ! local integers
[11613]	72	INTEGER :: ikbu, ikbv ! local integers
[9019]	73	REAL(wp) :: zzwi, ze3ua, zdt ! local scalars
	74	REAL(wp) :: zzws, ze3va ! - -
[10364]	75	REAL(wp) :: z1_e3ua, z1_e3va ! - -
	76	REAL(wp) :: zWu , zWv ! - -
	77	REAL(wp) :: zWui, zWvi ! - -
	78	REAL(wp) :: zWus, zWvs ! - -
[9019]	79	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace
	80	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -
[456]	81	!!---------------------------------------------------------------------
[3294]	82	!
[9019]	83	IF( ln_timing ) CALL timing_start('dyn_zdf')
[3294]	84	!
[9019]	85	IF( kt == nit000 ) THEN !* initialization
	86	IF(lwp) WRITE(numout,*)
	87	IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
	88	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
	89	!
	90	If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator
	91	ELSE ; r_vvl = 1._wp
	92	ENDIF
	93	ENDIF
	94	! !* set time step
[6140]	95	IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping)
	96	ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog)
[456]	97	ENDIF
[9250]	98	!
	99	! !* explicit top/bottom drag case
	100	IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, ub, vb, ua, va ) ! add top/bottom friction trend to (ua,va)
	101	!
	102	!
[9019]	103	! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in ua,va)
	104	!
	105	! ! time stepping except vertical diffusion
	106	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
	107	DO jk = 1, jpkm1
	108	ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
	109	va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)
	110	END DO
	111	ELSE ! applied on thickness weighted velocity
	112	DO jk = 1, jpkm1
	113	ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) &
	114	& + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)
	115	va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) &
	116	& + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)
	117	END DO
	118	ENDIF
	119	! ! add top/bottom friction
	120	! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
	121	! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
	122	! not lead to the effective stress seen over the whole barotropic loop.
	123	! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
	124	IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
	125	DO jk = 1, jpkm1 ! remove barotropic velocities
	126	ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
	127	va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
	128	END DO
[11613]	129	IF( l_trddyn ) THEN !* temporary save of ta and sa trends
	130	ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
	131	ztrdu(:,:,:) = ua(:,:,:)
	132	ztrdv(:,:,:) = va(:,:,:)
	133	ENDIF
[9019]	134	DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only
	135	DO ji = fs_2, fs_jpim1 ! vector opt.
	136	iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
	137	ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
	138	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
	139	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
	140	ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) ua_b(ji,jj) / ze3ua
	141	va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) va_b(ji,jj) / ze3va
	142	END DO
	143	END DO
[11613]	144	IF( l_trddyn ) THEN ! save explicit part of bottom friction trends
	145	ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / r2dt
	146	ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / r2dt
	147	CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfr, kt )
	148	ENDIF
[9019]	149	IF( ln_isfcav ) THEN ! Ocean cavities (ISF)
	150	DO jj = 2, jpjm1
	151	DO ji = fs_2, fs_jpim1 ! vector opt.
	152	iku = miku(ji,jj) ! top ocean level at u- and v-points
	153	ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
	154	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
	155	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
	156	ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) ua_b(ji,jj) / ze3ua
	157	va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) va_b(ji,jj) / ze3va
	158	END DO
	159	END DO
	160	END IF
	161	ENDIF
	162	!
[11613]	163	IF( l_trddyn ) THEN !* temporary save of ta and sa trends
	164	IF( .NOT. ALLOCATED(ztrdu) ) ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )
	165	ztrdu(:,:,:) = ua(:,:,:)
	166	ztrdv(:,:,:) = va(:,:,:)
	167	ENDIF
	168	!
[9019]	169	! !== Vertical diffusion on u ==!
	170	!
	171	! !* Matrix construction
	172	zdt = r2dt * 0.5
[10364]	173	IF( ln_zad_Aimp ) THEN !!
	174	SELECT CASE( nldf_dyn )
	175	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
	176	DO jk = 1, jpkm1
	177	DO jj = 2, jpjm1
	178	DO ji = fs_2, fs_jpim1 ! vector opt.
	179	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	180	zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
	181	& / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk )
	182	zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
	183	& / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	184	zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )
	185	zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
	186	zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
	187	zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
	188	zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
	189	END DO
[9019]	190	END DO
	191	END DO
[10364]	192	CASE DEFAULT ! iso-level lateral mixing
	193	DO jk = 1, jpkm1
	194	DO jj = 2, jpjm1
	195	DO ji = fs_2, fs_jpim1 ! vector opt.
	196	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	197	zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk )
	198	zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	199	zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )
	200	zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
	201	zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
	202	zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
	203	zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
	204	END DO
	205	END DO
	206	END DO
	207	END SELECT
	208	DO jj = 2, jpjm1 !* Surface boundary conditions
	209	DO ji = fs_2, fs_jpim1 ! vector opt.
	210	zwi(ji,jj,1) = 0._wp
	211	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
	212	zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2)
	213	zWus = 0.5_wp * ( wi(ji ,jj,2) + wi(ji+1,jj,2) )
	214	zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
	215	zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
	216	END DO
[9019]	217	END DO
[10364]	218	ELSE
	219	SELECT CASE( nldf_dyn )
	220	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
	221	DO jk = 1, jpkm1
	222	DO jj = 2, jpjm1
	223	DO ji = fs_2, fs_jpim1 ! vector opt.
	224	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	225	zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
	226	& / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk )
	227	zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
	228	& / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	229	zwi(ji,jj,jk) = zzwi
	230	zws(ji,jj,jk) = zzws
	231	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	232	END DO
[9019]	233	END DO
	234	END DO
[10364]	235	CASE DEFAULT ! iso-level lateral mixing
	236	DO jk = 1, jpkm1
	237	DO jj = 2, jpjm1
	238	DO ji = fs_2, fs_jpim1 ! vector opt.
	239	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	240	zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk )
	241	zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	242	zwi(ji,jj,jk) = zzwi
	243	zws(ji,jj,jk) = zzws
	244	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	245	END DO
	246	END DO
	247	END DO
	248	END SELECT
	249	DO jj = 2, jpjm1 !* Surface boundary conditions
	250	DO ji = fs_2, fs_jpim1 ! vector opt.
	251	zwi(ji,jj,1) = 0._wp
	252	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
	253	END DO
[9019]	254	END DO
[10364]	255	ENDIF
[9019]	256	!
	257	!
	258	! !== Apply semi-implicit bottom friction ==!
	259	!
	260	! Only needed for semi-implicit bottom friction setup. The explicit
	261	! bottom friction has been included in "u(v)a" which act as the R.H.S
	262	! column vector of the tri-diagonal matrix equation
	263	!
	264	IF ( ln_drgimp ) THEN ! implicit bottom friction
	265	DO jj = 2, jpjm1
	266	DO ji = 2, jpim1
	267	iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
	268	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point
	269	zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
	270	END DO
	271	END DO
	272	IF ( ln_isfcav ) THEN ! top friction (always implicit)
	273	DO jj = 2, jpjm1
	274	DO ji = 2, jpim1
	275	!!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed
	276	iku = miku(ji,jj) ! ocean top level at u- and v-points
	277	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point
	278	zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
	279	END DO
	280	END DO
	281	END IF
	282	ENDIF
	283	!
	284	! Matrix inversion starting from the first level
	285	!-----------------------------------------------------------------------
	286	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	287	!
	288	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	289	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	290	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	291	! ( ... )( ... ) ( ... )
	292	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	293	!
	294	! m is decomposed in the product of an upper and a lower triangular matrix
	295	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	296	! The solution (the after velocity) is in ua
	297	!-----------------------------------------------------------------------
	298	!
	299	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
	300	DO jj = 2, jpjm1
	301	DO ji = fs_2, fs_jpim1 ! vector opt.
	302	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	303	END DO
	304	END DO
	305	END DO
	306	!
	307	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
	308	DO ji = fs_2, fs_jpim1 ! vector opt.
	309	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
	310	ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
	311	& / ( ze3ua * rau0 ) * umask(ji,jj,1)
	312	END DO
	313	END DO
	314	DO jk = 2, jpkm1
	315	DO jj = 2, jpjm1
	316	DO ji = fs_2, fs_jpim1
	317	ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
	318	END DO
	319	END DO
	320	END DO
	321	!
	322	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
	323	DO ji = fs_2, fs_jpim1 ! vector opt.
	324	ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	325	END DO
	326	END DO
	327	DO jk = jpk-2, 1, -1
	328	DO jj = 2, jpjm1
	329	DO ji = fs_2, fs_jpim1
	330	ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	331	END DO
	332	END DO
	333	END DO
	334	!
	335	! !== Vertical diffusion on v ==!
	336	!
	337	! !* Matrix construction
	338	zdt = r2dt * 0.5
[10364]	339	IF( ln_zad_Aimp ) THEN !!
	340	SELECT CASE( nldf_dyn )
	341	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
	342	DO jk = 1, jpkm1
	343	DO jj = 2, jpjm1
	344	DO ji = fs_2, fs_jpim1 ! vector opt.
	345	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point
	346	zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
	347	& / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	348	zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
	349	& / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
	350	zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk )
	351	zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
	352	zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
	353	zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
	354	zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
	355	END DO
[9019]	356	END DO
	357	END DO
[10364]	358	CASE DEFAULT ! iso-level lateral mixing
	359	DO jk = 1, jpkm1
	360	DO jj = 2, jpjm1
	361	DO ji = fs_2, fs_jpim1 ! vector opt.
	362	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point
	363	zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	364	zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
	365	zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk )
	366	zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1)
	367	zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
	368	zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
	369	zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
	370	END DO
	371	END DO
	372	END DO
	373	END SELECT
	374	DO jj = 2, jpjm1 !* Surface boundary conditions
	375	DO ji = fs_2, fs_jpim1 ! vector opt.
	376	zwi(ji,jj,1) = 0._wp
	377	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
	378	zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2)
	379	zWvs = 0.5_wp * ( wi(ji,jj ,2) + wi(ji,jj+1,2) )
	380	zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
	381	zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
	382	END DO
[9019]	383	END DO
[10364]	384	ELSE
	385	SELECT CASE( nldf_dyn )
	386	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
	387	DO jk = 1, jpkm1
	388	DO jj = 2, jpjm1
	389	DO ji = fs_2, fs_jpim1 ! vector opt.
	390	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point
	391	zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
	392	& / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	393	zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
	394	& / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
	395	zwi(ji,jj,jk) = zzwi
	396	zws(ji,jj,jk) = zzws
	397	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	398	END DO
[9019]	399	END DO
	400	END DO
[10364]	401	CASE DEFAULT ! iso-level lateral mixing
	402	DO jk = 1, jpkm1
	403	DO jj = 2, jpjm1
	404	DO ji = fs_2, fs_jpim1 ! vector opt.
	405	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point
	406	zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	407	zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
	408	zwi(ji,jj,jk) = zzwi
	409	zws(ji,jj,jk) = zzws
	410	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	411	END DO
	412	END DO
	413	END DO
	414	END SELECT
	415	DO jj = 2, jpjm1 !* Surface boundary conditions
	416	DO ji = fs_2, fs_jpim1 ! vector opt.
	417	zwi(ji,jj,1) = 0._wp
	418	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
	419	END DO
[9019]	420	END DO
[10364]	421	ENDIF
[9019]	422	!
	423	! !== Apply semi-implicit top/bottom friction ==!
	424	!
	425	! Only needed for semi-implicit bottom friction setup. The explicit
	426	! bottom friction has been included in "u(v)a" which act as the R.H.S
	427	! column vector of the tri-diagonal matrix equation
	428	!
	429	IF( ln_drgimp ) THEN
	430	DO jj = 2, jpjm1
	431	DO ji = 2, jpim1
	432	ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
	433	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point
	434	zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va
	435	END DO
	436	END DO
	437	IF ( ln_isfcav ) THEN
	438	DO jj = 2, jpjm1
	439	DO ji = 2, jpim1
	440	ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
	441	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point
	442	zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
	443	END DO
	444	END DO
	445	ENDIF
	446	ENDIF
[456]	447
[9019]	448	! Matrix inversion
	449	!-----------------------------------------------------------------------
	450	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
[503]	451	!
[9019]	452	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	453	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	454	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	455	! ( ... )( ... ) ( ... )
	456	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
[503]	457	!
[9019]	458	! m is decomposed in the product of an upper and lower triangular matrix
	459	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	460	! The solution (after velocity) is in 2d array va
	461	!-----------------------------------------------------------------------
	462	!
	463	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
	464	DO jj = 2, jpjm1
	465	DO ji = fs_2, fs_jpim1 ! vector opt.
	466	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	467	END DO
	468	END DO
	469	END DO
	470	!
	471	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
	472	DO ji = fs_2, fs_jpim1 ! vector opt.
	473	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
	474	va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
	475	& / ( ze3va * rau0 ) * vmask(ji,jj,1)
	476	END DO
	477	END DO
	478	DO jk = 2, jpkm1
	479	DO jj = 2, jpjm1
	480	DO ji = fs_2, fs_jpim1 ! vector opt.
	481	va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
	482	END DO
	483	END DO
	484	END DO
	485	!
	486	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
	487	DO ji = fs_2, fs_jpim1 ! vector opt.
	488	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	489	END DO
	490	END DO
	491	DO jk = jpk-2, 1, -1
	492	DO jj = 2, jpjm1
	493	DO ji = fs_2, fs_jpim1
	494	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	495	END DO
	496	END DO
	497	END DO
	498	!
[503]	499	IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics
[11613]	500
	501	ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / r2dt
	502	ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / r2dt
[4990]	503	CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
[11613]	504	!
	505	IF( ln_drgimp ) THEN
	506	ztrdu(:,:,:) = 0._wp ; ztrdv(:,:,:) = 0._wp
	507	DO jk = 1, jpkm1
	508	DO jj = 2, jpjm1
	509	DO ji = 2, jpim1
	510	ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels
	511	ikbv = mbkv(ji,jj)
	512	ztrdu(ji,jj,jk) = 0.5 * ( rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj) ) &
	513	& * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
	514	ztrdv(ji,jj,jk) = 0.5 * ( rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj) ) &
	515	& * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
	516	END DO
	517	END DO
	518	END DO
	519	CALL trd_dyn( ztrdu, ztrdv, jpdyn_bfri, kt )
	520	ENDIF
	521	!
[9019]	522	DEALLOCATE( ztrdu, ztrdv )
[456]	523	ENDIF
	524	! ! print mean trends (used for debugging)
	525	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, &
[5836]	526	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
	527	!
[9019]	528	IF( ln_timing ) CALL timing_stop('dyn_zdf')
[503]	529	!
[456]	530	END SUBROUTINE dyn_zdf
	531
	532	!!==============================================================================
	533	END MODULE dynzdf

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/DYN/dynzdf.F90 @ 11613

Download in other formats: