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/UKMO/r5518_sst_landsea_cpl/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/UKMO/r5518_sst_landsea_cpl/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 7147

Last change on this file since 7147 was 7147, checked in by jcastill, 7 years ago

Remove svn keywords

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