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

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

source: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

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