Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

dynzdf_imp.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 16.0 KB

Rev	Line
[3]	1	MODULE dynzdf_imp
[2715]	2	!!======================================================================
[3]	3	!! * MODULE dynzdf_imp *
[6140]	4	!! Ocean dynamics: vertical component(s) of the momentum mixing trend, implicit scheme
[2715]	5	!!======================================================================
[2528]	6	!! History : OPA ! 1990-10 (B. Blanke) Original code
	7	!! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal
[2715]	8	!! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module
[2528]	9	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps
[3294]	10	!! 3.4 ! 2012-01 (H. Liu) Semi-implicit bottom friction
[503]	11	!!----------------------------------------------------------------------
[3]	12
	13	!!----------------------------------------------------------------------
[6140]	14	!! dyn_zdf_imp : compute the vertical diffusion using a implicit scheme
	15	!! together with the Leap-Frog time integration.
[3]	16	!!----------------------------------------------------------------------
[6140]	17	USE oce ! ocean dynamics and tracers
	18	USE phycst ! physical constants
	19	USE dom_oce ! ocean space and time domain
	20	USE domvvl ! variable volume
	21	USE sbc_oce ! surface boundary condition: ocean
	22	USE dynadv , ONLY: ln_dynadv_vec ! Momentum advection form
	23	USE zdf_oce ! ocean vertical physics
	24	USE zdfbfr ! Bottom friction setup
	25	!
	26	USE in_out_manager ! I/O manager
	27	USE lib_mpp ! MPP library
	28	USE timing ! Timing
[3]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
[2528]	33	PUBLIC dyn_zdf_imp ! called by step.F90
[3]	34
[6140]	35	REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
[4292]	36
[3]	37	!! * Substitutions
	38	# include "vectopt_loop_substitute.h90"
	39	!!----------------------------------------------------------------------
[2528]	40	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	41	!! $Id$
[2528]	42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	43	!!----------------------------------------------------------------------
	44	CONTAINS
	45
[503]	46	SUBROUTINE dyn_zdf_imp( kt, p2dt )
[3]	47	!!----------------------------------------------------------------------
	48	!! * ROUTINE dyn_zdf_imp *
	49	!!
	50	!! ** Purpose : Compute the trend due to the vert. momentum diffusion
[6140]	51	!! together with the Leap-Frog time stepping using an
	52	!! implicit scheme.
[3]	53	!!
[6140]	54	!! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing
	55	!! ua = ub + 2dt ua vector form or linear free surf.
	56	!! ua = ( e3u_bub + 2dt * e3u_n*ua ) / e3u_a otherwise
	57	!! - update the after velocity with the implicit vertical mixing.
	58	!! This requires to solver the following system:
	59	!! ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ]
	60	!! with the following surface/top/bottom boundary condition:
	61	!! surface: wind stress input (averaged over kt-1/2 & kt+1/2)
	62	!! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfbfr.F)
[3]	63	!!
[6140]	64	!! ** Action : (ua,va) after velocity
[3]	65	!!---------------------------------------------------------------------
[3294]	66	INTEGER , INTENT(in) :: kt ! ocean time-step index
[2715]	67	REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step
[6140]	68	!
	69	INTEGER :: ji, jj, jk ! dummy loop indices
	70	INTEGER :: ikbu, ikbv ! local integers
	71	REAL(wp) :: zzwi, ze3ua ! local scalars
	72	REAL(wp) :: zzws, ze3va ! - -
[7910]	73	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws
[3294]	74	!!----------------------------------------------------------------------
	75	!
	76	IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_imp')
	77	!
	78	!
[3]	79	IF( kt == nit000 ) THEN
	80	IF(lwp) WRITE(numout,*)
	81	IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
	82	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
[4292]	83	!
[6140]	84	If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator
	85	ELSE ; r_vvl = 1._wp
[4292]	86	ENDIF
[3]	87	ENDIF
[6140]	88	!
	89	! !== Time step dynamics ==!
	90	!
	91	IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
[5930]	92	DO jk = 1, jpkm1
[7753]	93	ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
	94	va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
[5930]	95	END DO
[6140]	96	ELSE ! applied on thickness weighted velocity
[5930]	97	DO jk = 1, jpkm1
[7753]	98	ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) &
	99	& + p2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)
	100	va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) &
	101	& + p2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)
[5930]	102	END DO
	103	ENDIF
[6140]	104	!
	105	! !== Apply semi-implicit bottom friction ==!
	106	!
[3294]	107	! Only needed for semi-implicit bottom friction setup. The explicit
	108	! bottom friction has been included in "u(v)a" which act as the R.H.S
	109	! column vector of the tri-diagonal matrix equation
	110	!
	111	IF( ln_bfrimp ) THEN
[4292]	112	DO jj = 2, jpjm1
	113	DO ji = 2, jpim1
	114	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	115	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
[6140]	116	avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1)
	117	avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1)
[4292]	118	END DO
[3294]	119	END DO
[5120]	120	IF ( ln_isfcav ) THEN
	121	DO jj = 2, jpjm1
	122	DO ji = 2, jpim1
	123	ikbu = miku(ji,jj) ! ocean top level at u- and v-points
	124	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
[6140]	125	IF( ikbu >= 2 ) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu)
	126	IF( ikbv >= 2 ) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv)
[5120]	127	END DO
	128	END DO
	129	END IF
[3294]	130	ENDIF
[6140]	131	!
[5930]	132	! With split-explicit free surface, barotropic stress is treated explicitly
	133	! Update velocities at the bottom.
	134	! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
	135	! not lead to the effective stress seen over the whole barotropic loop.
[6140]	136	! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
	137	IF( ln_bfrimp .AND. ln_dynspg_ts ) THEN
	138	DO jk = 1, jpkm1 ! remove barotropic velocities
[7753]	139	ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
	140	va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
[4990]	141	END DO
[6140]	142	DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only
[4292]	143	DO ji = fs_2, fs_jpim1 ! vector opt.
	144	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	145	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
[6140]	146	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
	147	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
[4292]	148	ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
	149	va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
	150	END DO
	151	END DO
[6140]	152	IF( ln_isfcav ) THEN ! Ocean cavities (ISF)
[5120]	153	DO jj = 2, jpjm1
	154	DO ji = fs_2, fs_jpim1 ! vector opt.
	155	ikbu = miku(ji,jj) ! top ocean level at u- and v-points
	156	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
[6140]	157	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
	158	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
[5120]	159	ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
	160	va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
	161	END DO
	162	END DO
	163	END IF
[4292]	164	ENDIF
[6140]	165	!
	166	! !== Vertical diffusion on u ==!
	167	!
[3]	168	! Matrix and second member construction
[1662]	169	! bottom boundary condition: both zwi and zws must be masked as avmu can take
[3294]	170	! non zero value at the ocean bottom depending on the bottom friction used.
[2528]	171	!
	172	DO jk = 1, jpkm1 ! Matrix
[3]	173	DO jj = 2, jpjm1
	174	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	175	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point
	176	zzwi = - p2dt * avmu(ji,jj,jk ) / ( ze3ua * e3uw_n(ji,jj,jk ) )
	177	zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) )
	178	zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk )
	179	zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1)
[5120]	180	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
[3]	181	END DO
	182	END DO
	183	END DO
[4292]	184	DO jj = 2, jpjm1 ! Surface boundary conditions
[3]	185	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	186	zwi(ji,jj,1) = 0._wp
	187	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	188	END DO
	189	END DO
	190
	191	! Matrix inversion starting from the first level
	192	!-----------------------------------------------------------------------
	193	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	194	!
	195	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	196	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	197	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	198	! ( ... )( ... ) ( ... )
	199	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	200	!
	201	! m is decomposed in the product of an upper and a lower triangular matrix
	202	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	203	! The solution (the after velocity) is in ua
	204	!-----------------------------------------------------------------------
[2528]	205	!
[5836]	206	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[5120]	207	DO jj = 2, jpjm1
	208	DO ji = fs_2, fs_jpim1 ! vector opt.
[3]	209	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	210	END DO
	211	END DO
	212	END DO
[2528]	213	!
[6140]	214	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
[3]	215	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	216	ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)
[5120]	217	ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
	218	& / ( ze3ua * rau0 ) * umask(ji,jj,1)
	219	END DO
	220	END DO
	221	DO jk = 2, jpkm1
	222	DO jj = 2, jpjm1
	223	DO ji = fs_2, fs_jpim1
[6140]	224	ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
[3]	225	END DO
	226	END DO
	227	END DO
[2528]	228	!
[6140]	229	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
[3]	230	DO ji = fs_2, fs_jpim1 ! vector opt.
	231	ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
[5120]	232	END DO
	233	END DO
	234	DO jk = jpk-2, 1, -1
	235	DO jj = 2, jpjm1
	236	DO ji = fs_2, fs_jpim1
[2528]	237	ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]	238	END DO
	239	END DO
	240	END DO
[6140]	241	!
	242	! !== Vertical diffusion on v ==!
	243	!
[3]	244	! Matrix and second member construction
[1662]	245	! bottom boundary condition: both zwi and zws must be masked as avmv can take
[3294]	246	! non zero value at the ocean bottom depending on the bottom friction used
[2528]	247	!
	248	DO jk = 1, jpkm1 ! Matrix
[3]	249	DO jj = 2, jpjm1
	250	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	251	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point
	252	zzwi = - p2dt * avmv (ji,jj,jk ) / ( ze3va * e3vw_n(ji,jj,jk ) )
	253	zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) )
	254	zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk )
	255	zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
[5120]	256	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
[3]	257	END DO
	258	END DO
	259	END DO
[4292]	260	DO jj = 2, jpjm1 ! Surface boundary conditions
[3]	261	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	262	zwi(ji,jj,1) = 0._wp
	263	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	264	END DO
	265	END DO
	266
	267	! Matrix inversion
	268	!-----------------------------------------------------------------------
	269	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	270	!
	271	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	272	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	273	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	274	! ( ... )( ... ) ( ... )
	275	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	276	!
[2528]	277	! m is decomposed in the product of an upper and lower triangular matrix
[3]	278	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	279	! The solution (after velocity) is in 2d array va
	280	!-----------------------------------------------------------------------
[2528]	281	!
[5836]	282	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[5120]	283	DO jj = 2, jpjm1
	284	DO ji = fs_2, fs_jpim1 ! vector opt.
[3]	285	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	286	END DO
	287	END DO
	288	END DO
[2528]	289	!
[6140]	290	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
[5930]	291	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	292	ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)
[5120]	293	va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
[6752]	294	& / ( ze3va * rau0 ) * vmask(ji,jj,1)
[5120]	295	END DO
	296	END DO
	297	DO jk = 2, jpkm1
	298	DO jj = 2, jpjm1
	299	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	300	va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
[3]	301	END DO
	302	END DO
	303	END DO
[2528]	304	!
[6140]	305	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
[3]	306	DO ji = fs_2, fs_jpim1 ! vector opt.
	307	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
[5120]	308	END DO
	309	END DO
	310	DO jk = jpk-2, 1, -1
	311	DO jj = 2, jpjm1
	312	DO ji = fs_2, fs_jpim1
[2528]	313	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]	314	END DO
	315	END DO
	316	END DO
[6140]	317
[4292]	318	! J. Chanut: Lines below are useless ?
[6140]	319	!! restore bottom layer avmu(v)
	320	!!gm I almost sure it is !!!!
[3294]	321	IF( ln_bfrimp ) THEN
[4990]	322	DO jj = 2, jpjm1
	323	DO ji = 2, jpim1
	324	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	325	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
[6140]	326	avmu(ji,jj,ikbu+1) = 0._wp
	327	avmv(ji,jj,ikbv+1) = 0._wp
[4990]	328	END DO
	329	END DO
[5120]	330	IF (ln_isfcav) THEN
	331	DO jj = 2, jpjm1
	332	DO ji = 2, jpim1
	333	ikbu = miku(ji,jj) ! ocean top level at u- and v-points
	334	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
[6140]	335	IF( ikbu > 1 ) avmu(ji,jj,ikbu) = 0._wp
	336	IF( ikbv > 1 ) avmv(ji,jj,ikbv) = 0._wp
[5120]	337	END DO
	338	END DO
[6140]	339	ENDIF
[3294]	340	ENDIF
[2528]	341	!
[2715]	342	!
[6140]	343	IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp')
[3294]	344	!
[3]	345	END SUBROUTINE dyn_zdf_imp
	346
	347	!!==============================================================================
	348	END MODULE dynzdf_imp

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 7910

Download in other formats: