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

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

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

Download in other formats: