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.

dynzdf.F90 in NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90 @ 10019

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

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

Download in other formats: