Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

trdmxl_trc.F90

Last change on this file was 15090, checked in by cetlod, 3 years ago
trunk : minor changes in PISCES to make it work in debug mode when nn_hls=2 ( still not working when enabling tiles )
Property svn:keywords set to `Id`
File size: 50.8 KB

Rev	Line
[4804]	1	MODULE trdmxl_trc
	2	!!======================================================================
	3	!! * MODULE trdmxl_trc *
	4	!! Ocean diagnostics: mixed layer passive tracer trends
	5	!!======================================================================
	6	!! History : 9.0 ! 06-08 (C. Deltel) Original code (from trdmxl.F90)
	7	!! ! 07-04 (C. Deltel) Bug fix : add trcrad trends
	8	!! ! 07-06 (C. Deltel) key_gyre : do not call lbc_lnk
	9	!!----------------------------------------------------------------------
[5836]	10	#if defined key_top && defined key_trdmxl_trc
[4804]	11	!!----------------------------------------------------------------------
	12	!! 'key_trdmxl_trc' mixed layer trend diagnostics
	13	!!----------------------------------------------------------------------
	14	!! trd_mxl_trc : passive tracer cumulated trends averaged over ML
	15	!! trd_mxl_trc_zint : passive tracer trends vertical integration
	16	!! trd_mxl_trc_init : initialization step
	17	!!----------------------------------------------------------------------
[12377]	18	USE trc ! tracer definitions (tr etc.)
[4804]	19	USE dom_oce ! domain definition
	20	USE zdfmxl , ONLY : nmln ! number of level in the mixed layer
[9019]	21	USE zdf_oce , ONLY : avs ! vert. diffusivity coef. at w-point for temp
[4804]	22	USE trdtrc_oce ! definition of main arrays used for trends computations
	23	USE in_out_manager ! I/O manager
	24	USE dianam ! build the name of file (routine)
	25	USE ldfslp ! iso-neutral slopes
	26	USE ioipsl ! NetCDF library
	27	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	28	USE lib_mpp ! MPP library
	29	USE trdmxl_trc_rst ! restart for diagnosing the ML trends
	30	USE prtctl ! print control
	31	USE sms_pisces ! PISCES bio-model
	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
	36	PUBLIC trd_mxl_trc
	37	PUBLIC trd_mxl_trc_alloc
	38	PUBLIC trd_mxl_trc_init
	39	PUBLIC trd_mxl_trc_zint
	40
	41	CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file
	42	INTEGER :: nmoymltrd
[7646]	43	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndextrd1, nidtrd, nh_t
[4804]	44	INTEGER :: ndimtrd1
	45	INTEGER, SAVE :: ionce, icount
	46	LOGICAL :: llwarn = .TRUE. ! this should always be .TRUE.
	47	LOGICAL :: lldebug = .TRUE.
	48
	49	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: ztmltrd2 !
	50
[12377]	51	!! * Substitutions
	52	# include "do_loop_substitute.h90"
[13237]	53	# include "domzgr_substitute.h90"
[4804]	54	!!----------------------------------------------------------------------
[10067]	55	!! NEMO/TOP 4.0 , NEMO Consortium (2018)
[5215]	56	!! $Id$
[10068]	57	!! Software governed by the CeCILL license (see ./LICENSE)
[4804]	58	!!----------------------------------------------------------------------
	59	CONTAINS
	60
	61	INTEGER FUNCTION trd_mxl_trc_alloc()
	62	!!----------------------------------------------------------------------
	63	!! * ROUTINE trd_mxl_trc_alloc *
	64	!!----------------------------------------------------------------------
	65	ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) , &
[7646]	66	& ndextrd1(jpi*jpj), nidtrd(jptra), nh_t(jptra), STAT=trd_mxl_trc_alloc)
[4804]	67	!
[10425]	68	CALL mpp_sum ( 'trdmxl_trc', trd_mxl_trc_alloc )
	69	IF( trd_mxl_trc_alloc /=0 ) CALL ctl_stop( 'STOP', 'trd_mxl_trc_alloc: failed to allocate arrays' )
[4804]	70	!
	71	END FUNCTION trd_mxl_trc_alloc
	72
	73
[12377]	74	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn, Kmm )
[4804]	75	!!----------------------------------------------------------------------
	76	!! * ROUTINE trd_mxl_trc_zint *
	77	!!
	78	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
	79	!! to the subroutine. This vertical average is performed from ocean
	80	!! surface down to a chosen control surface.
	81	!!
	82	!! ** Method/usage :
	83	!! The control surface can be either a mixed layer depth (time varying)
	84	!! or a fixed surface (jk level or bowl).
	85	!! Choose control surface with nctls_trc in namelist NAMTRD :
	86	!! nctls_trc = -2 : use isopycnal surface
	87	!! nctls_trc = -1 : use euphotic layer with light criterion
	88	!! nctls_trc = 0 : use mixed layer with density criterion
	89	!! nctls_trc = 1 : read index from file 'ctlsurf_idx'
	90	!! nctls_trc > 1 : use fixed level surface jk = nctls_trc
	91	!! Note: in the remainder of the routine, the volume between the
	92	!! surface and the control surface is called "mixed-layer"
	93	!!----------------------------------------------------------------------
	94	!!
	95	INTEGER, INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
[12377]	96	INTEGER, INTENT( in ) :: Kmm ! time level index
[4804]	97	CHARACTER(len=2), INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
	98	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: ptrc_trdmxl ! passive tracer trend
	99	!
	100	INTEGER :: ji, jj, jk, isum
[9125]	101	REAL(wp), DIMENSION(jpi,jpj) :: zvlmsk
[4804]	102	!!----------------------------------------------------------------------
	103
	104	! I. Definition of control surface and integration weights
	105	! --------------------------------------------------------
	106
	107	ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
	108	!
	109	tmltrd_trc(:,:,:,:) = 0.e0 ! <<< reset trend arrays to zero
	110
	111	! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
	112	SELECT CASE ( nn_ctls_trc ) ! choice of the control surface
[11536]	113	CASE ( -2 ) ; CALL ctl_stop( 'STOP', 'trdmxl_trc : not ready ' ) ! -> isopycnal surface (see ???)
[4804]	114	CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion
	115	CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl)
	116	CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file
	117	CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
	118	nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level
	119	END SELECT
	120
	121	! ... Compute ndextrd1 and ndimtrd1 ??? role de jpktrd_trc
	122	IF( ionce == 1 ) THEN
	123	!
	124	isum = 0 ; zvlmsk(:,:) = 0.e0
	125
[15090]	126	IF( jpktrd_trc < jpk ) THEN
	127	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
[12377]	128	IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
	129	zvlmsk(ji,jj) = tmask(ji,jj,1)
	130	ELSE
	131	isum = isum + 1
	132	zvlmsk(ji,jj) = 0.e0
	133	ENDIF
	134	END_2D
[4804]	135	ENDIF
	136
	137	IF( isum > 0 ) THEN ! index of ocean points (2D only)
	138	WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
	139	CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
	140	ELSE
	141	CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
	142	ENDIF
	143
	144	ionce = 0 ! no more pass here
	145	!
	146	ENDIF ! ionce == 1
	147
	148	! ... Weights for vertical averaging
	149	wkx_trc(:,:,:) = 0.e0
[15090]	150	DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpktrd_trc ) ! description ???
[12377]	151	IF( jk - nmld_trc(ji,jj) < 0 ) wkx_trc(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
	152	END_3D
[4804]	153
	154	rmld_trc(:,:) = 0.e0
	155	DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc
	156	rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
	157	END DO
	158
	159	DO jk = 1, jpktrd_trc ! compute integration weights
	160	wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
	161	END DO
	162
	163	icount = 0 ! <<< flag = off : control surface & integr. weights
	164	! ! computed only once per time step
	165	ENDIF ONCE_PER_TIME_STEP
	166
	167	! II. Vertical integration of trends in the mixed-layer
	168	! -----------------------------------------------------
	169
	170	SELECT CASE ( ctype )
	171	CASE ( '3D' ) ! mean passive tracer trends in the mixed-layer
	172	DO jk = 1, jpktrd_trc
	173	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)
	174	END DO
	175	CASE ( '2D' ) ! forcing at upper boundary of the mixed-layer
	176	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,1) * wkx_trc(:,:,1) ! non penetrative
	177	END SELECT
	178	!
	179	END SUBROUTINE trd_mxl_trc_zint
	180
	181
[12377]	182	SUBROUTINE trd_mxl_trc( kt, Kmm )
[4804]	183	!!----------------------------------------------------------------------
	184	!! * ROUTINE trd_mxl_trc *
	185	!!
	186	!! ** Purpose : Compute and cumulate the mixed layer trends over an analysis
[6140]	187	!! period, and write NetCDF outputs.
[4804]	188	!!
	189	!! ** Method/usage :
	190	!! The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
	191	!! logical namelist variable) :
	192	!! 1) to explain the difference between initial and final
	193	!! mixed-layer T & S (where initial and final relate to the
	194	!! current analysis window, defined by ntrc_trc in the namelist)
	195	!! 2) to explain the difference between the current and previous
	196	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
	197	!! performed over each analysis window).
	198	!!
	199	!! ** Consistency check :
	200	!! If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
	201	!! entrainment) should be zero, at machine accuracy. Note that in the case
	202	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
	203	!! over the first two analysis windows (except if restart).
[9873]	204	!! N.B. For ORCA2_ICE, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8
[4804]	205	!! for checking residuals.
	206	!! On a NEC-SX5 computer, this typically leads to:
	207	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
	208	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
	209	!!
	210	!! ** Action :
	211	!! At each time step, mixed-layer averaged trends are stored in the
	212	!! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
	213	!! This array is known when trd_mld is called, at the end of the stp subroutine,
	214	!! except for the purely vertical K_z diffusion term, which is embedded in the
	215	!! lateral diffusion trend.
	216	!!
	217	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
	218	!! from the lateral diffusion trend.
	219	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
	220	!! arrays are updated.
	221	!! In III), called only once per analysis window, we compute the total trends,
	222	!! along with the residuals and the Asselin correction terms.
	223	!! In IV), the appropriate trends are written in the trends NetCDF file.
	224	!!
	225	!! References :
	226	!! - Vialard & al.
	227	!! - See NEMO documentation (in preparation)
	228	!!----------------------------------------------------------------------
	229	!
	230	INTEGER, INTENT(in) :: kt ! ocean time-step index
[12377]	231	INTEGER, INTENT(in) :: Kmm ! time level index
[4804]	232	!
	233	INTEGER :: ji, jj, jk, jl, ik, it, itmod, jn
	234	REAL(wp) :: zavt, zfn, zfn2
	235	!
[9125]	236	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmltot ! d(trc)/dt over the anlysis window (incl. Asselin)
	237	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlres ! residual = dh/dt entrainment term
	238	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlatf ! for storage only
	239	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlrad ! for storage only (for trb<0 corr in trcrad)
[4804]	240	!
[9125]	241	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmltot2 ! -+
	242	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlres2 ! \| working arrays to diagnose the trends
	243	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmltrdm2 ! \| associated with the time meaned ML
	244	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlatf2 ! \| passive tracers
	245	REAL(wp), DIMENSION(jpi,jpj,jptra) :: ztmlrad2 ! \| (-> for trb<0 corr in trcrad)
[4804]	246	!
	247	CHARACTER (LEN=10) :: clvar
	248	!!----------------------------------------------------------------------
	249
	250
	251	! ======================================================================
	252	! I. Diagnose the purely vertical (K_z) diffusion trend
	253	! ======================================================================
	254
	255	! ... These terms can be estimated by flux computation at the lower boundary of the ML
	256	! (we compute (-1/h) * K_z * d_z( tracer ))
	257
	258	IF( ln_trcldf_iso ) THEN
	259	!
[9019]	260	DO jn = 1, jptra
[15090]	261	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! description ???
[12377]	262	ik = nmld_trc(ji,jj)
	263	IF( ln_trdtrc(jn) ) &
	264	tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - avs(ji,jj,ik) / e3w(ji,jj,ik,Kmm) * tmask(ji,jj,ik) &
	265	& * ( tr(ji,jj,ik-1,jn,Kmm) - tr(ji,jj,ik,jn,Kmm) ) &
	266	& / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
	267	END_2D
[4804]	268	END DO
	269
	270	DO jn = 1, jptra
	271	! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
	272	IF( ln_trdtrc(jn) ) &
	273	tmltrd_trc(:,:,jpmxl_trc_ldf,jn) = tmltrd_trc(:,:,jpmxl_trc_ldf,jn) - tmltrd_trc(:,:,jpmxl_trc_zdf,jn)
	274
	275	END DO
	276	!
	277	ENDIF
	278
[7646]	279	!!gm Test removed, nothing specific to a configuration should survive out of usrdef modules
	280	!!gm IF ( cn_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
	281	!!gm ! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm.
	282	!!gm ! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
[4804]	283	DO jn = 1, jptra
	284	IF( ln_trdtrc(jn) ) THEN
	285	DO jl = 1, jpltrd_trc
[10425]	286	CALL lbc_lnk( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions
[4804]	287	END DO
	288	ENDIF
	289	END DO
[7646]	290	!!gm ENDIF
	291
[4804]	292	! ======================================================================
	293	! II. Cumulate the trends over the analysis window
	294	! ======================================================================
	295
	296	ztmltrd2(:,:,:,:) = 0.e0 ; ztmltot2(:,:,:) = 0.e0 ! <<< reset arrays to zero
	297	ztmlres2(:,:,:) = 0.e0 ; ztmlatf2(:,:,:) = 0.e0
	298	ztmlrad2(:,:,:) = 0.e0
	299
	300	! II.1 Set before values of vertically averages passive tracers
	301	! -------------------------------------------------------------
	302	IF( kt > nittrc000 ) THEN
	303	DO jn = 1, jptra
	304	IF( ln_trdtrc(jn) ) THEN
	305	tmlb_trc (:,:,jn) = tml_trc (:,:,jn)
	306	tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_atf,jn)
	307	tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_radb,jn)
	308	ENDIF
	309	END DO
	310	ENDIF
	311
	312	! II.2 Vertically averaged passive tracers
	313	! ----------------------------------------
	314	tml_trc(:,:,:) = 0.e0
	315	DO jk = 1, jpktrd_trc ! - 1 ???
	316	DO jn = 1, jptra
	317	IF( ln_trdtrc(jn) ) &
[12377]	318	tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * tr(:,:,jk,jn,Kmm)
[4804]	319	END DO
	320	END DO
	321
	322	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
	323	! ------------------------------------------------------------------------
[12377]	324	IF( kt == nittrc000 + 1 ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1) ???
[4804]	325	!
	326	DO jn = 1, jptra
	327	IF( ln_trdtrc(jn) ) THEN
	328	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
	329	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
	330
	331	tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0 ; tmltrd_atf_sumb_trc (:,:,jn) = 0.e0
	332	tmltrd_rad_sumb_trc (:,:,jn) = 0.e0
	333	ENDIF
	334	END DO
	335
	336	rmldbn_trc(:,:) = rmld_trc(:,:)
	337	!
	338	ENDIF
	339
	340	! II.4 Cumulated trends over the analysis period
	341	! ----------------------------------------------
	342	!
	343	! [ 1rst analysis window ] [ 2nd analysis window ]
	344	!
	345	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
	346	! ntrd 2*ntrd etc.
	347	! 1 2 3 4 =5 e.g. =10
	348	!
	349	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN ! ???
	350	!
	351	nmoymltrd = nmoymltrd + 1
	352
	353
	354	! ... Cumulate over BOTH physical contributions AND over time steps
	355	DO jn = 1, jptra
	356	IF( ln_trdtrc(jn) ) THEN
	357	DO jl = 1, jpltrd_trc
	358	tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
	359	END DO
	360	ENDIF
	361	END DO
	362
	363	DO jn = 1, jptra
	364	IF( ln_trdtrc(jn) ) THEN
	365	! ... Special handling of the Asselin trend
	366	tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
	367	tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
	368
	369	! ... Trends associated with the time mean of the ML passive tracers
	370	tmltrd_sum_trc (:,:,:,jn) = tmltrd_sum_trc (:,:,:,jn) + tmltrd_trc (:,:,:,jn)
	371	tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
	372	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) + tml_trc (:,:,jn)
	373	ENDIF
	374	ENDDO
	375
	376	rmld_sum_trc (:,:) = rmld_sum_trc (:,:) + rmld_trc (:,:)
	377	!
	378	ENDIF
	379
	380	! ======================================================================
	381	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
	382	! ======================================================================
	383
	384	! Convert to appropriate physical units
	385	tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
	386
	387	itmod = kt - nittrc000 + 1
	388	it = kt
	389
	390	MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of nn_trd_trc
	391	!
	392	ztmltot (:,:,:) = 0.e0 ! reset arrays to zero
	393	ztmlres (:,:,:) = 0.e0
	394	ztmltot2(:,:,:) = 0.e0
	395	ztmlres2(:,:,:) = 0.e0
	396
	397	zfn = FLOAT( nmoymltrd ) ; zfn2 = zfn * zfn
	398
	399	! III.1 Prepare fields for output ("instantaneous" diagnostics)
	400	! -------------------------------------------------------------
	401
	402	DO jn = 1, jptra
	403	IF( ln_trdtrc(jn) ) THEN
[12489]	404	!-- Compute total trends
[4804]	405	IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes
[12489]	406	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rn_Dt
[4804]	407	ELSE ! LEAP-FROG schemes
[12489]	408	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rn_Dt)
[4804]	409	ENDIF
	410
	411	!-- Compute residuals
	412	ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
	413	& - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) )
	414
	415	!-- Diagnose Asselin trend over the analysis window
	416	ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
	417	ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
	418
	419	!-- Lateral boundary conditions
[7646]	420	IF ( cn_cfg .NE. 'gyre' ) THEN
[14433]	421	CALL lbc_lnk( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1. , ztmlres(:,:,jn) , 'T', 1., &
	422	& ztmlatf(:,:,jn) , 'T', 1. , ztmlrad(:,:,jn) , 'T', 1. )
[4804]	423	ENDIF
	424
	425
	426	#if defined key_diainstant
[11536]	427	CALL ctl_stop( 'STOP', 'tmltrd_trc : key_diainstant was never checked within trdmxl. Comment this to proceed.' )
[4804]	428	#endif
	429	ENDIF
	430	END DO
	431
	432	! III.2 Prepare fields for output ("mean" diagnostics)
	433	! ----------------------------------------------------
	434
	435	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
	436	rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
	437
	438	!-- Compute passive tracer total trends
	439	DO jn = 1, jptra
	440	IF( ln_trdtrc(jn) ) THEN
	441	tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
[12489]	442	ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rn_Dt ) ! now tracer unit is /sec
[4804]	443	ENDIF
	444	END DO
	445
	446	!-- Compute passive tracer residuals
	447	DO jn = 1, jptra
	448	IF( ln_trdtrc(jn) ) THEN
	449	!
	450	DO jl = 1, jpltrd_trc
	451	ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
	452	END DO
	453
	454	ztmltrdm2(:,:,jn) = 0.e0
	455	DO jl = 1, jpltrd_trc
	456	ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
	457	END DO
	458
	459	ztmlres2(:,:,jn) = ztmltot2(:,:,jn) - &
	460	& ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
	461	& - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
	462	!
	463
	464	!-- Diagnose Asselin trend over the analysis window
	465	ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) &
	466	& + tmltrd_atf_sumb_trc(:,:,jn)
	467	ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) &
	468	& + tmltrd_rad_sumb_trc(:,:,jn)
	469
	470	!-- Lateral boundary conditions
[7646]	471	IF ( cn_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
[14433]	472	CALL lbc_lnk( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1., ztmlres2(:,:,jn), 'T', 1. )
[4804]	473	DO jl = 1, jpltrd_trc
[10425]	474	CALL lbc_lnk( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file
[4804]	475	END DO
	476	ENDIF
	477
	478	ENDIF
	479	END DO
	480
	481	! * Debugging information *
	482	IF( lldebug ) THEN
	483	!
	484	WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
	485	WRITE(numout,*) '~~~~~~~~~~~ '
	486	WRITE(numout,*)
	487	WRITE(numout,*) 'TRC kt = ', kt, ' nmoymltrd = ', nmoymltrd
	488
	489	DO jn = 1, jptra
	490
	491	IF( ln_trdtrc(jn) ) THEN
	492	WRITE(numout, *)
	493	WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
	494
	495	WRITE(numout, *)
	496	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres : ', SUM2D(ztmlres(:,:,jn))
	497	!CD??? PREVOIR: z2d = ztmlres(:,:,jn) ; CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres - : ')
	498
	499	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
	500	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
	501	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot : ', SUM2D(+ztmltot (:,:,jn))
	502	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
	503	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
	504	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
	505	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
	506	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
	507	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
	508
	509	WRITE(numout, *)
	510	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2 : ', SUM2D(ztmlres2(:,:,jn))
	511	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
	512	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
	513	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2 : ', SUM2D(+ztmltot2(:,:,jn))
	514	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2 : ', SUM2D(+ztmltrdm2(:,:,jn))
	515	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn))
	516	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
	517	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
	518	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
	519	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
	520
	521	WRITE(numout, *)
	522	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot : ', SUM2D(ztmltot (:,:,jn))
	523	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
	524	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc : ', SUM2D(tml_trc (:,:,jn))
	525	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc : ', SUM2D(tmlbn_trc (:,:,jn))
	526	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc : ', SUM2D(tmlb_trc (:,:,jn))
	527	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc : ', SUM2D(tmlbb_trc (:,:,jn))
	528	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
	529
	530	WRITE(numout, *)
	531	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
	532	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
	533	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
	534	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
	535	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
	536
	537	WRITE(numout, *)
	538	DO jl = 1, jpltrd_trc
	539	WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
	540	& ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
	541	END DO
	542
	543	WRITE(numout,*)
	544	WRITE(numout,) ' ******************** ZTMLRES, ZTMLRES2 ********************* '
	545	WRITE(numout,*)
	546	WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
	547	WRITE(numout,*)
	548	WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
	549
	550	WRITE(numout,*)
	551	WRITE(numout,) ' ******************** ZTMLRES ********************* '
	552	WRITE(numout,*)
	553
	554	WRITE(numout,*) '...................................................'
	555	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
	556	DO jj = 5, 1, -1
	557	WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
	558	END DO
	559
	560	WRITE(numout,*)
	561	WRITE(numout,) ' ******************** ZTMLRES2 ********************* '
	562	WRITE(numout,*)
	563
	564	WRITE(numout,*) '...................................................'
	565	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
	566	DO jj = 5, 1, -1
	567	WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
	568	END DO
	569	!
	570	ENDIF
	571	!
	572	END DO
	573
	574
	575	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
	576	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
	577	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
	578	WRITE(numout,*)
	579	!
	580	ENDIF
	581
	582	! III.3 Time evolution array swap
	583	! -------------------------------
	584	! ML depth
	585	rmldbn_trc(:,:) = rmld_trc(:,:)
	586	rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
	587	DO jn = 1, jptra
	588	IF( ln_trdtrc(jn) ) THEN
	589	! For passive tracer instantaneous diagnostics
	590	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
	591	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
	592
	593	! For passive tracer mean diagnostics
	594	tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
	595	tml_sumb_trc (:,:,jn) = tml_sum_trc(:,:,jn)
	596	tmltrd_atf_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
	597	tmltrd_rad_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
	598
	599
	600	! III.4 Convert to appropriate physical units
	601	! -------------------------------------------
	602	ztmltot (:,:,jn) = ztmltot (:,:,jn) * rn_ucf_trc/zfn ! instant diags
	603	ztmlres (:,:,jn) = ztmlres (:,:,jn) * rn_ucf_trc/zfn
	604	ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * rn_ucf_trc/zfn
	605	ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * rn_ucf_trc/zfn
	606	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags
	607	ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * rn_ucf_trc/zfn2
	608	ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * rn_ucf_trc/zfn2
	609	ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * rn_ucf_trc/zfn2
	610	ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * rn_ucf_trc/zfn2
	611	ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * rn_ucf_trc/zfn2
	612	ENDIF
	613	END DO
	614	!
	615	ENDIF MODULO_NTRD
	616
	617	! ======================================================================
	618	! IV. Write trends in the NetCDF file
	619	! ======================================================================
	620
[6140]	621	! IV.1 Code for IOIPSL/NetCDF output
[4804]	622	! ----------------------------------
	623
	624	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
	625	WRITE(numout,*) ' '
	626	WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
	627	WRITE(numout,*) '~~~~~~~~~~~ '
	628	WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt
	629	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
	630	WRITE(numout,*) ' '
	631	ENDIF
	632
	633	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
	634	!
	635
	636	DO jn = 1, jptra
	637	!
	638	IF( ln_trdtrc(jn) ) THEN
	639	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
	640	!-- Output the fields
	641	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
	642	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 )
	643	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 )
	644	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 )
	645
	646	DO jl = 1, jpltrd_trc - 2
	647	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
	648	& it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
	649	END DO
	650
	651	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
	652	& it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
	653
	654	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
	655	& it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
	656
	657	ENDIF
	658	END DO
	659
	660	IF( kt == nitend ) THEN
	661	DO jn = 1, jptra
	662	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
	663	END DO
	664	ENDIF
	665
	666	ELSE ! <<< write the trends for passive tracer mean diagnostics
	667
	668	DO jn = 1, jptra
	669	!
	670	IF( ln_trdtrc(jn) ) THEN
	671	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 )
	672	!-- Output the fields
	673	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
	674
	675	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
	676	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 )
	677	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 )
	678
	679	DO jl = 1, jpltrd_trc - 2
	680	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
	681	& it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
	682	END DO
	683
	684	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
	685	& it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
	686
	687	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
	688	& it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
	689
	690	ENDIF
	691	!
	692	END DO
	693	IF( kt == nitend ) THEN
	694	DO jn = 1, jptra
	695	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
	696	END DO
	697	ENDIF
	698
	699	!
	700	ENDIF NETCDF_OUTPUT
	701
	702	! Compute the control surface (for next time step) : flag = on
	703	icount = 1
	704
	705	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
	706	!
	707	! Reset cumulative arrays to zero
	708	! -------------------------------
	709	nmoymltrd = 0
	710	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
	711	tmlradm_trc (:,:,:) = 0.e0 ; tml_sum_trc (:,:,:) = 0.e0
	712	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0
	713	rmld_sum_trc (:,:) = 0.e0
	714	!
	715	ENDIF
	716
	717	! ======================================================================
	718	! V. Write restart file
	719	! ======================================================================
	720
	721	IF( lrst_trc ) CALL trd_mxl_trc_rst_write( kt ) ! this must be after the array swap above (III.3)
	722	!
	723	END SUBROUTINE trd_mxl_trc
	724
	725	REAL FUNCTION sum2d( ztab )
	726	!!----------------------------------------------------------------------
	727	!! CD ??? prevoir d'utiliser plutot prtctl
	728	!!----------------------------------------------------------------------
	729	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: ztab
	730	!!----------------------------------------------------------------------
	731	sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
	732	END FUNCTION sum2d
	733
	734
	735	SUBROUTINE trd_mxl_trc_init
	736	!!----------------------------------------------------------------------
	737	!! * ROUTINE trd_mxl_init *
	738	!!
	739	!! ** Purpose : computation of vertically integrated T and S budgets
	740	!! from ocean surface down to control surface (NetCDF output)
	741	!!
	742	!! ** Method/usage :
	743	!!
	744	!!----------------------------------------------------------------------
	745	INTEGER :: inum ! logical unit
	746	INTEGER :: ilseq, jl, jn, iiter
	747	REAL(wp) :: zjulian, zsto, zout
	748	CHARACTER (LEN=40) :: clop
	749	CHARACTER (LEN=15) :: csuff
	750	CHARACTER (LEN=12) :: clmxl
	751	CHARACTER (LEN=16) :: cltrcu
	752	CHARACTER (LEN=10) :: clvar
	753
	754	!!----------------------------------------------------------------------
	755
	756	! ======================================================================
	757	! I. initialization
	758	! ======================================================================
	759
	760	IF(lwp) THEN
	761	WRITE(numout,*)
	762	WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers '
	763	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
	764	WRITE(numout,*)
	765	ENDIF
	766
	767
	768	! I.1 Check consistency of user defined preferences
	769	! -------------------------------------------------
	770
	771	IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
[10425]	772	WRITE(ctmp1,*) ' Your nitend parameter, nitend = ', nitend
	773	WRITE(ctmp2,*) ' is no multiple of the trends diagnostics frequency '
	774	WRITE(ctmp3,*) ' you defined, nn_trd_trc = ', nn_trd_trc
	775	WRITE(ctmp4,*) ' This will not allow you to restart from this simulation. '
	776	WRITE(ctmp5,*) ' You should reconsider this choice. '
	777	WRITE(ctmp6,*)
	778	WRITE(ctmp7,*) ' N.B. the nitend parameter is also constrained to be a '
	779	WRITE(ctmp8,*) ' multiple of the sea-ice frequency parameter (typically 5) '
	780	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4, ctmp5, ctmp6, ctmp7, ctmp8 )
[4804]	781	ENDIF
	782
	783	! * Debugging information *
	784	IF( lldebug ) THEN
	785	WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl
	786	WRITE(numout,*) ' ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
	787	ENDIF
	788
	789	IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
[10425]	790	WRITE(ctmp1,*) ' Currently, you can NOT use simultaneously tracer MUSCL '
	791	WRITE(ctmp2,*) ' advection and window averaged diagnostics of ML trends. '
	792	WRITE(ctmp3,*) ' WHY? Everything in trdmxl_trc is coded for leap-frog, and '
	793	WRITE(ctmp4,*) ' MUSCL scheme is Euler forward for passive tracers (note '
	794	WRITE(ctmp5,*) ' that MUSCL is leap-frog for active tracers T/S). '
	795	WRITE(ctmp6,*) ' In particuliar, entrainment trend would be FALSE. However '
	796	WRITE(ctmp7,*) ' this residual is correct for instantaneous ML diagnostics.'
	797	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4, ctmp5, ctmp6, ctmp7 )
[4804]	798	ENDIF
	799
	800	! I.2 Initialize arrays to zero or read a restart file
	801	! ----------------------------------------------------
	802	nmoymltrd = 0
	803
	804	rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst.
	805	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
	806	tmlradm_trc (:,:,:) = 0.e0
	807
	808	tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean
	809	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0
	810
	811	IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
	812	CALL trd_mxl_trc_rst_read
	813	ELSE
	814	tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst.
	815	tmlbn_trc (:,:,:) = 0.e0
	816
	817	tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean
	818	tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0
	819
	820	ENDIF
	821
	822	icount = 1 ; ionce = 1 ! open specifier
	823
	824
	825	! I.3 Read control surface from file ctlsurf_idx
	826	! ----------------------------------------------
	827	IF( nn_ctls_trc == 1 ) THEN
	828	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
	829	READ ( inum ) nbol_trc
	830	CLOSE( inum )
	831	ENDIF
	832
	833	! ======================================================================
	834	! II. netCDF output initialization
	835	! ======================================================================
	836
	837	! clmxl = legend root for netCDF output
	838	IF( nn_ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion
	839	clmxl = 'Mixed Layer '
	840	ELSE IF( nn_ctls_trc == 1 ) THEN ! control surface = read index from file
	841	clmxl = ' Bowl '
	842	ELSE IF( nn_ctls_trc >= 2 ) THEN ! control surface = model level
	843	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
	844	ENDIF
	845
	846	! II.1 Define frequency of output and means
	847	! -----------------------------------------
	848	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
	849	ELSE ; clop = "x" ! no use of the mask value (require less cp time)
	850	ENDIF
	851	# if defined key_diainstant
	852	IF( .NOT. ln_trdmxl_trc_instant ) THEN
[11536]	853	CALL ctl_stop( 'STOP', 'trd_mxl_trc : this was never checked. Comment this line to proceed...' )
[4804]	854	ENDIF
[12489]	855	zsto = nn_trd_trc * rn_Dt
[4804]	856	clop = "inst("//TRIM(clop)//")"
	857	# else
	858	IF( ln_trdmxl_trc_instant ) THEN
[12489]	859	zsto = rn_Dt ! inst. diags : we use IOIPSL time averaging
[4804]	860	ELSE
[12489]	861	zsto = nn_trd_trc * rn_Dt ! mean diags : we DO NOT use any IOIPSL time averaging
[4804]	862	ENDIF
	863	clop = "ave("//TRIM(clop)//")"
	864	# endif
[12489]	865	zout = nn_trd_trc * rn_Dt
[12377]	866	iiter = nittrc000 - 1
[4804]	867
	868	IF(lwp) WRITE (numout,*) ' netCDF initialization'
	869
	870	! II.2 Compute julian date from starting date of the run
	871	! ------------------------------------------------------
[12489]	872	CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
[4804]	873	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
	874	IF(lwp) WRITE(numout,*)' '
	875	IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000 &
	876	& ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday &
	877	& ,'Julian day : ', zjulian
	878
	879	! II.3 Define the T grid trend file (nidtrd)
	880	! ------------------------------------------
	881
	882	!-- Define long and short names for the NetCDF output variables
	883	! ==> choose them according to trdmxl_trc_oce.F90 <==
	884
	885	ctrd_trc(jpmxl_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmxl_trc_xad ,2) = "_xad"
	886	ctrd_trc(jpmxl_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmxl_trc_yad ,2) = "_yad"
	887	ctrd_trc(jpmxl_trc_zad ,1) = " Vertical advection" ; ctrd_trc(jpmxl_trc_zad ,2) = "_zad"
	888	ctrd_trc(jpmxl_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmxl_trc_ldf ,2) = "_ldf"
	889	ctrd_trc(jpmxl_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmxl_trc_zdf ,2) = "_zdf"
	890	ctrd_trc(jpmxl_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmxl_trc_bbl ,2) = "_bbl"
	891	ctrd_trc(jpmxl_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmxl_trc_dmp ,2) = "_dmp"
	892	ctrd_trc(jpmxl_trc_sbc ,1) = " Surface boundary cond." ; ctrd_trc(jpmxl_trc_sbc ,2) = "_sbc"
	893	ctrd_trc(jpmxl_trc_sms, 1) = " Sources minus sinks" ; ctrd_trc(jpmxl_trc_sms ,2) = "_sms"
	894	ctrd_trc(jpmxl_trc_radb ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radb ,2) = "_radb"
	895	ctrd_trc(jpmxl_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radn ,2) = "_radn"
	896	ctrd_trc(jpmxl_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmxl_trc_atf ,2) = "_atf"
	897
	898	DO jn = 1, jptra
	899	!-- Create a NetCDF file and enter the define mode
	900	IF( ln_trdtrc(jn) ) THEN
	901	csuff="ML_"//ctrcnm(jn)
	902	CALL dia_nam( clhstnam, nn_trd_trc, csuff )
	903	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
[12489]	904	& 1, jpi, 1, jpj, iiter, zjulian, rn_Dt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
[4804]	905
	906	!-- Define the ML depth variable
	907	CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m", &
	908	& jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	909
	910	ENDIF
	911	END DO
	912
	913	!-- Define physical units
	914	IF( rn_ucf_trc == 1. ) THEN
	915	cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit
	916	ELSEIF ( rn_ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3)
	917	cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf
	918	ELSE
	919	cltrcu = "unknown?"
	920	ENDIF
	921
	922	!-- Define miscellaneous passive tracer mixed-layer variables
	923	IF( jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
[11536]	924	CALL ctl_stop( 'STOP', 'Error : jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb' ) ! see below
[4804]	925	ENDIF
	926
	927	DO jn = 1, jptra
	928	!
	929	IF( ln_trdtrc(jn) ) THEN
	930	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc.
	931	CALL histdef(nidtrd(jn), trim(clvar), clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", &
	932	& "mmole-N/m3", jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	933	CALL histdef(nidtrd(jn), trim(clvar)//"_tot" , clmxl//" "//trim(ctrcnm(jn))//" Total trend ", &
	934	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
	935	CALL histdef(nidtrd(jn), trim(clvar)//"_res" , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)", &
	936	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
	937
	938	DO jl = 1, jpltrd_trc - 2 ! <== only true if jpltrd_trc == jpmxl_trc_atf
	939	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1), &
	940	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
[12489]	941	END DO ! if zsto=rn_Dt above
[4804]	942
	943	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), &
	944	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
	945
	946	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1), &
	947	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
	948	!
	949	ENDIF
	950	END DO
	951
	952	!-- Leave IOIPSL/NetCDF define mode
	953	DO jn = 1, jptra
	954	IF( ln_trdtrc(jn) ) CALL histend( nidtrd(jn), snc4set )
	955	END DO
	956
	957	IF(lwp) WRITE(numout,*)
	958
	959	END SUBROUTINE trd_mxl_trc_init
	960
	961	#else
	962	!!----------------------------------------------------------------------
	963	!! Default option : Empty module
	964	!!----------------------------------------------------------------------
	965	CONTAINS
[12377]	966	SUBROUTINE trd_mxl_trc( kt, Kmm ) ! Empty routine
[4804]	967	INTEGER, INTENT( in) :: kt
[12377]	968	INTEGER, INTENT( in) :: Kmm ! time level index
[4804]	969	WRITE(,) 'trd_mxl_trc: You should not have seen this print! error?', kt
	970	END SUBROUTINE trd_mxl_trc
[12377]	971	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn, Kmm )
[4804]	972	INTEGER , INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
[12377]	973	INTEGER , INTENT( in ) :: Kmm ! time level index
[4804]	974	CHARACTER(len=2) , INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
	975	REAL, DIMENSION(:,:,:), INTENT( in ) :: ptrc_trdmxl ! passive trc trend
	976	WRITE(,) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
	977	WRITE(,) ' " " : You should not have seen this print! error?', ctype
	978	WRITE(,) ' " " : You should not have seen this print! error?', ktrd
	979	WRITE(,) ' " " : You should not have seen this print! error?', kjn
	980	END SUBROUTINE trd_mxl_trc_zint
	981	SUBROUTINE trd_mxl_trc_init ! Empty routine
	982	WRITE(,) 'trd_mxl_trc_init: You should not have seen this print! error?'
	983	END SUBROUTINE trd_mxl_trc_init
	984	#endif
	985
	986	!!======================================================================
	987	END MODULE trdmxl_trc

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/TOP/TRP/trdmxl_trc.F90

Download in other formats: