New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

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

Context Navigation

source: trunk/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 1283

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

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

Download in other formats: