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

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

trdmld_trc.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 2636

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago
dynamic mem: #785 ; move ctl_stop & warn in lib_mpp to avoid a circular dependency + ctl_stop improvment
Property svn:keywords set to `Id`
File size: 70.1 KB

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

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

Download in other formats: