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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trdmxl_trc.F90 in branches/2015/dev_merge_2015/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 6060

Last change on this file since 6060 was 6060, checked in by timgraham, 8 years ago

Merged dev_r5836_noc2_VVL_BY_DEFAULT into branch

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