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 trunk/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 6137

Last change on this file since 6137 was 5836, checked in by cetlod, 9 years ago
merge the simplification branch onto the trunk, see ticket #1612
Property svn:keywords set to `Id`
File size: 68.6 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 (or dimg) 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	#if defined key_dimgout
393	INTEGER :: iyear,imon,iday
394	CHARACTER(LEN=80) :: cltext, clmode
395	#endif
396	!!----------------------------------------------------------------------
397
398	! Set-up pointers into sub-arrays of workspaces
399	CALL wrk_alloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
400	CALL wrk_alloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
401
402	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
403
404	! ======================================================================
405	! I. Diagnose the purely vertical (K_z) diffusion trend
406	! ======================================================================
407
408	! ... These terms can be estimated by flux computation at the lower boundary of the ML
409	! (we compute (-1/h) * K_z * d_z( tracer ))
410
411	IF( ln_trcldf_iso ) THEN
412	!
413	DO jj = 1,jpj
414	DO ji = 1,jpi
415	ik = nmld_trc(ji,jj)
416	zavt = fsavs(ji,jj,ik)
417	DO jn = 1, jptra
418	IF( ln_trdtrc(jn) ) &
419	tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) &
420	& * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) ) &
421	& / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
422	END DO
423	END DO
424	END DO
425
426	DO jn = 1, jptra
427	! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
428	IF( ln_trdtrc(jn) ) &
429	tmltrd_trc(:,:,jpmxl_trc_ldf,jn) = tmltrd_trc(:,:,jpmxl_trc_ldf,jn) - tmltrd_trc(:,:,jpmxl_trc_zdf,jn)
430
431	END DO
432	!
433	ENDIF
434
435	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
436	! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm.
437	! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
438	DO jn = 1, jptra
439	IF( ln_trdtrc(jn) ) THEN
440	DO jl = 1, jpltrd_trc
441	CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions
442	END DO
443	ENDIF
444	END DO
445	ENDIF
446	! ======================================================================
447	! II. Cumulate the trends over the analysis window
448	! ======================================================================
449
450	ztmltrd2(:,:,:,:) = 0.e0 ; ztmltot2(:,:,:) = 0.e0 ! <<< reset arrays to zero
451	ztmlres2(:,:,:) = 0.e0 ; ztmlatf2(:,:,:) = 0.e0
452	ztmlrad2(:,:,:) = 0.e0
453
454	! II.1 Set before values of vertically averages passive tracers
455	! -------------------------------------------------------------
456	IF( kt > nittrc000 ) THEN
457	DO jn = 1, jptra
458	IF( ln_trdtrc(jn) ) THEN
459	tmlb_trc (:,:,jn) = tml_trc (:,:,jn)
460	tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_atf,jn)
461	tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_radb,jn)
462	ENDIF
463	END DO
464	ENDIF
465
466	! II.2 Vertically averaged passive tracers
467	! ----------------------------------------
468	tml_trc(:,:,:) = 0.e0
469	DO jk = 1, jpktrd_trc ! - 1 ???
470	DO jn = 1, jptra
471	IF( ln_trdtrc(jn) ) &
472	tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
473	END DO
474	END DO
475
476	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
477	! ------------------------------------------------------------------------
478	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1) ???
479	!
480	DO jn = 1, jptra
481	IF( ln_trdtrc(jn) ) THEN
482	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
483	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
484
485	tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0 ; tmltrd_atf_sumb_trc (:,:,jn) = 0.e0
486	tmltrd_rad_sumb_trc (:,:,jn) = 0.e0
487	ENDIF
488	END DO
489
490	rmldbn_trc(:,:) = rmld_trc(:,:)
491	!
492	ENDIF
493
494	! II.4 Cumulated trends over the analysis period
495	! ----------------------------------------------
496	!
497	! [ 1rst analysis window ] [ 2nd analysis window ]
498	!
499	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
500	! ntrd 2*ntrd etc.
501	! 1 2 3 4 =5 e.g. =10
502	!
503	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN ! ???
504	!
505	nmoymltrd = nmoymltrd + 1
506
507
508	! ... Cumulate over BOTH physical contributions AND over time steps
509	DO jn = 1, jptra
510	IF( ln_trdtrc(jn) ) THEN
511	DO jl = 1, jpltrd_trc
512	tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
513	END DO
514	ENDIF
515	END DO
516
517	DO jn = 1, jptra
518	IF( ln_trdtrc(jn) ) THEN
519	! ... Special handling of the Asselin trend
520	tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
521	tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
522
523	! ... Trends associated with the time mean of the ML passive tracers
524	tmltrd_sum_trc (:,:,:,jn) = tmltrd_sum_trc (:,:,:,jn) + tmltrd_trc (:,:,:,jn)
525	tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
526	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) + tml_trc (:,:,jn)
527	ENDIF
528	ENDDO
529
530	rmld_sum_trc (:,:) = rmld_sum_trc (:,:) + rmld_trc (:,:)
531	!
532	ENDIF
533
534	! ======================================================================
535	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
536	! ======================================================================
537
538	! Convert to appropriate physical units
539	tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
540
541	itmod = kt - nittrc000 + 1
542	it = kt
543
544	MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of nn_trd_trc
545	!
546	ztmltot (:,:,:) = 0.e0 ! reset arrays to zero
547	ztmlres (:,:,:) = 0.e0
548	ztmltot2(:,:,:) = 0.e0
549	ztmlres2(:,:,:) = 0.e0
550
551	zfn = FLOAT( nmoymltrd ) ; zfn2 = zfn * zfn
552
553	! III.1 Prepare fields for output ("instantaneous" diagnostics)
554	! -------------------------------------------------------------
555
556	DO jn = 1, jptra
557	IF( ln_trdtrc(jn) ) THEN
558	!-- Compute total trends (use rdttrc instead of rdt ???)
559	IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes
560	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
561	ELSE ! LEAP-FROG schemes
562	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
563	ENDIF
564
565	!-- Compute residuals
566	ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
567	& - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) )
568
569	!-- Diagnose Asselin trend over the analysis window
570	ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
571	ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
572
573	!-- Lateral boundary conditions
574	IF ( cp_cfg .NE. 'gyre' ) THEN
575	CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
576	CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
577	ENDIF
578
579
580	#if defined key_diainstant
581	STOP 'tmltrd_trc : key_diainstant was never checked within trdmxl. Comment this to proceed.'
582	#endif
583	ENDIF
584	END DO
585
586	! III.2 Prepare fields for output ("mean" diagnostics)
587	! ----------------------------------------------------
588
589	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
590	rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
591
592	!-- Compute passive tracer total trends
593	DO jn = 1, jptra
594	IF( ln_trdtrc(jn) ) THEN
595	tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
596	ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rdt ) ! now tracer unit is /sec
597	ENDIF
598	END DO
599
600	!-- Compute passive tracer residuals
601	DO jn = 1, jptra
602	IF( ln_trdtrc(jn) ) THEN
603	!
604	DO jl = 1, jpltrd_trc
605	ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
606	END DO
607
608	ztmltrdm2(:,:,jn) = 0.e0
609	DO jl = 1, jpltrd_trc
610	ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
611	END DO
612
613	ztmlres2(:,:,jn) = ztmltot2(:,:,jn) - &
614	& ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
615	& - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
616	!
617
618	!-- Diagnose Asselin trend over the analysis window
619	ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) &
620	& + tmltrd_atf_sumb_trc(:,:,jn)
621	ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) &
622	& + tmltrd_rad_sumb_trc(:,:,jn)
623
624	!-- Lateral boundary conditions
625	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
626	CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
627	CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
628	DO jl = 1, jpltrd_trc
629	CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file
630	END DO
631	ENDIF
632
633	ENDIF
634	END DO
635
636	! * Debugging information *
637	IF( lldebug ) THEN
638	!
639	WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
640	WRITE(numout,*) '~~~~~~~~~~~ '
641	WRITE(numout,*)
642	WRITE(numout,*) 'TRC kt = ', kt, ' nmoymltrd = ', nmoymltrd
643
644	DO jn = 1, jptra
645
646	IF( ln_trdtrc(jn) ) THEN
647	WRITE(numout, *)
648	WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
649
650	WRITE(numout, *)
651	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres : ', SUM2D(ztmlres(:,:,jn))
652	!CD??? PREVOIR: z2d = ztmlres(:,:,jn) ; CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres - : ')
653
654	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
655	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
656	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot : ', SUM2D(+ztmltot (:,:,jn))
657	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
658	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
659	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
660	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
661	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
662	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
663
664	WRITE(numout, *)
665	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2 : ', SUM2D(ztmlres2(:,:,jn))
666	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
667	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
668	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2 : ', SUM2D(+ztmltot2(:,:,jn))
669	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2 : ', SUM2D(+ztmltrdm2(:,:,jn))
670	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn))
671	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
672	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
673	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
674	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
675
676	WRITE(numout, *)
677	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot : ', SUM2D(ztmltot (:,:,jn))
678	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
679	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc : ', SUM2D(tml_trc (:,:,jn))
680	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc : ', SUM2D(tmlbn_trc (:,:,jn))
681	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc : ', SUM2D(tmlb_trc (:,:,jn))
682	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc : ', SUM2D(tmlbb_trc (:,:,jn))
683	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
684
685	WRITE(numout, *)
686	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
687	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
688	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
689	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
690	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
691
692	WRITE(numout, *)
693	DO jl = 1, jpltrd_trc
694	WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
695	& ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
696	END DO
697
698	WRITE(numout,*)
699	WRITE(numout,) ' ******************** ZTMLRES, ZTMLRES2 ********************* '
700	WRITE(numout,*)
701	WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
702	WRITE(numout,*)
703	WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
704
705	WRITE(numout,*)
706	WRITE(numout,) ' ******************** ZTMLRES ********************* '
707	WRITE(numout,*)
708
709	WRITE(numout,*) '...................................................'
710	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
711	DO jj = 5, 1, -1
712	WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
713	END DO
714
715	WRITE(numout,*)
716	WRITE(numout,) ' ******************** ZTMLRES2 ********************* '
717	WRITE(numout,*)
718
719	WRITE(numout,*) '...................................................'
720	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
721	DO jj = 5, 1, -1
722	WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
723	END DO
724	!
725	ENDIF
726	!
727	END DO
728
729
730	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
731	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
732	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
733	WRITE(numout,*)
734	!
735	ENDIF
736
737	! III.3 Time evolution array swap
738	! -------------------------------
739	! ML depth
740	rmldbn_trc(:,:) = rmld_trc(:,:)
741	rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
742	DO jn = 1, jptra
743	IF( ln_trdtrc(jn) ) THEN
744	! For passive tracer instantaneous diagnostics
745	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
746	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
747
748	! For passive tracer mean diagnostics
749	tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
750	tml_sumb_trc (:,:,jn) = tml_sum_trc(:,:,jn)
751	tmltrd_atf_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
752	tmltrd_rad_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
753
754
755	! III.4 Convert to appropriate physical units
756	! -------------------------------------------
757	ztmltot (:,:,jn) = ztmltot (:,:,jn) * rn_ucf_trc/zfn ! instant diags
758	ztmlres (:,:,jn) = ztmlres (:,:,jn) * rn_ucf_trc/zfn
759	ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * rn_ucf_trc/zfn
760	ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * rn_ucf_trc/zfn
761	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags
762	ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * rn_ucf_trc/zfn2
763	ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * rn_ucf_trc/zfn2
764	ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * rn_ucf_trc/zfn2
765	ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * rn_ucf_trc/zfn2
766	ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * rn_ucf_trc/zfn2
767	ENDIF
768	END DO
769	!
770	ENDIF MODULO_NTRD
771
772	! ======================================================================
773	! IV. Write trends in the NetCDF file
774	! ======================================================================
775
776	! IV.1 Code for dimg mpp output
777	! -----------------------------
778
779	# if defined key_dimgout
780	STOP 'Not implemented'
781	# else
782
783	! IV.2 Code for IOIPSL/NetCDF output
784	! ----------------------------------
785
786	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
787	WRITE(numout,*) ' '
788	WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
789	WRITE(numout,*) '~~~~~~~~~~~ '
790	WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt
791	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
792	WRITE(numout,*) ' '
793	ENDIF
794
795	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
796	!
797
798	DO jn = 1, jptra
799	!
800	IF( ln_trdtrc(jn) ) THEN
801	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
802	!-- Output the fields
803	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
804	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 )
805	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 )
806	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 )
807
808	DO jl = 1, jpltrd_trc - 2
809	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
810	& it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
811	END DO
812
813	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
814	& it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
815
816	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
817	& it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
818
819	ENDIF
820	END DO
821
822	IF( kt == nitend ) THEN
823	DO jn = 1, jptra
824	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
825	END DO
826	ENDIF
827
828	ELSE ! <<< write the trends for passive tracer mean diagnostics
829
830	DO jn = 1, jptra
831	!
832	IF( ln_trdtrc(jn) ) THEN
833	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 )
834	!-- Output the fields
835	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
836
837	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
838	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 )
839	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 )
840
841	DO jl = 1, jpltrd_trc - 2
842	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
843	& it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
844	END DO
845
846	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
847	& it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
848
849	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
850	& it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
851
852	ENDIF
853	!
854	END DO
855	IF( kt == nitend ) THEN
856	DO jn = 1, jptra
857	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
858	END DO
859	ENDIF
860
861	!
862	ENDIF NETCDF_OUTPUT
863
864	! Compute the control surface (for next time step) : flag = on
865	icount = 1
866
867	# endif /* key_dimgout */
868
869	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
870	!
871	! Reset cumulative arrays to zero
872	! -------------------------------
873	nmoymltrd = 0
874	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
875	tmlradm_trc (:,:,:) = 0.e0 ; tml_sum_trc (:,:,:) = 0.e0
876	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0
877	rmld_sum_trc (:,:) = 0.e0
878	!
879	ENDIF
880
881	! ======================================================================
882	! V. Write restart file
883	! ======================================================================
884
885	IF( lrst_trc ) CALL trd_mxl_trc_rst_write( kt ) ! this must be after the array swap above (III.3)
886
887	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
888	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
889	!
890	END SUBROUTINE trd_mxl_trc
891
892
893	SUBROUTINE trd_mxl_bio( kt )
894	!!----------------------------------------------------------------------
895	!! * ROUTINE trd_mld *
896	!!
897	!! ** Purpose : Compute and cumulate the mixed layer biological trends over an analysis
898	!! period, and write NetCDF (or dimg) outputs.
899	!!
900	!! ** Method/usage :
901	!! The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
902	!! logical namelist variable) :
903	!! 1) to explain the difference between initial and final
904	!! mixed-layer T & S (where initial and final relate to the
905	!! current analysis window, defined by ntrd in the namelist)
906	!! 2) to explain the difference between the current and previous
907	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
908	!! performed over each analysis window).
909	!!
910	!! ** Consistency check :
911	!! If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt
912	!! entrainment) should be zero, at machine accuracy. Note that in the case
913	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
914	!! over the first two analysis windows (except if restart).
915	!! N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8
916	!! for checking residuals.
917	!! On a NEC-SX5 computer, this typically leads to:
918	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
919	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
920	!!
921	!! ** Action :
922	!! At each time step, mixed-layer averaged trends are stored in the
923	!! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
924	!! This array is known when trd_mld is called, at the end of the stp subroutine,
925	!! except for the purely vertical K_z diffusion term, which is embedded in the
926	!! lateral diffusion trend.
927	!!
928	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
929	!! from the lateral diffusion trend.
930	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
931	!! arrays are updated.
932	!! In III), called only once per analysis window, we compute the total trends,
933	!! along with the residuals and the Asselin correction terms.
934	!! In IV), the appropriate trends are written in the trends NetCDF file.
935	!!
936	!! References :
937	!! - Vialard & al.
938	!! - See NEMO documentation (in preparation)
939	!!----------------------------------------------------------------------
940	INTEGER, INTENT( in ) :: kt ! ocean time-step index
941	#if defined key_pisces_reduced
942	INTEGER :: jl, it, itmod
943	LOGICAL :: llwarn = .TRUE., lldebug = .TRUE.
944	REAL(wp) :: zfn, zfn2
945	#if defined key_dimgout
946	INTEGER :: iyear,imon,iday
947	CHARACTER(LEN=80) :: cltext, clmode
948	#endif
949	!!----------------------------------------------------------------------
950	! ... Warnings
951	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
952
953	! ======================================================================
954	! II. Cumulate the trends over the analysis window
955	! ======================================================================
956
957	ztmltrdbio2(:,:,:) = 0.e0 ! <<< reset arrays to zero
958
959	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
960	! ------------------------------------------------------------------------
961	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
962	!
963	tmltrd_csum_ub_bio (:,:,:) = 0.e0
964	!
965	END IF
966
967	! II.4 Cumulated trends over the analysis period
968	! ----------------------------------------------
969	!
970	! [ 1rst analysis window ] [ 2nd analysis window ]
971	!
972	!
973	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
974	! ntrd 2*ntrd etc.
975	! 1 2 3 4 =5 e.g. =10
976	!
977	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN
978	!
979	nmoymltrdbio = nmoymltrdbio + 1
980
981	! ... Trends associated with the time mean of the ML passive tracers
982	tmltrd_sum_bio (:,:,:) = tmltrd_sum_bio (:,:,:) + tmltrd_bio (:,:,:)
983	tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
984	!
985	END IF
986
987	! ======================================================================
988	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
989	! ======================================================================
990
991	! Convert to appropriate physical units
992	tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc
993
994	MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd
995	!
996	zfn = float(nmoymltrdbio) ; zfn2 = zfn * zfn
997
998	! III.1 Prepare fields for output ("instantaneous" diagnostics)
999	! -------------------------------------------------------------
1000
1001	#if defined key_diainstant
1002	STOP 'tmltrd_bio : key_diainstant was never checked within trdmxl. Comment this to proceed.'
1003	#endif
1004	! III.2 Prepare fields for output ("mean" diagnostics)
1005	! ----------------------------------------------------
1006
1007	ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)
1008
1009	!-- Lateral boundary conditions
1010	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
1011	! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut
1012	DO jn = 1, jpdiabio
1013	CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )
1014	ENDDO
1015	ENDIF
1016
1017	IF( lldebug ) THEN
1018	!
1019	WRITE(numout,*) 'trd_mxl_bio : write trends in the Mixed Layer for debugging process:'
1020	WRITE(numout,*) '~~~~~~~~~~~ '
1021	WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio
1022	WRITE(numout,*)
1023
1024	DO jl = 1, jpdiabio
1025	IF( ln_trdmxl_trc_instant ) THEN
1026	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1027	& ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))
1028	ELSE
1029	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1030	& ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))
1031	endif
1032	END DO
1033
1034	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
1035	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
1036	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
1037	WRITE(numout,*)
1038	!
1039	ENDIF
1040
1041	! III.3 Time evolution array swap
1042	! -------------------------------
1043
1044	! For passive tracer mean diagnostics
1045	tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)
1046
1047	! III.4 Convert to appropriate physical units
1048	! -------------------------------------------
1049	ztmltrdbio2 (:,:,:) = ztmltrdbio2 (:,:,:) * rn_ucf_trc/zfn2
1050
1051	END IF MODULO_NTRD
1052
1053	! ======================================================================
1054	! IV. Write trends in the NetCDF file
1055	! ======================================================================
1056
1057	! IV.1 Code for dimg mpp output
1058	! -----------------------------
1059
1060	# if defined key_dimgout
1061	STOP 'Not implemented'
1062	# else
1063
1064	! IV.2 Code for IOIPSL/NetCDF output
1065	! ----------------------------------
1066
1067	! define time axis
1068	itmod = kt - nittrc000 + 1
1069	it = kt
1070
1071	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
1072	WRITE(numout,*) ' '
1073	WRITE(numout,*) 'trd_mxl_bio : write ML bio trends in the NetCDF file :'
1074	WRITE(numout,*) '~~~~~~~~~~~ '
1075	WRITE(numout,*) ' ', TRIM(clhstnam), ' at kt = ', kt
1076	WRITE(numout,*) ' N.B. nmoymltrdbio = ', nmoymltrdbio
1077	WRITE(numout,*) ' '
1078	END IF
1079
1080
1081	! 2. Start writing data
1082	! ---------------------
1083
1084	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
1085	!
1086	DO jl = 1, jpdiabio
1087	CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) , &
1088	& it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
1089	END DO
1090
1091
1092	IF( kt == nitend ) CALL histclo( nidtrdbio )
1093
1094	ELSE ! <<< write the trends for passive tracer mean diagnostics
1095
1096	DO jl = 1, jpdiabio
1097	CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) , &
1098	& it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
1099	END DO
1100
1101	IF( kt == nitend ) CALL histclo( nidtrdbio )
1102	!
1103	END IF NETCDF_OUTPUT
1104
1105	! Compute the control surface (for next time step) : flag = on
1106	icountbio = 1
1107
1108
1109	# endif /* key_dimgout */
1110
1111	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
1112	!
1113	! III.5 Reset cumulative arrays to zero
1114	! -------------------------------------
1115	nmoymltrdbio = 0
1116	tmltrd_csum_ln_bio (:,:,:) = 0.e0
1117	tmltrd_sum_bio (:,:,:) = 0.e0
1118	END IF
1119
1120	! ======================================================================
1121	! Write restart file
1122	! ======================================================================
1123
1124	! restart write is done in trd_mxl_trc_write which is called by trd_mxl_bio (Marina)
1125	!
1126	#endif
1127	END SUBROUTINE trd_mxl_bio
1128
1129
1130	REAL FUNCTION sum2d( ztab )
1131	!!----------------------------------------------------------------------
1132	!! CD ??? prevoir d'utiliser plutot prtctl
1133	!!----------------------------------------------------------------------
1134	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: ztab
1135	!!----------------------------------------------------------------------
1136	sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
1137	END FUNCTION sum2d
1138
1139
1140	SUBROUTINE trd_mxl_trc_init
1141	!!----------------------------------------------------------------------
1142	!! * ROUTINE trd_mxl_init *
1143	!!
1144	!! ** Purpose : computation of vertically integrated T and S budgets
1145	!! from ocean surface down to control surface (NetCDF output)
1146	!!
1147	!! ** Method/usage :
1148	!!
1149	!!----------------------------------------------------------------------
1150	INTEGER :: inum ! logical unit
1151	INTEGER :: ilseq, jl, jn, iiter
1152	REAL(wp) :: zjulian, zsto, zout
1153	CHARACTER (LEN=40) :: clop
1154	CHARACTER (LEN=15) :: csuff
1155	CHARACTER (LEN=12) :: clmxl
1156	CHARACTER (LEN=16) :: cltrcu
1157	CHARACTER (LEN=10) :: clvar
1158
1159	!!----------------------------------------------------------------------
1160
1161	! ======================================================================
1162	! I. initialization
1163	! ======================================================================
1164
1165	IF(lwp) THEN
1166	WRITE(numout,*)
1167	WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers '
1168	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
1169	WRITE(numout,*)
1170	ENDIF
1171
1172
1173	! I.1 Check consistency of user defined preferences
1174	! -------------------------------------------------
1175
1176	IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
1177	WRITE(numout,cform_err)
1178	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
1179	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
1180	WRITE(numout,*) ' you defined, nn_trd_trc = ', nn_trd_trc
1181	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
1182	WRITE(numout,*) ' You should reconsider this choice. '
1183	WRITE(numout,*)
1184	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
1185	WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) '
1186	nstop = nstop + 1
1187	ENDIF
1188
1189	! * Debugging information *
1190	IF( lldebug ) THEN
1191	WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl
1192	WRITE(numout,*) ' ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
1193	ENDIF
1194
1195	IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
1196	WRITE(numout,cform_err)
1197	WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL '
1198	WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '
1199	WRITE(numout,*) ' WHY? Everything in trdmxl_trc is coded for leap-frog, and '
1200	WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note '
1201	WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). '
1202	WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '
1203	WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'
1204	WRITE(numout,*)
1205	nstop = nstop + 1
1206	ENDIF
1207
1208	! I.2 Initialize arrays to zero or read a restart file
1209	! ----------------------------------------------------
1210	nmoymltrd = 0
1211
1212	rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst.
1213	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
1214	tmlradm_trc (:,:,:) = 0.e0
1215
1216	tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean
1217	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0
1218
1219	#if defined key_pisces_reduced
1220	nmoymltrdbio = 0
1221	tmltrd_sum_bio (:,:,:) = 0.e0 ; tmltrd_csum_ln_bio (:,:,:) = 0.e0
1222	DO jl = 1, jp_pisces_trd
1223	ctrd_bio(jl,1) = ctrbil(jl) ! long name
1224	ctrd_bio(jl,2) = ctrbio(jl) ! short name
1225	ENDDO
1226	#endif
1227
1228	IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
1229	CALL trd_mxl_trc_rst_read
1230	ELSE
1231	tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst.
1232	tmlbn_trc (:,:,:) = 0.e0
1233
1234	tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean
1235	tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0
1236	#if defined key_pisces_reduced
1237	tmltrd_csum_ub_bio (:,:,:) = 0.e0
1238	#endif
1239
1240	ENDIF
1241
1242	icount = 1 ; ionce = 1 ! open specifier
1243
1244	#if defined key_pisces_reduced
1245	icountbio = 1 ; ioncebio = 1 ! open specifier
1246	#endif
1247
1248	! I.3 Read control surface from file ctlsurf_idx
1249	! ----------------------------------------------
1250	IF( nn_ctls_trc == 1 ) THEN
1251	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
1252	READ ( inum ) nbol_trc
1253	CLOSE( inum )
1254	ENDIF
1255
1256	! ======================================================================
1257	! II. netCDF output initialization
1258	! ======================================================================
1259
1260	#if defined key_dimgout
1261	???
1262	#else
1263	! clmxl = legend root for netCDF output
1264	IF( nn_ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion
1265	clmxl = 'Mixed Layer '
1266	ELSE IF( nn_ctls_trc == 1 ) THEN ! control surface = read index from file
1267	clmxl = ' Bowl '
1268	ELSE IF( nn_ctls_trc >= 2 ) THEN ! control surface = model level
1269	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
1270	ENDIF
1271
1272	! II.1 Define frequency of output and means
1273	! -----------------------------------------
1274	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
1275	ELSE ; clop = "x" ! no use of the mask value (require less cp time)
1276	ENDIF
1277	# if defined key_diainstant
1278	IF( .NOT. ln_trdmxl_trc_instant ) THEN
1279	STOP 'trd_mxl_trc : this was never checked. Comment this line to proceed...'
1280	ENDIF
1281	zsto = nn_trd_trc * rdt
1282	clop = "inst("//TRIM(clop)//")"
1283	# else
1284	IF( ln_trdmxl_trc_instant ) THEN
1285	zsto = rdt ! inst. diags : we use IOIPSL time averaging
1286	ELSE
1287	zsto = nn_trd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging
1288	ENDIF
1289	clop = "ave("//TRIM(clop)//")"
1290	# endif
1291	zout = nn_trd_trc * rdt
1292	iiter = ( nittrc000 - 1 ) / nn_dttrc
1293
1294	IF(lwp) WRITE (numout,*) ' netCDF initialization'
1295
1296	! II.2 Compute julian date from starting date of the run
1297	! ------------------------------------------------------
1298	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
1299	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
1300	IF(lwp) WRITE(numout,*)' '
1301	IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000 &
1302	& ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday &
1303	& ,'Julian day : ', zjulian
1304
1305	! II.3 Define the T grid trend file (nidtrd)
1306	! ------------------------------------------
1307
1308	!-- Define long and short names for the NetCDF output variables
1309	! ==> choose them according to trdmxl_trc_oce.F90 <==
1310
1311	ctrd_trc(jpmxl_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmxl_trc_xad ,2) = "_xad"
1312	ctrd_trc(jpmxl_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmxl_trc_yad ,2) = "_yad"
1313	ctrd_trc(jpmxl_trc_zad ,1) = " Vertical advection" ; ctrd_trc(jpmxl_trc_zad ,2) = "_zad"
1314	ctrd_trc(jpmxl_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmxl_trc_ldf ,2) = "_ldf"
1315	ctrd_trc(jpmxl_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmxl_trc_zdf ,2) = "_zdf"
1316	ctrd_trc(jpmxl_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmxl_trc_bbl ,2) = "_bbl"
1317	ctrd_trc(jpmxl_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmxl_trc_dmp ,2) = "_dmp"
1318	ctrd_trc(jpmxl_trc_sbc ,1) = " Surface boundary cond." ; ctrd_trc(jpmxl_trc_sbc ,2) = "_sbc"
1319	ctrd_trc(jpmxl_trc_sms, 1) = " Sources minus sinks" ; ctrd_trc(jpmxl_trc_sms ,2) = "_sms"
1320	ctrd_trc(jpmxl_trc_radb ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radb ,2) = "_radb"
1321	ctrd_trc(jpmxl_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radn ,2) = "_radn"
1322	ctrd_trc(jpmxl_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmxl_trc_atf ,2) = "_atf"
1323
1324	DO jn = 1, jptra
1325	!-- Create a NetCDF file and enter the define mode
1326	IF( ln_trdtrc(jn) ) THEN
1327	csuff="ML_"//ctrcnm(jn)
1328	CALL dia_nam( clhstnam, nn_trd_trc, csuff )
1329	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1330	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
1331
1332	!-- Define the ML depth variable
1333	CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m", &
1334	& jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1335
1336	ENDIF
1337	END DO
1338
1339	#if defined key_pisces_reduced
1340	!-- Create a NetCDF file and enter the define mode
1341	CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )
1342	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1343	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom, snc4chunks=snc4set )
1344	#endif
1345
1346	!-- Define physical units
1347	IF( rn_ucf_trc == 1. ) THEN
1348	cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit
1349	ELSEIF ( rn_ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3)
1350	cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf
1351	ELSE
1352	cltrcu = "unknown?"
1353	ENDIF
1354
1355	!-- Define miscellaneous passive tracer mixed-layer variables
1356	IF( jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
1357	STOP 'Error : jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb' ! see below
1358	ENDIF
1359
1360	DO jn = 1, jptra
1361	!
1362	IF( ln_trdtrc(jn) ) THEN
1363	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc.
1364	CALL histdef(nidtrd(jn), trim(clvar), clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", &
1365	& "mmole-N/m3", jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1366	CALL histdef(nidtrd(jn), trim(clvar)//"_tot" , clmxl//" "//trim(ctrcnm(jn))//" Total trend ", &
1367	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1368	CALL histdef(nidtrd(jn), trim(clvar)//"_res" , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)", &
1369	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1370
1371	DO jl = 1, jpltrd_trc - 2 ! <== only true if jpltrd_trc == jpmxl_trc_atf
1372	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1), &
1373	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1374	END DO ! if zsto=rdt above
1375
1376	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), &
1377	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1378
1379	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1), &
1380	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1381	!
1382	ENDIF
1383	END DO
1384
1385	#if defined key_pisces_reduced
1386	DO jl = 1, jp_pisces_trd
1387	CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1)) , &
1388	& cltrcu, jpi, jpj, nh_tb, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1389	END DO ! if zsto=rdt above
1390	#endif
1391
1392	!-- Leave IOIPSL/NetCDF define mode
1393	DO jn = 1, jptra
1394	IF( ln_trdtrc(jn) ) CALL histend( nidtrd(jn), snc4set )
1395	END DO
1396
1397	#if defined key_pisces_reduced
1398	!-- Leave IOIPSL/NetCDF define mode
1399	CALL histend( nidtrdbio, snc4set )
1400
1401	IF(lwp) WRITE(numout,*)
1402	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
1403	#endif
1404
1405	#endif /* key_dimgout */
1406	END SUBROUTINE trd_mxl_trc_init
1407
1408	#else
1409	!!----------------------------------------------------------------------
1410	!! Default option : Empty module
1411	!!----------------------------------------------------------------------
1412	CONTAINS
1413	SUBROUTINE trd_mxl_trc( kt ) ! Empty routine
1414	INTEGER, INTENT( in) :: kt
1415	WRITE(,) 'trd_mxl_trc: You should not have seen this print! error?', kt
1416	END SUBROUTINE trd_mxl_trc
1417	SUBROUTINE trd_mxl_bio( kt )
1418	INTEGER, INTENT( in) :: kt
1419	WRITE(,) 'trd_mxl_bio: You should not have seen this print! error?', kt
1420	END SUBROUTINE trd_mxl_bio
1421	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
1422	INTEGER , INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
1423	CHARACTER(len=2) , INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
1424	REAL, DIMENSION(:,:,:), INTENT( in ) :: ptrc_trdmxl ! passive trc trend
1425	WRITE(,) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
1426	WRITE(,) ' " " : You should not have seen this print! error?', ctype
1427	WRITE(,) ' " " : You should not have seen this print! error?', ktrd
1428	WRITE(,) ' " " : You should not have seen this print! error?', kjn
1429	END SUBROUTINE trd_mxl_trc_zint
1430	SUBROUTINE trd_mxl_trc_init ! Empty routine
1431	WRITE(,) 'trd_mxl_trc_init: You should not have seen this print! error?'
1432	END SUBROUTINE trd_mxl_trc_init
1433	#endif
1434
1435	!!======================================================================
1436	END MODULE trdmxl_trc

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

Download in other formats: