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

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

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

Download in other formats: