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 branches/2012/dev_LOCEAN_2012/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: branches/2012/dev_LOCEAN_2012/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 3584

Last change on this file since 3584 was 3584, checked in by cetlod, 11 years ago
Add in branch 2012/dev_LOCEAN_2012 changes from dev_r3438_LOCEAN15_PISLOB & dev_r3387_LOCEAN6_AGRIF_LIM, see ticket 1000
Property svn:keywords set to `Id`
File size: 68.4 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 wrk_nemo ! Memory allocation
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_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_mld_bio()
65	#endif
66
67	!! * Substitutions
68	# include "top_substitute.h90"
69	# include "zdfddm_substitute.h90"
70	!!----------------------------------------------------------------------
71	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
72	!! $Header: $
73	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
74	!!----------------------------------------------------------------------
75	CONTAINS
76
77	INTEGER FUNCTION trd_mld_trc_alloc()
78	!!----------------------------------------------------------------------
79	!! * ROUTINE trd_mld_trc_alloc *
80	!!----------------------------------------------------------------------
81	ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) , &
82	#if defined key_pisces_reduced
83	& ztmltrdbio2(jpi,jpj,jpdiabio) , &
84	#endif
85	& ndextrd1(jpi*jpj) , STAT=trd_mld_trc_alloc)
86	!
87	IF( lk_mpp ) CALL mpp_sum ( trd_mld_trc_alloc )
88	IF( trd_mld_trc_alloc /=0 ) CALL ctl_warn('trd_mld_trc_alloc: failed to allocate arrays')
89	!
90	END FUNCTION trd_mld_trc_alloc
91
92
93	SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
94	!!----------------------------------------------------------------------
95	!! * ROUTINE trd_mld_trc_zint *
96	!!
97	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
98	!! to the subroutine. This vertical average is performed from ocean
99	!! surface down to a chosen control surface.
100	!!
101	!! ** Method/usage :
102	!! The control surface can be either a mixed layer depth (time varying)
103	!! or a fixed surface (jk level or bowl).
104	!! Choose control surface with nctls_trc in namelist NAMTRD :
105	!! nctls_trc = -2 : use isopycnal surface
106	!! nctls_trc = -1 : use euphotic layer with light criterion
107	!! nctls_trc = 0 : use mixed layer with density criterion
108	!! nctls_trc = 1 : read index from file 'ctlsurf_idx'
109	!! nctls_trc > 1 : use fixed level surface jk = nctls_trc
110	!! Note: in the remainder of the routine, the volume between the
111	!! surface and the control surface is called "mixed-layer"
112	!!----------------------------------------------------------------------
113	!!
114	INTEGER, INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
115	CHARACTER(len=2), INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
116	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: ptrc_trdmld ! passive tracer trend
117	!
118	INTEGER :: ji, jj, jk, isum
119	REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
120	!!----------------------------------------------------------------------
121
122	CALL wrk_alloc( jpi, jpj, zvlmsk )
123
124	! I. Definition of control surface and integration weights
125	! --------------------------------------------------------
126
127	ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
128	!
129	tmltrd_trc(:,:,:,:) = 0.e0 ! <<< reset trend arrays to zero
130
131	! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
132	SELECT CASE ( nn_ctls_trc ) ! choice of the control surface
133	CASE ( -2 ) ; STOP 'trdmld_trc : not ready ' ! -> isopycnal surface (see ???)
134	#if defined key_pisces \|\| defined key_pisces_reduced
135	CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion
136	#endif
137	CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl)
138	CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file
139	CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
140	nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level
141	END SELECT
142
143	! ... Compute ndextrd1 and ndimtrd1 ??? role de jpktrd_trc
144	IF( ionce == 1 ) THEN
145	!
146	isum = 0 ; zvlmsk(:,:) = 0.e0
147
148	IF( jpktrd_trc < jpk ) THEN ! description ???
149	DO jj = 1, jpj
150	DO ji = 1, jpi
151	IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
152	zvlmsk(ji,jj) = tmask(ji,jj,1)
153	ELSE
154	isum = isum + 1
155	zvlmsk(ji,jj) = 0.e0
156	ENDIF
157	END DO
158	END DO
159	ENDIF
160
161	IF( isum > 0 ) THEN ! index of ocean points (2D only)
162	WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
163	CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
164	ELSE
165	CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
166	ENDIF
167
168	ionce = 0 ! no more pass here
169	!
170	ENDIF ! ionce == 1
171
172	! ... Weights for vertical averaging
173	wkx_trc(:,:,:) = 0.e0
174	DO jk = 1, jpktrd_trc ! initialize wkx_trc with vertical scale factor in mixed-layer
175	DO jj = 1, jpj
176	DO ji = 1, jpi
177	IF( jk - nmld_trc(ji,jj) < 0 ) wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
178	END DO
179	END DO
180	END DO
181
182	rmld_trc(:,:) = 0.e0
183	DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc
184	rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
185	END DO
186
187	DO jk = 1, jpktrd_trc ! compute integration weights
188	wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
189	END DO
190
191	icount = 0 ! <<< flag = off : control surface & integr. weights
192	! ! computed only once per time step
193	ENDIF ONCE_PER_TIME_STEP
194
195	! II. Vertical integration of trends in the mixed-layer
196	! -----------------------------------------------------
197
198	SELECT CASE ( ctype )
199	CASE ( '3D' ) ! mean passive tracer trends in the mixed-layer
200	DO jk = 1, jpktrd_trc
201	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,jk) * wkx_trc(:,:,jk)
202	END DO
203	CASE ( '2D' ) ! forcing at upper boundary of the mixed-layer
204	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,1) * wkx_trc(:,:,1) ! non penetrative
205	END SELECT
206	!
207	CALL wrk_dealloc( jpi, jpj, zvlmsk )
208	!
209	END SUBROUTINE trd_mld_trc_zint
210
211
212	SUBROUTINE trd_mld_bio_zint( ptrc_trdmld, ktrd )
213	!!----------------------------------------------------------------------
214	!! * ROUTINE trd_mld_bio_zint *
215	!!
216	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
217	!! to the subroutine. This vertical average is performed from ocean
218	!! surface down to a chosen control surface.
219	!!
220	!! ** Method/usage :
221	!! The control surface can be either a mixed layer depth (time varying)
222	!! or a fixed surface (jk level or bowl).
223	!! Choose control surface with nctls in namelist NAMTRD :
224	!! nctls_trc = 0 : use mixed layer with density criterion
225	!! nctls_trc = 1 : read index from file 'ctlsurf_idx'
226	!! nctls_trc > 1 : use fixed level surface jk = nctls_trc
227	!! Note: in the remainder of the routine, the volume between the
228	!! surface and the control surface is called "mixed-layer"
229	!!----------------------------------------------------------------------
230	!!
231	INTEGER , INTENT(in) :: ktrd ! bio trend index
232	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrc_trdmld ! passive trc trend
233	#if defined key_pisces_reduced
234	!
235	INTEGER :: ji, jj, jk, isum
236	REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
237	!!----------------------------------------------------------------------
238
239	CALL wrk_alloc( jpi, jpj, zvlmsk )
240
241	! I. Definition of control surface and integration weights
242	! --------------------------------------------------------
243	! ==> only once per time step <==
244
245	IF( icountbio == 1 ) THEN
246	!
247	tmltrd_bio(:,:,:) = 0.e0 ! <<< reset trend arrays to zero
248	! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
249	SELECT CASE ( nn_ctls_trc ) ! choice of the control surface
250	CASE ( -2 ) ; STOP 'trdmld_trc : not ready ' ! -> isopycnal surface (see ???)
251	CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion
252	CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl)
253	CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file
254	CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
255	nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level
256	END SELECT
257
258	! ... Compute ndextrd1 and ndimtrd1 only once
259	IF( ioncebio == 1 ) THEN
260	!
261	! Check of validity : nmld_trc(ji,jj) <= jpktrd_trc
262	isum = 0
263	zvlmsk(:,:) = 0.e0
264
265	IF( jpktrd_trc < jpk ) THEN
266	DO jj = 1, jpj
267	DO ji = 1, jpi
268	IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
269	zvlmsk(ji,jj) = tmask(ji,jj,1)
270	ELSE
271	isum = isum + 1
272	zvlmsk(ji,jj) = 0.
273	END IF
274	END DO
275	END DO
276	END IF
277
278	! Index of ocean points (2D only)
279	IF( isum > 0 ) THEN
280	WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
281	CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
282	ELSE
283	CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
284	END IF
285
286	ioncebio = 0 ! no more pass here
287	!
288	END IF ! ( ioncebio == 1 )
289
290	! ... Weights for vertical averaging
291	wkx_trc(:,:,:) = 0.e0
292	DO jk = 1, jpktrd_trc ! initialize wkx_trc with vertical scale factor in mixed-layer
293	DO jj = 1,jpj
294	DO ji = 1,jpi
295	IF( jk - nmld_trc(ji,jj) < 0. ) wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
296	END DO
297	END DO
298	END DO
299
300	rmld_trc(:,:) = 0.
301	DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc
302	rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
303	END DO
304
305	DO jk = 1, jpktrd_trc ! compute integration weights
306	wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
307	END DO
308
309	icountbio = 0 ! <<< flag = off : control surface & integr. weights
310	! ! computed only once per time step
311	END IF ! ( icountbio == 1 )
312
313	! II. Vertical integration of trends in the mixed-layer
314	! -----------------------------------------------------
315
316
317	DO jk = 1, jpktrd_trc
318	tmltrd_bio(:,:,ktrd) = tmltrd_bio(:,:,ktrd) + ptrc_trdmld(:,:,jk) * wkx_trc(:,:,jk)
319	END DO
320
321	CALL wrk_alloc( jpi, jpj, zvlmsk )
322	#endif
323	!
324	END SUBROUTINE trd_mld_bio_zint
325
326
327	SUBROUTINE trd_mld_trc( kt )
328	!!----------------------------------------------------------------------
329	!! * ROUTINE trd_mld_trc *
330	!!
331	!! ** Purpose : Compute and cumulate the mixed layer trends over an analysis
332	!! period, and write NetCDF (or dimg) outputs.
333	!!
334	!! ** Method/usage :
335	!! The stored trends can be chosen twofold (according to the ln_trdmld_trc_instant
336	!! logical namelist variable) :
337	!! 1) to explain the difference between initial and final
338	!! mixed-layer T & S (where initial and final relate to the
339	!! current analysis window, defined by ntrc_trc in the namelist)
340	!! 2) to explain the difference between the current and previous
341	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
342	!! performed over each analysis window).
343	!!
344	!! ** Consistency check :
345	!! If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
346	!! entrainment) should be zero, at machine accuracy. Note that in the case
347	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
348	!! over the first two analysis windows (except if restart).
349	!! N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8
350	!! for checking residuals.
351	!! On a NEC-SX5 computer, this typically leads to:
352	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmld_trc_instant=.false.
353	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmld_trc_instant=.true.
354	!!
355	!! ** Action :
356	!! At each time step, mixed-layer averaged trends are stored in the
357	!! tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
358	!! This array is known when trd_mld is called, at the end of the stp subroutine,
359	!! except for the purely vertical K_z diffusion term, which is embedded in the
360	!! lateral diffusion trend.
361	!!
362	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
363	!! from the lateral diffusion trend.
364	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
365	!! arrays are updated.
366	!! In III), called only once per analysis window, we compute the total trends,
367	!! along with the residuals and the Asselin correction terms.
368	!! In IV), the appropriate trends are written in the trends NetCDF file.
369	!!
370	!! References :
371	!! - Vialard & al.
372	!! - See NEMO documentation (in preparation)
373	!!----------------------------------------------------------------------
374	!
375	INTEGER, INTENT(in) :: kt ! ocean time-step index
376	!
377	INTEGER :: ji, jj, jk, jl, ik, it, itmod, jn
378	REAL(wp) :: zavt, zfn, zfn2
379	!
380	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltot ! d(trc)/dt over the anlysis window (incl. Asselin)
381	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlres ! residual = dh/dt entrainment term
382	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlatf ! for storage only
383	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlrad ! for storage only (for trb<0 corr in trcrad)
384	!
385	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltot2 ! -+
386	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlres2 ! \| working arrays to diagnose the trends
387	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltrdm2 ! \| associated with the time meaned ML
388	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlatf2 ! \| passive tracers
389	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlrad2 ! \| (-> for trb<0 corr in trcrad)
390	!
391	CHARACTER (LEN=10) :: clvar
392	#if defined key_dimgout
393	INTEGER :: iyear,imon,iday
394	CHARACTER(LEN=80) :: cltext, clmode
395	#endif
396	!!----------------------------------------------------------------------
397
398	! Set-up pointers into sub-arrays of workspaces
399	CALL wrk_alloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
400	CALL wrk_alloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
401
402	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
403
404	! ======================================================================
405	! I. Diagnose the purely vertical (K_z) diffusion trend
406	! ======================================================================
407
408	! ... These terms can be estimated by flux computation at the lower boundary of the ML
409	! (we compute (-1/h) * K_z * d_z( tracer ))
410
411	IF( ln_trcldf_iso ) THEN
412	!
413	DO jj = 1,jpj
414	DO ji = 1,jpi
415	ik = nmld_trc(ji,jj)
416	zavt = fsavs(ji,jj,ik)
417	DO jn = 1, jptra
418	IF( ln_trdtrc(jn) ) &
419	tmltrd_trc(ji,jj,jpmld_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) &
420	& * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) ) &
421	& / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
422	END DO
423	END DO
424	END DO
425
426	DO jn = 1, jptra
427	! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
428	IF( ln_trdtrc(jn) ) &
429	tmltrd_trc(:,:,jpmld_trc_ldf,jn) = tmltrd_trc(:,:,jpmld_trc_ldf,jn) - tmltrd_trc(:,:,jpmld_trc_zdf,jn)
430
431	END DO
432	!
433	ENDIF
434
435	#if ! defined key_gyre
436	! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm.
437	! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
438	DO jn = 1, jptra
439	IF( ln_trdtrc(jn) ) THEN
440	DO jl = 1, jpltrd_trc
441	CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions
442	END DO
443	ENDIF
444	END DO
445	#endif
446	! ======================================================================
447	! II. Cumulate the trends over the analysis window
448	! ======================================================================
449
450	ztmltrd2(:,:,:,:) = 0.e0 ; ztmltot2(:,:,:) = 0.e0 ! <<< reset arrays to zero
451	ztmlres2(:,:,:) = 0.e0 ; ztmlatf2(:,:,:) = 0.e0
452	ztmlrad2(:,:,:) = 0.e0
453
454	! II.1 Set before values of vertically averages passive tracers
455	! -------------------------------------------------------------
456	IF( kt > nittrc000 ) THEN
457	DO jn = 1, jptra
458	IF( ln_trdtrc(jn) ) THEN
459	tmlb_trc (:,:,jn) = tml_trc (:,:,jn)
460	tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_atf,jn)
461	tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_radb,jn)
462	ENDIF
463	END DO
464	ENDIF
465
466	! II.2 Vertically averaged passive tracers
467	! ----------------------------------------
468	tml_trc(:,:,:) = 0.e0
469	DO jk = 1, jpktrd_trc ! - 1 ???
470	DO jn = 1, jptra
471	IF( ln_trdtrc(jn) ) &
472	tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
473	END DO
474	END DO
475
476	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
477	! ------------------------------------------------------------------------
478	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1) ???
479	!
480	DO jn = 1, jptra
481	IF( ln_trdtrc(jn) ) THEN
482	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
483	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
484
485	tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0 ; tmltrd_atf_sumb_trc (:,:,jn) = 0.e0
486	tmltrd_rad_sumb_trc (:,:,jn) = 0.e0
487	ENDIF
488	END DO
489
490	rmldbn_trc(:,:) = rmld_trc(:,:)
491	!
492	ENDIF
493
494	! II.4 Cumulated trends over the analysis period
495	! ----------------------------------------------
496	!
497	! [ 1rst analysis window ] [ 2nd analysis window ]
498	!
499	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
500	! ntrd 2*ntrd etc.
501	! 1 2 3 4 =5 e.g. =10
502	!
503	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN ! ???
504	!
505	nmoymltrd = nmoymltrd + 1
506
507
508	! ... Cumulate over BOTH physical contributions AND over time steps
509	DO jn = 1, jptra
510	IF( ln_trdtrc(jn) ) THEN
511	DO jl = 1, jpltrd_trc
512	tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
513	END DO
514	ENDIF
515	END DO
516
517	DO jn = 1, jptra
518	IF( ln_trdtrc(jn) ) THEN
519	! ... Special handling of the Asselin trend
520	tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
521	tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
522
523	! ... Trends associated with the time mean of the ML passive tracers
524	tmltrd_sum_trc (:,:,:,jn) = tmltrd_sum_trc (:,:,:,jn) + tmltrd_trc (:,:,:,jn)
525	tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
526	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) + tml_trc (:,:,jn)
527	ENDIF
528	ENDDO
529
530	rmld_sum_trc (:,:) = rmld_sum_trc (:,:) + rmld_trc (:,:)
531	!
532	ENDIF
533
534	! ======================================================================
535	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
536	! ======================================================================
537
538	! Convert to appropriate physical units
539	tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
540
541	itmod = kt - nittrc000 + 1
542	it = kt
543
544	MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of nn_trd_trc
545	!
546	ztmltot (:,:,:) = 0.e0 ! reset arrays to zero
547	ztmlres (:,:,:) = 0.e0
548	ztmltot2(:,:,:) = 0.e0
549	ztmlres2(:,:,:) = 0.e0
550
551	zfn = FLOAT( nmoymltrd ) ; zfn2 = zfn * zfn
552
553	! III.1 Prepare fields for output ("instantaneous" diagnostics)
554	! -------------------------------------------------------------
555
556	DO jn = 1, jptra
557	IF( ln_trdtrc(jn) ) THEN
558	!-- Compute total trends (use rdttrc instead of rdt ???)
559	IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes
560	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
561	ELSE ! LEAP-FROG schemes
562	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
563	ENDIF
564
565	!-- Compute residuals
566	ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
567	& - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) )
568
569	!-- Diagnose Asselin trend over the analysis window
570	ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
571	ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
572
573	!-- Lateral boundary conditions
574	#if ! defined key_gyre
575
576	CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
577	CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
578
579	#endif
580
581	#if defined key_diainstant
582	STOP 'tmltrd_trc : key_diainstant was never checked within trdmld. Comment this to proceed.'
583	#endif
584	ENDIF
585	END DO
586
587	! III.2 Prepare fields for output ("mean" diagnostics)
588	! ----------------------------------------------------
589
590	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
591	rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
592
593	!-- Compute passive tracer total trends
594	DO jn = 1, jptra
595	IF( ln_trdtrc(jn) ) THEN
596	tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
597	ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rdt ) ! now tracer unit is /sec
598	ENDIF
599	END DO
600
601	!-- Compute passive tracer residuals
602	DO jn = 1, jptra
603	IF( ln_trdtrc(jn) ) THEN
604	!
605	DO jl = 1, jpltrd_trc
606	ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
607	END DO
608
609	ztmltrdm2(:,:,jn) = 0.e0
610	DO jl = 1, jpltrd_trc
611	ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
612	END DO
613
614	ztmlres2(:,:,jn) = ztmltot2(:,:,jn) - &
615	& ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
616	& - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
617	!
618
619	!-- Diagnose Asselin trend over the analysis window
620	ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) &
621	& + tmltrd_atf_sumb_trc(:,:,jn)
622	ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) &
623	& + tmltrd_rad_sumb_trc(:,:,jn)
624
625	!-- Lateral boundary conditions
626	#if ! defined key_gyre
627	CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
628	CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
629	DO jl = 1, jpltrd_trc
630	CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file
631	END DO
632	#endif
633	ENDIF
634	END DO
635
636	! * Debugging information *
637	IF( lldebug ) THEN
638	!
639	WRITE(numout,*) 'trd_mld_trc : write trends in the Mixed Layer for debugging process:'
640	WRITE(numout,*) '~~~~~~~~~~~ '
641	WRITE(numout,*)
642	WRITE(numout,*) 'TRC kt = ', kt, ' nmoymltrd = ', nmoymltrd
643
644	DO jn = 1, jptra
645
646	IF( ln_trdtrc(jn) ) THEN
647	WRITE(numout, *)
648	WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
649
650	WRITE(numout, *)
651	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres : ', SUM2D(ztmlres(:,:,jn))
652	!CD??? PREVOIR: z2d = ztmlres(:,:,jn) ; CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres - : ')
653
654	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
655	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
656	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot : ', SUM2D(+ztmltot (:,:,jn))
657	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
658	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
659	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
660	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
661	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
662	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
663
664	WRITE(numout, *)
665	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2 : ', SUM2D(ztmlres2(:,:,jn))
666	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
667	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
668	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2 : ', SUM2D(+ztmltot2(:,:,jn))
669	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2 : ', SUM2D(+ztmltrdm2(:,:,jn))
670	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_atf,jn))
671	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
672	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_radb,jn))
673	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
674	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
675
676	WRITE(numout, *)
677	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot : ', SUM2D(ztmltot (:,:,jn))
678	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
679	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc : ', SUM2D(tml_trc (:,:,jn))
680	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc : ', SUM2D(tmlbn_trc (:,:,jn))
681	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc : ', SUM2D(tmlb_trc (:,:,jn))
682	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc : ', SUM2D(tmlbb_trc (:,:,jn))
683	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
684
685	WRITE(numout, *)
686	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
687	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
688	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
689	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
690	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
691
692	WRITE(numout, *)
693	DO jl = 1, jpltrd_trc
694	WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmld_trc_xxx = ', jl, &
695	& ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
696	END DO
697
698	WRITE(numout,*)
699	WRITE(numout,) ' ******************** ZTMLRES, ZTMLRES2 ********************* '
700	WRITE(numout,*)
701	WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
702	WRITE(numout,*)
703	WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
704
705	WRITE(numout,*)
706	WRITE(numout,) ' ******************** ZTMLRES ********************* '
707	WRITE(numout,*)
708
709	WRITE(numout,*) '...................................................'
710	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
711	DO jj = 5, 1, -1
712	WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
713	END DO
714
715	WRITE(numout,*)
716	WRITE(numout,) ' ******************** ZTMLRES2 ********************* '
717	WRITE(numout,*)
718
719	WRITE(numout,*) '...................................................'
720	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
721	DO jj = 5, 1, -1
722	WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
723	END DO
724	!
725	ENDIF
726	!
727	END DO
728
729
730	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
731	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
732	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
733	WRITE(numout,*)
734	!
735	ENDIF
736
737	! III.3 Time evolution array swap
738	! -------------------------------
739	! ML depth
740	rmldbn_trc(:,:) = rmld_trc(:,:)
741	rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
742	DO jn = 1, jptra
743	IF( ln_trdtrc(jn) ) THEN
744	! For passive tracer instantaneous diagnostics
745	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
746	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
747
748	! For passive tracer mean diagnostics
749	tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
750	tml_sumb_trc (:,:,jn) = tml_sum_trc(:,:,jn)
751	tmltrd_atf_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn)
752	tmltrd_rad_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmld_trc_radb,jn)
753
754
755	! III.4 Convert to appropriate physical units
756	! -------------------------------------------
757	ztmltot (:,:,jn) = ztmltot (:,:,jn) * rn_ucf_trc/zfn ! instant diags
758	ztmlres (:,:,jn) = ztmlres (:,:,jn) * rn_ucf_trc/zfn
759	ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * rn_ucf_trc/zfn
760	ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * rn_ucf_trc/zfn
761	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags
762	ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * rn_ucf_trc/zfn2
763	ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * rn_ucf_trc/zfn2
764	ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * rn_ucf_trc/zfn2
765	ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * rn_ucf_trc/zfn2
766	ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * rn_ucf_trc/zfn2
767	ENDIF
768	END DO
769	!
770	ENDIF MODULO_NTRD
771
772	! ======================================================================
773	! IV. Write trends in the NetCDF file
774	! ======================================================================
775
776	! IV.1 Code for dimg mpp output
777	! -----------------------------
778
779	# if defined key_dimgout
780	STOP 'Not implemented'
781	# else
782
783	! IV.2 Code for IOIPSL/NetCDF output
784	! ----------------------------------
785
786	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
787	WRITE(numout,*) ' '
788	WRITE(numout,*) 'trd_mld_trc : write passive tracer trends in the NetCDF file :'
789	WRITE(numout,*) '~~~~~~~~~~~ '
790	WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt
791	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
792	WRITE(numout,*) ' '
793	ENDIF
794
795	NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
796	!
797
798	DO jn = 1, jptra
799	!
800	IF( ln_trdtrc(jn) ) THEN
801	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
802	!-- Output the fields
803	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
804	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 )
805	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 )
806	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 )
807
808	DO jl = 1, jpltrd_trc - 2
809	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
810	& it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
811	END DO
812
813	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
814	& it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
815
816	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_atf,2)), & ! now Asselin : jpltrd_trc
817	& it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
818
819	ENDIF
820	END DO
821
822	IF( kt == nitend ) THEN
823	DO jn = 1, jptra
824	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
825	END DO
826	ENDIF
827
828	ELSE ! <<< write the trends for passive tracer mean diagnostics
829
830	DO jn = 1, jptra
831	!
832	IF( ln_trdtrc(jn) ) THEN
833	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 )
834	!-- Output the fields
835	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
836
837	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
838	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 )
839	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 )
840
841	DO jl = 1, jpltrd_trc - 2
842	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
843	& it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
844	END DO
845
846	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
847	& it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
848
849	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(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_pisces_reduced
942	INTEGER :: jl, it, itmod
943	LOGICAL :: llwarn = .TRUE., lldebug = .TRUE.
944	REAL(wp) :: zfn, zfn2
945	#if defined key_dimgout
946	INTEGER :: iyear,imon,iday
947	CHARACTER(LEN=80) :: cltext, clmode
948	#endif
949	!!----------------------------------------------------------------------
950	! ... Warnings
951	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
952
953	! ======================================================================
954	! II. Cumulate the trends over the analysis window
955	! ======================================================================
956
957	ztmltrdbio2(:,:,:) = 0.e0 ! <<< reset arrays to zero
958
959	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
960	! ------------------------------------------------------------------------
961	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
962	!
963	tmltrd_csum_ub_bio (:,:,:) = 0.e0
964	!
965	END IF
966
967	! II.4 Cumulated trends over the analysis period
968	! ----------------------------------------------
969	!
970	! [ 1rst analysis window ] [ 2nd analysis window ]
971	!
972	!
973	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
974	! ntrd 2*ntrd etc.
975	! 1 2 3 4 =5 e.g. =10
976	!
977	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN
978	!
979	nmoymltrdbio = nmoymltrdbio + 1
980
981	! ... Trends associated with the time mean of the ML passive tracers
982	tmltrd_sum_bio (:,:,:) = tmltrd_sum_bio (:,:,:) + tmltrd_bio (:,:,:)
983	tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
984	!
985	END IF
986
987	! ======================================================================
988	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
989	! ======================================================================
990
991	! Convert to appropriate physical units
992	tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc
993
994	MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd
995	!
996	zfn = float(nmoymltrdbio) ; zfn2 = zfn * zfn
997
998	! III.1 Prepare fields for output ("instantaneous" diagnostics)
999	! -------------------------------------------------------------
1000
1001	#if defined key_diainstant
1002	STOP 'tmltrd_bio : key_diainstant was never checked within 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, iiter
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=10) :: 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 .OR. ln_trcadv_muscl2 ) .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	! I.2 Initialize arrays to zero or read a restart file
1208	! ----------------------------------------------------
1209	nmoymltrd = 0
1210
1211	rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst.
1212	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
1213	tmlradm_trc (:,:,:) = 0.e0
1214
1215	tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean
1216	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0
1217
1218	#if defined key_pisces_reduced
1219	nmoymltrdbio = 0
1220	tmltrd_sum_bio (:,:,:) = 0.e0 ; tmltrd_csum_ln_bio (:,:,:) = 0.e0
1221	DO jl = 1, jp_pisces_trd
1222	ctrd_bio(jl,1) = ctrbil(jl) ! long name
1223	ctrd_bio(jl,2) = ctrbio(jl) ! short name
1224	ENDDO
1225	#endif
1226
1227	IF( ln_rsttr .AND. ln_trdmld_trc_restart ) THEN
1228	CALL trd_mld_trc_rst_read
1229	ELSE
1230	tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst.
1231	tmlbn_trc (:,:,:) = 0.e0
1232
1233	tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean
1234	tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0
1235	#if defined key_pisces_reduced
1236	tmltrd_csum_ub_bio (:,:,:) = 0.e0
1237	#endif
1238
1239	ENDIF
1240
1241	icount = 1 ; ionce = 1 ! open specifier
1242
1243	#if defined key_pisces_reduced
1244	icountbio = 1 ; ioncebio = 1 ! open specifier
1245	#endif
1246
1247	! I.3 Read control surface from file ctlsurf_idx
1248	! ----------------------------------------------
1249	IF( nn_ctls_trc == 1 ) THEN
1250	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
1251	READ ( inum ) nbol_trc
1252	CLOSE( inum )
1253	ENDIF
1254
1255	! ======================================================================
1256	! II. netCDF output initialization
1257	! ======================================================================
1258
1259	#if defined key_dimgout
1260	???
1261	#else
1262	! clmxl = legend root for netCDF output
1263	IF( nn_ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion
1264	clmxl = 'Mixed Layer '
1265	ELSE IF( nn_ctls_trc == 1 ) THEN ! control surface = read index from file
1266	clmxl = ' Bowl '
1267	ELSE IF( nn_ctls_trc >= 2 ) THEN ! control surface = model level
1268	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
1269	ENDIF
1270
1271	! II.1 Define frequency of output and means
1272	! -----------------------------------------
1273	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
1274	ELSE ; clop = "x" ! no use of the mask value (require less cpu time)
1275	ENDIF
1276	# if defined key_diainstant
1277	IF( .NOT. ln_trdmld_trc_instant ) THEN
1278	STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...'
1279	ENDIF
1280	zsto = nn_trd_trc * rdt
1281	clop = "inst("//TRIM(clop)//")"
1282	# else
1283	IF( ln_trdmld_trc_instant ) THEN
1284	zsto = rdt ! inst. diags : we use IOIPSL time averaging
1285	ELSE
1286	zsto = nn_trd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging
1287	ENDIF
1288	clop = "ave("//TRIM(clop)//")"
1289	# endif
1290	zout = nn_trd_trc * rdt
1291	iiter = ( nittrc000 - 1 ) / nn_dttrc
1292
1293	IF(lwp) WRITE (numout,*) ' netCDF initialization'
1294
1295	! II.2 Compute julian date from starting date of the run
1296	! ------------------------------------------------------
1297	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
1298	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
1299	IF(lwp) WRITE(numout,*)' '
1300	IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000 &
1301	& ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday &
1302	& ,'Julian day : ', zjulian
1303
1304	! II.3 Define the T grid trend file (nidtrd)
1305	! ------------------------------------------
1306
1307	!-- Define long and short names for the NetCDF output variables
1308	! ==> choose them according to trdmld_trc_oce.F90 <==
1309
1310	ctrd_trc(jpmld_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmld_trc_xad ,2) = "_xad"
1311	ctrd_trc(jpmld_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmld_trc_yad ,2) = "_yad"
1312	ctrd_trc(jpmld_trc_zad ,1) = " Vertical advection" ; ctrd_trc(jpmld_trc_zad ,2) = "_zad"
1313	ctrd_trc(jpmld_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmld_trc_ldf ,2) = "_ldf"
1314	ctrd_trc(jpmld_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmld_trc_zdf ,2) = "_zdf"
1315	ctrd_trc(jpmld_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmld_trc_bbl ,2) = "_bbl"
1316	ctrd_trc(jpmld_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmld_trc_dmp ,2) = "_dmp"
1317	ctrd_trc(jpmld_trc_sbc ,1) = " Surface boundary cond." ; ctrd_trc(jpmld_trc_sbc ,2) = "_sbc"
1318	ctrd_trc(jpmld_trc_sms, 1) = " Sources minus sinks" ; ctrd_trc(jpmld_trc_sms ,2) = "_sms"
1319	ctrd_trc(jpmld_trc_radb ,1) = " Correct negative concentrations" ; ctrd_trc(jpmld_trc_radb ,2) = "_radb"
1320	ctrd_trc(jpmld_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmld_trc_radn ,2) = "_radn"
1321	ctrd_trc(jpmld_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmld_trc_atf ,2) = "_atf"
1322
1323	DO jn = 1, jptra
1324	!-- Create a NetCDF file and enter the define mode
1325	IF( ln_trdtrc(jn) ) THEN
1326	csuff="ML_"//ctrcnm(jn)
1327	CALL dia_nam( clhstnam, nn_trd_trc, csuff )
1328	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1329	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
1330
1331	!-- Define the ML depth variable
1332	CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m", &
1333	& jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1334
1335	ENDIF
1336	END DO
1337
1338	#if defined key_pisces_reduced
1339	!-- Create a NetCDF file and enter the define mode
1340	CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )
1341	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1342	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom, snc4chunks=snc4set )
1343	#endif
1344
1345	!-- Define physical units
1346	IF( rn_ucf_trc == 1. ) THEN
1347	cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit
1348	ELSEIF ( rn_ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3)
1349	cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf
1350	ELSE
1351	cltrcu = "unknown?"
1352	ENDIF
1353
1354	!-- Define miscellaneous passive tracer mixed-layer variables
1355	IF( jpltrd_trc /= jpmld_trc_atf .OR. jpltrd_trc - 1 /= jpmld_trc_radb ) THEN
1356	STOP 'Error : jpltrd_trc /= jpmld_trc_atf .OR. jpltrd_trc - 1 /= jpmld_trc_radb' ! see below
1357	ENDIF
1358
1359	DO jn = 1, jptra
1360	!
1361	IF( ln_trdtrc(jn) ) THEN
1362	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc.
1363	CALL histdef(nidtrd(jn), trim(clvar), clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", &
1364	& "mmole-N/m3", jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1365	CALL histdef(nidtrd(jn), trim(clvar)//"_tot" , clmxl//" "//trim(ctrcnm(jn))//" Total trend ", &
1366	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1367	CALL histdef(nidtrd(jn), trim(clvar)//"_res" , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)", &
1368	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1369
1370	DO jl = 1, jpltrd_trc - 2 ! <== only true if jpltrd_trc == jpmld_trc_atf
1371	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1), &
1372	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1373	END DO ! if zsto=rdt above
1374
1375	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_radb,1), &
1376	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1377
1378	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmld_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_atf,1), &
1379	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1380	!
1381	ENDIF
1382	END DO
1383
1384	#if defined key_pisces_reduced
1385	DO jl = 1, jp_pisces_trd
1386	CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1)) , &
1387	& cltrcu, jpi, jpj, nh_tb, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1388	END DO ! if zsto=rdt above
1389	#endif
1390
1391	!-- Leave IOIPSL/NetCDF define mode
1392	DO jn = 1, jptra
1393	IF( ln_trdtrc(jn) ) CALL histend( nidtrd(jn), snc4set )
1394	END DO
1395
1396	#if defined key_pisces_reduced
1397	!-- Leave IOIPSL/NetCDF define mode
1398	CALL histend( nidtrdbio, snc4set )
1399
1400	IF(lwp) WRITE(numout,*)
1401	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
1402	#endif
1403
1404	#endif /* key_dimgout */
1405	END SUBROUTINE trd_mld_trc_init
1406
1407	#else
1408	!!----------------------------------------------------------------------
1409	!! Default option : Empty module
1410	!!----------------------------------------------------------------------
1411	CONTAINS
1412	SUBROUTINE trd_mld_trc( kt ) ! Empty routine
1413	INTEGER, INTENT( in) :: kt
1414	WRITE(,) 'trd_mld_trc: You should not have seen this print! error?', kt
1415	END SUBROUTINE trd_mld_trc
1416	SUBROUTINE trd_mld_bio( kt )
1417	INTEGER, INTENT( in) :: kt
1418	WRITE(,) 'trd_mld_bio: You should not have seen this print! error?', kt
1419	END SUBROUTINE trd_mld_bio
1420	SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
1421	INTEGER , INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
1422	CHARACTER(len=2) , INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
1423	REAL, DIMENSION(:,:,:), INTENT( in ) :: ptrc_trdmld ! passive trc trend
1424	WRITE(,) 'trd_mld_trc_zint: You should not have seen this print! error?', ptrc_trdmld(1,1,1)
1425	WRITE(,) ' " " : You should not have seen this print! error?', ctype
1426	WRITE(,) ' " " : You should not have seen this print! error?', ktrd
1427	WRITE(,) ' " " : You should not have seen this print! error?', kjn
1428	END SUBROUTINE trd_mld_trc_zint
1429	SUBROUTINE trd_mld_trc_init ! Empty routine
1430	WRITE(,) 'trd_mld_trc_init: You should not have seen this print! error?'
1431	END SUBROUTINE trd_mld_trc_init
1432	#endif
1433
1434	!!======================================================================
1435	END MODULE trdmld_trc

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

Download in other formats: