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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago
update licence of all NEMO files...
Property svn:keywords set to `Id`
File size: 14.7 KB

Rev	Line
[3]	1	MODULE trazdf_imp
[1438]	2	!!======================================================================
[457]	3	!! * MODULE trazdf_imp *
[2024]	4	!! Ocean 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
[2024]	16	!! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC
[3]	17	!!----------------------------------------------------------------------
[1438]	18
	19	!!----------------------------------------------------------------------
[457]	20	!! tra_zdf_imp : Update the tracer trend with the diagonal vertical
	21	!! part of the mixing tensor.
[3]	22	!!----------------------------------------------------------------------
[457]	23	USE oce ! ocean dynamics and tracers variables
	24	USE dom_oce ! ocean space and time domain variables
	25	USE zdf_oce ! ocean vertical physics variables
[216]	26	USE ldftra_oce ! ocean active tracers: lateral physics
[457]	27	USE ldfslp ! lateral physics: slope of diffusion
	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)
[592]	31	USE domvvl ! variable volume
[915]	32	USE ldftra ! lateral mixing type
[2082]	33	USE trc_oce ! share passive tracers/Ocean variables
[3]	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
[1438]	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	!!----------------------------------------------------------------------
[2287]	46	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1156]	47	!! $Id$
[2287]	48	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[457]	49	!!----------------------------------------------------------------------
[3]	50	CONTAINS
[2024]	51
[2034]	52	SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt )
[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)
[2024]	72	!! The vertical diffusion of the tracer t is given by:
[457]	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) )
[2024]	79	!! if lk_zdfddm=T, use avs for salinity or for passive tracers
	80	!! (sa = sa + dz( avs dz(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).
[2024]	85	!! if lk_zdfddm=T, use avt = zavt
[457]	86	!! (avs = zavs if lk_zdfddm=T )
[3]	87	!!
[2024]	88	!! ** Action : - Update (ta) with before vertical diffusion trend
[457]	89	!!
[3]	90	!!---------------------------------------------------------------------
[1438]	91	USE oce , ONLY : zwd => ua ! ua used as workspace
	92	USE oce , ONLY : zws => va ! va - -
[2034]	93	!!
[2104]	94	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	95	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	96	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	97	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	98	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields
	99	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[1438]	100	!!
[2024]	101	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2104]	102	REAL(wp) :: zavi, zrhs, znvvl ! local scalars
[2024]	103	REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors
	104	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt ! workspace arrays
[3]	105	!!---------------------------------------------------------------------
	106
[2104]	107	IF( kt == nit000 ) THEN
[457]	108	IF(lwp)WRITE(numout,*)
[2082]	109	IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
[457]	110	IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
	111	zavi = 0.e0 ! avoid warning at compilation phase when lk_ldfslp=F
	112	ENDIF
[2024]	113	!
[457]	114	! I. Local initialization
	115	! -----------------------
[2024]	116	zwd(1,:, : ) = 0.e0 ; zwd(jpi,:,:) = 0.e0
	117	zws(1,:, : ) = 0.e0 ; zws(jpi,:,:) = 0.e0
	118	zwi(1,:, : ) = 0.e0 ; zwi(jpi,:,:) = 0.e0
	119	zwt(1,:, : ) = 0.e0 ; zwt(jpi,:,:) = 0.e0
	120	zwt(:,:,jpk) = 0.e0 ; zwt( : ,:,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
[2024]	134	!
[457]	135	! II.0 Matrix construction
	136	! ------------------------
[2024]	137	DO jn = 1, kjpt
	138	!
	139	! Matrix construction
	140	! ------------------------
	141	IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
[457]	142	#if defined key_ldfslp
[2236]	143	IF( ln_traldf_grif ) THEN
	144	DO jk = 2, jpkm1
	145	DO jj = 2, jpjm1
	146	DO ji = fs_2, fs_jpim1 ! vector opt.
	147	zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing
	148	zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)
	149	END DO
	150	END DO
	151	END DO
[2024]	152	! update and save of avt (and avs if double diffusive mixing)
[2236]	153	ELSE IF( l_traldf_rot ) THEN
[2024]	154	DO jk = 2, jpkm1
	155	DO jj = 2, jpjm1
	156	DO ji = fs_2, fs_jpim1 ! vector opt.
	157	zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing
	158	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	159	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
	160	zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)
	161	END DO
	162	END DO
	163	END DO
	164	ELSE ! no rotation but key_ldfslp defined
	165	zwt (:,:,:) = avt(:,:,:)
	166	ENDIF
	167	#else
	168	! No isopycnal diffusion
	169	zwt(:,:,:) = avt(:,:,:)
	170	#endif
	171	! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)
	172	DO jk = 1, jpkm1
	173	DO jj = 2, jpjm1
	174	DO ji = fs_2, fs_jpim1 ! vector opt.
	175	ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point
	176	ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point
	177	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	178	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	179	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
	180	END DO
	181	END DO
	182	END DO
[2052]	183	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
	184	DO jj = 2, jpjm1
	185	DO ji = fs_2, fs_jpim1
	186	zwt(ji,jj,1) = zwd(ji,jj,1)
	187	END DO
	188	END DO
	189	DO jk = 2, jpkm1
	190	DO jj = 2, jpjm1
	191	DO ji = fs_2, fs_jpim1
	192	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	193	END DO
	194	END DO
	195	END DO
[2024]	196	!
	197	ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
	198	#if defined key_ldfslp
[2236]	199	IF( ln_traldf_grif ) THEN
[2024]	200	DO jk = 2, jpkm1
	201	DO jj = 2, jpjm1
	202	DO ji = fs_2, fs_jpim1 ! vector opt.
[2236]	203	zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing
	204	zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)
	205	END DO
	206	END DO
	207	END DO
	208	ELSE IF( l_traldf_rot ) THEN
	209	DO jk = 2, jpkm1
	210	DO jj = 2, jpjm1
	211	DO ji = fs_2, fs_jpim1 ! vector opt.
[2024]	212	zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing
	213	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
	214	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
[2052]	215	zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on salinity)
[2024]	216	END DO
	217	END DO
	218	END DO
	219	ELSE ! no rotation but key_ldfslp defined
	220	zwt (:,:,:) = fsavs(:,:,:)
	221	ENDIF
[255]	222	#else
[2024]	223	! No isopycnal diffusion
	224	zwt(:,:,:) = fsavs(:,:,:)
[592]	225	#endif
[2024]	226	! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)
	227	DO jk = 1, jpkm1
	228	DO jj = 2, jpjm1
	229	DO ji = fs_2, fs_jpim1 ! vector opt.
	230	ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point
	231	ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point
	232	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
	233	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
	234	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
	235	END DO
	236	END DO
[3]	237	END DO
[2024]	238	! Surface boudary conditions
	239	DO jj = 2, jpjm1
	240	DO ji = fs_2, fs_jpim1 ! vector opt.
	241	ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) ! after scale factor at T-point
	242	zwi(ji,jj,1) = 0.e0
	243	zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
	244	END DO
	245	END DO
	246	!
[2052]	247	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
[2024]	248	DO jj = 2, jpjm1
	249	DO ji = fs_2, fs_jpim1
[2052]	250	zwt(ji,jj,1) = zwd(ji,jj,1)
[2024]	251	END DO
[457]	252	END DO
[2052]	253	DO jk = 2, jpkm1
	254	DO jj = 2, jpjm1
	255	DO ji = fs_2, fs_jpim1
	256	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	257	END DO
	258	END DO
	259	END DO
	260	!
	261	END IF
	262
	263	! II.1. Vertical diffusion on tracer
	264	! ---------------------------
[2024]	265
	266	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
[457]	267	DO jj = 2, jpjm1
	268	DO ji = fs_2, fs_jpim1
[2024]	269	ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
	270	ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
[2034]	271	pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
[457]	272	END DO
	273	END DO
[2024]	274	DO jk = 2, jpkm1
	275	DO jj = 2, jpjm1
	276	DO ji = fs_2, fs_jpim1
	277	ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
	278	ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk)
[2034]	279	zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side
	280	pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
[2024]	281	END DO
[3]	282	END DO
	283	END DO
[2236]	284
	285	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
[457]	286	DO jj = 2, jpjm1
	287	DO ji = fs_2, fs_jpim1
[2034]	288	pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
[3]	289	END DO
	290	END DO
[2024]	291	DO jk = jpk-2, 1, -1
	292	DO jj = 2, jpjm1
	293	DO ji = fs_2, fs_jpim1
[2034]	294	pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
[2104]	295	& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
[2024]	296	END DO
[457]	297	END DO
	298	END DO
[2024]	299	!
[457]	300	END DO
[1438]	301	!
[3]	302	END SUBROUTINE tra_zdf_imp
	303
	304	!!==============================================================================
	305	END MODULE trazdf_imp

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

Download in other formats: