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.

trdmxl_trc.F90 in branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 5282

Last change on this file since 5282 was 5282, checked in by diovino, 9 years ago
Dev. branch CMCC4_simplification ticket #1456
File size: 68.2 KB

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

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

Download in other formats: