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.F90 in NEMO/trunk/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/TRA/trazdf.F90 @ 13291

Last change on this file since 13291 was 13237, checked in by smasson, 4 years ago
trunk: Mid-year merge, merge back KERNEL-06_techene_e3
Property svn:keywords set to `Id`
File size: 12.8 KB

Rev	Line
[458]	1	MODULE trazdf
	2	!!==============================================================================
	3	!! * MODULE trazdf *
	4	!! Ocean active tracers: vertical component of the tracer mixing trend
	5	!!==============================================================================
[9019]	6	!! History : 1.0 ! 2005-11 (G. Madec) Original code
	7	!! 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA
	8	!! 4.0 ! 2017-06 (G. Madec) remove explict time-stepping option
[458]	9	!!----------------------------------------------------------------------
[503]	10
	11	!!----------------------------------------------------------------------
[6140]	12	!! tra_zdf : Update the tracer trend with the vertical diffusion
[458]	13	!!----------------------------------------------------------------------
[6140]	14	USE oce ! ocean dynamics and tracers variables
	15	USE dom_oce ! ocean space and time domain variables
	16	USE domvvl ! variable volume
	17	USE phycst ! physical constant
	18	USE zdf_oce ! ocean vertical physics variables
	19	USE sbc_oce ! surface boundary condition: ocean
	20	USE ldftra ! lateral diffusion: eddy diffusivity
	21	USE ldfslp ! lateral diffusion: iso-neutral slope
	22	USE trd_oce ! trends: ocean variables
	23	USE trdtra ! trends: tracer trend manager
[5836]	24	!
[6140]	25	USE in_out_manager ! I/O manager
	26	USE prtctl ! Print control
	27	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	28	USE lib_mpp ! MPP library
	29	USE timing ! Timing
[592]	30
[458]	31	IMPLICIT NONE
	32	PRIVATE
	33
[9019]	34	PUBLIC tra_zdf ! called by step.F90
	35	PUBLIC tra_zdf_imp ! called by trczdf.F90
[458]	36
	37	!! * Substitutions
[12377]	38	# include "do_loop_substitute.h90"
[13237]	39	# include "domzgr_substitute.h90"
[458]	40	!!----------------------------------------------------------------------
[9598]	41	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[888]	42	!! $Id$
[10068]	43	!! Software governed by the CeCILL license (see ./LICENSE)
[458]	44	!!----------------------------------------------------------------------
[2715]	45	CONTAINS
[458]	46
[12377]	47	SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa )
[458]	48	!!----------------------------------------------------------------------
	49	!! * ROUTINE tra_zdf *
	50	!!
	51	!! ** Purpose : compute the vertical ocean tracer physics.
	52	!!---------------------------------------------------------------------
[12377]	53	INTEGER , INTENT(in) :: kt ! ocean time-step index
	54	INTEGER , INTENT(in) :: Kbb, Kmm, Krhs, Kaa ! time level indices
	55	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
[6140]	56	!
[9019]	57	INTEGER :: jk ! Dummy loop indices
	58	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace
[458]	59	!!---------------------------------------------------------------------
[3294]	60	!
[9019]	61	IF( ln_timing ) CALL timing_start('tra_zdf')
[3294]	62	!
[9124]	63	IF( kt == nit000 ) THEN
	64	IF(lwp)WRITE(numout,*)
	65	IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S'
	66	IF(lwp)WRITE(numout,*) '~~~~~~~ '
	67	ENDIF
	68	!
[9019]	69	IF( l_trdtra ) THEN !* Save ta and sa trends
	70	ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
[12377]	71	ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa)
	72	ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
[458]	73	ENDIF
[6140]	74	!
[9019]	75	! !* compute lateral mixing trend and add it to the general trend
[12489]	76	CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts )
[9019]	77
[5836]	78	!!gm WHY here ! and I don't like that !
[5385]	79	! DRAKKAR SSS control {
	80	! JMM avoid negative salinities near river outlet ! Ugly fix
	81	! JMM : restore negative salinities to small salinities:
[12377]	82	WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp
[5836]	83	!!gm
[458]	84
[1110]	85	IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics
	86	DO jk = 1, jpkm1
[13237]	87	ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) &
	88	& - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) &
	89	& / ( e3t(:,:,jk,Kmm)*rDt ) ) &
	90	& - ztrdt(:,:,jk)
	91	ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) &
	92	& - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) &
	93	& / ( e3t(:,:,jk,Kmm)*rDt ) ) &
	94	& - ztrds(:,:,jk)
[1110]	95	END DO
[5836]	96	!!gm this should be moved in trdtra.F90 and done on all trends
[13226]	97	CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1.0_wp , ztrds, 'T', 1.0_wp )
[5836]	98	!!gm
[12377]	99	CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt )
	100	CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds )
[9019]	101	DEALLOCATE( ztrdt , ztrds )
[1110]	102	ENDIF
	103	! ! print mean trends (used for debugging)
[12377]	104	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf - Ta: ', mask1=tmask, &
	105	& tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
[2715]	106	!
[9019]	107	IF( ln_timing ) CALL timing_stop('tra_zdf')
[3294]	108	!
[458]	109	END SUBROUTINE tra_zdf
	110
[9019]	111
[12377]	112	SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
[458]	113	!!----------------------------------------------------------------------
[9019]	114	!! * ROUTINE tra_zdf_imp *
[458]	115	!!
[9019]	116	!! ** Purpose : Compute the after tracer through a implicit computation
	117	!! of the vertical tracer diffusion (including the vertical component
	118	!! of lateral mixing (only for 2nd order operator, for fourth order
	119	!! it is already computed and add to the general trend in traldf)
[458]	120	!!
[9019]	121	!! ** Method : The vertical diffusion of a tracer ,t , is given by:
	122	!! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
	123	!! It is computed using a backward time scheme (t=after field)
	124	!! which provide directly the after tracer field.
	125	!! If ln_zdfddm=T, use avs for salinity or for passive tracers
	126	!! Surface and bottom boundary conditions: no diffusive flux on
	127	!! both tracers (bottom, applied through the masked field avt).
	128	!! If iso-neutral mixing, add to avt the contribution due to lateral mixing.
	129	!!
[12377]	130	!! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer
[9019]	131	!!---------------------------------------------------------------------
[12377]	132	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	133	INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
	134	INTEGER , INTENT(in ) :: kit000 ! first time step index
	135	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	136	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	137	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
	138	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
[6140]	139	!
[9019]	140	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[10364]	141	REAL(wp) :: zrhs, zzwi, zzws ! local scalars
[9019]	142	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt, zwd, zws
	143	!!---------------------------------------------------------------------
	144	!
	145	! ! ============= !
	146	DO jn = 1, kjpt ! tracer loop !
	147	! ! ============= !
	148	! Matrix construction
	149	! --------------------
	150	! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
	151	!
	152	IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
	153	& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
[6140]	154	!
[9019]	155	! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
	156	IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt(:,:,2:jpk)
	157	ELSE ; zwt(:,:,2:jpk) = avs(:,:,2:jpk)
	158	ENDIF
	159	zwt(:,:,1) = 0._wp
	160	!
	161	IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution
	162	IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
[12377]	163	DO_3D_00_00( 2, jpkm1 )
	164	zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)
	165	END_3D
[9019]	166	ELSE ! standard or triad iso-neutral operator
[12377]	167	DO_3D_00_00( 2, jpkm1 )
	168	zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
	169	END_3D
[9019]	170	ENDIF
	171	ENDIF
	172	!
	173	! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
[10364]	174	IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
[12377]	175	DO_3D_00_00( 1, jpkm1 )
	176	zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm)
	177	zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
	178	zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
	179	& + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
	180	zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp )
	181	zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp )
	182	END_3D
[10364]	183	ELSE
[12377]	184	DO_3D_00_00( 1, jpkm1 )
	185	zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm)
	186	zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
	187	zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
	188	END_3D
[10364]	189	ENDIF
[9019]	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 3d arrays: zwd, zws, zwi.
	203	! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
	204	! and "superior" (above diagonal) components of the tridiagonal system.
	205	! The solution will be in the 4d array pta.
	206	! The 3d array zwt is used as a work space array.
[12377]	207	! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
[9019]	208	! used as a work space array: its value is modified.
	209	!
[12377]	210	DO_2D_00_00
	211	zwt(ji,jj,1) = zwd(ji,jj,1)
	212	END_2D
	213	DO_3D_00_00( 2, jpkm1 )
	214	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
	215	END_3D
[9019]	216	!
	217	ENDIF
	218	!
[12377]	219	DO_2D_00_00
[13237]	220	pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) &
	221	& + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
[12377]	222	END_2D
	223	DO_3D_00_00( 2, jpkm1 )
[13237]	224	zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) &
	225	& + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
[12377]	226	pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
	227	END_3D
[9019]	228	!
[12377]	229	DO_2D_00_00
	230	pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
	231	END_2D
	232	DO_3DS_00_00( jpk-2, 1, -1 )
	233	pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
	234	& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
	235	END_3D
[9019]	236	! ! ================= !
	237	END DO ! end tracer loop !
	238	! ! ================= !
	239	END SUBROUTINE tra_zdf_imp
[458]	240
	241	!!==============================================================================
	242	END MODULE trazdf

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

Download in other formats: