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

Context Navigation

source: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 7351

Last change on this file since 7351 was 7351, checked in by emanuelaclementi, 7 years ago
ticket #1805 step 3: /2016/dev_INGV_UKMO_2016 aligned to the trunk at revision 7161
Property svn:keywords set to `Id`
File size: 68.0 KB

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

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

Download in other formats: