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.

trdmld_trc.F90 in trunk/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:keywords set to `Id`
File size: 69.0 KB

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

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

Download in other formats: