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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 4409

Last change on this file since 4409 was 4406, checked in by trackstand2, 10 years ago
Move from jpk to jpkf - trim sub-domains in z
Property svn:keywords set to `Id`
File size: 15.9 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
[503]	10	!!----------------------------------------------------------------------
[3]	11
	12	!!----------------------------------------------------------------------
[2715]	13	!! dyn_zdf_imp : update the momentum trend with the vertical diffusion using a implicit time-stepping
[3]	14	!!----------------------------------------------------------------------
	15	USE oce ! ocean dynamics and tracers
	16	USE dom_oce ! ocean space and time domain
[888]	17	USE sbc_oce ! surface boundary condition: ocean
	18	USE zdf_oce ! ocean vertical physics
[719]	19	USE phycst ! physical constants
[3]	20	USE in_out_manager ! I/O manager
[2715]	21	USE lib_mpp ! MPP library
[3]	22
	23	IMPLICIT NONE
	24	PRIVATE
	25
[2528]	26	PUBLIC dyn_zdf_imp ! called by step.F90
[3]	27
[3211]	28	!! * Control permutation of array indices
	29	# include "oce_ftrans.h90"
	30	# include "dom_oce_ftrans.h90"
	31	# include "sbc_oce_ftrans.h90"
	32	# include "zdf_oce_ftrans.h90"
	33
[3]	34	!! * Substitutions
	35	# include "domzgr_substitute.h90"
	36	# include "vectopt_loop_substitute.h90"
	37	!!----------------------------------------------------------------------
[2528]	38	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	39	!! $Id$
[2528]	40	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	41	!!----------------------------------------------------------------------
	42	CONTAINS
	43
[503]	44	SUBROUTINE dyn_zdf_imp( kt, p2dt )
[3]	45	!!----------------------------------------------------------------------
	46	!! * ROUTINE dyn_zdf_imp *
	47	!!
	48	!! ** Purpose : Compute the trend due to the vert. momentum diffusion
	49	!! and the surface forcing, and add it to the general trend of
	50	!! the momentum equations.
	51	!!
	52	!! ** Method : The vertical momentum mixing trend is given by :
	53	!! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
	54	!! backward time stepping
[2528]	55	!! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
[3]	56	!! Bottom boundary conditions : bottom stress (cf zdfbfr.F)
	57	!! Add this trend to the general trend ua :
	58	!! ua = ua + dz( avmu dz(u) )
	59	!!
[2528]	60	!! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
[3]	61	!!---------------------------------------------------------------------
[2715]	62	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	63	USE oce , ONLY: zwd => ta , zws => sa ! (ta,sa) used as 3D workspace
	64	USE wrk_nemo, ONLY: zwi => wrk_3d_3 ! 3D workspace
[3211]	65	!! DCSE_NEMO: need additional directives for renamed module variables
	66	!FTRANS zwd :I :I :z
	67	!FTRANS zws :I :I :z
	68	!FTRANS zwi :I :I :z
[2528]	69	!!
[2715]	70	INTEGER , INTENT(in) :: kt ! ocean time-step index
	71	REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step
[2528]	72	!!
[2715]	73	INTEGER :: ji, jj, jk ! dummy loop indices
[3211]	74	REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs, zzwibd ! local scalars
[3]	75	!!----------------------------------------------------------------------
	76
[2715]	77	IF( wrk_in_use(3, 3) ) THEN
	78	CALL ctl_stop('dyn_zdf_imp: requested workspace array unavailable') ; RETURN
	79	END IF
	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,*) '~~~~~~~~~~~ '
	85	ENDIF
	86
	87	! 0. Local constant initialization
	88	! --------------------------------
[2528]	89	z1_p2dt = 1._wp / p2dt ! inverse of the timestep
[455]	90
[3837]	91
[3]	92	! 1. Vertical diffusion on u
	93	! ---------------------------
	94	! Matrix and second member construction
[1662]	95	! bottom boundary condition: both zwi and zws must be masked as avmu can take
[3]	96	! non zero value at the ocean bottom depending on the bottom friction
[1662]	97	! used but the bottom velocities have already been updated with the bottom
	98	! friction velocity in dyn_bfr using values from the previous timestep. There
	99	! is no need to include these in the implicit calculation.
[2528]	100	!
[3211]	101	#if defined key_z_first
	102	DO jj = 2, jpjm1
	103	DO ji = 2, jpim1
[4406]	104	DO jk = 1, jpkfm1
[3211]	105	zcoef = - p2dt / fse3u(ji,jj,jk)
	106	zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )
	107	zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
	108	zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
	109	zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
	110	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
	111	END DO
	112	! Surface boundary conditions
	113	zwi(ji,jj,1) = 0._wp
	114	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
	115	END DO
	116	END DO
	117	#else
[4406]	118	DO jk = 1, jpkfm1 ! Matrix
[3]	119	DO jj = 2, jpjm1
	120	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	121	zcoef = - p2dt / fse3u(ji,jj,jk)
[2528]	122	zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )
	123	zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
	124	zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
	125	zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
	126	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
[3]	127	END DO
	128	END DO
	129	END DO
[2528]	130	DO jj = 2, jpjm1 ! Surface boudary conditions
[3]	131	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	132	zwi(ji,jj,1) = 0._wp
	133	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	134	END DO
	135	END DO
[3211]	136	#endif
[3]	137
	138	! Matrix inversion starting from the first level
	139	!-----------------------------------------------------------------------
	140	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	141	!
	142	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	143	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	144	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	145	! ( ... )( ... ) ( ... )
	146	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	147	!
	148	! m is decomposed in the product of an upper and a lower triangular matrix
	149	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	150	! The solution (the after velocity) is in ua
	151	!-----------------------------------------------------------------------
[2528]	152	!
[3211]	153	#if defined key_z_first
	154	DO jj = 2, jpjm1
	155	DO ji = 2, jpim1
	156	!== Do first and second recurrences in the same loop
	157	ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
	158	& / ( fse3u(ji,jj,1) * rau0 ) )
[4406]	159	DO jk = 2, jpkfm1
[3211]	160	zzwibd = zwi(ji,jj,jk) / zwd(ji,jj,jk-1)
	161	!== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
	162	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zzwibd * zws(ji,jj,jk-1)
	163	!== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
	164	zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side
	165	ua(ji,jj,jk) = zrhs - zzwibd * ua(ji,jj,jk-1)
	166	END DO
	167	!== third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==
[4406]	168	ua(ji,jj,jpkfm1) = ua(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
	169	DO jk = jpkf-2, 1, -1
[3211]	170	ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	171	END DO
	172	! Normalization to obtain the general momentum trend ua
[4406]	173	DO jk = 1, jpkfm1
[3211]	174	ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
	175	END DO
	176	END DO
	177	END DO
	178	#else
[4406]	179	DO jk = 2, jpkfm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	180	DO jj = 2, jpjm1
	181	DO ji = fs_2, fs_jpim1 ! vector opt.
	182	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	183	END DO
	184	END DO
	185	END DO
[2528]	186	!
	187	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	188	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	189	ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
	190	& / ( fse3u(ji,jj,1) * rau0 ) )
[3]	191	END DO
	192	END DO
[4406]	193	DO jk = 2, jpkfm1
[3]	194	DO jj = 2, jpjm1
	195	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	196	zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side
[3]	197	ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
	198	END DO
	199	END DO
	200	END DO
[2528]	201	!
[3211]	202	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==
[3]	203	DO ji = fs_2, fs_jpim1 ! vector opt.
[4406]	204	ua(ji,jj,jpkfm1) = ua(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
[3]	205	END DO
	206	END DO
[4406]	207	DO jk = jpkf-2, 1, -1
[3]	208	DO jj = 2, jpjm1
	209	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	210	ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]	211	END DO
	212	END DO
	213	END DO
	214	! Normalization to obtain the general momentum trend ua
[4406]	215	DO jk = 1, jpkfm1
[3]	216	DO jj = 2, jpjm1
	217	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	218	ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
[3]	219	END DO
	220	END DO
	221	END DO
[3211]	222	#endif
[3]	223
	224	! 2. Vertical diffusion on v
	225	! ---------------------------
	226	! Matrix and second member construction
[1662]	227	! bottom boundary condition: both zwi and zws must be masked as avmv can take
[3]	228	! non zero value at the ocean bottom depending on the bottom friction
[1662]	229	! used but the bottom velocities have already been updated with the bottom
	230	! friction velocity in dyn_bfr using values from the previous timestep. There
	231	! is no need to include these in the implicit calculation.
[2528]	232	!
[3211]	233	#if defined key_z_first
	234	DO jj = 2, jpjm1
	235	DO ji = 2, jpim1
[4406]	236	DO jk = 1, jpkfm1 ! Matrix
[3211]	237	zcoef = -p2dt / fse3v(ji,jj,jk)
	238	zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )
	239	zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk)
	240	zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
	241	zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1)
	242	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
	243	END DO
	244	! Surface boundary conditions
	245	zwi(ji,jj,1) = 0._wp
	246	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
	247	END DO
	248	END DO
	249	#else
[4406]	250	DO jk = 1, jpkfm1 ! Matrix
[3]	251	DO jj = 2, jpjm1
	252	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	253	zcoef = -p2dt / fse3v(ji,jj,jk)
[2528]	254	zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )
[1662]	255	zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk)
[2528]	256	zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
[3]	257	zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1)
[2528]	258	zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
[3]	259	END DO
	260	END DO
	261	END DO
[2528]	262	DO jj = 2, jpjm1 ! Surface boudary conditions
[3]	263	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	264	zwi(ji,jj,1) = 0._wp
	265	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]	266	END DO
	267	END DO
[3211]	268	#endif
[3]	269
	270	! Matrix inversion
	271	!-----------------------------------------------------------------------
	272	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	273	!
	274	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	275	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	276	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	277	! ( ... )( ... ) ( ... )
	278	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	279	!
[2528]	280	! m is decomposed in the product of an upper and lower triangular matrix
[3]	281	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	282	! The solution (after velocity) is in 2d array va
	283	!-----------------------------------------------------------------------
[2528]	284	!
[3211]	285	#if defined key_z_first
	286	DO jj = 2, jpjm1
	287	DO ji = 2, jpim1
	288	!== Do first and second recurrences in the same loop
	289	va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
	290	& / ( fse3v(ji,jj,1) * rau0 ) )
[4406]	291	DO jk = 2, jpkfm1
[3211]	292	zzwibd = zwi(ji,jj,jk) / zwd(ji,jj,jk-1)
	293	!== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
	294	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zzwibd * zws(ji,jj,jk-1)
	295	!== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
	296	zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side
	297	va(ji,jj,jk) = zrhs - zzwibd * va(ji,jj,jk-1)
	298	END DO
	299	!== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==
[4406]	300	va(ji,jj,jpkfm1) = va(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
	301	DO jk = jpkf-2, 1, -1
[3211]	302	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	303	END DO
	304	! Normalization to obtain the general momentum trend va
[4406]	305	DO jk = 1, jpkfm1
[3211]	306	va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
	307	END DO
	308	END DO
	309	END DO
	310	#else
[4406]	311	DO jk = 2, jpkfm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	312	DO jj = 2, jpjm1
	313	DO ji = fs_2, fs_jpim1 ! vector opt.
	314	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	315	END DO
	316	END DO
	317	END DO
[2528]	318	!
	319	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	320	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	321	va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
	322	& / ( fse3v(ji,jj,1) * rau0 ) )
[3]	323	END DO
	324	END DO
[4406]	325	DO jk = 2, jpkfm1
[3]	326	DO jj = 2, jpjm1
	327	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	328	zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side
[3]	329	va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
	330	END DO
	331	END DO
	332	END DO
[2528]	333	!
[3211]	334	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==
[3]	335	DO ji = fs_2, fs_jpim1 ! vector opt.
[4406]	336	va(ji,jj,jpkfm1) = va(ji,jj,jpkfm1) / zwd(ji,jj,jpkfm1)
[3]	337	END DO
	338	END DO
[4406]	339	DO jk = jpkf-2, 1, -1
[3]	340	DO jj = 2, jpjm1
	341	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	342	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]	343	END DO
	344	END DO
	345	END DO
	346
	347	! Normalization to obtain the general momentum trend va
[4406]	348	DO jk = 1, jpkfm1
[3]	349	DO jj = 2, jpjm1
	350	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	351	va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
[3]	352	END DO
	353	END DO
	354	END DO
[3211]	355	#endif
[2528]	356	!
[2715]	357	IF( wrk_not_released(3, 3) ) CALL ctl_stop('dyn_zdf_imp: failed to release workspace array')
	358	!
[3]	359	END SUBROUTINE dyn_zdf_imp
	360
	361	!!==============================================================================
	362	END MODULE dynzdf_imp

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

Download in other formats: