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/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2148

Last change on this file since 2148 was 2148, checked in by cetlod, 14 years ago
merge LOCEAN 2010 developments branches
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 11.3 KB

Rev	Line
[3]	1	MODULE dynzdf_imp
	2	!!==============================================================================
	3	!! * MODULE dynzdf_imp *
	4	!! Ocean dynamics: vertical component(s) of the momentum mixing trend
	5	!!==============================================================================
[2148]	6	!! History : OPA ! 1990-10 (B. Blanke) Original code
	7	!! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal
	8	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	9	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps
[503]	10	!!----------------------------------------------------------------------
[3]	11
	12	!!----------------------------------------------------------------------
	13	!! dyn_zdf_imp : update the momentum trend with the vertical diffu-
	14	!! sion using a implicit time-stepping.
	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and tracers
	17	USE dom_oce ! ocean space and time domain
[888]	18	USE sbc_oce ! surface boundary condition: ocean
	19	USE zdf_oce ! ocean vertical physics
[719]	20	USE phycst ! physical constants
[3]	21	USE in_out_manager ! I/O manager
	22
	23	IMPLICIT NONE
	24	PRIVATE
	25
[2148]	26	PUBLIC dyn_zdf_imp ! called by step.F90
[3]	27
	28	!! * Substitutions
	29	# include "domzgr_substitute.h90"
	30	# include "vectopt_loop_substitute.h90"
	31	!!----------------------------------------------------------------------
[2148]	32	!! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
[888]	33	!! $Id$
[2148]	34	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	35	!!----------------------------------------------------------------------
	36
	37	CONTAINS
	38
[503]	39	SUBROUTINE dyn_zdf_imp( kt, p2dt )
[3]	40	!!----------------------------------------------------------------------
	41	!! * ROUTINE dyn_zdf_imp *
	42	!!
	43	!! ** Purpose : Compute the trend due to the vert. momentum diffusion
	44	!! and the surface forcing, and add it to the general trend of
	45	!! the momentum equations.
	46	!!
	47	!! ** Method : The vertical momentum mixing trend is given by :
	48	!! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
	49	!! backward time stepping
[2148]	50	!! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
[3]	51	!! Bottom boundary conditions : bottom stress (cf zdfbfr.F)
	52	!! Add this trend to the general trend ua :
	53	!! ua = ua + dz( avmu dz(u) )
	54	!!
[2148]	55	!! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
[3]	56	!!---------------------------------------------------------------------
[2148]	57	USE oce, ONLY : zwd => ta ! use ta as workspace
	58	USE oce, ONLY : zws => sa ! use sa as workspace
	59	!!
	60	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	61	REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step
	62	!!
	63	INTEGER :: ji, jj, jk ! dummy loop indices
	64	REAL(wp) :: zrau0r, zcoef ! temporary scalars
	65	REAL(wp) :: zzwi, zzws, zrhs ! temporary scalars
	66	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! 3D workspace
[3]	67	!!----------------------------------------------------------------------
	68
	69	IF( kt == nit000 ) THEN
	70	IF(lwp) WRITE(numout,*)
	71	IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
	72	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
	73	ENDIF
	74
	75	! 0. Local constant initialization
	76	! --------------------------------
	77	zrau0r = 1. / rau0 ! inverse of the reference density
[455]	78
[3]	79	! 1. Vertical diffusion on u
	80	! ---------------------------
	81	! Matrix and second member construction
[1662]	82	! bottom boundary condition: both zwi and zws must be masked as avmu can take
[3]	83	! non zero value at the ocean bottom depending on the bottom friction
[1662]	84	! used but the bottom velocities have already been updated with the bottom
	85	! friction velocity in dyn_bfr using values from the previous timestep. There
	86	! is no need to include these in the implicit calculation.
[2148]	87	!
	88	DO jk = 1, jpkm1 ! Matrix
[3]	89	DO jj = 2, jpjm1
	90	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	91	zcoef = - p2dt / fse3u(ji,jj,jk)
[2148]	92	zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )
	93	zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
	94	zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
	95	zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
[3]	96	zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
	97	END DO
	98	END DO
	99	END DO
[2148]	100	DO jj = 2, jpjm1 ! Surface boudary conditions
[3]	101	DO ji = fs_2, fs_jpim1 ! vector opt.
	102	zwi(ji,jj,1) = 0.
	103	zwd(ji,jj,1) = 1. - zws(ji,jj,1)
	104	END DO
	105	END DO
	106
	107	! Matrix inversion starting from the first level
	108	!-----------------------------------------------------------------------
	109	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	110	!
	111	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	112	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	113	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	114	! ( ... )( ... ) ( ... )
	115	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	116	!
	117	! m is decomposed in the product of an upper and a lower triangular matrix
	118	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	119	! The solution (the after velocity) is in ua
	120	!-----------------------------------------------------------------------
[2148]	121	!
	122	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	123	DO jj = 2, jpjm1
	124	DO ji = fs_2, fs_jpim1 ! vector opt.
	125	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	126	END DO
	127	END DO
	128	END DO
[2148]	129	!
	130	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	131	DO ji = fs_2, fs_jpim1 ! vector opt.
[2148]	132	ua(ji,jj,1) = ub(ji,jj,1) &
	133	& + p2dt * ( ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 ) )
[3]	134	END DO
	135	END DO
	136	DO jk = 2, jpkm1
	137	DO jj = 2, jpjm1
	138	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	139	zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side
[3]	140	ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
	141	END DO
	142	END DO
	143	END DO
[2148]	144	!
	145	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==
[3]	146	DO ji = fs_2, fs_jpim1 ! vector opt.
	147	ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	148	END DO
	149	END DO
	150	DO jk = jpk-2, 1, -1
	151	DO jj = 2, jpjm1
	152	DO ji = fs_2, fs_jpim1 ! vector opt.
	153	ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	154	END DO
	155	END DO
	156	END DO
[2148]	157	!
	158	DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend ua ==
[3]	159	DO jj = 2, jpjm1
	160	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	161	ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / p2dt
[3]	162	END DO
	163	END DO
	164	END DO
	165
	166
	167	! 2. Vertical diffusion on v
	168	! ---------------------------
	169	! Matrix and second member construction
[1662]	170	! bottom boundary condition: both zwi and zws must be masked as avmv can take
[3]	171	! non zero value at the ocean bottom depending on the bottom friction
[1662]	172	! used but the bottom velocities have already been updated with the bottom
	173	! friction velocity in dyn_bfr using values from the previous timestep. There
	174	! is no need to include these in the implicit calculation.
[2148]	175	!
	176	DO jk = 1, jpkm1 ! Matrix
[3]	177	DO jj = 2, jpjm1
	178	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	179	zcoef = -p2dt / fse3v(ji,jj,jk)
[2148]	180	zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )
[1662]	181	zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk)
[2148]	182	zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
[3]	183	zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1)
	184	zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
	185	END DO
	186	END DO
	187	END DO
[2148]	188	DO jj = 2, jpjm1 ! Surface boudary conditions
[3]	189	DO ji = fs_2, fs_jpim1 ! vector opt.
	190	zwi(ji,jj,1) = 0.e0
	191	zwd(ji,jj,1) = 1. - zws(ji,jj,1)
	192	END DO
	193	END DO
	194
	195	! Matrix inversion
	196	!-----------------------------------------------------------------------
	197	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	198	!
	199	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	200	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	201	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	202	! ( ... )( ... ) ( ... )
	203	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	204	!
[2148]	205	! m is decomposed in the product of an upper and lower triangular matrix
[3]	206	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	207	! The solution (after velocity) is in 2d array va
	208	!-----------------------------------------------------------------------
[2148]	209	!
	210	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	211	DO jj = 2, jpjm1
	212	DO ji = fs_2, fs_jpim1 ! vector opt.
	213	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	214	END DO
	215	END DO
	216	END DO
[2148]	217	!
	218	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	219	DO ji = fs_2, fs_jpim1 ! vector opt.
[2148]	220	va(ji,jj,1) = vb(ji,jj,1) &
	221	& + p2dt * ( va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 ) )
[3]	222	END DO
	223	END DO
	224	DO jk = 2, jpkm1
	225	DO jj = 2, jpjm1
	226	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	227	zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side
[3]	228	va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
	229	END DO
	230	END DO
	231	END DO
[2148]	232	!
	233	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==
[3]	234	DO ji = fs_2, fs_jpim1 ! vector opt.
	235	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	236	END DO
	237	END DO
	238	DO jk = jpk-2, 1, -1
	239	DO jj = 2, jpjm1
	240	DO ji = fs_2, fs_jpim1 ! vector opt.
	241	va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	242	END DO
	243	END DO
	244	END DO
[2148]	245	!
	246	DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend va ==
[3]	247	DO jj = 2, jpjm1
	248	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	249	va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / p2dt
[3]	250	END DO
	251	END DO
	252	END DO
[2148]	253	!
[3]	254	END SUBROUTINE dyn_zdf_imp
	255
	256	!!==============================================================================
	257	END MODULE dynzdf_imp

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

Download in other formats: