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.

trdmxl_trc.F90 in branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 6051

Last change on this file since 6051 was 6051, checked in by lovato, 8 years ago
Merge branches/2015/dev_r5056_CMCC4_simplification (see ticket #1456)
Property svn:keywords set to `Id`
File size: 68.0 KB

Line
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	!!----------------------------------------------------------------------
10	#if defined key_top && defined key_trdmxl_trc
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	!!----------------------------------------------------------------------
18	USE trc ! tracer definitions (trn, trb, tra, etc.)
19	USE trc_oce, ONLY : nn_dttrc ! frequency of step on passive tracers
20	USE dom_oce ! domain definition
21	USE zdfmxl , ONLY : nmln ! number of level in the mixed layer
22	USE zdf_oce , ONLY : avt ! vert. diffusivity coef. at w-point for temp
23	# if defined key_zdfddm
24	USE zdfddm , ONLY : avs ! salinity vertical diffusivity coeff. at w-point
25	# endif
26	USE trdtrc_oce ! definition of main arrays used for trends computations
27	USE in_out_manager ! I/O manager
28	USE dianam ! build the name of file (routine)
29	USE ldfslp ! iso-neutral slopes
30	USE ioipsl ! NetCDF library
31	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
32	USE lib_mpp ! MPP library
33	USE trdmxl_trc_rst ! restart for diagnosing the ML trends
34	USE prtctl ! print control
35	USE sms_pisces ! PISCES bio-model
36	USE wrk_nemo ! Memory allocation
37
38	IMPLICIT NONE
39	PRIVATE
40
41	PUBLIC trd_mxl_trc
42	PUBLIC trd_mxl_trc_alloc
43	PUBLIC trd_mxl_bio
44	PUBLIC trd_mxl_trc_init
45	PUBLIC trd_mxl_trc_zint
46	PUBLIC trd_mxl_bio_zint
47
48	CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file
49	INTEGER :: nmoymltrd
50	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndextrd1
51	INTEGER, DIMENSION(jptra) :: nidtrd, nh_t
52	INTEGER :: ndimtrd1
53	INTEGER, SAVE :: ionce, icount
54	#if defined key_pisces_reduced
55	INTEGER :: nidtrdbio, nh_tb
56	INTEGER, SAVE :: ioncebio, icountbio
57	INTEGER, SAVE :: nmoymltrdbio
58	#endif
59	LOGICAL :: llwarn = .TRUE. ! this should always be .TRUE.
60	LOGICAL :: lldebug = .TRUE.
61
62	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: ztmltrd2 !
63	#if defined key_pisces_reduced
64	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ztmltrdbio2 ! only needed for mean diagnostics in trd_mxl_bio()
65	#endif
66
67	!! * Substitutions
68	# include "domzgr_substitute.h90"
69	# include "zdfddm_substitute.h90"
70	!!----------------------------------------------------------------------
71	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
72	!! $Id$
73	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
74	!!----------------------------------------------------------------------
75	CONTAINS
76
77	INTEGER FUNCTION trd_mxl_trc_alloc()
78	!!----------------------------------------------------------------------
79	!! * ROUTINE trd_mxl_trc_alloc *
80	!!----------------------------------------------------------------------
81	ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) , &
82	#if defined key_pisces_reduced
83	& ztmltrdbio2(jpi,jpj,jpdiabio) , &
84	#endif
85	& ndextrd1(jpi*jpj) , STAT=trd_mxl_trc_alloc)
86	!
87	IF( lk_mpp ) CALL mpp_sum ( trd_mxl_trc_alloc )
88	IF( trd_mxl_trc_alloc /=0 ) CALL ctl_warn('trd_mxl_trc_alloc: failed to allocate arrays')
89	!
90	END FUNCTION trd_mxl_trc_alloc
91
92
93	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
94	!!----------------------------------------------------------------------
95	!! * ROUTINE trd_mxl_trc_zint *
96	!!
97	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
98	!! to the subroutine. This vertical average is performed from ocean
99	!! surface down to a chosen control surface.
100	!!
101	!! ** Method/usage :
102	!! The control surface can be either a mixed layer depth (time varying)
103	!! or a fixed surface (jk level or bowl).
104	!! Choose control surface with nctls_trc in namelist NAMTRD :
105	!! nctls_trc = -2 : use isopycnal surface
106	!! nctls_trc = -1 : use euphotic layer with light criterion
107	!! nctls_trc = 0 : use mixed layer with density criterion
108	!! nctls_trc = 1 : read index from file 'ctlsurf_idx'
109	!! nctls_trc > 1 : use fixed level surface jk = nctls_trc
110	!! Note: in the remainder of the routine, the volume between the
111	!! surface and the control surface is called "mixed-layer"
112	!!----------------------------------------------------------------------
113	!!
114	INTEGER, INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
115	CHARACTER(len=2), INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
116	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: ptrc_trdmxl ! passive tracer trend
117	!
118	INTEGER :: ji, jj, jk, isum
119	REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
120	!!----------------------------------------------------------------------
121
122	CALL wrk_alloc( jpi, jpj, zvlmsk )
123
124	! I. Definition of control surface and integration weights
125	! --------------------------------------------------------
126
127	ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
128	!
129	tmltrd_trc(:,:,:,:) = 0.e0 ! <<< reset trend arrays to zero
130
131	! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
132	SELECT CASE ( nn_ctls_trc ) ! choice of the control surface
133	CASE ( -2 ) ; STOP 'trdmxl_trc : not ready ' ! -> isopycnal surface (see ???)
134	#if defined key_pisces \|\| defined key_pisces_reduced
135	CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion
136	#endif
137	CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl)
138	CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file
139	CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
140	nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level
141	END SELECT
142
143	! ... Compute ndextrd1 and ndimtrd1 ??? role de jpktrd_trc
144	IF( ionce == 1 ) THEN
145	!
146	isum = 0 ; zvlmsk(:,:) = 0.e0
147
148	IF( jpktrd_trc < jpk ) THEN ! description ???
149	DO jj = 1, jpj
150	DO ji = 1, jpi
151	IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
152	zvlmsk(ji,jj) = tmask(ji,jj,1)
153	ELSE
154	isum = isum + 1
155	zvlmsk(ji,jj) = 0.e0
156	ENDIF
157	END DO
158	END DO
159	ENDIF
160
161	IF( isum > 0 ) THEN ! index of ocean points (2D only)
162	WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
163	CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
164	ELSE
165	CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
166	ENDIF
167
168	ionce = 0 ! no more pass here
169	!
170	ENDIF ! ionce == 1
171
172	! ... Weights for vertical averaging
173	wkx_trc(:,:,:) = 0.e0
174	DO jk = 1, jpktrd_trc ! initialize wkx_trc with vertical scale factor in mixed-layer
175	DO jj = 1, jpj
176	DO ji = 1, jpi
177	IF( jk - nmld_trc(ji,jj) < 0 ) wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
178	END DO
179	END DO
180	END DO
181
182	rmld_trc(:,:) = 0.e0
183	DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc
184	rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
185	END DO
186
187	DO jk = 1, jpktrd_trc ! compute integration weights
188	wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
189	END DO
190
191	icount = 0 ! <<< flag = off : control surface & integr. weights
192	! ! computed only once per time step
193	ENDIF ONCE_PER_TIME_STEP
194
195	! II. Vertical integration of trends in the mixed-layer
196	! -----------------------------------------------------
197
198	SELECT CASE ( ctype )
199	CASE ( '3D' ) ! mean passive tracer trends in the mixed-layer
200	DO jk = 1, jpktrd_trc
201	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)
202	END DO
203	CASE ( '2D' ) ! forcing at upper boundary of the mixed-layer
204	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,1) * wkx_trc(:,:,1) ! non penetrative
205	END SELECT
206	!
207	CALL wrk_dealloc( jpi, jpj, zvlmsk )
208	!
209	END SUBROUTINE trd_mxl_trc_zint
210
211
212	SUBROUTINE trd_mxl_bio_zint( ptrc_trdmxl, ktrd )
213	!!----------------------------------------------------------------------
214	!! * ROUTINE trd_mxl_bio_zint *
215	!!
216	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
217	!! to the subroutine. This vertical average is performed from ocean
218	!! surface down to a chosen control surface.
219	!!
220	!! ** Method/usage :
221	!! The control surface can be either a mixed layer depth (time varying)
222	!! or a fixed surface (jk level or bowl).
223	!! Choose control surface with nctls in namelist NAMTRD :
224	!! nctls_trc = 0 : use mixed layer with density criterion
225	!! nctls_trc = 1 : read index from file 'ctlsurf_idx'
226	!! nctls_trc > 1 : use fixed level surface jk = nctls_trc
227	!! Note: in the remainder of the routine, the volume between the
228	!! surface and the control surface is called "mixed-layer"
229	!!----------------------------------------------------------------------
230	!!
231	INTEGER , INTENT(in) :: ktrd ! bio trend index
232	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrc_trdmxl ! passive trc trend
233	#if defined key_pisces_reduced
234	!
235	INTEGER :: ji, jj, jk, isum
236	REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
237	!!----------------------------------------------------------------------
238
239	CALL wrk_alloc( jpi, jpj, zvlmsk )
240
241	! I. Definition of control surface and integration weights
242	! --------------------------------------------------------
243	! ==> only once per time step <==
244
245	IF( icountbio == 1 ) THEN
246	!
247	tmltrd_bio(:,:,:) = 0.e0 ! <<< reset trend arrays to zero
248	! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
249	SELECT CASE ( nn_ctls_trc ) ! choice of the control surface
250	CASE ( -2 ) ; STOP 'trdmxl_trc : not ready ' ! -> isopycnal surface (see ???)
251	CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion
252	CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl)
253	CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file
254	CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
255	nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level
256	END SELECT
257
258	! ... Compute ndextrd1 and ndimtrd1 only once
259	IF( ioncebio == 1 ) THEN
260	!
261	! Check of validity : nmld_trc(ji,jj) <= jpktrd_trc
262	isum = 0
263	zvlmsk(:,:) = 0.e0
264
265	IF( jpktrd_trc < jpk ) THEN
266	DO jj = 1, jpj
267	DO ji = 1, jpi
268	IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
269	zvlmsk(ji,jj) = tmask(ji,jj,1)
270	ELSE
271	isum = isum + 1
272	zvlmsk(ji,jj) = 0.
273	END IF
274	END DO
275	END DO
276	END IF
277
278	! Index of ocean points (2D only)
279	IF( isum > 0 ) THEN
280	WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
281	CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
282	ELSE
283	CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
284	END IF
285
286	ioncebio = 0 ! no more pass here
287	!
288	END IF ! ( ioncebio == 1 )
289
290	! ... Weights for vertical averaging
291	wkx_trc(:,:,:) = 0.e0
292	DO jk = 1, jpktrd_trc ! initialize wkx_trc with vertical scale factor in mixed-layer
293	DO jj = 1,jpj
294	DO ji = 1,jpi
295	IF( jk - nmld_trc(ji,jj) < 0. ) wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
296	END DO
297	END DO
298	END DO
299
300	rmld_trc(:,:) = 0.
301	DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc
302	rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
303	END DO
304
305	DO jk = 1, jpktrd_trc ! compute integration weights
306	wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
307	END DO
308
309	icountbio = 0 ! <<< flag = off : control surface & integr. weights
310	! ! computed only once per time step
311	END IF ! ( icountbio == 1 )
312
313	! II. Vertical integration of trends in the mixed-layer
314	! -----------------------------------------------------
315
316
317	DO jk = 1, jpktrd_trc
318	tmltrd_bio(:,:,ktrd) = tmltrd_bio(:,:,ktrd) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)
319	END DO
320
321	CALL wrk_dealloc( jpi, jpj, zvlmsk )
322	#endif
323	!
324	END SUBROUTINE trd_mxl_bio_zint
325
326
327	SUBROUTINE trd_mxl_trc( kt )
328	!!----------------------------------------------------------------------
329	!! * ROUTINE trd_mxl_trc *
330	!!
331	!! ** Purpose : Compute and cumulate the mixed layer trends over an analysis
332	!! period, and write NetCDF outputs.
333	!!
334	!! ** Method/usage :
335	!! The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
336	!! logical namelist variable) :
337	!! 1) to explain the difference between initial and final
338	!! mixed-layer T & S (where initial and final relate to the
339	!! current analysis window, defined by ntrc_trc in the namelist)
340	!! 2) to explain the difference between the current and previous
341	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
342	!! performed over each analysis window).
343	!!
344	!! ** Consistency check :
345	!! If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
346	!! entrainment) should be zero, at machine accuracy. Note that in the case
347	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
348	!! over the first two analysis windows (except if restart).
349	!! N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8
350	!! for checking residuals.
351	!! On a NEC-SX5 computer, this typically leads to:
352	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
353	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
354	!!
355	!! ** Action :
356	!! At each time step, mixed-layer averaged trends are stored in the
357	!! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
358	!! This array is known when trd_mld is called, at the end of the stp subroutine,
359	!! except for the purely vertical K_z diffusion term, which is embedded in the
360	!! lateral diffusion trend.
361	!!
362	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
363	!! from the lateral diffusion trend.
364	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
365	!! arrays are updated.
366	!! In III), called only once per analysis window, we compute the total trends,
367	!! along with the residuals and the Asselin correction terms.
368	!! In IV), the appropriate trends are written in the trends NetCDF file.
369	!!
370	!! References :
371	!! - Vialard & al.
372	!! - See NEMO documentation (in preparation)
373	!!----------------------------------------------------------------------
374	!
375	INTEGER, INTENT(in) :: kt ! ocean time-step index
376	!
377	INTEGER :: ji, jj, jk, jl, ik, it, itmod, jn
378	REAL(wp) :: zavt, zfn, zfn2
379	!
380	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltot ! d(trc)/dt over the anlysis window (incl. Asselin)
381	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlres ! residual = dh/dt entrainment term
382	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlatf ! for storage only
383	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlrad ! for storage only (for trb<0 corr in trcrad)
384	!
385	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltot2 ! -+
386	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlres2 ! \| working arrays to diagnose the trends
387	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltrdm2 ! \| associated with the time meaned ML
388	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlatf2 ! \| passive tracers
389	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlrad2 ! \| (-> for trb<0 corr in trcrad)
390	!
391	CHARACTER (LEN=10) :: clvar
392	!!----------------------------------------------------------------------
393
394	! Set-up pointers into sub-arrays of workspaces
395	CALL wrk_alloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
396	CALL wrk_alloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
397
398	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
399
400	! ======================================================================
401	! I. Diagnose the purely vertical (K_z) diffusion trend
402	! ======================================================================
403
404	! ... These terms can be estimated by flux computation at the lower boundary of the ML
405	! (we compute (-1/h) * K_z * d_z( tracer ))
406
407	IF( ln_trcldf_iso ) THEN
408	!
409	DO jj = 1,jpj
410	DO ji = 1,jpi
411	ik = nmld_trc(ji,jj)
412	zavt = fsavs(ji,jj,ik)
413	DO jn = 1, jptra
414	IF( ln_trdtrc(jn) ) &
415	tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) &
416	& * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) ) &
417	& / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
418	END DO
419	END DO
420	END DO
421
422	DO jn = 1, jptra
423	! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
424	IF( ln_trdtrc(jn) ) &
425	tmltrd_trc(:,:,jpmxl_trc_ldf,jn) = tmltrd_trc(:,:,jpmxl_trc_ldf,jn) - tmltrd_trc(:,:,jpmxl_trc_zdf,jn)
426
427	END DO
428	!
429	ENDIF
430
431	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
432	! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm.
433	! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
434	DO jn = 1, jptra
435	IF( ln_trdtrc(jn) ) THEN
436	DO jl = 1, jpltrd_trc
437	CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions
438	END DO
439	ENDIF
440	END DO
441	ENDIF
442	! ======================================================================
443	! II. Cumulate the trends over the analysis window
444	! ======================================================================
445
446	ztmltrd2(:,:,:,:) = 0.e0 ; ztmltot2(:,:,:) = 0.e0 ! <<< reset arrays to zero
447	ztmlres2(:,:,:) = 0.e0 ; ztmlatf2(:,:,:) = 0.e0
448	ztmlrad2(:,:,:) = 0.e0
449
450	! II.1 Set before values of vertically averages passive tracers
451	! -------------------------------------------------------------
452	IF( kt > nittrc000 ) THEN
453	DO jn = 1, jptra
454	IF( ln_trdtrc(jn) ) THEN
455	tmlb_trc (:,:,jn) = tml_trc (:,:,jn)
456	tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_atf,jn)
457	tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_radb,jn)
458	ENDIF
459	END DO
460	ENDIF
461
462	! II.2 Vertically averaged passive tracers
463	! ----------------------------------------
464	tml_trc(:,:,:) = 0.e0
465	DO jk = 1, jpktrd_trc ! - 1 ???
466	DO jn = 1, jptra
467	IF( ln_trdtrc(jn) ) &
468	tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
469	END DO
470	END DO
471
472	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
473	! ------------------------------------------------------------------------
474	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1) ???
475	!
476	DO jn = 1, jptra
477	IF( ln_trdtrc(jn) ) THEN
478	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
479	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
480
481	tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0 ; tmltrd_atf_sumb_trc (:,:,jn) = 0.e0
482	tmltrd_rad_sumb_trc (:,:,jn) = 0.e0
483	ENDIF
484	END DO
485
486	rmldbn_trc(:,:) = rmld_trc(:,:)
487	!
488	ENDIF
489
490	! II.4 Cumulated trends over the analysis period
491	! ----------------------------------------------
492	!
493	! [ 1rst analysis window ] [ 2nd analysis window ]
494	!
495	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
496	! ntrd 2*ntrd etc.
497	! 1 2 3 4 =5 e.g. =10
498	!
499	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN ! ???
500	!
501	nmoymltrd = nmoymltrd + 1
502
503
504	! ... Cumulate over BOTH physical contributions AND over time steps
505	DO jn = 1, jptra
506	IF( ln_trdtrc(jn) ) THEN
507	DO jl = 1, jpltrd_trc
508	tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
509	END DO
510	ENDIF
511	END DO
512
513	DO jn = 1, jptra
514	IF( ln_trdtrc(jn) ) THEN
515	! ... Special handling of the Asselin trend
516	tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
517	tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
518
519	! ... Trends associated with the time mean of the ML passive tracers
520	tmltrd_sum_trc (:,:,:,jn) = tmltrd_sum_trc (:,:,:,jn) + tmltrd_trc (:,:,:,jn)
521	tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
522	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) + tml_trc (:,:,jn)
523	ENDIF
524	ENDDO
525
526	rmld_sum_trc (:,:) = rmld_sum_trc (:,:) + rmld_trc (:,:)
527	!
528	ENDIF
529
530	! ======================================================================
531	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
532	! ======================================================================
533
534	! Convert to appropriate physical units
535	tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
536
537	itmod = kt - nittrc000 + 1
538	it = kt
539
540	MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of nn_trd_trc
541	!
542	ztmltot (:,:,:) = 0.e0 ! reset arrays to zero
543	ztmlres (:,:,:) = 0.e0
544	ztmltot2(:,:,:) = 0.e0
545	ztmlres2(:,:,:) = 0.e0
546
547	zfn = FLOAT( nmoymltrd ) ; zfn2 = zfn * zfn
548
549	! III.1 Prepare fields for output ("instantaneous" diagnostics)
550	! -------------------------------------------------------------
551
552	DO jn = 1, jptra
553	IF( ln_trdtrc(jn) ) THEN
554	!-- Compute total trends (use rdttrc instead of rdt ???)
555	IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes
556	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
557	ELSE ! LEAP-FROG schemes
558	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
559	ENDIF
560
561	!-- Compute residuals
562	ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
563	& - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) )
564
565	!-- Diagnose Asselin trend over the analysis window
566	ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
567	ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
568
569	!-- Lateral boundary conditions
570	IF ( cp_cfg .NE. 'gyre' ) THEN
571	CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
572	CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
573	ENDIF
574
575
576	#if defined key_diainstant
577	STOP 'tmltrd_trc : key_diainstant was never checked within trdmxl. Comment this to proceed.'
578	#endif
579	ENDIF
580	END DO
581
582	! III.2 Prepare fields for output ("mean" diagnostics)
583	! ----------------------------------------------------
584
585	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
586	rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
587
588	!-- Compute passive tracer total trends
589	DO jn = 1, jptra
590	IF( ln_trdtrc(jn) ) THEN
591	tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
592	ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rdt ) ! now tracer unit is /sec
593	ENDIF
594	END DO
595
596	!-- Compute passive tracer residuals
597	DO jn = 1, jptra
598	IF( ln_trdtrc(jn) ) THEN
599	!
600	DO jl = 1, jpltrd_trc
601	ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
602	END DO
603
604	ztmltrdm2(:,:,jn) = 0.e0
605	DO jl = 1, jpltrd_trc
606	ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
607	END DO
608
609	ztmlres2(:,:,jn) = ztmltot2(:,:,jn) - &
610	& ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
611	& - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
612	!
613
614	!-- Diagnose Asselin trend over the analysis window
615	ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) &
616	& + tmltrd_atf_sumb_trc(:,:,jn)
617	ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) &
618	& + tmltrd_rad_sumb_trc(:,:,jn)
619
620	!-- Lateral boundary conditions
621	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
622	CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
623	CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
624	DO jl = 1, jpltrd_trc
625	CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file
626	END DO
627	ENDIF
628
629	ENDIF
630	END DO
631
632	! * Debugging information *
633	IF( lldebug ) THEN
634	!
635	WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
636	WRITE(numout,*) '~~~~~~~~~~~ '
637	WRITE(numout,*)
638	WRITE(numout,*) 'TRC kt = ', kt, ' nmoymltrd = ', nmoymltrd
639
640	DO jn = 1, jptra
641
642	IF( ln_trdtrc(jn) ) THEN
643	WRITE(numout, *)
644	WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
645
646	WRITE(numout, *)
647	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres : ', SUM2D(ztmlres(:,:,jn))
648	!CD??? PREVOIR: z2d = ztmlres(:,:,jn) ; CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres - : ')
649
650	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
651	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
652	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot : ', SUM2D(+ztmltot (:,:,jn))
653	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
654	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
655	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
656	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
657	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
658	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
659
660	WRITE(numout, *)
661	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2 : ', SUM2D(ztmlres2(:,:,jn))
662	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
663	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
664	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2 : ', SUM2D(+ztmltot2(:,:,jn))
665	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2 : ', SUM2D(+ztmltrdm2(:,:,jn))
666	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn))
667	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
668	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
669	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
670	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
671
672	WRITE(numout, *)
673	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot : ', SUM2D(ztmltot (:,:,jn))
674	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
675	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc : ', SUM2D(tml_trc (:,:,jn))
676	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc : ', SUM2D(tmlbn_trc (:,:,jn))
677	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc : ', SUM2D(tmlb_trc (:,:,jn))
678	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc : ', SUM2D(tmlbb_trc (:,:,jn))
679	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
680
681	WRITE(numout, *)
682	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
683	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
684	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
685	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
686	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
687
688	WRITE(numout, *)
689	DO jl = 1, jpltrd_trc
690	WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
691	& ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
692	END DO
693
694	WRITE(numout,*)
695	WRITE(numout,) ' ******************** ZTMLRES, ZTMLRES2 ********************* '
696	WRITE(numout,*)
697	WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
698	WRITE(numout,*)
699	WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
700
701	WRITE(numout,*)
702	WRITE(numout,) ' ******************** ZTMLRES ********************* '
703	WRITE(numout,*)
704
705	WRITE(numout,*) '...................................................'
706	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
707	DO jj = 5, 1, -1
708	WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
709	END DO
710
711	WRITE(numout,*)
712	WRITE(numout,) ' ******************** ZTMLRES2 ********************* '
713	WRITE(numout,*)
714
715	WRITE(numout,*) '...................................................'
716	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
717	DO jj = 5, 1, -1
718	WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
719	END DO
720	!
721	ENDIF
722	!
723	END DO
724
725
726	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
727	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
728	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
729	WRITE(numout,*)
730	!
731	ENDIF
732
733	! III.3 Time evolution array swap
734	! -------------------------------
735	! ML depth
736	rmldbn_trc(:,:) = rmld_trc(:,:)
737	rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
738	DO jn = 1, jptra
739	IF( ln_trdtrc(jn) ) THEN
740	! For passive tracer instantaneous diagnostics
741	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
742	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
743
744	! For passive tracer mean diagnostics
745	tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
746	tml_sumb_trc (:,:,jn) = tml_sum_trc(:,:,jn)
747	tmltrd_atf_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
748	tmltrd_rad_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
749
750
751	! III.4 Convert to appropriate physical units
752	! -------------------------------------------
753	ztmltot (:,:,jn) = ztmltot (:,:,jn) * rn_ucf_trc/zfn ! instant diags
754	ztmlres (:,:,jn) = ztmlres (:,:,jn) * rn_ucf_trc/zfn
755	ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * rn_ucf_trc/zfn
756	ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * rn_ucf_trc/zfn
757	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags
758	ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * rn_ucf_trc/zfn2
759	ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * rn_ucf_trc/zfn2
760	ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * rn_ucf_trc/zfn2
761	ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * rn_ucf_trc/zfn2
762	ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * rn_ucf_trc/zfn2
763	ENDIF
764	END DO
765	!
766	ENDIF MODULO_NTRD
767
768	! ======================================================================
769	! IV. Write trends in the NetCDF file
770	! ======================================================================
771
772	! IV.1 Code for IOIPSL/NetCDF output
773	! ----------------------------------
774
775	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
776	WRITE(numout,*) ' '
777	WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
778	WRITE(numout,*) '~~~~~~~~~~~ '
779	WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt
780	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
781	WRITE(numout,*) ' '
782	ENDIF
783
784	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
785	!
786
787	DO jn = 1, jptra
788	!
789	IF( ln_trdtrc(jn) ) THEN
790	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
791	!-- Output the fields
792	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
793	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 )
794	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 )
795	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 )
796
797	DO jl = 1, jpltrd_trc - 2
798	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
799	& it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
800	END DO
801
802	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
803	& it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
804
805	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
806	& it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
807
808	ENDIF
809	END DO
810
811	IF( kt == nitend ) THEN
812	DO jn = 1, jptra
813	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
814	END DO
815	ENDIF
816
817	ELSE ! <<< write the trends for passive tracer mean diagnostics
818
819	DO jn = 1, jptra
820	!
821	IF( ln_trdtrc(jn) ) THEN
822	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 )
823	!-- Output the fields
824	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
825
826	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
827	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 )
828	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 )
829
830	DO jl = 1, jpltrd_trc - 2
831	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
832	& it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
833	END DO
834
835	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
836	& it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
837
838	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
839	& it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
840
841	ENDIF
842	!
843	END DO
844	IF( kt == nitend ) THEN
845	DO jn = 1, jptra
846	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
847	END DO
848	ENDIF
849
850	!
851	ENDIF NETCDF_OUTPUT
852
853	! Compute the control surface (for next time step) : flag = on
854	icount = 1
855
856	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
857	!
858	! Reset cumulative arrays to zero
859	! -------------------------------
860	nmoymltrd = 0
861	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
862	tmlradm_trc (:,:,:) = 0.e0 ; tml_sum_trc (:,:,:) = 0.e0
863	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0
864	rmld_sum_trc (:,:) = 0.e0
865	!
866	ENDIF
867
868	! ======================================================================
869	! V. Write restart file
870	! ======================================================================
871
872	IF( lrst_trc ) CALL trd_mxl_trc_rst_write( kt ) ! this must be after the array swap above (III.3)
873
874	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
875	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
876	!
877	END SUBROUTINE trd_mxl_trc
878
879
880	SUBROUTINE trd_mxl_bio( kt )
881	!!----------------------------------------------------------------------
882	!! * ROUTINE trd_mld *
883	!!
884	!! ** Purpose : Compute and cumulate the mixed layer biological trends over an analysis
885	!! period, and write NetCDF outputs.
886	!!
887	!! ** Method/usage :
888	!! The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
889	!! logical namelist variable) :
890	!! 1) to explain the difference between initial and final
891	!! mixed-layer T & S (where initial and final relate to the
892	!! current analysis window, defined by ntrd in the namelist)
893	!! 2) to explain the difference between the current and previous
894	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
895	!! performed over each analysis window).
896	!!
897	!! ** Consistency check :
898	!! If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt
899	!! entrainment) should be zero, at machine accuracy. Note that in the case
900	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
901	!! over the first two analysis windows (except if restart).
902	!! N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8
903	!! for checking residuals.
904	!! On a NEC-SX5 computer, this typically leads to:
905	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
906	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
907	!!
908	!! ** Action :
909	!! At each time step, mixed-layer averaged trends are stored in the
910	!! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
911	!! This array is known when trd_mld is called, at the end of the stp subroutine,
912	!! except for the purely vertical K_z diffusion term, which is embedded in the
913	!! lateral diffusion trend.
914	!!
915	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
916	!! from the lateral diffusion trend.
917	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
918	!! arrays are updated.
919	!! In III), called only once per analysis window, we compute the total trends,
920	!! along with the residuals and the Asselin correction terms.
921	!! In IV), the appropriate trends are written in the trends NetCDF file.
922	!!
923	!! References :
924	!! - Vialard & al.
925	!! - See NEMO documentation (in preparation)
926	!!----------------------------------------------------------------------
927	INTEGER, INTENT( in ) :: kt ! ocean time-step index
928	#if defined key_pisces_reduced
929	INTEGER :: jl, it, itmod
930	LOGICAL :: llwarn = .TRUE., lldebug = .TRUE.
931	REAL(wp) :: zfn, zfn2
932	!!----------------------------------------------------------------------
933	! ... Warnings
934	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
935
936	! ======================================================================
937	! II. Cumulate the trends over the analysis window
938	! ======================================================================
939
940	ztmltrdbio2(:,:,:) = 0.e0 ! <<< reset arrays to zero
941
942	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
943	! ------------------------------------------------------------------------
944	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
945	!
946	tmltrd_csum_ub_bio (:,:,:) = 0.e0
947	!
948	END IF
949
950	! II.4 Cumulated trends over the analysis period
951	! ----------------------------------------------
952	!
953	! [ 1rst analysis window ] [ 2nd analysis window ]
954	!
955	!
956	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
957	! ntrd 2*ntrd etc.
958	! 1 2 3 4 =5 e.g. =10
959	!
960	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN
961	!
962	nmoymltrdbio = nmoymltrdbio + 1
963
964	! ... Trends associated with the time mean of the ML passive tracers
965	tmltrd_sum_bio (:,:,:) = tmltrd_sum_bio (:,:,:) + tmltrd_bio (:,:,:)
966	tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
967	!
968	END IF
969
970	! ======================================================================
971	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
972	! ======================================================================
973
974	! Convert to appropriate physical units
975	tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc
976
977	MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd
978	!
979	zfn = float(nmoymltrdbio) ; zfn2 = zfn * zfn
980
981	! III.1 Prepare fields for output ("instantaneous" diagnostics)
982	! -------------------------------------------------------------
983
984	#if defined key_diainstant
985	STOP 'tmltrd_bio : key_diainstant was never checked within trdmxl. Comment this to proceed.'
986	#endif
987	! III.2 Prepare fields for output ("mean" diagnostics)
988	! ----------------------------------------------------
989
990	ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)
991
992	!-- Lateral boundary conditions
993	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
994	! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut
995	DO jn = 1, jpdiabio
996	CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )
997	ENDDO
998	ENDIF
999
1000	IF( lldebug ) THEN
1001	!
1002	WRITE(numout,*) 'trd_mxl_bio : write trends in the Mixed Layer for debugging process:'
1003	WRITE(numout,*) '~~~~~~~~~~~ '
1004	WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio
1005	WRITE(numout,*)
1006
1007	DO jl = 1, jpdiabio
1008	IF( ln_trdmxl_trc_instant ) THEN
1009	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1010	& ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))
1011	ELSE
1012	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1013	& ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))
1014	endif
1015	END DO
1016
1017	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
1018	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
1019	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
1020	WRITE(numout,*)
1021	!
1022	ENDIF
1023
1024	! III.3 Time evolution array swap
1025	! -------------------------------
1026
1027	! For passive tracer mean diagnostics
1028	tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)
1029
1030	! III.4 Convert to appropriate physical units
1031	! -------------------------------------------
1032	ztmltrdbio2 (:,:,:) = ztmltrdbio2 (:,:,:) * rn_ucf_trc/zfn2
1033
1034	END IF MODULO_NTRD
1035
1036	! ======================================================================
1037	! IV. Write trends in the NetCDF file
1038	! ======================================================================
1039
1040	! IV.1 Code for IOIPSL/NetCDF output
1041	! ----------------------------------
1042
1043	! define time axis
1044	itmod = kt - nittrc000 + 1
1045	it = kt
1046
1047	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
1048	WRITE(numout,*) ' '
1049	WRITE(numout,*) 'trd_mxl_bio : write ML bio trends in the NetCDF file :'
1050	WRITE(numout,*) '~~~~~~~~~~~ '
1051	WRITE(numout,*) ' ', TRIM(clhstnam), ' at kt = ', kt
1052	WRITE(numout,*) ' N.B. nmoymltrdbio = ', nmoymltrdbio
1053	WRITE(numout,*) ' '
1054	END IF
1055
1056
1057	! 2. Start writing data
1058	! ---------------------
1059
1060	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
1061	!
1062	DO jl = 1, jpdiabio
1063	CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) , &
1064	& it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
1065	END DO
1066
1067
1068	IF( kt == nitend ) CALL histclo( nidtrdbio )
1069
1070	ELSE ! <<< write the trends for passive tracer mean diagnostics
1071
1072	DO jl = 1, jpdiabio
1073	CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) , &
1074	& it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
1075	END DO
1076
1077	IF( kt == nitend ) CALL histclo( nidtrdbio )
1078	!
1079	END IF NETCDF_OUTPUT
1080
1081	! Compute the control surface (for next time step) : flag = on
1082	icountbio = 1
1083
1084
1085
1086	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
1087	!
1088	! III.5 Reset cumulative arrays to zero
1089	! -------------------------------------
1090	nmoymltrdbio = 0
1091	tmltrd_csum_ln_bio (:,:,:) = 0.e0
1092	tmltrd_sum_bio (:,:,:) = 0.e0
1093	END IF
1094
1095	! ======================================================================
1096	! Write restart file
1097	! ======================================================================
1098
1099	! restart write is done in trd_mxl_trc_write which is called by trd_mxl_bio (Marina)
1100	!
1101	#endif
1102	END SUBROUTINE trd_mxl_bio
1103
1104
1105	REAL FUNCTION sum2d( ztab )
1106	!!----------------------------------------------------------------------
1107	!! CD ??? prevoir d'utiliser plutot prtctl
1108	!!----------------------------------------------------------------------
1109	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: ztab
1110	!!----------------------------------------------------------------------
1111	sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
1112	END FUNCTION sum2d
1113
1114
1115	SUBROUTINE trd_mxl_trc_init
1116	!!----------------------------------------------------------------------
1117	!! * ROUTINE trd_mxl_init *
1118	!!
1119	!! ** Purpose : computation of vertically integrated T and S budgets
1120	!! from ocean surface down to control surface (NetCDF output)
1121	!!
1122	!! ** Method/usage :
1123	!!
1124	!!----------------------------------------------------------------------
1125	INTEGER :: inum ! logical unit
1126	INTEGER :: ilseq, jl, jn, iiter
1127	REAL(wp) :: zjulian, zsto, zout
1128	CHARACTER (LEN=40) :: clop
1129	CHARACTER (LEN=15) :: csuff
1130	CHARACTER (LEN=12) :: clmxl
1131	CHARACTER (LEN=16) :: cltrcu
1132	CHARACTER (LEN=10) :: clvar
1133
1134	!!----------------------------------------------------------------------
1135
1136	! ======================================================================
1137	! I. initialization
1138	! ======================================================================
1139
1140	IF(lwp) THEN
1141	WRITE(numout,*)
1142	WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers '
1143	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
1144	WRITE(numout,*)
1145	ENDIF
1146
1147
1148	! I.1 Check consistency of user defined preferences
1149	! -------------------------------------------------
1150
1151	IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
1152	WRITE(numout,cform_err)
1153	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
1154	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
1155	WRITE(numout,*) ' you defined, nn_trd_trc = ', nn_trd_trc
1156	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
1157	WRITE(numout,*) ' You should reconsider this choice. '
1158	WRITE(numout,*)
1159	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
1160	WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) '
1161	nstop = nstop + 1
1162	ENDIF
1163
1164	! * Debugging information *
1165	IF( lldebug ) THEN
1166	WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl
1167	WRITE(numout,*) ' ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
1168	ENDIF
1169
1170	IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
1171	WRITE(numout,cform_err)
1172	WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL '
1173	WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '
1174	WRITE(numout,*) ' WHY? Everything in trdmxl_trc is coded for leap-frog, and '
1175	WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note '
1176	WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). '
1177	WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '
1178	WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'
1179	WRITE(numout,*)
1180	nstop = nstop + 1
1181	ENDIF
1182
1183	! I.2 Initialize arrays to zero or read a restart file
1184	! ----------------------------------------------------
1185	nmoymltrd = 0
1186
1187	rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst.
1188	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
1189	tmlradm_trc (:,:,:) = 0.e0
1190
1191	tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean
1192	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0
1193
1194	#if defined key_pisces_reduced
1195	nmoymltrdbio = 0
1196	tmltrd_sum_bio (:,:,:) = 0.e0 ; tmltrd_csum_ln_bio (:,:,:) = 0.e0
1197	DO jl = 1, jp_pisces_trd
1198	ctrd_bio(jl,1) = ctrbil(jl) ! long name
1199	ctrd_bio(jl,2) = ctrbio(jl) ! short name
1200	ENDDO
1201	#endif
1202
1203	IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
1204	CALL trd_mxl_trc_rst_read
1205	ELSE
1206	tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst.
1207	tmlbn_trc (:,:,:) = 0.e0
1208
1209	tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean
1210	tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0
1211	#if defined key_pisces_reduced
1212	tmltrd_csum_ub_bio (:,:,:) = 0.e0
1213	#endif
1214
1215	ENDIF
1216
1217	icount = 1 ; ionce = 1 ! open specifier
1218
1219	#if defined key_pisces_reduced
1220	icountbio = 1 ; ioncebio = 1 ! open specifier
1221	#endif
1222
1223	! I.3 Read control surface from file ctlsurf_idx
1224	! ----------------------------------------------
1225	IF( nn_ctls_trc == 1 ) THEN
1226	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
1227	READ ( inum ) nbol_trc
1228	CLOSE( inum )
1229	ENDIF
1230
1231	! ======================================================================
1232	! II. netCDF output initialization
1233	! ======================================================================
1234
1235	! clmxl = legend root for netCDF output
1236	IF( nn_ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion
1237	clmxl = 'Mixed Layer '
1238	ELSE IF( nn_ctls_trc == 1 ) THEN ! control surface = read index from file
1239	clmxl = ' Bowl '
1240	ELSE IF( nn_ctls_trc >= 2 ) THEN ! control surface = model level
1241	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
1242	ENDIF
1243
1244	! II.1 Define frequency of output and means
1245	! -----------------------------------------
1246	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
1247	ELSE ; clop = "x" ! no use of the mask value (require less cp time)
1248	ENDIF
1249	# if defined key_diainstant
1250	IF( .NOT. ln_trdmxl_trc_instant ) THEN
1251	STOP 'trd_mxl_trc : this was never checked. Comment this line to proceed...'
1252	ENDIF
1253	zsto = nn_trd_trc * rdt
1254	clop = "inst("//TRIM(clop)//")"
1255	# else
1256	IF( ln_trdmxl_trc_instant ) THEN
1257	zsto = rdt ! inst. diags : we use IOIPSL time averaging
1258	ELSE
1259	zsto = nn_trd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging
1260	ENDIF
1261	clop = "ave("//TRIM(clop)//")"
1262	# endif
1263	zout = nn_trd_trc * rdt
1264	iiter = ( nittrc000 - 1 ) / nn_dttrc
1265
1266	IF(lwp) WRITE (numout,*) ' netCDF initialization'
1267
1268	! II.2 Compute julian date from starting date of the run
1269	! ------------------------------------------------------
1270	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
1271	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
1272	IF(lwp) WRITE(numout,*)' '
1273	IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000 &
1274	& ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday &
1275	& ,'Julian day : ', zjulian
1276
1277	! II.3 Define the T grid trend file (nidtrd)
1278	! ------------------------------------------
1279
1280	!-- Define long and short names for the NetCDF output variables
1281	! ==> choose them according to trdmxl_trc_oce.F90 <==
1282
1283	ctrd_trc(jpmxl_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmxl_trc_xad ,2) = "_xad"
1284	ctrd_trc(jpmxl_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmxl_trc_yad ,2) = "_yad"
1285	ctrd_trc(jpmxl_trc_zad ,1) = " Vertical advection" ; ctrd_trc(jpmxl_trc_zad ,2) = "_zad"
1286	ctrd_trc(jpmxl_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmxl_trc_ldf ,2) = "_ldf"
1287	ctrd_trc(jpmxl_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmxl_trc_zdf ,2) = "_zdf"
1288	ctrd_trc(jpmxl_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmxl_trc_bbl ,2) = "_bbl"
1289	ctrd_trc(jpmxl_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmxl_trc_dmp ,2) = "_dmp"
1290	ctrd_trc(jpmxl_trc_sbc ,1) = " Surface boundary cond." ; ctrd_trc(jpmxl_trc_sbc ,2) = "_sbc"
1291	ctrd_trc(jpmxl_trc_sms, 1) = " Sources minus sinks" ; ctrd_trc(jpmxl_trc_sms ,2) = "_sms"
1292	ctrd_trc(jpmxl_trc_radb ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radb ,2) = "_radb"
1293	ctrd_trc(jpmxl_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radn ,2) = "_radn"
1294	ctrd_trc(jpmxl_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmxl_trc_atf ,2) = "_atf"
1295
1296	DO jn = 1, jptra
1297	!-- Create a NetCDF file and enter the define mode
1298	IF( ln_trdtrc(jn) ) THEN
1299	csuff="ML_"//ctrcnm(jn)
1300	CALL dia_nam( clhstnam, nn_trd_trc, csuff )
1301	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1302	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
1303
1304	!-- Define the ML depth variable
1305	CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m", &
1306	& jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1307
1308	ENDIF
1309	END DO
1310
1311	#if defined key_pisces_reduced
1312	!-- Create a NetCDF file and enter the define mode
1313	CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )
1314	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1315	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom, snc4chunks=snc4set )
1316	#endif
1317
1318	!-- Define physical units
1319	IF( rn_ucf_trc == 1. ) THEN
1320	cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit
1321	ELSEIF ( rn_ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3)
1322	cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf
1323	ELSE
1324	cltrcu = "unknown?"
1325	ENDIF
1326
1327	!-- Define miscellaneous passive tracer mixed-layer variables
1328	IF( jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
1329	STOP 'Error : jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb' ! see below
1330	ENDIF
1331
1332	DO jn = 1, jptra
1333	!
1334	IF( ln_trdtrc(jn) ) THEN
1335	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc.
1336	CALL histdef(nidtrd(jn), trim(clvar), clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", &
1337	& "mmole-N/m3", jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1338	CALL histdef(nidtrd(jn), trim(clvar)//"_tot" , clmxl//" "//trim(ctrcnm(jn))//" Total trend ", &
1339	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1340	CALL histdef(nidtrd(jn), trim(clvar)//"_res" , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)", &
1341	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1342
1343	DO jl = 1, jpltrd_trc - 2 ! <== only true if jpltrd_trc == jpmxl_trc_atf
1344	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1), &
1345	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1346	END DO ! if zsto=rdt above
1347
1348	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), &
1349	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1350
1351	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1), &
1352	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1353	!
1354	ENDIF
1355	END DO
1356
1357	#if defined key_pisces_reduced
1358	DO jl = 1, jp_pisces_trd
1359	CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1)) , &
1360	& cltrcu, jpi, jpj, nh_tb, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1361	END DO ! if zsto=rdt above
1362	#endif
1363
1364	!-- Leave IOIPSL/NetCDF define mode
1365	DO jn = 1, jptra
1366	IF( ln_trdtrc(jn) ) CALL histend( nidtrd(jn), snc4set )
1367	END DO
1368
1369	#if defined key_pisces_reduced
1370	!-- Leave IOIPSL/NetCDF define mode
1371	CALL histend( nidtrdbio, snc4set )
1372
1373	IF(lwp) WRITE(numout,*)
1374	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
1375	#endif
1376
1377	END SUBROUTINE trd_mxl_trc_init
1378
1379	#else
1380	!!----------------------------------------------------------------------
1381	!! Default option : Empty module
1382	!!----------------------------------------------------------------------
1383	CONTAINS
1384	SUBROUTINE trd_mxl_trc( kt ) ! Empty routine
1385	INTEGER, INTENT( in) :: kt
1386	WRITE(,) 'trd_mxl_trc: You should not have seen this print! error?', kt
1387	END SUBROUTINE trd_mxl_trc
1388	SUBROUTINE trd_mxl_bio( kt )
1389	INTEGER, INTENT( in) :: kt
1390	WRITE(,) 'trd_mxl_bio: You should not have seen this print! error?', kt
1391	END SUBROUTINE trd_mxl_bio
1392	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
1393	INTEGER , INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
1394	CHARACTER(len=2) , INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
1395	REAL, DIMENSION(:,:,:), INTENT( in ) :: ptrc_trdmxl ! passive trc trend
1396	WRITE(,) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
1397	WRITE(,) ' " " : You should not have seen this print! error?', ctype
1398	WRITE(,) ' " " : You should not have seen this print! error?', ktrd
1399	WRITE(,) ' " " : You should not have seen this print! error?', kjn
1400	END SUBROUTINE trd_mxl_trc_zint
1401	SUBROUTINE trd_mxl_trc_init ! Empty routine
1402	WRITE(,) 'trd_mxl_trc_init: You should not have seen this print! error?'
1403	END SUBROUTINE trd_mxl_trc_init
1404	#endif
1405
1406	!!======================================================================
1407	END MODULE trdmxl_trc

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

Download in other formats: