source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

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