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

Last change on this file since 9976 was 9976, checked in by acc, 6 years ago
Branch: dev_r9956_ENHANCE05_ZAD_AIMP. First set of changes to implement an adaptive-implicit vertical advection option (see ticket #2042). This code compiles and runs but has issues when the new option is activated.
Property svn:keywords set to `Id`
File size: 30.1 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 ! - -
[9976]	73	REAL(wp) :: z1_e3un, z1_e3vn ! - -
	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.
	168	z1_e3un = 1._wp / e3u_n(ji,jj,jk)
	169	zzwi = ( ( avm (ji+1,jj,jk ) + avm (ji,jj,jk ) + akzu(ji,jj,jk ) ) &
	170	& / e3uw_a(ji ,jj,jk ) ) * z1_e3un * wumask(ji,jj,jk )
	171	zzws = ( ( avm (ji+1,jj,jk+1) + avm (ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
	172	& / e3uw_a(ji ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
	173	zWu = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) &
	174	& + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
	175	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
	176	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
	177	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
	178	END DO
	179	END DO
[9019]	180	END DO
[9976]	181	CASE DEFAULT ! iso-level lateral mixing
	182	DO jk = 1, jpkm1
	183	DO jj = 2, jpjm1
	184	DO ji = fs_2, fs_jpim1 ! vector opt.
	185	z1_e3un = 1._wp / e3u_n(ji,jj,jk)
	186	zzwi = ( ( avm (ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
	187	& / e3uw_a(ji ,jj,jk ) ) * z1_e3un * wumask(ji,jj,jk )
	188	zzws = ( ( avm (ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
	189	& / e3uw_a(ji ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1)
	190	zWu = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) &
	191	& + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
	192	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un )
	193	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un )
	194	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un )
	195	END DO
	196	END DO
	197	END DO
	198	END SELECT
	199	ELSE !== Flux form advection ==!
	200	SELECT CASE( nldf_dyn )
	201	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
	202	DO jk = 1, jpkm1
	203	DO jj = 2, jpjm1
	204	DO ji = fs_2, fs_jpim1 ! vector opt.
	205	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	206	zzwi = ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
	207	& / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk )
	208	zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
	209	& / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	210	zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )
	211	zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
	212	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWui, 0._wp ) )
	213	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWus, 0._wp ) )
	214	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) )
	215	END DO
	216	END DO
	217	END DO
	218	CASE DEFAULT ! iso-level lateral mixing
	219	DO jk = 1, jpkm1
	220	DO jj = 2, jpjm1
	221	DO ji = fs_2, fs_jpim1 ! vector opt.
	222	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	223	zzwi = ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk )
	224	zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	225	zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) )
	226	zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) )
	227	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWui, 0._wp ) )
	228	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWus, 0._wp ) )
	229	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) )
	230	END DO
	231	END DO
	232	END DO
	233	END SELECT
	234	ENDIF
	235	ELSE
	236	SELECT CASE( nldf_dyn )
	237	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
	238	DO jk = 1, jpkm1
	239	DO jj = 2, jpjm1
	240	DO ji = fs_2, fs_jpim1 ! vector opt.
	241	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	242	zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
	243	& / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk )
	244	zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
	245	& / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	246	zwi(ji,jj,jk) = zzwi
	247	zws(ji,jj,jk) = zzws
	248	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	249	END DO
	250	END DO
[9019]	251	END DO
[9976]	252	CASE DEFAULT ! iso-level lateral mixing
	253	DO jk = 1, jpkm1
	254	DO jj = 2, jpjm1
	255	DO ji = fs_2, fs_jpim1 ! vector opt.
	256	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point
	257	zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk )
	258	zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
	259	zwi(ji,jj,jk) = zzwi
	260	zws(ji,jj,jk) = zzws
	261	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	262	END DO
[9019]	263	END DO
	264	END DO
[9976]	265	END SELECT
	266	ENDIF
[9019]	267	!
	268	DO jj = 2, jpjm1 !* Surface boundary conditions
	269	DO ji = fs_2, fs_jpim1 ! vector opt.
	270	zwi(ji,jj,1) = 0._wp
	271	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
	272	END DO
	273	END DO
	274	!
	275	! !== Apply semi-implicit bottom friction ==!
	276	!
	277	! Only needed for semi-implicit bottom friction setup. The explicit
	278	! bottom friction has been included in "u(v)a" which act as the R.H.S
	279	! column vector of the tri-diagonal matrix equation
	280	!
	281	IF ( ln_drgimp ) THEN ! implicit bottom friction
	282	DO jj = 2, jpjm1
	283	DO ji = 2, jpim1
	284	iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
	285	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point
	286	zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
	287	END DO
	288	END DO
	289	IF ( ln_isfcav ) THEN ! top friction (always implicit)
	290	DO jj = 2, jpjm1
	291	DO ji = 2, jpim1
	292	!!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed
	293	iku = miku(ji,jj) ! ocean top level at u- and v-points
	294	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point
	295	zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
	296	END DO
	297	END DO
	298	END IF
	299	ENDIF
	300	!
	301	! Matrix inversion starting from the first level
	302	!-----------------------------------------------------------------------
	303	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	304	!
	305	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	306	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	307	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	308	! ( ... )( ... ) ( ... )
	309	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	310	!
	311	! m is decomposed in the product of an upper and a lower triangular matrix
	312	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	313	! The solution (the after velocity) is in ua
	314	!-----------------------------------------------------------------------
	315	!
	316	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
	317	DO jj = 2, jpjm1
	318	DO ji = fs_2, fs_jpim1 ! vector opt.
	319	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	320	END DO
	321	END DO
	322	END DO
	323	!
	324	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
	325	DO ji = fs_2, fs_jpim1 ! vector opt.
	326	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
	327	ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
	328	& / ( ze3ua * rau0 ) * umask(ji,jj,1)
	329	END DO
	330	END DO
	331	DO jk = 2, jpkm1
	332	DO jj = 2, jpjm1
	333	DO ji = fs_2, fs_jpim1
	334	ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
	335	END DO
	336	END DO
	337	END DO
	338	!
	339	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
	340	DO ji = fs_2, fs_jpim1 ! vector opt.
	341	ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	342	END DO
	343	END DO
	344	DO jk = jpk-2, 1, -1
	345	DO jj = 2, jpjm1
	346	DO ji = fs_2, fs_jpim1
	347	ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	348	END DO
	349	END DO
	350	END DO
	351	!
	352	! !== Vertical diffusion on v ==!
	353	!
	354	! !* Matrix construction
	355	zdt = r2dt * 0.5
[9976]	356	IF( ln_zad_Aimp ) THEN !!
	357	IF( ln_dynadv_vec ) THEN !== Vector invariant advection ==!
	358	SELECT CASE( nldf_dyn )
	359	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
	360	DO jk = 1, jpkm1
	361	DO jj = 2, jpjm1
	362	DO ji = fs_2, fs_jpim1 ! vector opt.
	363	z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
	364	zzwi = ( ( avm (ji,jj+1,jk ) + avm (ji,jj,jk ) + akzv(ji,jj,jk ) ) &
	365	& / e3vw_a(ji,jj ,jk ) ) * z1_e3vn * wvmask(ji,jj,jk )
	366	zzws = ( ( avm (ji,jj+1,jk+1) + avm (ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
	367	& / e3vw_a(ji,jj ,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
	368	zWv = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) &
	369	& + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
	370	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
	371	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
	372	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
	373	END DO
	374	END DO
[9019]	375	END DO
[9976]	376	CASE DEFAULT ! iso-level lateral mixing
	377	DO jk = 1, jpkm1
	378	DO jj = 2, jpjm1
	379	DO ji = fs_2, fs_jpim1 ! vector opt.
	380	z1_e3vn = 1._wp / e3v_n(ji,jj,jk)
	381	zzwi = ( ( avm (ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
	382	& / e3vw_a(ji,jj ,jk ) ) * z1_e3vn * wvmask(ji,jj,jk )
	383	zzws = ( ( avm (ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
	384	& / e3vw_a(ji ,jj,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1)
	385	zWv = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) &
	386	& + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
	387	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn )
	388	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn )
	389	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn )
	390	END DO
	391	END DO
	392	END DO
	393	END SELECT
	394	ELSE !== Flux form advection ==!
	395	SELECT CASE( nldf_dyn )
	396	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
	397	DO jk = 1, jpkm1
	398	DO jj = 2, jpjm1
	399	DO ji = fs_2, fs_jpim1 ! vector opt.
	400	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at U-point
	401	zzwi = ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
	402	& / ( ze3va * e3vw_a(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	403	zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
	404	& / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
	405	zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) )
	406	zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
	407	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWvi, 0._wp ) )
	408	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWvs, 0._wp ) )
	409	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
	410	END DO
	411	END DO
	412	END DO
	413	CASE DEFAULT ! iso-level lateral mixing
	414	DO jk = 1, jpkm1
	415	DO jj = 2, jpjm1
	416	DO ji = fs_2, fs_jpim1 ! vector opt.
	417	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at U-point
	418	zzwi = ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_a(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	419	zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
	420	zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) )
	421	zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) )
	422	zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWvi, 0._wp ) )
	423	zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWvs, 0._wp ) )
	424	zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
	425	END DO
	426	END DO
	427	END DO
	428	END SELECT
	429	ENDIF
	430	ELSE
	431	SELECT CASE( nldf_dyn )
	432	CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
	433	DO jk = 1, jpkm1
	434	DO jj = 2, jpjm1
	435	DO ji = fs_2, fs_jpim1 ! vector opt.
	436	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point
	437	zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
	438	& / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	439	zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
	440	& / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
	441	zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk )
	442	zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
	443	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	444	END DO
	445	END DO
[9019]	446	END DO
[9976]	447	CASE DEFAULT ! iso-level lateral mixing
	448	DO jk = 1, jpkm1
	449	DO jj = 2, jpjm1
	450	DO ji = fs_2, fs_jpim1 ! vector opt.
	451	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point
	452	zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )
	453	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)
	454	zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk )
	455	zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
	456	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
	457	END DO
[9019]	458	END DO
	459	END DO
[9976]	460	END SELECT
	461	ENDIF
[9019]	462	!
	463	DO jj = 2, jpjm1 !* Surface boundary conditions
	464	DO ji = fs_2, fs_jpim1 ! vector opt.
	465	zwi(ji,jj,1) = 0._wp
	466	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
	467	END DO
	468	END DO
	469	! !== Apply semi-implicit top/bottom friction ==!
	470	!
	471	! Only needed for semi-implicit bottom friction setup. The explicit
	472	! bottom friction has been included in "u(v)a" which act as the R.H.S
	473	! column vector of the tri-diagonal matrix equation
	474	!
	475	IF( ln_drgimp ) THEN
	476	DO jj = 2, jpjm1
	477	DO ji = 2, jpim1
	478	ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
	479	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point
	480	zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va
	481	END DO
	482	END DO
	483	IF ( ln_isfcav ) THEN
	484	DO jj = 2, jpjm1
	485	DO ji = 2, jpim1
	486	ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
	487	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point
	488	zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
	489	END DO
	490	END DO
	491	ENDIF
	492	ENDIF
[456]	493
[9019]	494	! Matrix inversion
	495	!-----------------------------------------------------------------------
	496	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
[503]	497	!
[9019]	498	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	499	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	500	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	501	! ( ... )( ... ) ( ... )
	502	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
[503]	503	!
[9019]	504	! m is decomposed in the product of an upper and lower triangular matrix
	505	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	506	! The solution (after velocity) is in 2d array va
	507	!-----------------------------------------------------------------------
	508	!
	509	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
	510	DO jj = 2, jpjm1
	511	DO ji = fs_2, fs_jpim1 ! vector opt.
	512	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	513	END DO
	514	END DO
	515	END DO
	516	!
	517	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
	518	DO ji = fs_2, fs_jpim1 ! vector opt.
	519	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
	520	va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
	521	& / ( ze3va * rau0 ) * vmask(ji,jj,1)
	522	END DO
	523	END DO
	524	DO jk = 2, jpkm1
	525	DO jj = 2, jpjm1
	526	DO ji = fs_2, fs_jpim1 ! vector opt.
	527	va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
	528	END DO
	529	END DO
	530	END DO
	531	!
	532	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
	533	DO ji = fs_2, fs_jpim1 ! vector opt.
	534	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	535	END DO
	536	END DO
	537	DO jk = jpk-2, 1, -1
	538	DO jj = 2, jpjm1
	539	DO ji = fs_2, fs_jpim1
	540	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	541	END DO
	542	END DO
	543	END DO
	544	!
[503]	545	IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics
[7753]	546	ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
	547	ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
[4990]	548	CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
[9019]	549	DEALLOCATE( ztrdu, ztrdv )
[456]	550	ENDIF
	551	! ! print mean trends (used for debugging)
	552	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, &
[5836]	553	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
	554	!
[9019]	555	IF( ln_timing ) CALL timing_stop('dyn_zdf')
[503]	556	!
[456]	557	END SUBROUTINE dyn_zdf
	558
	559	!!==============================================================================
	560	END MODULE dynzdf

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

Download in other formats: