Context Navigation

dynzdf_imp.F90 @ 2243

Last change on this file since 2243 was 2236, checked in by cetlod, 14 years ago
First guess of NEMO_v3.3
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 12.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
[2236]	22	#if defined key_zdfgls
	23	USE zdfbfr, ONLY : bfrua, bfrva, wbotu, wbotv ! bottom stresses
	24	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	25	#endif
[3]	26
	27	IMPLICIT NONE
	28	PRIVATE
	29
[2148]	30	PUBLIC dyn_zdf_imp ! called by step.F90
[3]	31
	32	!! * Substitutions
	33	# include "domzgr_substitute.h90"
	34	# include "vectopt_loop_substitute.h90"
	35	!!----------------------------------------------------------------------
[2148]	36	!! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
[888]	37	!! $Id$
[2148]	38	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	39	!!----------------------------------------------------------------------
	40
	41	CONTAINS
	42
[503]	43	SUBROUTINE dyn_zdf_imp( kt, p2dt )
[3]	44	!!----------------------------------------------------------------------
	45	!! * ROUTINE dyn_zdf_imp *
	46	!!
	47	!! ** Purpose : Compute the trend due to the vert. momentum diffusion
	48	!! and the surface forcing, and add it to the general trend of
	49	!! the momentum equations.
	50	!!
	51	!! ** Method : The vertical momentum mixing trend is given by :
	52	!! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
	53	!! backward time stepping
[2148]	54	!! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
[3]	55	!! Bottom boundary conditions : bottom stress (cf zdfbfr.F)
	56	!! Add this trend to the general trend ua :
	57	!! ua = ua + dz( avmu dz(u) )
	58	!!
[2148]	59	!! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
[3]	60	!!---------------------------------------------------------------------
[2148]	61	USE oce, ONLY : zwd => ta ! use ta as workspace
	62	USE oce, ONLY : zws => sa ! use sa as workspace
	63	!!
	64	INTEGER , INTENT( in ) :: kt ! ocean time-step index
	65	REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step
	66	!!
	67	INTEGER :: ji, jj, jk ! dummy loop indices
	68	REAL(wp) :: zrau0r, zcoef ! temporary scalars
	69	REAL(wp) :: zzwi, zzws, zrhs ! temporary scalars
	70	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! 3D workspace
[2236]	71	#if defined key_zdfgls
	72	INTEGER :: ikbu, ikbv, ikbum1, ikbvm1
	73	REAL(wp) :: zcbcu, zcbcv
	74	#endif
[3]	75	!!----------------------------------------------------------------------
	76
	77	IF( kt == nit000 ) THEN
	78	IF(lwp) WRITE(numout,*)
	79	IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
	80	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
	81	ENDIF
	82
	83	! 0. Local constant initialization
	84	! --------------------------------
	85	zrau0r = 1. / rau0 ! inverse of the reference density
[455]	86
[3]	87	! 1. Vertical diffusion on u
	88	! ---------------------------
	89	! Matrix and second member construction
[1662]	90	! bottom boundary condition: both zwi and zws must be masked as avmu can take
[3]	91	! non zero value at the ocean bottom depending on the bottom friction
[1662]	92	! used but the bottom velocities have already been updated with the bottom
	93	! friction velocity in dyn_bfr using values from the previous timestep. There
	94	! is no need to include these in the implicit calculation.
[2148]	95	!
	96	DO jk = 1, jpkm1 ! Matrix
[3]	97	DO jj = 2, jpjm1
	98	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	99	zcoef = - p2dt / fse3u(ji,jj,jk)
[2148]	100	zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )
	101	zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)
	102	zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
	103	zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)
[3]	104	zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
	105	END DO
	106	END DO
	107	END DO
[2148]	108	DO jj = 2, jpjm1 ! Surface boudary conditions
[3]	109	DO ji = fs_2, fs_jpim1 ! vector opt.
	110	zwi(ji,jj,1) = 0.
	111	zwd(ji,jj,1) = 1. - zws(ji,jj,1)
	112	END DO
	113	END DO
	114
	115	! Matrix inversion starting from the first level
	116	!-----------------------------------------------------------------------
	117	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	118	!
	119	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	120	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	121	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	122	! ( ... )( ... ) ( ... )
	123	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	124	!
	125	! m is decomposed in the product of an upper and a lower triangular matrix
	126	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	127	! The solution (the after velocity) is in ua
	128	!-----------------------------------------------------------------------
[2148]	129	!
	130	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	131	DO jj = 2, jpjm1
	132	DO ji = fs_2, fs_jpim1 ! vector opt.
	133	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	134	END DO
	135	END DO
	136	END DO
[2148]	137	!
	138	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	139	DO ji = fs_2, fs_jpim1 ! vector opt.
[2148]	140	ua(ji,jj,1) = ub(ji,jj,1) &
	141	& + p2dt * ( ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 ) )
[3]	142	END DO
	143	END DO
	144	DO jk = 2, jpkm1
	145	DO jj = 2, jpjm1
	146	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	147	zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side
[3]	148	ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
	149	END DO
	150	END DO
	151	END DO
[2148]	152	!
	153	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==
[3]	154	DO ji = fs_2, fs_jpim1 ! vector opt.
	155	ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	156	END DO
	157	END DO
	158	DO jk = jpk-2, 1, -1
	159	DO jj = 2, jpjm1
	160	DO ji = fs_2, fs_jpim1 ! vector opt.
	161	ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	162	END DO
	163	END DO
	164	END DO
[2236]	165
	166	#if defined key_zdfgls
	167	! Save bottom stress for next time step
	168	DO jj = 2, jpjm1
	169	DO ji = fs_2, fs_jpim1 ! vector opt.
	170	ikbu = MIN( mbathy(ji+1,jj ), mbathy(ji,jj) )
	171	ikbum1 = MAX( ikbu-1, 1 )
	172	wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,ikbum1) * umask(ji,jj,ikbum1)
	173	END DO
	174	END DO
	175	CALL lbc_lnk( wbotu(:,:), 'U', -1. )
	176	#endif
	177
	178	! Normalization to obtain the general momentum trend ua
	179	DO jk = 1, jpkm1
[3]	180	DO jj = 2, jpjm1
	181	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	182	ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / p2dt
[3]	183	END DO
	184	END DO
	185	END DO
	186
	187
	188	! 2. Vertical diffusion on v
	189	! ---------------------------
	190	! Matrix and second member construction
[1662]	191	! bottom boundary condition: both zwi and zws must be masked as avmv can take
[3]	192	! non zero value at the ocean bottom depending on the bottom friction
[1662]	193	! used but the bottom velocities have already been updated with the bottom
	194	! friction velocity in dyn_bfr using values from the previous timestep. There
	195	! is no need to include these in the implicit calculation.
[2148]	196	!
	197	DO jk = 1, jpkm1 ! Matrix
[3]	198	DO jj = 2, jpjm1
	199	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	200	zcoef = -p2dt / fse3v(ji,jj,jk)
[2148]	201	zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )
[1662]	202	zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk)
[2148]	203	zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
[3]	204	zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1)
	205	zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
	206	END DO
	207	END DO
	208	END DO
[2148]	209	DO jj = 2, jpjm1 ! Surface boudary conditions
[3]	210	DO ji = fs_2, fs_jpim1 ! vector opt.
	211	zwi(ji,jj,1) = 0.e0
	212	zwd(ji,jj,1) = 1. - zws(ji,jj,1)
	213	END DO
	214	END DO
	215
	216	! Matrix inversion
	217	!-----------------------------------------------------------------------
	218	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	219	!
	220	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	221	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	222	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	223	! ( ... )( ... ) ( ... )
	224	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	225	!
[2148]	226	! m is decomposed in the product of an upper and lower triangular matrix
[3]	227	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	228	! The solution (after velocity) is in 2d array va
	229	!-----------------------------------------------------------------------
[2148]	230	!
	231	DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
[3]	232	DO jj = 2, jpjm1
	233	DO ji = fs_2, fs_jpim1 ! vector opt.
	234	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
	235	END DO
	236	END DO
	237	END DO
[2148]	238	!
	239	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
[3]	240	DO ji = fs_2, fs_jpim1 ! vector opt.
[2148]	241	va(ji,jj,1) = vb(ji,jj,1) &
	242	& + p2dt * ( va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 ) )
[3]	243	END DO
	244	END DO
	245	DO jk = 2, jpkm1
	246	DO jj = 2, jpjm1
	247	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	248	zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side
[3]	249	va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
	250	END DO
	251	END DO
	252	END DO
[2148]	253	!
	254	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==
[3]	255	DO ji = fs_2, fs_jpim1 ! vector opt.
	256	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
	257	END DO
	258	END DO
	259	DO jk = jpk-2, 1, -1
	260	DO jj = 2, jpjm1
	261	DO ji = fs_2, fs_jpim1 ! vector opt.
	262	va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
	263	END DO
	264	END DO
	265	END DO
[2236]	266
	267	#if defined key_zdfgls
	268	! Save bottom stress for next time step
	269	DO jj = 2, jpjm1
	270	DO ji = fs_2, fs_jpim1 ! vector opt.
	271	ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
	272	ikbvm1 = MAX( ikbv-1, 1 )
	273	wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1)
	274	END DO
	275	END DO
	276	CALL lbc_lnk( wbotv(:,:), 'V', -1. )
	277	#endif
	278
	279	! Normalization to obtain the general momentum trend va
	280	DO jk = 1, jpkm1
[3]	281	DO jj = 2, jpjm1
	282	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	283	va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / p2dt
[3]	284	END DO
	285	END DO
	286	END DO
[2148]	287	!
[3]	288	END SUBROUTINE dyn_zdf_imp
	289
	290	!!==============================================================================
	291	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/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2243

Download in other formats: