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 @ 5807

Last change on this file since 5807 was 5807, checked in by diovino, 9 years ago

File size: 68.1 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	! IV.1 Code for IOIPSL/NetCDF output
774	! ----------------------------------
775
776	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
777	WRITE(numout,*) ' '
778	WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
779	WRITE(numout,*) '~~~~~~~~~~~ '
780	WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt
781	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
782	WRITE(numout,*) ' '
783	ENDIF
784
785	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
786	!
787
788	DO jn = 1, jptra
789	!
790	IF( ln_trdtrc(jn) ) THEN
791	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
792	!-- Output the fields
793	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
794	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 )
795	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 )
796	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 )
797
798	DO jl = 1, jpltrd_trc - 2
799	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
800	& it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
801	END DO
802
803	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
804	& it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
805
806	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
807	& it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
808
809	ENDIF
810	END DO
811
812	IF( kt == nitend ) THEN
813	DO jn = 1, jptra
814	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
815	END DO
816	ENDIF
817
818	ELSE ! <<< write the trends for passive tracer mean diagnostics
819
820	DO jn = 1, jptra
821	!
822	IF( ln_trdtrc(jn) ) THEN
823	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 )
824	!-- Output the fields
825	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
826
827	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
828	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 )
829	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 )
830
831	DO jl = 1, jpltrd_trc - 2
832	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
833	& it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
834	END DO
835
836	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
837	& it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
838
839	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
840	& it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
841
842	ENDIF
843	!
844	END DO
845	IF( kt == nitend ) THEN
846	DO jn = 1, jptra
847	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
848	END DO
849	ENDIF
850
851	!
852	ENDIF NETCDF_OUTPUT
853
854	! Compute the control surface (for next time step) : flag = on
855	icount = 1
856
857	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
858	!
859	! Reset cumulative arrays to zero
860	! -------------------------------
861	nmoymltrd = 0
862	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
863	tmlradm_trc (:,:,:) = 0.e0 ; tml_sum_trc (:,:,:) = 0.e0
864	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0
865	rmld_sum_trc (:,:) = 0.e0
866	!
867	ENDIF
868
869	! ======================================================================
870	! V. Write restart file
871	! ======================================================================
872
873	IF( lrst_trc ) CALL trd_mxl_trc_rst_write( kt ) ! this must be after the array swap above (III.3)
874
875	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
876	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
877	!
878	END SUBROUTINE trd_mxl_trc
879
880
881	SUBROUTINE trd_mxl_bio( kt )
882	!!----------------------------------------------------------------------
883	!! * ROUTINE trd_mld *
884	!!
885	!! ** Purpose : Compute and cumulate the mixed layer biological trends over an analysis
886	!! period, and write NetCDF outputs.
887	!!
888	!! ** Method/usage :
889	!! The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
890	!! logical namelist variable) :
891	!! 1) to explain the difference between initial and final
892	!! mixed-layer T & S (where initial and final relate to the
893	!! current analysis window, defined by ntrd in the namelist)
894	!! 2) to explain the difference between the current and previous
895	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
896	!! performed over each analysis window).
897	!!
898	!! ** Consistency check :
899	!! If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt
900	!! entrainment) should be zero, at machine accuracy. Note that in the case
901	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
902	!! over the first two analysis windows (except if restart).
903	!! N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8
904	!! for checking residuals.
905	!! On a NEC-SX5 computer, this typically leads to:
906	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
907	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
908	!!
909	!! ** Action :
910	!! At each time step, mixed-layer averaged trends are stored in the
911	!! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
912	!! This array is known when trd_mld is called, at the end of the stp subroutine,
913	!! except for the purely vertical K_z diffusion term, which is embedded in the
914	!! lateral diffusion trend.
915	!!
916	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
917	!! from the lateral diffusion trend.
918	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
919	!! arrays are updated.
920	!! In III), called only once per analysis window, we compute the total trends,
921	!! along with the residuals and the Asselin correction terms.
922	!! In IV), the appropriate trends are written in the trends NetCDF file.
923	!!
924	!! References :
925	!! - Vialard & al.
926	!! - See NEMO documentation (in preparation)
927	!!----------------------------------------------------------------------
928	INTEGER, INTENT( in ) :: kt ! ocean time-step index
929	#if defined key_pisces_reduced
930	INTEGER :: jl, it, itmod
931	LOGICAL :: llwarn = .TRUE., lldebug = .TRUE.
932	REAL(wp) :: zfn, zfn2
933	!!----------------------------------------------------------------------
934	! ... Warnings
935	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
936
937	! ======================================================================
938	! II. Cumulate the trends over the analysis window
939	! ======================================================================
940
941	ztmltrdbio2(:,:,:) = 0.e0 ! <<< reset arrays to zero
942
943	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
944	! ------------------------------------------------------------------------
945	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
946	!
947	tmltrd_csum_ub_bio (:,:,:) = 0.e0
948	!
949	END IF
950
951	! II.4 Cumulated trends over the analysis period
952	! ----------------------------------------------
953	!
954	! [ 1rst analysis window ] [ 2nd analysis window ]
955	!
956	!
957	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
958	! ntrd 2*ntrd etc.
959	! 1 2 3 4 =5 e.g. =10
960	!
961	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN
962	!
963	nmoymltrdbio = nmoymltrdbio + 1
964
965	! ... Trends associated with the time mean of the ML passive tracers
966	tmltrd_sum_bio (:,:,:) = tmltrd_sum_bio (:,:,:) + tmltrd_bio (:,:,:)
967	tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
968	!
969	END IF
970
971	! ======================================================================
972	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
973	! ======================================================================
974
975	! Convert to appropriate physical units
976	tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc
977
978	MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd
979	!
980	zfn = float(nmoymltrdbio) ; zfn2 = zfn * zfn
981
982	! III.1 Prepare fields for output ("instantaneous" diagnostics)
983	! -------------------------------------------------------------
984
985	#if defined key_diainstant
986	STOP 'tmltrd_bio : key_diainstant was never checked within trdmxl. Comment this to proceed.'
987	#endif
988	! III.2 Prepare fields for output ("mean" diagnostics)
989	! ----------------------------------------------------
990
991	ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)
992
993	!-- Lateral boundary conditions
994	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
995	! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut
996	DO jn = 1, jpdiabio
997	CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )
998	ENDDO
999	ENDIF
1000
1001	IF( lldebug ) THEN
1002	!
1003	WRITE(numout,*) 'trd_mxl_bio : write trends in the Mixed Layer for debugging process:'
1004	WRITE(numout,*) '~~~~~~~~~~~ '
1005	WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio
1006	WRITE(numout,*)
1007
1008	DO jl = 1, jpdiabio
1009	IF( ln_trdmxl_trc_instant ) THEN
1010	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1011	& ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))
1012	ELSE
1013	WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX = ', jl, &
1014	& ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))
1015	endif
1016	END DO
1017
1018	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
1019	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
1020	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
1021	WRITE(numout,*)
1022	!
1023	ENDIF
1024
1025	! III.3 Time evolution array swap
1026	! -------------------------------
1027
1028	! For passive tracer mean diagnostics
1029	tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)
1030
1031	! III.4 Convert to appropriate physical units
1032	! -------------------------------------------
1033	ztmltrdbio2 (:,:,:) = ztmltrdbio2 (:,:,:) * rn_ucf_trc/zfn2
1034
1035	END IF MODULO_NTRD
1036
1037	! ======================================================================
1038	! IV. Write trends in the NetCDF file
1039	! ======================================================================
1040
1041	! IV.1 Code for IOIPSL/NetCDF output
1042	! ----------------------------------
1043
1044	! define time axis
1045	itmod = kt - nittrc000 + 1
1046	it = kt
1047
1048	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
1049	WRITE(numout,*) ' '
1050	WRITE(numout,*) 'trd_mxl_bio : write ML bio trends in the NetCDF file :'
1051	WRITE(numout,*) '~~~~~~~~~~~ '
1052	WRITE(numout,*) ' ', TRIM(clhstnam), ' at kt = ', kt
1053	WRITE(numout,*) ' N.B. nmoymltrdbio = ', nmoymltrdbio
1054	WRITE(numout,*) ' '
1055	END IF
1056
1057
1058	! 2. Start writing data
1059	! ---------------------
1060
1061	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
1062	!
1063	DO jl = 1, jpdiabio
1064	CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) , &
1065	& it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
1066	END DO
1067
1068
1069	IF( kt == nitend ) CALL histclo( nidtrdbio )
1070
1071	ELSE ! <<< write the trends for passive tracer mean diagnostics
1072
1073	DO jl = 1, jpdiabio
1074	CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) , &
1075	& it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
1076	END DO
1077
1078	IF( kt == nitend ) CALL histclo( nidtrdbio )
1079	!
1080	END IF NETCDF_OUTPUT
1081
1082	! Compute the control surface (for next time step) : flag = on
1083	icountbio = 1
1084
1085
1086
1087	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
1088	!
1089	! III.5 Reset cumulative arrays to zero
1090	! -------------------------------------
1091	nmoymltrdbio = 0
1092	tmltrd_csum_ln_bio (:,:,:) = 0.e0
1093	tmltrd_sum_bio (:,:,:) = 0.e0
1094	END IF
1095
1096	! ======================================================================
1097	! Write restart file
1098	! ======================================================================
1099
1100	! restart write is done in trd_mxl_trc_write which is called by trd_mxl_bio (Marina)
1101	!
1102	#endif
1103	END SUBROUTINE trd_mxl_bio
1104
1105
1106	REAL FUNCTION sum2d( ztab )
1107	!!----------------------------------------------------------------------
1108	!! CD ??? prevoir d'utiliser plutot prtctl
1109	!!----------------------------------------------------------------------
1110	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: ztab
1111	!!----------------------------------------------------------------------
1112	sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
1113	END FUNCTION sum2d
1114
1115
1116	SUBROUTINE trd_mxl_trc_init
1117	!!----------------------------------------------------------------------
1118	!! * ROUTINE trd_mxl_init *
1119	!!
1120	!! ** Purpose : computation of vertically integrated T and S budgets
1121	!! from ocean surface down to control surface (NetCDF output)
1122	!!
1123	!! ** Method/usage :
1124	!!
1125	!!----------------------------------------------------------------------
1126	INTEGER :: inum ! logical unit
1127	INTEGER :: ilseq, jl, jn, iiter
1128	REAL(wp) :: zjulian, zsto, zout
1129	CHARACTER (LEN=40) :: clop
1130	CHARACTER (LEN=15) :: csuff
1131	CHARACTER (LEN=12) :: clmxl
1132	CHARACTER (LEN=16) :: cltrcu
1133	CHARACTER (LEN=10) :: clvar
1134
1135	!!----------------------------------------------------------------------
1136
1137	! ======================================================================
1138	! I. initialization
1139	! ======================================================================
1140
1141	IF(lwp) THEN
1142	WRITE(numout,*)
1143	WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers '
1144	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
1145	WRITE(numout,*)
1146	ENDIF
1147
1148
1149	! I.1 Check consistency of user defined preferences
1150	! -------------------------------------------------
1151
1152	IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
1153	WRITE(numout,cform_err)
1154	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
1155	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
1156	WRITE(numout,*) ' you defined, nn_trd_trc = ', nn_trd_trc
1157	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
1158	WRITE(numout,*) ' You should reconsider this choice. '
1159	WRITE(numout,*)
1160	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
1161	WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) '
1162	nstop = nstop + 1
1163	ENDIF
1164
1165	! * Debugging information *
1166	IF( lldebug ) THEN
1167	WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl
1168	WRITE(numout,*) ' ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
1169	ENDIF
1170
1171	IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
1172	WRITE(numout,cform_err)
1173	WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL '
1174	WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '
1175	WRITE(numout,*) ' WHY? Everything in trdmxl_trc is coded for leap-frog, and '
1176	WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note '
1177	WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). '
1178	WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '
1179	WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'
1180	WRITE(numout,*)
1181	nstop = nstop + 1
1182	ENDIF
1183
1184	! I.2 Initialize arrays to zero or read a restart file
1185	! ----------------------------------------------------
1186	nmoymltrd = 0
1187
1188	rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst.
1189	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
1190	tmlradm_trc (:,:,:) = 0.e0
1191
1192	tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean
1193	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0
1194
1195	#if defined key_pisces_reduced
1196	nmoymltrdbio = 0
1197	tmltrd_sum_bio (:,:,:) = 0.e0 ; tmltrd_csum_ln_bio (:,:,:) = 0.e0
1198	DO jl = 1, jp_pisces_trd
1199	ctrd_bio(jl,1) = ctrbil(jl) ! long name
1200	ctrd_bio(jl,2) = ctrbio(jl) ! short name
1201	ENDDO
1202	#endif
1203
1204	IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
1205	CALL trd_mxl_trc_rst_read
1206	ELSE
1207	tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst.
1208	tmlbn_trc (:,:,:) = 0.e0
1209
1210	tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean
1211	tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0
1212	#if defined key_pisces_reduced
1213	tmltrd_csum_ub_bio (:,:,:) = 0.e0
1214	#endif
1215
1216	ENDIF
1217
1218	icount = 1 ; ionce = 1 ! open specifier
1219
1220	#if defined key_pisces_reduced
1221	icountbio = 1 ; ioncebio = 1 ! open specifier
1222	#endif
1223
1224	! I.3 Read control surface from file ctlsurf_idx
1225	! ----------------------------------------------
1226	IF( nn_ctls_trc == 1 ) THEN
1227	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
1228	READ ( inum ) nbol_trc
1229	CLOSE( inum )
1230	ENDIF
1231
1232	! ======================================================================
1233	! II. netCDF output initialization
1234	! ======================================================================
1235
1236	! clmxl = legend root for netCDF output
1237	IF( nn_ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion
1238	clmxl = 'Mixed Layer '
1239	ELSE IF( nn_ctls_trc == 1 ) THEN ! control surface = read index from file
1240	clmxl = ' Bowl '
1241	ELSE IF( nn_ctls_trc >= 2 ) THEN ! control surface = model level
1242	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
1243	ENDIF
1244
1245	! II.1 Define frequency of output and means
1246	! -----------------------------------------
1247	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
1248	ELSE ; clop = "x" ! no use of the mask value (require less cp time)
1249	ENDIF
1250	# if defined key_diainstant
1251	IF( .NOT. ln_trdmxl_trc_instant ) THEN
1252	STOP 'trd_mxl_trc : this was never checked. Comment this line to proceed...'
1253	ENDIF
1254	zsto = nn_trd_trc * rdt
1255	clop = "inst("//TRIM(clop)//")"
1256	# else
1257	IF( ln_trdmxl_trc_instant ) THEN
1258	zsto = rdt ! inst. diags : we use IOIPSL time averaging
1259	ELSE
1260	zsto = nn_trd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging
1261	ENDIF
1262	clop = "ave("//TRIM(clop)//")"
1263	# endif
1264	zout = nn_trd_trc * rdt
1265	iiter = ( nittrc000 - 1 ) / nn_dttrc
1266
1267	IF(lwp) WRITE (numout,*) ' netCDF initialization'
1268
1269	! II.2 Compute julian date from starting date of the run
1270	! ------------------------------------------------------
1271	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
1272	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
1273	IF(lwp) WRITE(numout,*)' '
1274	IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000 &
1275	& ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday &
1276	& ,'Julian day : ', zjulian
1277
1278	! II.3 Define the T grid trend file (nidtrd)
1279	! ------------------------------------------
1280
1281	!-- Define long and short names for the NetCDF output variables
1282	! ==> choose them according to trdmxl_trc_oce.F90 <==
1283
1284	ctrd_trc(jpmxl_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmxl_trc_xad ,2) = "_xad"
1285	ctrd_trc(jpmxl_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmxl_trc_yad ,2) = "_yad"
1286	ctrd_trc(jpmxl_trc_zad ,1) = " Vertical advection" ; ctrd_trc(jpmxl_trc_zad ,2) = "_zad"
1287	ctrd_trc(jpmxl_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmxl_trc_ldf ,2) = "_ldf"
1288	ctrd_trc(jpmxl_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmxl_trc_zdf ,2) = "_zdf"
1289	ctrd_trc(jpmxl_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmxl_trc_bbl ,2) = "_bbl"
1290	ctrd_trc(jpmxl_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmxl_trc_dmp ,2) = "_dmp"
1291	ctrd_trc(jpmxl_trc_sbc ,1) = " Surface boundary cond." ; ctrd_trc(jpmxl_trc_sbc ,2) = "_sbc"
1292	ctrd_trc(jpmxl_trc_sms, 1) = " Sources minus sinks" ; ctrd_trc(jpmxl_trc_sms ,2) = "_sms"
1293	ctrd_trc(jpmxl_trc_radb ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radb ,2) = "_radb"
1294	ctrd_trc(jpmxl_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radn ,2) = "_radn"
1295	ctrd_trc(jpmxl_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmxl_trc_atf ,2) = "_atf"
1296
1297	DO jn = 1, jptra
1298	!-- Create a NetCDF file and enter the define mode
1299	IF( ln_trdtrc(jn) ) THEN
1300	csuff="ML_"//ctrcnm(jn)
1301	CALL dia_nam( clhstnam, nn_trd_trc, csuff )
1302	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1303	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
1304
1305	!-- Define the ML depth variable
1306	CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m", &
1307	& jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1308
1309	ENDIF
1310	END DO
1311
1312	#if defined key_pisces_reduced
1313	!-- Create a NetCDF file and enter the define mode
1314	CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' )
1315	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
1316	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom, snc4chunks=snc4set )
1317	#endif
1318
1319	!-- Define physical units
1320	IF( rn_ucf_trc == 1. ) THEN
1321	cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit
1322	ELSEIF ( rn_ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3)
1323	cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf
1324	ELSE
1325	cltrcu = "unknown?"
1326	ENDIF
1327
1328	!-- Define miscellaneous passive tracer mixed-layer variables
1329	IF( jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
1330	STOP 'Error : jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb' ! see below
1331	ENDIF
1332
1333	DO jn = 1, jptra
1334	!
1335	IF( ln_trdtrc(jn) ) THEN
1336	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc.
1337	CALL histdef(nidtrd(jn), trim(clvar), clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", &
1338	& "mmole-N/m3", jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
1339	CALL histdef(nidtrd(jn), trim(clvar)//"_tot" , clmxl//" "//trim(ctrcnm(jn))//" Total trend ", &
1340	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1341	CALL histdef(nidtrd(jn), trim(clvar)//"_res" , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)", &
1342	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
1343
1344	DO jl = 1, jpltrd_trc - 2 ! <== only true if jpltrd_trc == jpmxl_trc_atf
1345	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1), &
1346	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1347	END DO ! if zsto=rdt above
1348
1349	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), &
1350	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1351
1352	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1), &
1353	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1354	!
1355	ENDIF
1356	END DO
1357
1358	#if defined key_pisces_reduced
1359	DO jl = 1, jp_pisces_trd
1360	CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1)) , &
1361	& cltrcu, jpi, jpj, nh_tb, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1362	END DO ! if zsto=rdt above
1363	#endif
1364
1365	!-- Leave IOIPSL/NetCDF define mode
1366	DO jn = 1, jptra
1367	IF( ln_trdtrc(jn) ) CALL histend( nidtrd(jn), snc4set )
1368	END DO
1369
1370	#if defined key_pisces_reduced
1371	!-- Leave IOIPSL/NetCDF define mode
1372	CALL histend( nidtrdbio, snc4set )
1373
1374	IF(lwp) WRITE(numout,*)
1375	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
1376	#endif
1377
1378	END SUBROUTINE trd_mxl_trc_init
1379
1380	#else
1381	!!----------------------------------------------------------------------
1382	!! Default option : Empty module
1383	!!----------------------------------------------------------------------
1384	CONTAINS
1385	SUBROUTINE trd_mxl_trc( kt ) ! Empty routine
1386	INTEGER, INTENT( in) :: kt
1387	WRITE(,) 'trd_mxl_trc: You should not have seen this print! error?', kt
1388	END SUBROUTINE trd_mxl_trc
1389	SUBROUTINE trd_mxl_bio( kt )
1390	INTEGER, INTENT( in) :: kt
1391	WRITE(,) 'trd_mxl_bio: You should not have seen this print! error?', kt
1392	END SUBROUTINE trd_mxl_bio
1393	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
1394	INTEGER , INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
1395	CHARACTER(len=2) , INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
1396	REAL, DIMENSION(:,:,:), INTENT( in ) :: ptrc_trdmxl ! passive trc trend
1397	WRITE(,) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
1398	WRITE(,) ' " " : You should not have seen this print! error?', ctype
1399	WRITE(,) ' " " : You should not have seen this print! error?', ktrd
1400	WRITE(,) ' " " : You should not have seen this print! error?', kjn
1401	END SUBROUTINE trd_mxl_trc_zint
1402	SUBROUTINE trd_mxl_trc_init ! Empty routine
1403	WRITE(,) 'trd_mxl_trc_init: You should not have seen this print! error?'
1404	END SUBROUTINE trd_mxl_trc_init
1405	#endif
1406
1407	!!======================================================================
1408	END MODULE trdmxl_trc

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

Download in other formats: