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/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 7068

Last change on this file since 7068 was 7068, checked in by cetlod, 7 years ago
ROBUST5_CNRS : implementation of part I of new TOP interface - 1st step -, see ticket #1782
Property svn:keywords set to `Id`
File size: 51.4 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
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 trdtrc_oce ! definition of main arrays used for trends computations
27	USE in_out_manager ! I/O manager
28	USE dianam ! build the name of file (routine)
29	USE ldfslp ! iso-neutral slopes
30	USE ioipsl ! NetCDF library
31	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
32	USE lib_mpp ! MPP library
33	USE trdmxl_trc_rst ! restart for diagnosing the ML trends
34	USE prtctl ! print control
35	USE sms_pisces ! PISCES bio-model
36	USE wrk_nemo ! Memory allocation
37
38	IMPLICIT NONE
39	PRIVATE
40
41	PUBLIC trd_mxl_trc
42	PUBLIC trd_mxl_trc_alloc
43	PUBLIC trd_mxl_trc_init
44	PUBLIC trd_mxl_trc_zint
45
46	CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file
47	INTEGER :: nmoymltrd
48	INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndextrd1, nidtrd, nh_t
49	INTEGER :: ndimtrd1
50	INTEGER, SAVE :: ionce, icount
51	LOGICAL :: llwarn = .TRUE. ! this should always be .TRUE.
52	LOGICAL :: lldebug = .TRUE.
53
54	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: ztmltrd2 !
55
56	!! * Substitutions
57	# include "zdfddm_substitute.h90"
58	!!----------------------------------------------------------------------
59	!! NEMO/TOP 3.3 , NEMO Consortium (2010)
60	!! $Id$
61	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
62	!!----------------------------------------------------------------------
63	CONTAINS
64
65	INTEGER FUNCTION trd_mxl_trc_alloc()
66	!!----------------------------------------------------------------------
67	!! * ROUTINE trd_mxl_trc_alloc *
68	!!----------------------------------------------------------------------
69	ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) , &
70	& ndextrd1(jpi*jpj), nidtrd(jptra), nh_t(jptra), STAT=trd_mxl_trc_alloc)
71	!
72	IF( lk_mpp ) CALL mpp_sum ( trd_mxl_trc_alloc )
73	IF( trd_mxl_trc_alloc /=0 ) CALL ctl_warn('trd_mxl_trc_alloc: failed to allocate arrays')
74	!
75	END FUNCTION trd_mxl_trc_alloc
76
77
78	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
79	!!----------------------------------------------------------------------
80	!! * ROUTINE trd_mxl_trc_zint *
81	!!
82	!! ** Purpose : Compute the vertical average of the 3D fields given as arguments
83	!! to the subroutine. This vertical average is performed from ocean
84	!! surface down to a chosen control surface.
85	!!
86	!! ** Method/usage :
87	!! The control surface can be either a mixed layer depth (time varying)
88	!! or a fixed surface (jk level or bowl).
89	!! Choose control surface with nctls_trc in namelist NAMTRD :
90	!! nctls_trc = -2 : use isopycnal surface
91	!! nctls_trc = -1 : use euphotic layer with light criterion
92	!! nctls_trc = 0 : use mixed layer with density criterion
93	!! nctls_trc = 1 : read index from file 'ctlsurf_idx'
94	!! nctls_trc > 1 : use fixed level surface jk = nctls_trc
95	!! Note: in the remainder of the routine, the volume between the
96	!! surface and the control surface is called "mixed-layer"
97	!!----------------------------------------------------------------------
98	!!
99	INTEGER, INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
100	CHARACTER(len=2), INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
101	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: ptrc_trdmxl ! passive tracer trend
102	!
103	INTEGER :: ji, jj, jk, isum
104	REAL(wp), POINTER, DIMENSION(:,:) :: zvlmsk
105	!!----------------------------------------------------------------------
106
107	CALL wrk_alloc( jpi, jpj, zvlmsk )
108
109	! I. Definition of control surface and integration weights
110	! --------------------------------------------------------
111
112	ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
113	!
114	tmltrd_trc(:,:,:,:) = 0.e0 ! <<< reset trend arrays to zero
115
116	! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
117	SELECT CASE ( nn_ctls_trc ) ! choice of the control surface
118	CASE ( -2 ) ; STOP 'trdmxl_trc : not ready ' ! -> isopycnal surface (see ???)
119	CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion
120	CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl)
121	CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file
122	CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
123	nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level
124	END SELECT
125
126	! ... Compute ndextrd1 and ndimtrd1 ??? role de jpktrd_trc
127	IF( ionce == 1 ) THEN
128	!
129	isum = 0 ; zvlmsk(:,:) = 0.e0
130
131	IF( jpktrd_trc < jpk ) THEN ! description ???
132	DO jj = 1, jpj
133	DO ji = 1, jpi
134	IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
135	zvlmsk(ji,jj) = tmask(ji,jj,1)
136	ELSE
137	isum = isum + 1
138	zvlmsk(ji,jj) = 0.e0
139	ENDIF
140	END DO
141	END DO
142	ENDIF
143
144	IF( isum > 0 ) THEN ! index of ocean points (2D only)
145	WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum
146	CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
147	ELSE
148	CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
149	ENDIF
150
151	ionce = 0 ! no more pass here
152	!
153	ENDIF ! ionce == 1
154
155	! ... Weights for vertical averaging
156	wkx_trc(:,:,:) = 0.e0
157	DO jk = 1, jpktrd_trc ! initialize wkx_trc with vertical scale factor in mixed-layer
158	DO jj = 1, jpj
159	DO ji = 1, jpi
160	IF( jk - nmld_trc(ji,jj) < 0 ) wkx_trc(ji,jj,jk) = e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
161	END DO
162	END DO
163	END DO
164
165	rmld_trc(:,:) = 0.e0
166	DO jk = 1, jpktrd_trc ! compute mixed-layer depth : rmld_trc
167	rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
168	END DO
169
170	DO jk = 1, jpktrd_trc ! compute integration weights
171	wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
172	END DO
173
174	icount = 0 ! <<< flag = off : control surface & integr. weights
175	! ! computed only once per time step
176	ENDIF ONCE_PER_TIME_STEP
177
178	! II. Vertical integration of trends in the mixed-layer
179	! -----------------------------------------------------
180
181	SELECT CASE ( ctype )
182	CASE ( '3D' ) ! mean passive tracer trends in the mixed-layer
183	DO jk = 1, jpktrd_trc
184	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)
185	END DO
186	CASE ( '2D' ) ! forcing at upper boundary of the mixed-layer
187	tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,1) * wkx_trc(:,:,1) ! non penetrative
188	END SELECT
189	!
190	CALL wrk_dealloc( jpi, jpj, zvlmsk )
191	!
192	END SUBROUTINE trd_mxl_trc_zint
193
194
195	SUBROUTINE trd_mxl_trc( kt )
196	!!----------------------------------------------------------------------
197	!! * ROUTINE trd_mxl_trc *
198	!!
199	!! ** Purpose : Compute and cumulate the mixed layer trends over an analysis
200	!! period, and write NetCDF outputs.
201	!!
202	!! ** Method/usage :
203	!! The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
204	!! logical namelist variable) :
205	!! 1) to explain the difference between initial and final
206	!! mixed-layer T & S (where initial and final relate to the
207	!! current analysis window, defined by ntrc_trc in the namelist)
208	!! 2) to explain the difference between the current and previous
209	!! TIME-AVERAGED mixed-layer T & S (where time-averaging is
210	!! performed over each analysis window).
211	!!
212	!! ** Consistency check :
213	!! If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
214	!! entrainment) should be zero, at machine accuracy. Note that in the case
215	!! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
216	!! over the first two analysis windows (except if restart).
217	!! N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8
218	!! for checking residuals.
219	!! On a NEC-SX5 computer, this typically leads to:
220	!! O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
221	!! O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
222	!!
223	!! ** Action :
224	!! At each time step, mixed-layer averaged trends are stored in the
225	!! tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
226	!! This array is known when trd_mld is called, at the end of the stp subroutine,
227	!! except for the purely vertical K_z diffusion term, which is embedded in the
228	!! lateral diffusion trend.
229	!!
230	!! In I), this K_z term is diagnosed and stored, thus its contribution is removed
231	!! from the lateral diffusion trend.
232	!! In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
233	!! arrays are updated.
234	!! In III), called only once per analysis window, we compute the total trends,
235	!! along with the residuals and the Asselin correction terms.
236	!! In IV), the appropriate trends are written in the trends NetCDF file.
237	!!
238	!! References :
239	!! - Vialard & al.
240	!! - See NEMO documentation (in preparation)
241	!!----------------------------------------------------------------------
242	!
243	INTEGER, INTENT(in) :: kt ! ocean time-step index
244	!
245	INTEGER :: ji, jj, jk, jl, ik, it, itmod, jn
246	REAL(wp) :: zavt, zfn, zfn2
247	!
248	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltot ! d(trc)/dt over the anlysis window (incl. Asselin)
249	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlres ! residual = dh/dt entrainment term
250	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlatf ! for storage only
251	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlrad ! for storage only (for trb<0 corr in trcrad)
252	!
253	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltot2 ! -+
254	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlres2 ! \| working arrays to diagnose the trends
255	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmltrdm2 ! \| associated with the time meaned ML
256	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlatf2 ! \| passive tracers
257	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmlrad2 ! \| (-> for trb<0 corr in trcrad)
258	!
259	CHARACTER (LEN=10) :: clvar
260	!!----------------------------------------------------------------------
261
262	! Set-up pointers into sub-arrays of workspaces
263	CALL wrk_alloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
264	CALL wrk_alloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
265
266	IF( nn_dttrc /= 1 ) CALL ctl_stop( " Be careful, trends diags never validated " )
267
268	! ======================================================================
269	! I. Diagnose the purely vertical (K_z) diffusion trend
270	! ======================================================================
271
272	! ... These terms can be estimated by flux computation at the lower boundary of the ML
273	! (we compute (-1/h) * K_z * d_z( tracer ))
274
275	IF( ln_trcldf_iso ) THEN
276	!
277	DO jj = 1,jpj
278	DO ji = 1,jpi
279	ik = nmld_trc(ji,jj)
280	zavt = fsavs(ji,jj,ik)
281	DO jn = 1, jptra
282	IF( ln_trdtrc(jn) ) &
283	tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - zavt / e3w_n(ji,jj,ik) * tmask(ji,jj,ik) &
284	& * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) ) &
285	& / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
286	END DO
287	END DO
288	END DO
289
290	DO jn = 1, jptra
291	! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
292	IF( ln_trdtrc(jn) ) &
293	tmltrd_trc(:,:,jpmxl_trc_ldf,jn) = tmltrd_trc(:,:,jpmxl_trc_ldf,jn) - tmltrd_trc(:,:,jpmxl_trc_zdf,jn)
294
295	END DO
296	!
297	ENDIF
298
299	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
300	! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm.
301	! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
302	DO jn = 1, jptra
303	IF( ln_trdtrc(jn) ) THEN
304	DO jl = 1, jpltrd_trc
305	CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions
306	END DO
307	ENDIF
308	END DO
309	ENDIF
310	! ======================================================================
311	! II. Cumulate the trends over the analysis window
312	! ======================================================================
313
314	ztmltrd2(:,:,:,:) = 0.e0 ; ztmltot2(:,:,:) = 0.e0 ! <<< reset arrays to zero
315	ztmlres2(:,:,:) = 0.e0 ; ztmlatf2(:,:,:) = 0.e0
316	ztmlrad2(:,:,:) = 0.e0
317
318	! II.1 Set before values of vertically averages passive tracers
319	! -------------------------------------------------------------
320	IF( kt > nittrc000 ) THEN
321	DO jn = 1, jptra
322	IF( ln_trdtrc(jn) ) THEN
323	tmlb_trc (:,:,jn) = tml_trc (:,:,jn)
324	tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_atf,jn)
325	tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_radb,jn)
326	ENDIF
327	END DO
328	ENDIF
329
330	! II.2 Vertically averaged passive tracers
331	! ----------------------------------------
332	tml_trc(:,:,:) = 0.e0
333	DO jk = 1, jpktrd_trc ! - 1 ???
334	DO jn = 1, jptra
335	IF( ln_trdtrc(jn) ) &
336	tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
337	END DO
338	END DO
339
340	! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
341	! ------------------------------------------------------------------------
342	IF( kt == nittrc000 + nn_dttrc ) THEN ! i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1) ???
343	!
344	DO jn = 1, jptra
345	IF( ln_trdtrc(jn) ) THEN
346	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
347	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
348
349	tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0 ; tmltrd_atf_sumb_trc (:,:,jn) = 0.e0
350	tmltrd_rad_sumb_trc (:,:,jn) = 0.e0
351	ENDIF
352	END DO
353
354	rmldbn_trc(:,:) = rmld_trc(:,:)
355	!
356	ENDIF
357
358	! II.4 Cumulated trends over the analysis period
359	! ----------------------------------------------
360	!
361	! [ 1rst analysis window ] [ 2nd analysis window ]
362	!
363	! o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
364	! ntrd 2*ntrd etc.
365	! 1 2 3 4 =5 e.g. =10
366	!
367	IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN ! ???
368	!
369	nmoymltrd = nmoymltrd + 1
370
371
372	! ... Cumulate over BOTH physical contributions AND over time steps
373	DO jn = 1, jptra
374	IF( ln_trdtrc(jn) ) THEN
375	DO jl = 1, jpltrd_trc
376	tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
377	END DO
378	ENDIF
379	END DO
380
381	DO jn = 1, jptra
382	IF( ln_trdtrc(jn) ) THEN
383	! ... Special handling of the Asselin trend
384	tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
385	tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
386
387	! ... Trends associated with the time mean of the ML passive tracers
388	tmltrd_sum_trc (:,:,:,jn) = tmltrd_sum_trc (:,:,:,jn) + tmltrd_trc (:,:,:,jn)
389	tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
390	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) + tml_trc (:,:,jn)
391	ENDIF
392	ENDDO
393
394	rmld_sum_trc (:,:) = rmld_sum_trc (:,:) + rmld_trc (:,:)
395	!
396	ENDIF
397
398	! ======================================================================
399	! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
400	! ======================================================================
401
402	! Convert to appropriate physical units
403	tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
404
405	itmod = kt - nittrc000 + 1
406	it = kt
407
408	MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of nn_trd_trc
409	!
410	ztmltot (:,:,:) = 0.e0 ! reset arrays to zero
411	ztmlres (:,:,:) = 0.e0
412	ztmltot2(:,:,:) = 0.e0
413	ztmlres2(:,:,:) = 0.e0
414
415	zfn = FLOAT( nmoymltrd ) ; zfn2 = zfn * zfn
416
417	! III.1 Prepare fields for output ("instantaneous" diagnostics)
418	! -------------------------------------------------------------
419
420	DO jn = 1, jptra
421	IF( ln_trdtrc(jn) ) THEN
422	!-- Compute total trends (use rdttrc instead of rdt ???)
423	IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes
424	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
425	ELSE ! LEAP-FROG schemes
426	ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
427	ENDIF
428
429	!-- Compute residuals
430	ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
431	& - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn) )
432
433	!-- Diagnose Asselin trend over the analysis window
434	ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
435	ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
436
437	!-- Lateral boundary conditions
438	IF ( cp_cfg .NE. 'gyre' ) THEN
439	CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
440	CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. ) ; CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
441	ENDIF
442
443
444	#if defined key_diainstant
445	STOP 'tmltrd_trc : key_diainstant was never checked within trdmxl. Comment this to proceed.'
446	#endif
447	ENDIF
448	END DO
449
450	! III.2 Prepare fields for output ("mean" diagnostics)
451	! ----------------------------------------------------
452
453	!-- Update the ML depth time sum (to build the Leap-Frog time mean)
454	rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
455
456	!-- Compute passive tracer total trends
457	DO jn = 1, jptra
458	IF( ln_trdtrc(jn) ) THEN
459	tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
460	ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rdt ) ! now tracer unit is /sec
461	ENDIF
462	END DO
463
464	!-- Compute passive tracer residuals
465	DO jn = 1, jptra
466	IF( ln_trdtrc(jn) ) THEN
467	!
468	DO jl = 1, jpltrd_trc
469	ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
470	END DO
471
472	ztmltrdm2(:,:,jn) = 0.e0
473	DO jl = 1, jpltrd_trc
474	ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
475	END DO
476
477	ztmlres2(:,:,jn) = ztmltot2(:,:,jn) - &
478	& ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
479	& - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
480	!
481
482	!-- Diagnose Asselin trend over the analysis window
483	ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) &
484	& + tmltrd_atf_sumb_trc(:,:,jn)
485	ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) &
486	& + tmltrd_rad_sumb_trc(:,:,jn)
487
488	!-- Lateral boundary conditions
489	IF ( cp_cfg .NE. 'gyre' ) THEN ! other than GYRE configuration
490	CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
491	CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
492	DO jl = 1, jpltrd_trc
493	CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. ) ! will be output in the NetCDF trends file
494	END DO
495	ENDIF
496
497	ENDIF
498	END DO
499
500	! * Debugging information *
501	IF( lldebug ) THEN
502	!
503	WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
504	WRITE(numout,*) '~~~~~~~~~~~ '
505	WRITE(numout,*)
506	WRITE(numout,*) 'TRC kt = ', kt, ' nmoymltrd = ', nmoymltrd
507
508	DO jn = 1, jptra
509
510	IF( ln_trdtrc(jn) ) THEN
511	WRITE(numout, *)
512	WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
513
514	WRITE(numout, *)
515	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres : ', SUM2D(ztmlres(:,:,jn))
516	!CD??? PREVOIR: z2d = ztmlres(:,:,jn) ; CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres - : ')
517
518	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
519	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
520	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot : ', SUM2D(+ztmltot (:,:,jn))
521	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
522	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
523	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
524	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
525	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
526	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
527
528	WRITE(numout, *)
529	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2 : ', SUM2D(ztmlres2(:,:,jn))
530	WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
531	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
532	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2 : ', SUM2D(+ztmltot2(:,:,jn))
533	WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2 : ', SUM2D(+ztmltrdm2(:,:,jn))
534	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn))
535	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
536	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
537	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
538	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
539
540	WRITE(numout, *)
541	WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot : ', SUM2D(ztmltot (:,:,jn))
542	WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
543	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc : ', SUM2D(tml_trc (:,:,jn))
544	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc : ', SUM2D(tmlbn_trc (:,:,jn))
545	WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc : ', SUM2D(tmlb_trc (:,:,jn))
546	WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc : ', SUM2D(tmlbb_trc (:,:,jn))
547	WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
548
549	WRITE(numout, *)
550	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
551	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
552	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
553	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
554	WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
555
556	WRITE(numout, *)
557	DO jl = 1, jpltrd_trc
558	WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
559	& ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
560	END DO
561
562	WRITE(numout,*)
563	WRITE(numout,) ' ******************** ZTMLRES, ZTMLRES2 ********************* '
564	WRITE(numout,*)
565	WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
566	WRITE(numout,*)
567	WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
568
569	WRITE(numout,*)
570	WRITE(numout,) ' ******************** ZTMLRES ********************* '
571	WRITE(numout,*)
572
573	WRITE(numout,*) '...................................................'
574	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
575	DO jj = 5, 1, -1
576	WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
577	END DO
578
579	WRITE(numout,*)
580	WRITE(numout,) ' ******************** ZTMLRES2 ********************* '
581	WRITE(numout,*)
582
583	WRITE(numout,*) '...................................................'
584	WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
585	DO jj = 5, 1, -1
586	WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
587	END DO
588	!
589	ENDIF
590	!
591	END DO
592
593
594	97 FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
595	98 FORMAT(a10, i3, 2x, a30, 2x, g20.10)
596	99 FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
597	WRITE(numout,*)
598	!
599	ENDIF
600
601	! III.3 Time evolution array swap
602	! -------------------------------
603	! ML depth
604	rmldbn_trc(:,:) = rmld_trc(:,:)
605	rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum
606	DO jn = 1, jptra
607	IF( ln_trdtrc(jn) ) THEN
608	! For passive tracer instantaneous diagnostics
609	tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn)
610	tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
611
612	! For passive tracer mean diagnostics
613	tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
614	tml_sumb_trc (:,:,jn) = tml_sum_trc(:,:,jn)
615	tmltrd_atf_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
616	tmltrd_rad_sumb_trc(:,:,jn) = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
617
618
619	! III.4 Convert to appropriate physical units
620	! -------------------------------------------
621	ztmltot (:,:,jn) = ztmltot (:,:,jn) * rn_ucf_trc/zfn ! instant diags
622	ztmlres (:,:,jn) = ztmlres (:,:,jn) * rn_ucf_trc/zfn
623	ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * rn_ucf_trc/zfn
624	ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * rn_ucf_trc/zfn
625	tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags
626	ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * rn_ucf_trc/zfn2
627	ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * rn_ucf_trc/zfn2
628	ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * rn_ucf_trc/zfn2
629	ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * rn_ucf_trc/zfn2
630	ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * rn_ucf_trc/zfn2
631	ENDIF
632	END DO
633	!
634	ENDIF MODULO_NTRD
635
636	! ======================================================================
637	! IV. Write trends in the NetCDF file
638	! ======================================================================
639
640	! IV.1 Code for IOIPSL/NetCDF output
641	! ----------------------------------
642
643	IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
644	WRITE(numout,*) ' '
645	WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
646	WRITE(numout,*) '~~~~~~~~~~~ '
647	WRITE(numout,*) ' ', trim(clhstnam), ' at kt = ', kt
648	WRITE(numout,*) ' N.B. nmoymltrd = ', nmoymltrd
649	WRITE(numout,*) ' '
650	ENDIF
651
652	NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN ! <<< write the trends for passive tracer instant. diags
653	!
654
655	DO jn = 1, jptra
656	!
657	IF( ln_trdtrc(jn) ) THEN
658	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
659	!-- Output the fields
660	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
661	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 )
662	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 )
663	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 )
664
665	DO jl = 1, jpltrd_trc - 2
666	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
667	& it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
668	END DO
669
670	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
671	& it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
672
673	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
674	& it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
675
676	ENDIF
677	END DO
678
679	IF( kt == nitend ) THEN
680	DO jn = 1, jptra
681	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
682	END DO
683	ENDIF
684
685	ELSE ! <<< write the trends for passive tracer mean diagnostics
686
687	DO jn = 1, jptra
688	!
689	IF( ln_trdtrc(jn) ) THEN
690	CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 )
691	!-- Output the fields
692	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, nh4ml, etc.
693
694	CALL histwrite( nidtrd(jn), trim(clvar) , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
695	CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot2(:,:,jn), ndimtrd1, ndextrd1 )
696	CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres2(:,:,jn), ndimtrd1, ndextrd1 )
697
698	DO jl = 1, jpltrd_trc - 2
699	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), &
700	& it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
701	END DO
702
703	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), & ! now trcrad : jpltrd_trc - 1
704	& it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
705
706	CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), & ! now Asselin : jpltrd_trc
707	& it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
708
709	ENDIF
710	!
711	END DO
712	IF( kt == nitend ) THEN
713	DO jn = 1, jptra
714	IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) )
715	END DO
716	ENDIF
717
718	!
719	ENDIF NETCDF_OUTPUT
720
721	! Compute the control surface (for next time step) : flag = on
722	icount = 1
723
724	IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
725	!
726	! Reset cumulative arrays to zero
727	! -------------------------------
728	nmoymltrd = 0
729	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
730	tmlradm_trc (:,:,:) = 0.e0 ; tml_sum_trc (:,:,:) = 0.e0
731	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0
732	rmld_sum_trc (:,:) = 0.e0
733	!
734	ENDIF
735
736	! ======================================================================
737	! V. Write restart file
738	! ======================================================================
739
740	IF( lrst_trc ) CALL trd_mxl_trc_rst_write( kt ) ! this must be after the array swap above (III.3)
741
742	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad )
743	CALL wrk_dealloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
744	!
745	END SUBROUTINE trd_mxl_trc
746
747
748	REAL FUNCTION sum2d( ztab )
749	!!----------------------------------------------------------------------
750	!! CD ??? prevoir d'utiliser plutot prtctl
751	!!----------------------------------------------------------------------
752	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: ztab
753	!!----------------------------------------------------------------------
754	sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
755	END FUNCTION sum2d
756
757
758	SUBROUTINE trd_mxl_trc_init
759	!!----------------------------------------------------------------------
760	!! * ROUTINE trd_mxl_init *
761	!!
762	!! ** Purpose : computation of vertically integrated T and S budgets
763	!! from ocean surface down to control surface (NetCDF output)
764	!!
765	!! ** Method/usage :
766	!!
767	!!----------------------------------------------------------------------
768	INTEGER :: inum ! logical unit
769	INTEGER :: ilseq, jl, jn, iiter
770	REAL(wp) :: zjulian, zsto, zout
771	CHARACTER (LEN=40) :: clop
772	CHARACTER (LEN=15) :: csuff
773	CHARACTER (LEN=12) :: clmxl
774	CHARACTER (LEN=16) :: cltrcu
775	CHARACTER (LEN=10) :: clvar
776
777	!!----------------------------------------------------------------------
778
779	! ======================================================================
780	! I. initialization
781	! ======================================================================
782
783	IF(lwp) THEN
784	WRITE(numout,*)
785	WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers '
786	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
787	WRITE(numout,*)
788	ENDIF
789
790
791	! I.1 Check consistency of user defined preferences
792	! -------------------------------------------------
793
794	IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
795	WRITE(numout,cform_err)
796	WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend
797	WRITE(numout,*) ' is no multiple of the trends diagnostics frequency '
798	WRITE(numout,*) ' you defined, nn_trd_trc = ', nn_trd_trc
799	WRITE(numout,*) ' This will not allow you to restart from this simulation. '
800	WRITE(numout,*) ' You should reconsider this choice. '
801	WRITE(numout,*)
802	WRITE(numout,*) ' N.B. the nitend parameter is also constrained to be a '
803	WRITE(numout,*) ' multiple of the sea-ice frequency parameter (typically 5) '
804	nstop = nstop + 1
805	ENDIF
806
807	! * Debugging information *
808	IF( lldebug ) THEN
809	WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl
810	WRITE(numout,*) ' ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
811	ENDIF
812
813	IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
814	WRITE(numout,cform_err)
815	WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer MUSCL '
816	WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '
817	WRITE(numout,*) ' WHY? Everything in trdmxl_trc is coded for leap-frog, and '
818	WRITE(numout,*) ' MUSCL scheme is Euler forward for passive tracers (note '
819	WRITE(numout,*) ' that MUSCL is leap-frog for active tracers T/S). '
820	WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '
821	WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'
822	WRITE(numout,*)
823	nstop = nstop + 1
824	ENDIF
825
826	! I.2 Initialize arrays to zero or read a restart file
827	! ----------------------------------------------------
828	nmoymltrd = 0
829
830	rmld_trc (:,:) = 0.e0 ; tml_trc (:,:,:) = 0.e0 ! inst.
831	tmltrdm_trc (:,:,:) = 0.e0 ; tmlatfm_trc (:,:,:) = 0.e0
832	tmlradm_trc (:,:,:) = 0.e0
833
834	tml_sum_trc (:,:,:) = 0.e0 ; tmltrd_sum_trc (:,:,:,:) = 0.e0 ! mean
835	tmltrd_csum_ln_trc (:,:,:,:) = 0.e0 ; rmld_sum_trc (:,:) = 0.e0
836
837	IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
838	CALL trd_mxl_trc_rst_read
839	ELSE
840	tmlb_trc (:,:,:) = 0.e0 ; tmlbb_trc (:,:,:) = 0.e0 ! inst.
841	tmlbn_trc (:,:,:) = 0.e0
842
843	tml_sumb_trc (:,:,:) = 0.e0 ; tmltrd_csum_ub_trc (:,:,:,:) = 0.e0 ! mean
844	tmltrd_atf_sumb_trc(:,:,:) = 0.e0 ; tmltrd_rad_sumb_trc(:,:,:) = 0.e0
845
846	ENDIF
847
848	icount = 1 ; ionce = 1 ! open specifier
849
850
851	! I.3 Read control surface from file ctlsurf_idx
852	! ----------------------------------------------
853	IF( nn_ctls_trc == 1 ) THEN
854	CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
855	READ ( inum ) nbol_trc
856	CLOSE( inum )
857	ENDIF
858
859	! ======================================================================
860	! II. netCDF output initialization
861	! ======================================================================
862
863	! clmxl = legend root for netCDF output
864	IF( nn_ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion
865	clmxl = 'Mixed Layer '
866	ELSE IF( nn_ctls_trc == 1 ) THEN ! control surface = read index from file
867	clmxl = ' Bowl '
868	ELSE IF( nn_ctls_trc >= 2 ) THEN ! control surface = model level
869	WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
870	ENDIF
871
872	! II.1 Define frequency of output and means
873	! -----------------------------------------
874	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
875	ELSE ; clop = "x" ! no use of the mask value (require less cp time)
876	ENDIF
877	# if defined key_diainstant
878	IF( .NOT. ln_trdmxl_trc_instant ) THEN
879	STOP 'trd_mxl_trc : this was never checked. Comment this line to proceed...'
880	ENDIF
881	zsto = nn_trd_trc * rdt
882	clop = "inst("//TRIM(clop)//")"
883	# else
884	IF( ln_trdmxl_trc_instant ) THEN
885	zsto = rdt ! inst. diags : we use IOIPSL time averaging
886	ELSE
887	zsto = nn_trd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging
888	ENDIF
889	clop = "ave("//TRIM(clop)//")"
890	# endif
891	zout = nn_trd_trc * rdt
892	iiter = ( nittrc000 - 1 ) / nn_dttrc
893
894	IF(lwp) WRITE (numout,*) ' netCDF initialization'
895
896	! II.2 Compute julian date from starting date of the run
897	! ------------------------------------------------------
898	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
899	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
900	IF(lwp) WRITE(numout,*)' '
901	IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000 &
902	& ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday &
903	& ,'Julian day : ', zjulian
904
905	! II.3 Define the T grid trend file (nidtrd)
906	! ------------------------------------------
907
908	!-- Define long and short names for the NetCDF output variables
909	! ==> choose them according to trdmxl_trc_oce.F90 <==
910
911	ctrd_trc(jpmxl_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmxl_trc_xad ,2) = "_xad"
912	ctrd_trc(jpmxl_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmxl_trc_yad ,2) = "_yad"
913	ctrd_trc(jpmxl_trc_zad ,1) = " Vertical advection" ; ctrd_trc(jpmxl_trc_zad ,2) = "_zad"
914	ctrd_trc(jpmxl_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmxl_trc_ldf ,2) = "_ldf"
915	ctrd_trc(jpmxl_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmxl_trc_zdf ,2) = "_zdf"
916	ctrd_trc(jpmxl_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmxl_trc_bbl ,2) = "_bbl"
917	ctrd_trc(jpmxl_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmxl_trc_dmp ,2) = "_dmp"
918	ctrd_trc(jpmxl_trc_sbc ,1) = " Surface boundary cond." ; ctrd_trc(jpmxl_trc_sbc ,2) = "_sbc"
919	ctrd_trc(jpmxl_trc_sms, 1) = " Sources minus sinks" ; ctrd_trc(jpmxl_trc_sms ,2) = "_sms"
920	ctrd_trc(jpmxl_trc_radb ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radb ,2) = "_radb"
921	ctrd_trc(jpmxl_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmxl_trc_radn ,2) = "_radn"
922	ctrd_trc(jpmxl_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmxl_trc_atf ,2) = "_atf"
923
924	DO jn = 1, jptra
925	!-- Create a NetCDF file and enter the define mode
926	IF( ln_trdtrc(jn) ) THEN
927	csuff="ML_"//ctrcnm(jn)
928	CALL dia_nam( clhstnam, nn_trd_trc, csuff )
929	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, &
930	& 1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
931
932	!-- Define the ML depth variable
933	CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m", &
934	& jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
935
936	ENDIF
937	END DO
938
939	!-- Define physical units
940	IF( rn_ucf_trc == 1. ) THEN
941	cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit
942	ELSEIF ( rn_ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3)
943	cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf
944	ELSE
945	cltrcu = "unknown?"
946	ENDIF
947
948	!-- Define miscellaneous passive tracer mixed-layer variables
949	IF( jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
950	STOP 'Error : jpltrd_trc /= jpmxl_trc_atf .OR. jpltrd_trc - 1 /= jpmxl_trc_radb' ! see below
951	ENDIF
952
953	DO jn = 1, jptra
954	!
955	IF( ln_trdtrc(jn) ) THEN
956	clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc.
957	CALL histdef(nidtrd(jn), trim(clvar), clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", &
958	& "mmole-N/m3", jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout )
959	CALL histdef(nidtrd(jn), trim(clvar)//"_tot" , clmxl//" "//trim(ctrcnm(jn))//" Total trend ", &
960	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
961	CALL histdef(nidtrd(jn), trim(clvar)//"_res" , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)", &
962	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout )
963
964	DO jl = 1, jpltrd_trc - 2 ! <== only true if jpltrd_trc == jpmxl_trc_atf
965	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1), &
966	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
967	END DO ! if zsto=rdt above
968
969	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), &
970	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
971
972	CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1), &
973	& cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
974	!
975	ENDIF
976	END DO
977
978	!-- Leave IOIPSL/NetCDF define mode
979	DO jn = 1, jptra
980	IF( ln_trdtrc(jn) ) CALL histend( nidtrd(jn), snc4set )
981	END DO
982
983	IF(lwp) WRITE(numout,*)
984
985	END SUBROUTINE trd_mxl_trc_init
986
987	#else
988	!!----------------------------------------------------------------------
989	!! Default option : Empty module
990	!!----------------------------------------------------------------------
991	CONTAINS
992	SUBROUTINE trd_mxl_trc( kt ) ! Empty routine
993	INTEGER, INTENT( in) :: kt
994	WRITE(,) 'trd_mxl_trc: You should not have seen this print! error?', kt
995	END SUBROUTINE trd_mxl_trc
996	SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
997	INTEGER , INTENT( in ) :: ktrd, kjn ! ocean trend index and passive tracer rank
998	CHARACTER(len=2) , INTENT( in ) :: ctype ! surface/bottom (2D) or interior (3D) physics
999	REAL, DIMENSION(:,:,:), INTENT( in ) :: ptrc_trdmxl ! passive trc trend
1000	WRITE(,) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
1001	WRITE(,) ' " " : You should not have seen this print! error?', ctype
1002	WRITE(,) ' " " : You should not have seen this print! error?', ktrd
1003	WRITE(,) ' " " : You should not have seen this print! error?', kjn
1004	END SUBROUTINE trd_mxl_trc_zint
1005	SUBROUTINE trd_mxl_trc_init ! Empty routine
1006	WRITE(,) 'trd_mxl_trc_init: You should not have seen this print! error?'
1007	END SUBROUTINE trd_mxl_trc_init
1008	#endif
1009
1010	!!======================================================================
1011	END MODULE trdmxl_trc

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

Download in other formats: