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

Context Navigation

source: tags/nemo_v3_2_beta/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 7039

Last change on this file since 7039 was 1581, checked in by smasson, 15 years ago
ctlopn cleanup, see ticket:515 and ticket:237
File size: 77.7 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, itmod, 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.( ln_rsttr ) ) 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	itmod = kt - nittrc000 + 1
578	it = kt
579
580	MODULO_NTRD : IF( MOD( itmod, ntrd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd_trc
581	!
582	ztmltot (:,:,:) = 0.e0 ! reset arrays to zero
583	ztmlres (:,:,:) = 0.e0
584	ztmltot2(:,:,:) = 0.e0
585	ztmlres2(:,:,:) = 0.e0
586
587	zfn = FLOAT( nmoymltrd ) ; zfn2 = zfn * zfn
588
589	! III.1 Prepare fields for output ("instantaneous" diagnostics)
590	! -------------------------------------------------------------
591
592	DO jn = 1, jptra
593	IF( luttrd(jn) ) THEN
594	!-- Compute total trends (use rdttrc instead of rdt ???)
595	IF ( ln_trcadv_smolar .OR. ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes
596	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
597	ELSE ! LEAP-FROG schemes
598	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
599	ENDIF
600
601	!-- Compute residuals
602	ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
603	& - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) )
604
605	!-- Diagnose Asselin trend over the analysis window
606	ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
607	ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
608
609	!-- Lateral boundary conditions
610	#if ! defined key_gyre
611
612	CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
613	CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
614
615	#endif
616
617	#if defined key_diainstant
618	STOP 'tmltrd_trc : key_diainstant was never checked within trdmld. Comment this to proceed.'
619	#endif
620	ENDIF
621	END DO
622
623	! III.2 Prepare fields for output ("mean" diagnostics)
624	! ----------------------------------------------------
625
626	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
627	rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
628
629	!-- Compute passive tracer total trends
630	DO jn = 1, jptra
631	IF( luttrd(jn) ) THEN
632	tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
633	ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rdt ) ! now tracer unit is /sec
634	ENDIF
635	END DO
636
637	!-- Compute passive tracer residuals
638	DO jn = 1, jptra
639	IF( luttrd(jn) ) THEN
640	!
641	DO jl = 1, jpltrd_trc
642	ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
643	END DO
644
645	ztmltrdm2(:,:,jn) = 0.e0
646	DO jl = 1, jpltrd_trc
647	ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
648	END DO
649
650	ztmlres2(:,:,jn) = ztmltot2(:,:,jn) - &
651	& ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
652	& - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
653	!
654
655	!-- Diagnose Asselin trend over the analysis window
656	ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) &
657	& + tmltrd_atf_sumb_trc(:,:,jn)
658	ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) &
659	& + tmltrd_rad_sumb_trc(:,:,jn)
660
661	!-- Lateral boundary conditions
662	#if ! defined key_gyre
663	CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
664	CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
665	DO jl = 1, jpltrd_trc
666	CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file
667	END DO
668	#endif
669	ENDIF
670	END DO
671
672	! * Debugging information *
673	IF( lldebug ) THEN
674	!
675	WRITE(numout,*) 'trd_mld_trc : write trends in the Mixed Layer for debugging process:'
676	WRITE(numout,*) '~~~~~~~~~~~ '
677	WRITE(numout,*)
678	WRITE(numout,*) 'TRC kt = ', kt, ' nmoymltrd = ', nmoymltrd
679
680	DO jn = 1, jptra
681
682	IF( luttrd(jn) ) THEN
683	WRITE(numout, *)
684	WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
685
686	WRITE(numout, *)
687	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres : ', SUM2D(ztmlres(:,:,jn))
688	!CD??? PREVOIR: z2d = ztmlres(:,:,jn) ; CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres - : ')
689
690	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
691	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
692	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot : ', SUM2D(+ztmltot (:,:,jn))
693	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
694	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
695	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
696	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
697	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
698	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
699
700	WRITE(numout, *)
701	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2 : ', SUM2D(ztmlres2(:,:,jn))
702	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
703	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
704	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2 : ', SUM2D(+ztmltot2(:,:,jn))
705	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2 : ', SUM2D(+ztmltrdm2(:,:,jn))
706	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_atf,jn))
707	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
708	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_radb,jn))
709	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
710	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
711
712	WRITE(numout, *)
713	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot : ', SUM2D(ztmltot (:,:,jn))
714	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
715	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc : ', SUM2D(tml_trc (:,:,jn))
716	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc : ', SUM2D(tmlbn_trc (:,:,jn))
717	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc : ', SUM2D(tmlb_trc (:,:,jn))
718	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc : ', SUM2D(tmlbb_trc (:,:,jn))
719	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
720
721	WRITE(numout, *)
722	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
723	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
724	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
725	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
726	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
727
728	WRITE(numout, *)
729	DO jl = 1, jpltrd_trc
730	WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmld_trc_xxx = ', jl, &
731	& ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
732	END DO
733
734	WRITE(numout,*)
735	WRITE(numout,) ' ******************** ZTMLRES, ZTMLRES2 ********************* '
736	WRITE(numout,*)
737	WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
738	WRITE(numout,*)
739	WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
740
741	WRITE(numout,*)
742	WRITE(numout,) ' ******************** ZTMLRES ********************* '
743	WRITE(numout,*)
744
745	WRITE(numout,*) '...................................................'
746	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
747	DO jj = 5, 1, -1
748	WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
749	END DO
750
751	WRITE(numout,*)
752	WRITE(numout,) ' ******************** ZTMLRES2 ********************* '
753	WRITE(numout,*)
754
755	WRITE(numout,*) '...................................................'
756	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
757	DO jj = 5, 1, -1
758	WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
759	END DO
760	!
761	ENDIF
762	!
763	END DO
764
765
766	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
767	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
768	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
769	WRITE(numout,*)
770	!
771	ENDIF
772
773	! III.3 Time evolution array swap
774	! -------------------------------
775	! ML depth
776	rmldbn_trc(:,:) = rmld_trc(:,:)
777	rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
778	DO jn = 1, jptra
779	IF( luttrd(jn) ) THEN
780	! For passive tracer instantaneous diagnostics
781	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
782	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
783
784	! For passive tracer mean diagnostics
785	tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
786	tml_sumb_trc (:,:,jn) = tml_sum_trc(:,:,jn)
787	tmltrd_atf_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn)
788	tmltrd_rad_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmld_trc_radb,jn)
789
790
791	! III.4 Convert to appropriate physical units
792	! -------------------------------------------
793	ztmltot (:,:,jn) = ztmltot (:,:,jn) * ucf_trc/zfn ! instant diags
794	ztmlres (:,:,jn) = ztmlres (:,:,jn) * ucf_trc/zfn
795	ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * ucf_trc/zfn
796	ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * ucf_trc/zfn
797	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags
798	ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * ucf_trc/zfn2
799	ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * ucf_trc/zfn2
800	ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * ucf_trc/zfn2
801	ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * ucf_trc/zfn2
802	ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * ucf_trc/zfn2
803	ENDIF
804	END DO
805	!
806	ENDIF MODULO_NTRD
807
808	! ======================================================================
809	! IV. Write trends in the NetCDF file
810	! ======================================================================
811
812	! IV.1 Code for dimg mpp output
813	! -----------------------------
814
815	# if defined key_dimgout
816	STOP 'Not implemented'
817	# else
818
819	! IV.2 Code for IOIPSL/NetCDF output
820	! ----------------------------------
821
822	IF( lwp .AND. MOD( itmod , ntrd_trc ) == 0 ) THEN
823	WRITE(numout,*) ' '
824	WRITE(numout,*) 'trd_mld_trc : write passive tracer trends in the NetCDF file :'
825	WRITE(numout,*) '~~~~~~~~~~~ '
826	WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt
827	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
828	WRITE(numout,*) ' '
829	ENDIF
830
831	NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
832	!
833
834	DO jn = 1, jptra
835	!
836	IF( luttrd(jn) ) THEN
837	!-- Specific treatment for EIV trends
838	! WARNING : When eiv is switched on but key_diaeiv is not, we do NOT diagnose
839	! u_eiv, v_eiv, and w_eiv : the exact eiv advective trends thus cannot be computed,
840	! only their sum makes sense => mask directional contrib. to avoid confusion
841	z2d(:,:) = tmltrd_trc(:,:,jpmld_trc_xei,jn) + tmltrd_trc(:,:,jpmld_trc_yei,jn) &
842	& + tmltrd_trc(:,:,jpmld_trc_zei,jn)
843	#if ( defined key_trcldf_eiv && defined key_diaeiv )
844	tmltrd_trc(:,:,jpmld_trc_xei,jn) = -999.
845	tmltrd_trc(:,:,jpmld_trc_yei,jn) = -999.
846	tmltrd_trc(:,:,jpmld_trc_zei,jn) = -999.
847	#endif
848	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
849	!-- Output the fields
850	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
851	CALL histwrite( nidtrd(jn), clvar , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 )
852	CALL histwrite( nidtrd(jn), clvar//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 )
853	CALL histwrite( nidtrd(jn), clvar//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 )
854
855	DO jl = 1, jpltrd_trc - 2
856	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), &
857	& it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
858	END DO
859
860	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
861	& it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
862
863	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), & ! now Asselin : jpltrd_trc
864	& it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
865
866	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)), & ! now total EIV : jpltrd_trc + 1
867	& it, z2d(:,:), ndimtrd1, ndextrd1 )
868	!
869	ENDIF
870	END DO
871
872	IF( kt == nitend ) THEN
873	DO jn = 1, jptra
874	IF( luttrd(jn) ) CALL histclo( nidtrd(jn) )
875	END DO
876	ENDIF
877
878	ELSE ! <<< write the trends for passive tracer mean diagnostics
879
880
881	DO jn = 1, jptra
882	!
883	IF( luttrd(jn) ) THEN
884	!-- Specific treatment for EIV trends
885	! WARNING : see above
886	z2d(:,:) = ztmltrd2(:,:,jpmld_trc_xei,jn) + ztmltrd2(:,:,jpmld_trc_yei,jn) &
887	& + ztmltrd2(:,:,jpmld_trc_zei,jn)
888
889	#if ( defined key_trcldf_eiv && defined key_diaeiv )
890	ztmltrd2(:,:,jpmld_trc_xei,jn) = -999.
891	ztmltrd2(:,:,jpmld_trc_yei,jn) = -999.
892	ztmltrd2(:,:,jpmld_trc_zei,jn) = -999.
893	#endif
894	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 )
895	!-- Output the fields
896	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
897
898	CALL histwrite( nidtrd(jn), clvar , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
899	CALL histwrite( nidtrd(jn), clvar//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 )
900	CALL histwrite( nidtrd(jn), clvar//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 )
901
902	DO jl = 1, jpltrd_trc - 2
903	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), &
904	& it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
905	END DO
906
907	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
908	& it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
909
910	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), & ! now Asselin : jpltrd_trc
911	& it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
912
913	CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)), & ! now total EIV : jpltrd_trc + 1
914	& it, z2d(:,:), ndimtrd1, ndextrd1 )
915
916	ENDIF
917	!
918	END DO
919	IF( kt == nitend ) THEN
920	DO jn = 1, jptra
921	IF( luttrd(jn) ) CALL histclo( nidtrd(jn) )
922	END DO
923	ENDIF
924
925	!
926	ENDIF NETCDF_OUTPUT
927
928	! Compute the control surface (for next time step) : flag = on
929	icount = 1
930
931	# endif /* key_dimgout */
932
933	IF( MOD( itmod, ntrd_trc ) == 0 ) THEN
934	!
935	! Reset cumulative arrays to zero
936	! -------------------------------
937	nmoymltrd = 0
938	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
939	tmlradm_trc (:,:,:) = 0.e0 ; tml_sum_trc (:,:,:) = 0.e0
940	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0
941	rmld_sum_trc (:,:) = 0.e0
942	!
943	ENDIF
944
945	! ======================================================================
946	! V. Write restart file
947	! ======================================================================
948
949	IF( lrst_trc ) CALL trd_mld_trc_rst_write( kt ) ! this must be after the array swap above (III.3)
950
951	END SUBROUTINE trd_mld_trc
952
953	SUBROUTINE trd_mld_bio( kt )
954	!!----------------------------------------------------------------------
955	!! * ROUTINE trd_mld *
956	!!
957	!! ** Purpose : Compute and cumulate the mixed layer biological trends over an analysis
958	!! period, and write NetCDF (or dimg) outputs.
959	!!
960	!! ** Method/usage :
961	!! The stored trends can be chosen twofold (according to the ln_trdmld_trc_instant
962	!! logical namelist variable) :
963	!! 1) to explain the difference between initial and final
964	!! mixed-layer T & S (where initial and final relate to the
965	!! current analysis window, defined by ntrd in the namelist)
966	!! 2) to explain the difference between the current and previous
967	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
968	!! performed over each analysis window).
969	!!
970	!! ** Consistency check :
971	!! If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt
972	!! entrainment) should be zero, at machine accuracy. Note that in the case
973	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
974	!! over the first two analysis windows (except if restart).
975	!! N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8
976	!! for checking residuals.
977	!! On a NEC-SX5 computer, this typically leads to:
978	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmld_trc_instant=.false.
979	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmld_trc_instant=.true.
980	!!
981	!! ** Action :
982	!! At each time step, mixed-layer averaged trends are stored in the
983	!! tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
984	!! This array is known when trd_mld is called, at the end of the stp subroutine,
985	!! except for the purely vertical K_z diffusion term, which is embedded in the
986	!! lateral diffusion trend.
987	!!
988	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
989	!! from the lateral diffusion trend.
990	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
991	!! arrays are updated.
992	!! In III), called only once per analysis window, we compute the total trends,
993	!! along with the residuals and the Asselin correction terms.
994	!! In IV), the appropriate trends are written in the trends NetCDF file.
995	!!
996	!! References :
997	!! - Vialard & al.
998	!! - See NEMO documentation (in preparation)
999	!!----------------------------------------------------------------------
1000	INTEGER, INTENT( in ) :: kt ! ocean time-step index
1001	#if defined key_lobster
1002	INTEGER :: jl, it, itmod
1003	LOGICAL :: llwarn = .TRUE., lldebug = .TRUE.
1004	REAL(wp), DIMENSION(jpi,jpj,jpdiabio) :: ztmltrdbio2 ! only needed for mean diagnostics
1005	REAL(wp) :: zfn, zfn2
1006	#if defined key_dimgout
1007	INTEGER :: iyear,imon,iday
1008	CHARACTER(LEN=80) :: cltext, clmode
1009	#endif
1010	!!----------------------------------------------------------------------
1011	! ... Warnings
1012	IF( llwarn ) THEN
1013	IF( ( nittrc000 /= nit000 ) &
1014	.OR.( ndttrc /= 1 ) ) THEN
1015
1016	WRITE(numout,*) 'Be careful, trends diags never validated'
1017	STOP 'Uncomment this line to proceed'
1018	END IF
1019	END IF
1020
1021	! ======================================================================
1022	! II. Cumulate the trends over the analysis window
1023	! ======================================================================
1024
1025	ztmltrdbio2(:,:,:) = 0.e0 ! <<< reset arrays to zero
1026
1027	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
1028	! ------------------------------------------------------------------------
1029	IF( kt == 2 ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
1030	!
1031	tmltrd_csum_ub_bio (:,:,:) = 0.e0
1032	!
1033	END IF
1034
1035	! II.4 Cumulated trends over the analysis period
1036	! ----------------------------------------------
1037	!
1038	! [ 1rst analysis window ] [ 2nd analysis window ]
1039	!
1040	!
1041	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
1042	! ntrd 2*ntrd etc.
1043	! 1 2 3 4 =5 e.g. =10
1044	!
1045	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN
1046	!
1047	nmoymltrdbio = nmoymltrdbio + 1
1048
1049	! ... Trends associated with the time mean of the ML passive tracers
1050	tmltrd_sum_bio (:,:,:) = tmltrd_sum_bio (:,:,:) + tmltrd_bio (:,:,:)
1051	tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
1052	!
1053	END IF
1054
1055	! ======================================================================
1056	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
1057	! ======================================================================
1058
1059	! Convert to appropriate physical units
1060	tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * ucf_trc
1061
1062	MODULO_NTRD : IF( MOD( kt, ntrd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd
1063	!
1064	zfn = float(nmoymltrdbio) ; zfn2 = zfn * zfn
1065
1066	! III.1 Prepare fields for output ("instantaneous" diagnostics)
1067	! -------------------------------------------------------------
1068
1069	#if defined key_diainstant
1070	STOP 'tmltrd_bio : key_diainstant was never checked within trdmld. Comment this to proceed.'
1071	#endif
1072	! III.2 Prepare fields for output ("mean" diagnostics)
1073	! ----------------------------------------------------
1074
1075	ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)
1076
1077	!-- Lateral boundary conditions
1078	#if ! defined key_gyre
1079	! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut
1080	DO jn = 1, jpdiabio
1081	CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )
1082	ENDDO
1083	#endif
1084	IF( lldebug ) THEN
1085	!
1086	WRITE(numout,*) 'trd_mld_bio : write trends in the Mixed Layer for debugging process:'
1087	WRITE(numout,*) '~~~~~~~~~~~ '
1088	WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio
1089	WRITE(numout,*)
1090
1091	DO jl = 1, jpdiabio
1092	IF( ln_trdmld_trc_instant ) THEN
1093	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1094	& ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))
1095	ELSE
1096	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1097	& ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))
1098	endif
1099	END DO
1100
1101	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
1102	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
1103	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
1104	WRITE(numout,*)
1105	!
1106	ENDIF
1107
1108	! III.3 Time evolution array swap
1109	! -------------------------------
1110
1111	! For passive tracer mean diagnostics
1112	tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)
1113
1114	! III.4 Convert to appropriate physical units
1115	! -------------------------------------------
1116	ztmltrdbio2 (:,:,:) = ztmltrdbio2 (:,:,:) * ucf_trc/zfn2
1117
1118	END IF MODULO_NTRD
1119
1120	! ======================================================================
1121	! IV. Write trends in the NetCDF file
1122	! ======================================================================
1123
1124	! IV.1 Code for dimg mpp output
1125	! -----------------------------
1126
1127	# if defined key_dimgout
1128	STOP 'Not implemented'
1129	# else
1130
1131	! IV.2 Code for IOIPSL/NetCDF output
1132	! ----------------------------------
1133
1134	! define time axis
1135	itmod = kt - nittrc000 + 1
1136	it = kt
1137
1138	IF( lwp .AND. MOD( itmod , ntrd_trc ) == 0 ) THEN
1139	WRITE(numout,*) ' '
1140	WRITE(numout,*) 'trd_mld_bio : write ML bio trends in the NetCDF file :'
1141	WRITE(numout,*) '~~~~~~~~~~~ '
1142	WRITE(numout,*) ' ', TRIM(clhstnam), ' at kt = ', kt
1143	WRITE(numout,*) ' N.B. nmoymltrdbio = ', nmoymltrdbio
1144	WRITE(numout,*) ' '
1145	END IF
1146
1147
1148	! 2. Start writing data
1149	! ---------------------
1150
1151	NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
1152	!
1153	DO jl = 1, jpdiabio
1154	CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) , &
1155	& it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
1156	END DO
1157
1158
1159	IF( kt == nitend ) CALL histclo( nidtrdbio )
1160
1161	ELSE ! <<< write the trends for passive tracer mean diagnostics
1162
1163	DO jl = 1, jpdiabio
1164	CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) , &
1165	& it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
1166	END DO
1167
1168	IF( kt == nitend ) CALL histclo( nidtrdbio )
1169	!
1170	END IF NETCDF_OUTPUT
1171
1172	! Compute the control surface (for next time step) : flag = on
1173	icountbio = 1
1174
1175
1176	# endif /* key_dimgout */
1177
1178	IF( MOD( itmod, ntrd_trc ) == 0 ) THEN
1179	!
1180	! III.5 Reset cumulative arrays to zero
1181	! -------------------------------------
1182	nmoymltrdbio = 0
1183	tmltrd_csum_ln_bio (:,:,:) = 0.e0
1184	tmltrd_sum_bio (:,:,:) = 0.e0
1185	END IF
1186
1187	! ======================================================================
1188	! Write restart file
1189	! ======================================================================
1190
1191	! restart write is done in trd_mld_trc_write which is called by trd_mld_bio (Marina)
1192	!
1193	#endif
1194	END SUBROUTINE trd_mld_bio
1195
1196	REAL FUNCTION sum2d( ztab )
1197	!!----------------------------------------------------------------------
1198	!! CD ??? prevoir d'utiliser plutot prtctl
1199	!!----------------------------------------------------------------------
1200	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: ztab
1201	!!----------------------------------------------------------------------
1202	sum2d = SUM(ztab(2:jpi-1,2:jpj-1))
1203	END FUNCTION sum2d
1204
1205	SUBROUTINE trd_mld_trc_init
1206	!!----------------------------------------------------------------------
1207	!! * ROUTINE trd_mld_init *
1208	!!
1209	!! ** Purpose : computation of vertically integrated T and S budgets
1210	!! from ocean surface down to control surface (NetCDF output)
1211	!!
1212	!! ** Method/usage :
1213	!!
1214	!!----------------------------------------------------------------------
1215	INTEGER :: ilseq, jl, jn
1216	REAL(wp) :: zjulian, zsto, zout
1217	CHARACTER (LEN=40) :: clop, cleiv
1218	CHARACTER (LEN=15) :: csuff
1219	CHARACTER (LEN=12) :: clmxl
1220	CHARACTER (LEN=16) :: cltrcu
1221	CHARACTER (LEN= 5) :: clvar
1222
1223	NAMELIST/namtoptrd/ ntrd_trc, nctls_trc, ucf_trc, &
1224	ln_trdmld_trc_restart, ln_trdmld_trc_instant, luttrd
1225
1226	!!----------------------------------------------------------------------
1227
1228	! ======================================================================
1229	! I. initialization
1230	! ======================================================================
1231
1232	IF(lwp) THEN
1233	WRITE(numout,*)
1234	WRITE(numout,*) ' trd_mld_trc_init : Mixed-layer trends for passive tracers '
1235	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
1236	WRITE(numout,*)
1237	ENDIF
1238
1239
1240	! I.1 Check consistency of user defined preferences
1241	! -------------------------------------------------
1242	#if defined key_trcldf_eiv
1243	IF( lk_trdmld_trc .AND. ln_trcldf_iso ) THEN
1244	WRITE(numout,cform_war)
1245	WRITE(numout,*) ' You asked for ML diagnostics with iso-neutral diffusion '
1246	WRITE(numout,*) ' and eiv physics. '
1247	WRITE(numout,*) ' Yet, key_diaeiv is NOT switched on, so the eddy induced '
1248	WRITE(numout,*) ' velocity is not diagnosed. '
1249	WRITE(numout,*) ' Therefore, we cannot deduce the eiv advective trends. '
1250	WRITE(numout,*) ' Only THE SUM of the i,j,k directional contributions then '
1251	WRITE(numout,*) ' makes sense => To avoid any confusion, we choosed to mask '
1252	WRITE(numout,*) ' these i,j,k directional contributions (with -999.) '
1253	nwarn = nwarn + 1
1254	ENDIF
1255	# endif
1256
1257	IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, ntrd_trc ) /= 0 ) ) THEN
1258	WRITE(numout,cform_err)
1259	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
1260	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
1261	WRITE(numout,*) ' you defined, ntrd_trc = ', ntrd_trc
1262	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
1263	WRITE(numout,*) ' You should reconsider this choice. '
1264	WRITE(numout,*)
1265	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
1266	WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) '
1267	nstop = nstop + 1
1268	ENDIF
1269
1270	IF( ( lk_trdmld_trc ) .AND. ( n_cla == 1 ) ) THEN
1271	WRITE(numout,cform_war)
1272	WRITE(numout,*) ' You set n_cla = 1. Note that the Mixed-Layer diagnostics '
1273	WRITE(numout,*) ' are not exact along the corresponding straits. '
1274	nwarn = nwarn + 1
1275	ENDIF
1276
1277
1278	! * Debugging information *
1279	IF( lldebug ) THEN
1280	WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl
1281	WRITE(numout,*) ' ln_trcadv_smolar = ' , ln_trcadv_smolar
1282	WRITE(numout,*) ' ln_trdmld_trc_instant = ', ln_trdmld_trc_instant
1283	ENDIF
1284
1285	IF( ln_trcadv_smolar .AND. .NOT. ln_trdmld_trc_instant ) THEN
1286	WRITE(numout,cform_err)
1287	WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer Smolark. '
1288	WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '
1289	WRITE(numout,*) ' WHY? Everything in trdmld_trc is coded for leap-frog, and '
1290	WRITE(numout,*) ' Smolarkiewicz scheme is Euler forward. '
1291	WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '
1292	WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'
1293	WRITE(numout,*)
1294	nstop = nstop + 1
1295	ENDIF
1296
1297	IF( ln_trcadv_muscl .AND. .NOT. ln_trdmld_trc_instant ) THEN
1298	WRITE(numout,cform_err)
1299	WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL '
1300	WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '
1301	WRITE(numout,*) ' WHY? Everything in trdmld_trc is coded for leap-frog, and '
1302	WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note '
1303	WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). '
1304	WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '
1305	WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'
1306	WRITE(numout,*)
1307	nstop = nstop + 1
1308	ENDIF
1309
1310	IF( ln_trcadv_muscl2 .AND. .NOT. ln_trdmld_trc_instant ) THEN
1311	WRITE(numout,cform_err)
1312	WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL2 '
1313	WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '
1314	WRITE(numout,*) ' WHY? Everything in trdmld_trc is coded for leap-frog, and '
1315	WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note '
1316	WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). '
1317	WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '
1318	WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'
1319	WRITE(numout,*)
1320	nstop = nstop + 1
1321	ENDIF
1322
1323	! I.2 Initialize arrays to zero or read a restart file
1324	! ----------------------------------------------------
1325	nmoymltrd = 0
1326
1327	rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst.
1328	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
1329	tmlradm_trc (:,:,:) = 0.e0
1330
1331	tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean
1332	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0
1333
1334	#if defined key_lobster
1335	nmoymltrdbio = 0
1336	tmltrd_sum_bio (:,:,:) = 0.e0 ; tmltrd_csum_ln_bio (:,:,:) = 0.e0
1337	DO jl = 1, jp_lobster_trd
1338	ctrd_bio(jl,1) = ctrbil(jl) ! long name
1339	ctrd_bio(jl,2) = ctrbio(jl) ! short name
1340	ENDDO
1341	#endif
1342
1343	IF( ln_rsttr .AND. ln_trdmld_trc_restart ) THEN
1344	CALL trd_mld_trc_rst_read
1345	ELSE
1346	tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst.
1347	tmlbn_trc (:,:,:) = 0.e0
1348
1349	tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean
1350	tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0
1351	#if defined key_lobster
1352	tmltrd_csum_ub_bio (:,:,:) = 0.e0
1353	#endif
1354
1355	ENDIF
1356
1357	icount = 1 ; ionce = 1 ! open specifier
1358
1359	#if defined key_lobster
1360	icountbio = 1 ; ioncebio = 1 ! open specifier
1361	#endif
1362
1363	! I.3 Read control surface from file ctlsurf_idx
1364	! ----------------------------------------------
1365	IF( nctls_trc == 1 ) THEN
1366	CALL ctl_opn( numbol, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
1367	READ ( numbol ) nbol_trc
1368	CLOSE( numbol )
1369	ENDIF
1370
1371	! ======================================================================
1372	! II. netCDF output initialization
1373	! ======================================================================
1374
1375	#if defined key_dimgout
1376	???
1377	#else
1378	! clmxl = legend root for netCDF output
1379	IF( nctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion
1380	clmxl = 'Mixed Layer '
1381	ELSE IF( nctls_trc == 1 ) THEN ! control surface = read index from file
1382	clmxl = ' Bowl '
1383	ELSE IF( nctls_trc >= 2 ) THEN ! control surface = model level
1384	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls_trc
1385	ENDIF
1386
1387	! II.1 Define frequency of output and means
1388	! -----------------------------------------
1389	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
1390	ELSE ; clop = "x" ! no use of the mask value (require less cpu time)
1391	ENDIF
1392	# if defined key_diainstant
1393	IF( .NOT. ln_trdmld_trc_instant ) THEN
1394	STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...'
1395	ENDIF
1396	zsto = ntrd_trc * rdt
1397	clop = "inst("//TRIM(clop)//")"
1398	# else
1399	IF( ln_trdmld_trc_instant ) THEN
1400	zsto = rdt ! inst. diags : we use IOIPSL time averaging
1401	ELSE
1402	zsto = ntrd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging
1403	ENDIF
1404	clop = "ave("//TRIM(clop)//")"
1405	# endif
1406	zout = ntrd_trc * rdt
1407
1408	IF(lwp) WRITE (numout,*) ' netCDF initialization'
1409
1410	! II.2 Compute julian date from starting date of the run
1411	! ------------------------------------------------------
1412	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
1413	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
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, nittrc000-ndttrc, 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, nittrc000-ndttrc, 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: