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.

trazdf_imp.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 915

Last change on this file since 915 was 915, checked in by rblod, 16 years ago
Correct the ldf_slp/trazdf_imp bug in a proper way, see ticket #118
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.2 KB

Rev	Line
[3]	1	MODULE trazdf_imp
	2	!!==============================================================================
[457]	3	!! * MODULE trazdf_imp *
[3]	4	!! Ocean active tracers: vertical component of the tracer mixing trend
	5	!!==============================================================================
[457]	6	!! History :
	7	!! 6.0 ! 90-10 (B. Blanke) Original code
	8	!! 7.0 ! 91-11 (G. Madec)
	9	!! ! 92-06 (M. Imbard) correction on tracer trend loops
	10	!! ! 96-01 (G. Madec) statement function for e3
	11	!! ! 97-05 (G. Madec) vertical component of isopycnal
	12	!! ! 97-07 (G. Madec) geopotential diffusion in s-coord
	13	!! ! 00-08 (G. Madec) double diffusive mixing
	14	!! 8.5 ! 02-08 (G. Madec) F90: Free form and module
[503]	15	!! 9.0 ! 06-11 (G. Madec) New step reorganisation
[3]	16	!!----------------------------------------------------------------------
[457]	17	!! tra_zdf_imp : Update the tracer trend with the diagonal vertical
	18	!! part of the mixing tensor.
[3]	19	!!----------------------------------------------------------------------
	20	!! * Modules used
[457]	21	USE oce ! ocean dynamics and tracers variables
	22	USE dom_oce ! ocean space and time domain variables
	23	USE zdf_oce ! ocean vertical physics variables
[216]	24	USE ldftra_oce ! ocean active tracers: lateral physics
[457]	25	USE ldfslp ! lateral physics: slope of diffusion
[216]	26	USE trdmod ! ocean active tracers trends
	27	USE trdmod_oce ! ocean variables trends
[457]	28	USE zdfddm ! ocean vertical physics: double diffusion
[3]	29	USE in_out_manager ! I/O manager
[457]	30	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	31	USE prtctl ! Print control
[592]	32	USE domvvl ! variable volume
[915]	33	USE ldftra ! lateral mixing type
[3]	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
	38	!! * Routine accessibility
[457]	39	PUBLIC tra_zdf_imp ! routine called by step.F90
[3]	40
	41	!! * Substitutions
	42	# include "domzgr_substitute.h90"
[457]	43	# include "ldftra_substitute.h90"
[3]	44	# include "zdfddm_substitute.h90"
[457]	45	# include "vectopt_loop_substitute.h90"
[3]	46	!!----------------------------------------------------------------------
[719]	47	!!----------------------------------------------------------------------
[457]	48	!! OPA 9.0 , LOCEAN-IPSL (2005)
	49	!!----------------------------------------------------------------------
[3]	50	CONTAINS
[457]	51
	52	SUBROUTINE tra_zdf_imp( kt, p2dt )
[3]	53	!!----------------------------------------------------------------------
	54	!! * ROUTINE tra_zdf_imp *
	55	!!
[457]	56	!! ** Purpose : Compute the trend due to the vertical tracer diffusion
	57	!! including the vertical component of lateral mixing (only for 2nd
	58	!! order operator, for fourth order it is already computed and add
	59	!! to the general trend in traldf.F) and add it to the general trend
	60	!! of the tracer equations.
[3]	61	!!
[457]	62	!! ** Method : The vertical component of the lateral diffusive trends
[503]	63	!! is provided by a 2nd order operator rotated along neutral or geo-
[457]	64	!! potential surfaces to which an eddy induced advection can be
	65	!! added. It is computed using before fields (forward in time) and
	66	!! isopycnal or geopotential slopes computed in routine ldfslp.
	67	!!
	68	!! Second part: vertical trend associated with the vertical physics
	69	!! =========== (including the vertical flux proportional to dk[t]
	70	!! associated with the lateral mixing, through the
	71	!! update of avt)
	72	!! The vertical diffusion of tracers (t & s) is given by:
	73	!! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
	74	!! It is computed using a backward time scheme (t=ta).
[3]	75	!! Surface and bottom boundary conditions: no diffusive flux on
	76	!! both tracers (bottom, applied through the masked field avt).
	77	!! Add this trend to the general trend ta,sa :
[457]	78	!! ta = ta + dz( avt dz(t) )
	79	!! (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
[3]	80	!!
[457]	81	!! Third part: recover avt resulting from the vertical physics
	82	!! ========== alone, for further diagnostics (for example to
	83	!! compute the turbocline depth in zdfmxl.F90).
	84	!! avt = zavt
	85	!! (avs = zavs if lk_zdfddm=T )
[3]	86	!!
[457]	87	!! ** Action : - Update (ta,sa) with before vertical diffusion trend
	88	!!
[3]	89	!!---------------------------------------------------------------------
[457]	90	!! * Modules used
	91	USE oce , ONLY : zwd => ua, & ! ua used as workspace
	92	zws => va ! va " "
[3]	93	!! * Arguments
[457]	94	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	95	REAL(wp), DIMENSION(jpk), INTENT( in ) :: &
	96	p2dt ! vertical profile of tracer time-step
[3]	97
	98	!! * Local declarations
[457]	99	INTEGER :: ji, jj, jk ! dummy loop indices
[592]	100	REAL(wp) :: zavi, zrhs, znvvl, & ! temporary scalars
	101	ze3tb, ze3tn, ze3ta, zvsfvvl ! variable vertical scale factors
[457]	102	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
	103	zwi, zwt, zavsi ! workspace arrays
[3]	104	!!---------------------------------------------------------------------
	105
[457]	106	IF( kt == nit000 ) THEN
	107	IF(lwp)WRITE(numout,*)
[789]	108	IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing'
[457]	109	IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
	110	zavi = 0.e0 ! avoid warning at compilation phase when lk_ldfslp=F
	111	ENDIF
[3]	112
[457]	113	! I. Local initialization
	114	! -----------------------
	115	zwd (1,:, : ) = 0.e0 ; zwd (jpi,:,:) = 0.e0
	116	zws (1,:, : ) = 0.e0 ; zws (jpi,:,:) = 0.e0
	117	zwi (1,:, : ) = 0.e0 ; zwi (jpi,:,:) = 0.e0
	118	zwt (1,:, : ) = 0.e0 ; zwt (jpi,:,:) = 0.e0
	119	zavsi(1,:, : ) = 0.e0 ; zavsi(jpi,:,:) = 0.e0
	120	zwt (:,:,jpk) = 0.e0 ; zwt ( : ,:,1) = 0.e0
	121	zavsi(:,:,jpk) = 0.e0 ; zavsi( : ,:,1) = 0.e0
[3]	122
[592]	123	! I.1 Variable volume : to take into account vertical variable vertical scale factors
	124	! -------------------
	125	IF( lk_vvl ) THEN ; znvvl = 1.
	126	ELSE ; znvvl = 0.e0
	127	ENDIF
	128
[457]	129	! II. Vertical trend associated with the vertical physics
	130	! =======================================================
	131	! (including the vertical flux proportional to dk[t] associated
	132	! with the lateral mixing, through the avt update)
	133	! dk[ avt dk[ (t,s) ] ] diffusive trends
[3]	134
[216]	135
[457]	136	! II.0 Matrix construction
	137	! ------------------------
[3]	138
[457]	139	#if defined key_ldfslp
	140	! update and save of avt (and avs if double diffusive mixing)
[915]	141	IF( l_traldf_rot ) THEN
[902]	142	DO jk = 2, jpkm1
	143	DO jj = 2, jpjm1
	144	DO ji = fs_2, fs_jpim1 ! vector opt.
	145	zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing
	146	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	147	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
	148	zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)
[457]	149	# if defined key_zdfddm
[902]	150	zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! dd mixing: zavsi = total vertical mixing coef. on salinity
[457]	151	# endif
[902]	152	END DO
[3]	153	END DO
	154	END DO
[902]	155	ENDIF
[255]	156	#else
[592]	157	! No isopycnal diffusion
	158	zwt(:,:,:) = avt(:,:,:)
	159	# if defined key_zdfddm
	160	zavsi(:,:,:) = avs(:,:,:)
	161	# endif
	162
	163	#endif
	164
[457]	165	! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)
	166	DO jk = 1, jpkm1
	167	DO jj = 2, jpjm1
	168	DO ji = fs_2, fs_jpim1 ! vector opt.
[592]	169	zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
	170	ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point
	171	ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl ! now scale factor at T-point
	172	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	173	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	174	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
[3]	175	END DO
	176	END DO
[457]	177	END DO
[3]	178
[457]	179	! Surface boudary conditions
	180	DO jj = 2, jpjm1
	181	DO ji = fs_2, fs_jpim1 ! vector opt.
[592]	182	zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
	183	ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point
[457]	184	zwi(ji,jj,1) = 0.e0
[592]	185	zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
[457]	186	END DO
	187	END DO
[3]	188
	189
[457]	190	! II.1. Vertical diffusion on t
	191	! ---------------------------
	192
	193	!! Matrix inversion from the first level
	194	!!----------------------------------------------------------------------
	195	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	196	!
	197	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	198	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	199	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	200	! ( ... )( ... ) ( ... )
	201	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	202	!
	203	! m is decomposed in the product of an upper and lower triangular matrix
	204	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	205	! The second member is in 2d array zwy
	206	! The solution is in 2d array zwx
	207	! The 3d arry zwt is a work space array
	208	! zwy is used and then used as a work space array : its value is modified!
	209
	210	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
	211	DO jj = 2, jpjm1
	212	DO ji = fs_2, fs_jpim1
	213	zwt(ji,jj,1) = zwd(ji,jj,1)
	214	END DO
	215	END DO
	216	DO jk = 2, jpkm1
	217	DO jj = 2, jpjm1
	218	DO ji = fs_2, fs_jpim1
	219	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
[3]	220	END DO
	221	END DO
[457]	222	END DO
[3]	223
[457]	224	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
	225	DO jj = 2, jpjm1
	226	DO ji = fs_2, fs_jpim1
[592]	227	zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )
	228	ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
	229	ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,1)
	230	ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1)
[457]	231	END DO
	232	END DO
	233	DO jk = 2, jpkm1
	234	DO jj = 2, jpjm1
	235	DO ji = fs_2, fs_jpim1
[592]	236	zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )
	237	ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
	238	ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk)
	239	zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk) ! zrhs=right hand side
[457]	240	ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1)
	241	END DO
	242	END DO
	243	END DO
[3]	244
[457]	245	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
	246	! Save the masked temperature after in ta
	247	! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
	248	DO jj = 2, jpjm1
	249	DO ji = fs_2, fs_jpim1
	250	ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
	251	END DO
	252	END DO
	253	DO jk = jpk-2, 1, -1
	254	DO jj = 2, jpjm1
	255	DO ji = fs_2, fs_jpim1
	256	ta(ji,jj,jk) = ( ta(ji,jj,jk) - zws(ji,jj,jk) * ta(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
	257	END DO
	258	END DO
	259	END DO
[3]	260
[457]	261	! II.2 Vertical diffusion on salinity
	262	! -----------------------------------
	263
[3]	264	#if defined key_zdfddm
[457]	265	! Rebuild the Matrix as avt /= avs
[3]	266
[457]	267	! Diagonal, inferior, superior (including the bottom boundary condition via avs masked)
	268	DO jk = 1, jpkm1
	269	DO jj = 2, jpjm1
	270	DO ji = fs_2, fs_jpim1 ! vector opt.
[592]	271	zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
	272	ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point
	273	ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl ! now scale factor at T-point
	274	zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	275	zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	276	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
[3]	277	END DO
	278	END DO
[457]	279	END DO
	280
	281	! Surface boudary conditions
	282	DO jj = 2, jpjm1
	283	DO ji = fs_2, fs_jpim1 ! vector opt.
[592]	284	zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
	285	ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point
[457]	286	zwi(ji,jj,1) = 0.e0
[592]	287	zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
[3]	288	END DO
[457]	289	END DO
[255]	290	#endif
[3]	291
	292
[457]	293	!! Matrix inversion from the first level
	294	!!----------------------------------------------------------------------
	295	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
	296	!
	297	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
	298	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
	299	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
	300	! ( ... )( ... ) ( ... )
	301	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
	302	!
	303	! m is decomposed in the product of an upper and lower triangular
	304	! matrix
	305	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
	306	! The second member is in 2d array zwy
	307	! The solution is in 2d array zwx
	308	! The 3d arry zwt is a work space array
	309	! zwy is used and then used as a work space array : its value is modified!
[3]	310
[457]	311	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
	312	DO jj = 2, jpjm1
	313	DO ji = fs_2, fs_jpim1
	314	zwt(ji,jj,1) = zwd(ji,jj,1)
	315	END DO
	316	END DO
	317	DO jk = 2, jpkm1
	318	DO jj = 2, jpjm1
	319	DO ji = fs_2, fs_jpim1
	320	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
[3]	321	END DO
	322	END DO
[457]	323	END DO
[3]	324
[457]	325	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
	326	DO jj = 2, jpjm1
	327	DO ji = fs_2, fs_jpim1
[592]	328	zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )
	329	ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl ! before scale factor at T-point
	330	ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,1) ! now scale factor at T-point
	331	sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1)
[457]	332	END DO
	333	END DO
	334	DO jk = 2, jpkm1
	335	DO jj = 2, jpjm1
	336	DO ji = fs_2, fs_jpim1
[592]	337	zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )
	338	ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl ! before scale factor at T-point
	339	ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,jk) ! now scale factor at T-point
	340	zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk) ! zrhs=right hand side
[457]	341	sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1)
	342	END DO
	343	END DO
	344	END DO
[3]	345
[457]	346	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
	347	! Save the masked temperature after in ta
	348	! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
	349	DO jj = 2, jpjm1
	350	DO ji = fs_2, fs_jpim1
	351	sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
[216]	352	END DO
[457]	353	END DO
	354	DO jk = jpk-2, 1, -1
	355	DO jj = 2, jpjm1
	356	DO ji = fs_2, fs_jpim1
	357	sa(ji,jj,jk) = ( sa(ji,jj,jk) - zws(ji,jj,jk) * sa(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
	358	END DO
	359	END DO
	360	END DO
[216]	361
[3]	362	END SUBROUTINE tra_zdf_imp
	363
	364	!!==============================================================================
	365	END MODULE trazdf_imp

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

Download in other formats: