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_imp.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 4487

Last change on this file since 4487 was 4370, checked in by jchanut, 10 years ago
Time split cleaning / AMM12 restart / BDY and Agrif. See tickets #1228, #1227, #1225 and #1133
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 *
	4	!! Ocean dynamics: vertical component(s) of the momentum mixing trend
[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	!!----------------------------------------------------------------------
[2715]	14	!! dyn_zdf_imp : update the momentum trend with the vertical diffusion using a implicit time-stepping
[3]	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and tracers
	17	USE dom_oce ! ocean space and time domain
[4292]	18	USE domvvl ! variable volume
[888]	19	USE sbc_oce ! surface boundary condition: ocean
	20	USE zdf_oce ! ocean vertical physics
[719]	21	USE phycst ! physical constants
[3]	22	USE in_out_manager ! I/O manager
[2715]	23	USE lib_mpp ! MPP library
[3294]	24	USE zdfbfr ! Bottom friction setup
	25	USE wrk_nemo ! Memory Allocation
	26	USE timing ! Timing
[4292]	27	USE dynadv ! dynamics: vector invariant versus flux form
[4354]	28	USE dynspg_oce, ONLY: lk_dynspg_ts
[3]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
[2528]	33	PUBLIC dyn_zdf_imp ! called by step.F90
[3]	34
[4292]	35	REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
	36
[3]	37	!! * Substitutions
	38	# include "domzgr_substitute.h90"
	39	# include "vectopt_loop_substitute.h90"
	40	!!----------------------------------------------------------------------
[2528]	41	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	42	!! $Id$
[2528]	43	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	44	!!----------------------------------------------------------------------
	45	CONTAINS
	46
[503]	47	SUBROUTINE dyn_zdf_imp( kt, p2dt )
[3]	48	!!----------------------------------------------------------------------
	49	!! * ROUTINE dyn_zdf_imp *
	50	!!
	51	!! ** Purpose : Compute the trend due to the vert. momentum diffusion
	52	!! and the surface forcing, and add it to the general trend of
	53	!! the momentum equations.
	54	!!
	55	!! ** Method : The vertical momentum mixing trend is given by :
	56	!! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
	57	!! backward time stepping
[2528]	58	!! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
[3]	59	!! Bottom boundary conditions : bottom stress (cf zdfbfr.F)
	60	!! Add this trend to the general trend ua :
	61	!! ua = ua + dz( avmu dz(u) )
	62	!!
[2528]	63	!! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
[3]	64	!!---------------------------------------------------------------------
[3294]	65	INTEGER , INTENT(in) :: kt ! ocean time-step index
[2715]	66	REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step
[2528]	67	!!
[2715]	68	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	69	INTEGER :: ikbu, ikbv ! local integers
[2715]	70	REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars
[4292]	71	REAL(wp) :: ze3ua, ze3va
[3]	72	!!----------------------------------------------------------------------
	73
[3294]	74	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws
	75	!!----------------------------------------------------------------------
	76	!
	77	IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_imp')
	78	!
	79	CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws )
	80	!
[3]	81	IF( kt == nit000 ) THEN
	82	IF(lwp) WRITE(numout,*)
	83	IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
	84	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
[4292]	85	!
	86	IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator
	87	ELSE ; r_vvl = 0._wp
	88	ENDIF
[3]	89	ENDIF
	90
	91	! 0. Local constant initialization
	92	! --------------------------------
[2528]	93	z1_p2dt = 1._wp / p2dt ! inverse of the timestep
[455]	94
[3294]	95	! 1. Apply semi-implicit bottom friction
	96	! --------------------------------------
	97	! Only needed for semi-implicit bottom friction setup. The explicit
	98	! bottom friction has been included in "u(v)a" which act as the R.H.S
	99	! column vector of the tri-diagonal matrix equation
	100	!
	101
	102	IF( ln_bfrimp ) THEN
	103	# if defined key_vectopt_loop
[4292]	104	DO jj = 1, 1
	105	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
[3294]	106	# else
[4292]	107	DO jj = 2, jpjm1
	108	DO ji = 2, jpim1
[3294]	109	# endif
[4292]	110	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	111	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
	112	avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
	113	avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
	114	END DO
[3294]	115	END DO
	116	ENDIF
	117
[4292]	118	#if defined key_dynspg_ts
	119	IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity
	120	DO jk = 1, jpkm1
	121	ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
	122	va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
	123	END DO
	124	ELSE ! applied on thickness weighted velocity
	125	DO jk = 1, jpkm1
	126	ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &
	127	& + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &
	128	& / fse3u_a(:,:,jk) * umask(:,:,jk)
	129	va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &
	130	& + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &
	131	& / fse3v_a(:,:,jk) * vmask(:,:,jk)
	132	END DO
	133	ENDIF
	134
	135	IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
	136	! remove barotropic velocities:
	137	DO jk = 1, jpkm1
	138	ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
	139	va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
	140	ENDDO
	141	! Add bottom stress due to barotropic component only:
	142	DO jj = 2, jpjm1
	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)
	146	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)
	147	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)
	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
	152	ENDIF
	153	#endif
	154
[3294]	155	! 2. Vertical diffusion on u
[3]	156	! ---------------------------
	157	! Matrix and second member construction
[1662]	158	! bottom boundary condition: both zwi and zws must be masked as avmu can take
[3294]	159	! non zero value at the ocean bottom depending on the bottom friction used.
[2528]	160	!
	161	DO jk = 1, jpkm1 ! Matrix
[3]	162	DO jj = 2, jpjm1
	163	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	164	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl * fse3u_a(ji,jj,jk) ! after scale factor at T-point
	165	zcoef = - p2dt / ze3ua
[2528]	166	zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )
	167	zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
	168	zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
	169	zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
	170	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
[3]	171	END DO
	172	END DO
	173	END DO
[4292]	174	DO jj = 2, jpjm1 ! Surface boundary conditions
[3]	175	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	176	zwi(ji,jj,1) = 0._wp
	177	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	178	END DO
	179	END DO
	180
	181	! Matrix inversion starting from the first level
	182	!-----------------------------------------------------------------------
	183	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	184	!
	185	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	186	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	187	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	188	! ( ... )( ... ) ( ... )
	189	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	190	!
	191	! m is decomposed in the product of an upper and a lower triangular matrix
	192	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	193	! The solution (the after velocity) is in ua
	194	!-----------------------------------------------------------------------
[2528]	195	!
	196	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	197	DO jj = 2, jpjm1
	198	DO ji = fs_2, fs_jpim1 ! vector opt.
	199	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	200	END DO
	201	END DO
	202	END DO
[2528]	203	!
	204	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	205	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	206	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1)
	207	#if defined key_dynspg_ts
	208	ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
	209	& / ( ze3ua * rau0 )
	210	#else
	211	ua(ji,jj,1) = ub(ji,jj,1) + p2dt (ua(ji,jj,1) + 0.5_wp ( utau_b(ji,jj) + utau(ji,jj) ) &
	212	& / ( fse3u(ji,jj,1) * rau0 ) )
	213	#endif
[3]	214	END DO
	215	END DO
	216	DO jk = 2, jpkm1
	217	DO jj = 2, jpjm1
	218	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	219	#if defined key_dynspg_ts
	220	zrhs = ua(ji,jj,jk) ! zrhs=right hand side
	221	#else
	222	zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
	223	#endif
[3]	224	ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
	225	END DO
	226	END DO
	227	END DO
[2528]	228	!
	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)
	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 ! vector opt.
[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
	241
[4292]	242	#if ! defined key_dynspg_ts
[3]	243	! Normalization to obtain the general momentum trend ua
	244	DO jk = 1, jpkm1
	245	DO jj = 2, jpjm1
	246	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	247	ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
[3]	248	END DO
	249	END DO
	250	END DO
[4292]	251	#endif
[3]	252
[3294]	253	! 3. Vertical diffusion on v
[3]	254	! ---------------------------
	255	! Matrix and second member construction
[1662]	256	! bottom boundary condition: both zwi and zws must be masked as avmv can take
[3294]	257	! non zero value at the ocean bottom depending on the bottom friction used
[2528]	258	!
	259	DO jk = 1, jpkm1 ! Matrix
[3]	260	DO jj = 2, jpjm1
	261	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	262	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk) + r_vvl * fse3v_a(ji,jj,jk) ! after scale factor at T-point
	263	zcoef = - p2dt / ze3va
[2528]	264	zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )
[1662]	265	zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk)
[2528]	266	zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
[3]	267	zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1)
[2528]	268	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
[3]	269	END DO
	270	END DO
	271	END DO
[4292]	272	DO jj = 2, jpjm1 ! Surface boundary conditions
[3]	273	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	274	zwi(ji,jj,1) = 0._wp
	275	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	276	END DO
	277	END DO
	278
	279	! Matrix inversion
	280	!-----------------------------------------------------------------------
	281	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	282	!
	283	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	284	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	285	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	286	! ( ... )( ... ) ( ... )
	287	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	288	!
[2528]	289	! m is decomposed in the product of an upper and lower triangular matrix
[3]	290	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	291	! The solution (after velocity) is in 2d array va
	292	!-----------------------------------------------------------------------
[2528]	293	!
	294	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	295	DO jj = 2, jpjm1
	296	DO ji = fs_2, fs_jpim1 ! vector opt.
	297	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	298	END DO
	299	END DO
	300	END DO
[2528]	301	!
	302	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	303	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	304	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1)
	305	#if defined key_dynspg_ts
	306	va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
	307	& / ( ze3va * rau0 )
	308	#else
	309	va(ji,jj,1) = vb(ji,jj,1) + p2dt (va(ji,jj,1) + 0.5_wp ( vtau_b(ji,jj) + vtau(ji,jj) ) &
	310	& / ( fse3v(ji,jj,1) * rau0 ) )
	311	#endif
[3]	312	END DO
	313	END DO
	314	DO jk = 2, jpkm1
	315	DO jj = 2, jpjm1
	316	DO ji = fs_2, fs_jpim1 ! vector opt.
[4292]	317	#if defined key_dynspg_ts
	318	zrhs = va(ji,jj,jk) ! zrhs=right hand side
	319	#else
	320	zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
	321	#endif
[3]	322	va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
	323	END DO
	324	END DO
	325	END DO
[2528]	326	!
[4292]	327	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==
[3]	328	DO ji = fs_2, fs_jpim1 ! vector opt.
	329	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	330	END DO
	331	END DO
	332	DO jk = jpk-2, 1, -1
	333	DO jj = 2, jpjm1
	334	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	335	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]	336	END DO
	337	END DO
	338	END DO
	339
	340	! Normalization to obtain the general momentum trend va
[4292]	341	#if ! defined key_dynspg_ts
[3]	342	DO jk = 1, jpkm1
	343	DO jj = 2, jpjm1
	344	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	345	va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
[3]	346	END DO
	347	END DO
	348	END DO
[4292]	349	#endif
[3294]	350
[4292]	351	! J. Chanut: Lines below are useless ?
[3294]	352	!! restore bottom layer avmu(v)
	353	IF( ln_bfrimp ) THEN
	354	# if defined key_vectopt_loop
	355	DO jj = 1, 1
	356	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
	357	# else
	358	DO jj = 2, jpjm1
	359	DO ji = 2, jpim1
	360	# endif
	361	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
	362	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
[4292]	363	avmu(ji,jj,ikbu+1) = 0.e0
	364	avmv(ji,jj,ikbv+1) = 0.e0
[3294]	365	END DO
	366	END DO
	367	ENDIF
[2528]	368	!
[3294]	369	CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)
[2715]	370	!
[3294]	371	IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp')
	372	!
[3]	373	END SUBROUTINE dyn_zdf_imp
	374
	375	!!==============================================================================
	376	END MODULE dynzdf_imp

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

Download in other formats: