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

Context Navigation

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 6596

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

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

Download in other formats: