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

Last change on this file since 1517 was 1517, checked in by ctlod, 15 years ago
see tickets: #363 & #433
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 16.1 KB

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

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

Download in other formats: