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

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

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

Context Navigation

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

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago
Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
Property svn:keywords set to `Id`
File size: 69.9 KB

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

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

Download in other formats: