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

Context Navigation

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 2082

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

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

Download in other formats: