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/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 2830

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

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

Download in other formats: