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

source: trunk/NEMO/TOP_SRC/TRP/trdmld_trc.F90 @ 1174

Last change on this file since 1174 was 1174, checked in by cetlod, 16 years ago

New trends diagnostics organization in TOP, see ticket:248

File size: 60.2 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 trp_trc               ! tracer definitions (trn, trb, tra, etc.)
19   USE oce_trc           ! needed for namelist logicals, and euphotic layer arrays
20   USE trctrp_lec
21   USE trdmld_trc_oce        ! definition of main arrays used for trends computations
22   USE in_out_manager    ! I/O manager
23   USE dianam            ! build the name of file (routine)
24   USE ldfslp            ! iso-neutral slopes
25   USE ioipsl            ! NetCDF library
26   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
27   USE trdmld_trc_rst    ! restart for diagnosing the ML trends
28   USE prtctl            ! print control
29   USE sms_pisces       
30   USE sms_lobster
31   USE trcsms_cfc
32   USE trc
33   USE trcrst            ! for lrst_trc -> circ. dep. ??? we put lrst_trc in trc_oce
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC trd_mod_trc                                             ! routine called by step.F90
39   PUBLIC trd_mld_trc
40   PUBLIC trd_mld_trc_init
41
42   CHARACTER (LEN=40) ::  clhstnam                                ! name of the trends NetCDF file
43   INTEGER ::   nmoymltrd
44   INTEGER ::   ndextrd1(jpi*jpj)
45   INTEGER, DIMENSION(jptra) ::   nidtrd, nh_t
46   INTEGER ::   ndimtrd1                       
47   INTEGER, SAVE ::  ionce, icount
48   LOGICAL :: llwarn  = .TRUE.                                    ! this should always be .TRUE.
49   LOGICAL :: lldebug = .TRUE.
50
51   !! * Substitutions
52#  include "top_substitute.h90"
53   !!----------------------------------------------------------------------
54   !!   TOP 1.0 , LOCEAN-IPSL (2007)
55   !! $Header:  $
56   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58
59CONTAINS
60
61   SUBROUTINE trd_mod_trc( ptrtrd, kjn, ktrd, kt )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE trd_mod_trc  ***
64      !!----------------------------------------------------------------------
65#if defined key_trcbbl_adv
66      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zun, zvn                   ! temporary arrays
67#else
68      USE oce_trc,   zun => un                                            ! When no bbl, zun == un
69      USE oce_trc,   zvn => vn                                            ! When no bbl, zvn == vn
70#endif
71      INTEGER, INTENT( in )  ::   kt                                  ! time step
72      INTEGER, INTENT( in )  ::   kjn                                 ! tracer index
73      INTEGER, INTENT( in )  ::   ktrd                                ! tracer trend index
74      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
75      !!----------------------------------------------------------------------
76
77      IF( kt == nittrc000 ) THEN
78!         IF(lwp)WRITE(numout,*)
79!         IF(lwp)WRITE(numout,*) 'trd_mod_trc:'
80!         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~'
81      ENDIF
82
83      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
84      ! Mixed layer trends for passive tracers
85      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
86
87      SELECT CASE ( ktrd )
88         CASE ( jptrc_trd_xad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xad    , '3D', kjn )
89         CASE ( jptrc_trd_yad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yad    , '3D', kjn )
90         CASE ( jptrc_trd_zad     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zad    , '3D', kjn )
91         CASE ( jptrc_trd_ldf     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf    , '3D', kjn )
92         CASE ( jptrc_trd_xei     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xei    , '3D', kjn )
93         CASE ( jptrc_trd_yei     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yei    , '3D', kjn )
94         CASE ( jptrc_trd_bbl     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbl    , '3D', kjn )
95         CASE ( jptrc_trd_zdf     )
96            IF( ln_trcldf_iso ) THEN
97               CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf, '3D', kjn )
98            ELSE
99               CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zdf, '3D', kjn )
100            ENDIF
101         CASE ( jptrc_trd_zei     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zei    , '3D', kjn )
102         CASE ( jptrc_trd_dmp     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_dmp    , '3D', kjn )
103         CASE ( jptrc_trd_sbc     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sbc    , '2D', kjn )
104#if defined key_lobster
105         CASE ( jptrc_trd_sms_sed )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_sed, '3D', kjn )
106         CASE ( jptrc_trd_sms_bio )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_bio, '3D', kjn )
107         CASE ( jptrc_trd_sms_exp )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms_exp, '3D', kjn )
108#else
109         CASE ( jptrc_trd_sms     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms    , '3D', kjn )
110#endif
111         CASE ( jptrc_trd_bbc     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbc    , '3D', kjn )
112         CASE ( jptrc_trd_radb    )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radb   , '3D', kjn )
113         CASE ( jptrc_trd_radn    )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radn   , '3D', kjn )
114         CASE ( jptrc_trd_atf     )   ;   CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_atf    , '3D', kjn )
115      END SELECT
116
117
118   END SUBROUTINE trd_mod_trc
119
120   SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
121      !!----------------------------------------------------------------------
122      !!                  ***  ROUTINE trd_mld_trc_zint  ***
123      !!
124      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments
125      !!                to the subroutine. This vertical average is performed from ocean
126      !!                surface down to a chosen control surface.
127      !!
128      !! ** Method/usage :
129      !!      The control surface can be either a mixed layer depth (time varying)
130      !!      or a fixed surface (jk level or bowl).
131      !!      Choose control surface with nctls_trc in namelist NAMTRD :
132      !!        nctls_trc = -2  : use isopycnal surface
133      !!        nctls_trc = -1  : use euphotic layer with light criterion
134      !!        nctls_trc =  0  : use mixed layer with density criterion
135      !!        nctls_trc =  1  : read index from file 'ctlsurf_idx'
136      !!        nctls_trc >  1  : use fixed level surface jk = nctls_trc
137      !!      Note: in the remainder of the routine, the volume between the
138      !!            surface and the control surface is called "mixed-layer"
139      !!----------------------------------------------------------------------
140      INTEGER, INTENT( in ) ::   ktrd, kjn                        ! ocean trend index and passive tracer rank
141      CHARACTER(len=2), INTENT( in ) ::  ctype                    ! surface/bottom (2D) or interior (3D) physics
142      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::  ptrc_trdmld ! passive tracer trend
143      INTEGER ::   ji, jj, jk, isum
144      REAL(wp), DIMENSION(jpi,jpj) ::  zvlmsk
145      !!----------------------------------------------------------------------
146
147      ! I. Definition of control surface and integration weights
148      ! --------------------------------------------------------
149
150      ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
151         !
152         tmltrd_trc(:,:,:,:) = 0.e0                               ! <<< reset trend arrays to zero
153         
154         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
155         SELECT CASE ( nctls_trc )                                ! choice of the control surface
156            CASE ( -2  )   ;   STOP 'trdmld_trc : not ready '     !     -> isopycnal surface (see ???)
157#if defined key_pisces || defined key_lobster
158            CASE ( -1  )   ;   nmld_trc(:,:) = neln(:,:)          !     -> euphotic layer with light criterion
159#endif
160            CASE (  0  )   ;   nmld_trc(:,:) = nmln(:,:)          !     -> ML with density criterion (see zdfmxl)
161            CASE (  1  )   ;   nmld_trc(:,:) = nbol_trc(:,:)          !     -> read index from file
162            CASE (  2: )   ;   nctls_trc = MIN( nctls_trc, jpktrd_trc - 1 )
163                               nmld_trc(:,:) = nctls_trc + 1      !     -> model level
164         END SELECT
165
166         ! ... Compute ndextrd1 and ndimtrd1  ??? role de jpktrd_trc
167         IF( ionce == 1 ) THEN
168            !
169            isum  = 0   ;   zvlmsk(:,:) = 0.e0
170
171            IF( jpktrd_trc < jpk ) THEN                           ! description ???
172               DO jj = 1, jpj
173                  DO ji = 1, jpi
174                     IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
175                        zvlmsk(ji,jj) = tmask(ji,jj,1)
176                     ELSE
177                        isum = isum + 1
178                        zvlmsk(ji,jj) = 0.e0
179                     ENDIF
180                  END DO
181               END DO
182            ENDIF
183
184            IF( isum > 0 ) THEN                                   ! index of ocean points (2D only)
185               WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum 
186               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
187            ELSE
188               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
189            ENDIF                               
190
191            ionce = 0                                             ! no more pass here
192            !
193         ENDIF   ! ionce == 1
194         
195         ! ... Weights for vertical averaging
196         wkx_trc(:,:,:) = 0.e0
197         DO jk = 1, jpktrd_trc                                    ! initialize wkx_trc with vertical scale factor in mixed-layer
198            DO jj = 1, jpj
199               DO ji = 1, jpi
200                  IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
201               END DO
202            END DO
203         END DO
204         
205         rmld_trc(:,:) = 0.e0
206         DO jk = 1, jpktrd_trc                                    ! compute mixed-layer depth : rmld_trc
207            rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
208         END DO
209         
210         DO jk = 1, jpktrd_trc                                    ! compute integration weights
211            wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
212         END DO
213
214         icount = 0                                               ! <<< flag = off : control surface & integr. weights
215         !                                                        !     computed only once per time step
216      ENDIF ONCE_PER_TIME_STEP
217
218      ! II. Vertical integration of trends in the mixed-layer
219      ! -----------------------------------------------------
220
221      SELECT CASE ( ctype )
222         CASE ( '3D' )                                            ! mean passive tracer trends in the mixed-layer
223            DO jk = 1, jpktrd_trc
224               tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,jk) * wkx_trc(:,:,jk)   
225            END DO
226         CASE ( '2D' )                                            ! forcing at upper boundary of the mixed-layer
227            tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmld(:,:,1) * wkx_trc(:,:,1)  ! non penetrative
228      END SELECT
229
230    END SUBROUTINE trd_mld_trc_zint
231   
232
233    SUBROUTINE trd_mld_trc( kt )
234      !!----------------------------------------------------------------------
235      !!                  ***  ROUTINE trd_mld_trc  ***
236      !!
237      !! ** Purpose :  Compute and cumulate the mixed layer trends over an analysis
238      !!               period, and write NetCDF (or dimg) outputs.
239      !!
240      !! ** Method/usage :
241      !!          The stored trends can be chosen twofold (according to the ln_trdmld_trc_instant
242      !!          logical namelist variable) :
243      !!          1) to explain the difference between initial and final
244      !!             mixed-layer T & S (where initial and final relate to the
245      !!             current analysis window, defined by ntrc_trc in the namelist)
246      !!          2) to explain the difference between the current and previous
247      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
248      !!             performed over each analysis window).
249      !!
250      !! ** Consistency check :
251      !!        If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
252      !!        entrainment) should be zero, at machine accuracy. Note that in the case
253      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
254      !!        over the first two analysis windows (except if restart).
255      !!        N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, ucf_trc=1., nctls_trc=8
256      !!             for checking residuals.
257      !!             On a NEC-SX5 computer, this typically leads to:
258      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmld_trc_instant=.false.
259      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmld_trc_instant=.true.
260      !!
261      !! ** Action :
262      !!       At each time step, mixed-layer averaged trends are stored in the
263      !!       tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
264      !!       This array is known when trd_mld is called, at the end of the stp subroutine,
265      !!       except for the purely vertical K_z diffusion term, which is embedded in the
266      !!       lateral diffusion trend.
267      !!
268      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
269      !!       from the lateral diffusion trend.
270      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
271      !!       arrays are updated.
272      !!       In III), called only once per analysis window, we compute the total trends,
273      !!       along with the residuals and the Asselin correction terms.
274      !!       In IV), the appropriate trends are written in the trends NetCDF file.
275      !!
276      !! References :
277      !!       - Vialard & al.
278      !!       - See NEMO documentation (in preparation)
279      !!----------------------------------------------------------------------
280      INTEGER, INTENT( in ) ::   kt                               ! ocean time-step index
281      INTEGER ::   ji, jj, jk, jl, ik, it, jn
282      REAL(wp) ::   zavt, zfn, zfn2
283      !!
284      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltot             ! d(trc)/dt over the anlysis window (incl. Asselin)
285      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlres             ! residual = dh/dt entrainment term
286      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlatf             ! for storage only
287      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlrad             ! for storage only (for trb<0 corr in trcrad)
288      !!
289      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltot2            ! -+
290      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlres2            !  | working arrays to diagnose the trends
291      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltrdm2           !  | associated with the time meaned ML
292      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlatf2            !  | passive tracers
293      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlrad2            !  | (-> for trb<0 corr in trcrad)
294      REAL(wp), DIMENSION(jpi,jpj,jpltrd_trc,jptra) ::  ztmltrd2  ! -+
295      !!
296      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                       ! temporary array, used for eiv arrays
297      CHARACTER (LEN= 5) ::   clvar
298#if defined key_dimgout
299      INTEGER ::   iyear,imon,iday
300      CHARACTER(LEN=80) ::   cltext, clmode
301#endif
302      !!----------------------------------------------------------------------
303
304      IF( llwarn ) THEN                                           ! warnings
305         IF(      ( nittrc000 /= nit000   ) &
306              .OR.( ndttrc    /= 1        )    ) THEN
307
308            WRITE(numout,*) 'Be careful, trends diags never validated'
309            STOP 'Uncomment this line to proceed'
310         ENDIF
311      ENDIF
312
313      ! ======================================================================
314      ! I. Diagnose the purely vertical (K_z) diffusion trend
315      ! ======================================================================
316
317      ! ... These terms can be estimated by flux computation at the lower boundary of the ML
318      !     (we compute (-1/h) * K_z * d_z( tracer ))
319
320      IF( ln_trcldf_iso ) THEN
321         !
322         DO jj = 1,jpj
323            DO ji = 1,jpi
324               ik = nmld_trc(ji,jj)
325               zavt = avt(ji,jj,ik)
326               DO jn = 1, jptra
327                  IF( luttrd(jn) )    &
328                  tmltrd_trc(ji,jj,jpmld_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
329                       &                    * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) )            &
330                       &                    / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
331               END DO
332            END DO
333         END DO
334
335         DO jn = 1, jptra
336         ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
337            IF( luttrd(jn) ) &
338                 tmltrd_trc(:,:,jpmld_trc_ldf,jn) = tmltrd_trc(:,:,jpmld_trc_ldf,jn) - tmltrd_trc(:,:,jpmld_trc_zdf,jn)
339   
340         END DO
341         !     
342      ENDIF
343
344#if ! defined key_gyre
345      ! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm.
346      ! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
347      DO jn = 1, jptra
348         IF( luttrd(jn) ) THEN
349            DO jl = 1, jpltrd_trc
350               CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. )        ! lateral boundary conditions
351            END DO
352         ENDIF
353      END DO
354#endif
355      ! ======================================================================
356      ! II. Cumulate the trends over the analysis window
357      ! ======================================================================
358
359      ztmltrd2(:,:,:,:) = 0.e0   ;   ztmltot2(:,:,:)   = 0.e0     ! <<< reset arrays to zero
360      ztmlres2(:,:,:)   = 0.e0   ;   ztmlatf2(:,:,:)   = 0.e0
361      ztmlrad2(:,:,:)   = 0.e0
362
363      ! II.1 Set before values of vertically averages passive tracers
364      ! -------------------------------------------------------------
365      IF( kt > nittrc000 ) THEN
366         DO jn = 1, jptra
367            IF( luttrd(jn) ) THEN
368               tmlb_trc   (:,:,jn) = tml_trc   (:,:,jn)
369               tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_atf,jn)
370               tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_radb,jn)
371            ENDIF
372         END DO
373      ENDIF
374
375      ! II.2 Vertically averaged passive tracers
376      ! ----------------------------------------
377      tml_trc(:,:,:) = 0.e0
378      DO jk = 1, jpktrd_trc ! - 1 ???
379         DO jn = 1, jptra
380            IF( luttrd(jn) ) &
381               tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn)
382         END DO
383      END DO
384
385      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window   
386      ! ------------------------------------------------------------------------
387      IF( kt == 2 ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)    ???
388         !
389         DO jn = 1, jptra
390            IF( luttrd(jn) ) THEN
391               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
392               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
393               
394               tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0   ;   tmltrd_atf_sumb_trc  (:,:,jn) = 0.e0
395               tmltrd_rad_sumb_trc  (:,:,jn) = 0.e0
396            ENDIF
397         END DO
398         
399         rmldbn_trc(:,:) = rmld_trc(:,:)
400         !
401      ENDIF
402
403      ! II.4 Cumulated trends over the analysis period
404      ! ----------------------------------------------
405      !
406      !         [  1rst analysis window ] [     2nd analysis window     ]                       
407      !
408      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
409      !                            ntrd                             2*ntrd       etc.
410      !     1      2     3     4    =5 e.g.                          =10
411      !
412      IF( ( kt >= 2 ).OR.( lrsttr ) ) THEN                        ! ???
413         !
414         nmoymltrd = nmoymltrd + 1
415
416
417         ! ... Cumulate over BOTH physical contributions AND over time steps
418         DO jn = 1, jptra
419            IF( luttrd(jn) ) THEN
420               DO jl = 1, jpltrd_trc
421                  tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
422               END DO
423            ENDIF
424         END DO
425
426         DO jn = 1, jptra
427            IF( luttrd(jn) ) THEN
428               ! ... Special handling of the Asselin trend
429               tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
430               tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
431
432               ! ... Trends associated with the time mean of the ML passive tracers
433               tmltrd_sum_trc    (:,:,:,jn) = tmltrd_sum_trc    (:,:,:,jn) + tmltrd_trc    (:,:,:,jn)
434               tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
435               tml_sum_trc       (:,:,jn)   = tml_sum_trc       (:,:,jn)   + tml_trc       (:,:,jn)
436            ENDIF
437         ENDDO
438
439         rmld_sum_trc      (:,:)     = rmld_sum_trc      (:,:)     + rmld_trc      (:,:)
440         !
441      ENDIF
442
443      ! ======================================================================
444      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
445      ! ======================================================================
446
447      ! Convert to appropriate physical units
448      tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * ucf_trc
449
450      MODULO_NTRD : IF( MOD( kt, ntrd_trc ) == 0 ) THEN           ! nitend MUST be multiple of ntrd_trc
451         !
452         ztmltot (:,:,:) = 0.e0                                   ! reset arrays to zero
453         ztmlres (:,:,:) = 0.e0
454         ztmltot2(:,:,:) = 0.e0
455         ztmlres2(:,:,:) = 0.e0
456     
457         zfn  = FLOAT( nmoymltrd )    ;    zfn2 = zfn * zfn
458         
459         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
460         ! -------------------------------------------------------------
461
462         DO jn = 1, jptra
463            IF( luttrd(jn) ) THEN
464               !-- Compute total trends    (use rdttrc instead of rdt ???)
465               IF ( ln_trcadv_smolar .OR. ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN  ! EULER-FORWARD schemes
466                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
467               ELSE                                                                     ! LEAP-FROG schemes
468                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
469               ENDIF
470               
471               !-- Compute residuals
472               ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
473                  &                                                 - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)  )
474               
475               !-- Diagnose Asselin trend over the analysis window
476               ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
477               ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
478               
479         !-- Lateral boundary conditions
480#if ! defined key_gyre
481
482               CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
483               CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
484
485#endif
486
487#if defined key_diainstant
488               STOP 'tmltrd_trc : key_diainstant was never checked within trdmld. Comment this to proceed.'
489#endif
490            ENDIF
491         END DO
492
493         ! III.2 Prepare fields for output ("mean" diagnostics)
494         ! ----------------------------------------------------
495         
496         !-- Update the ML depth time sum (to build the Leap-Frog time mean)
497         rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
498
499               !-- Compute passive tracer total trends
500         DO jn = 1, jptra
501            IF( luttrd(jn) ) THEN
502               tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
503               ztmltot2   (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) /  ( 2.*rdt )    ! now tracer unit is /sec
504            ENDIF
505         END DO
506
507         !-- Compute passive tracer residuals
508         DO jn = 1, jptra
509            IF( luttrd(jn) ) THEN
510               !
511               DO jl = 1, jpltrd_trc
512                  ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
513               END DO
514               
515               ztmltrdm2(:,:,jn) = 0.e0
516               DO jl = 1, jpltrd_trc
517                  ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
518               END DO
519               
520               ztmlres2(:,:,jn) =  ztmltot2(:,:,jn)  -       &
521                  & ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
522                  &                     - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
523               !
524
525               !-- Diagnose Asselin trend over the analysis window
526               ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) &
527                  &                                               + tmltrd_atf_sumb_trc(:,:,jn)
528               ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) &
529                  &                                               + tmltrd_rad_sumb_trc(:,:,jn)
530
531         !-- Lateral boundary conditions
532#if ! defined key_gyre         
533               CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
534               CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
535               DO jl = 1, jpltrd_trc
536                  CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
537               END DO
538#endif
539            ENDIF
540         END DO
541
542         ! * Debugging information *
543         IF( lldebug ) THEN
544            !
545            WRITE(numout,*) 'trd_mld_trc : write trends in the Mixed Layer for debugging process:'
546            WRITE(numout,*) '~~~~~~~~~~~  '
547            WRITE(numout,*)
548            WRITE(numout,*) 'TRC kt = ', kt, '    nmoymltrd = ', nmoymltrd
549
550            DO jn = 1, jptra
551
552               IF( luttrd(jn) ) THEN
553                  WRITE(numout, *)
554                  WRITE(numout, *) '>>>>>>>>>>>>>>>>>>  TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
555                 
556                  WRITE(numout, *)
557                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres     : ', SUM2D(ztmlres(:,:,jn))
558                  !CD??? PREVOIR: z2d = ztmlres(:,:,jn)   ;   CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres   -   : ')
559                 
560                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
561                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
562                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot    : ', SUM2D(+ztmltot    (:,:,jn))
563                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
564                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
565                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
566                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
567                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
568                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
569                 
570                  WRITE(numout, *)
571                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2    : ', SUM2D(ztmlres2(:,:,jn))
572                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
573                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
574                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2            : ', SUM2D(+ztmltot2(:,:,jn))
575                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2           : ', SUM2D(+ztmltrdm2(:,:,jn)) 
576                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_atf,jn)) 
577                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
578                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_radb,jn))
579                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
580                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
581                 
582                  WRITE(numout, *)
583                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot     : ', SUM2D(ztmltot    (:,:,jn))
584                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
585                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc    : ', SUM2D(tml_trc    (:,:,jn))
586                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc  : ', SUM2D(tmlbn_trc  (:,:,jn))
587                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc   : ', SUM2D(tmlb_trc   (:,:,jn))
588                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc  : ', SUM2D(tmlbb_trc  (:,:,jn))
589                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
590                 
591                  WRITE(numout, *)
592                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
593                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
594                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
595                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
596                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
597                 
598                  WRITE(numout, *)
599                  DO jl = 1, jpltrd_trc
600                     WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmld_trc_xxx = ', jl, &
601                        & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
602                  END DO
603               
604                  WRITE(numout,*) 
605                  WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** '
606                  WRITE(numout,*)
607                  WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
608                  WRITE(numout,*)
609                  WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
610                 
611                  WRITE(numout,*) 
612                  WRITE(numout,*) ' *********************** ZTMLRES *********************** '
613                  WRITE(numout,*)
614                 
615                  WRITE(numout,*) '...................................................'
616                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
617                  DO jj = 5, 1, -1
618                     WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
619                  END DO
620                 
621                  WRITE(numout,*) 
622                  WRITE(numout,*) ' *********************** ZTMLRES2 *********************** '
623                  WRITE(numout,*)
624                 
625                  WRITE(numout,*) '...................................................'
626                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
627                  DO jj = 5, 1, -1
628                     WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
629                  END DO
630                  !
631               ENDIF
632               !
633            END DO
634
635
63697            FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
63798            FORMAT(a10, i3, 2x, a30, 2x, g20.10)
63899            FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
639              WRITE(numout,*)
640            !
641         ENDIF
642
643         ! III.3 Time evolution array swap
644         ! -------------------------------
645         ! ML depth
646         rmldbn_trc(:,:)   = rmld_trc(:,:)
647         rmld_sum_trc(:,:)     = rmld_sum_trc(:,:)     /      (2*zfn)  ! similar to tml_sum and sml_sum
648         DO jn = 1, jptra
649            IF( luttrd(jn) ) THEN       
650               ! For passive tracer instantaneous diagnostics
651               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
652               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
653               
654               ! For passive tracer mean diagnostics
655               tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
656               tml_sumb_trc       (:,:,jn)   = tml_sum_trc(:,:,jn)
657               tmltrd_atf_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn)
658               tmltrd_rad_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmld_trc_radb,jn)
659               
660               
661               ! III.4 Convert to appropriate physical units
662               ! -------------------------------------------
663               ztmltot     (:,:,jn)   = ztmltot     (:,:,jn)   * ucf_trc/zfn   ! instant diags
664               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * ucf_trc/zfn
665               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * ucf_trc/zfn
666               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * ucf_trc/zfn
667               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
668               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * ucf_trc/zfn2
669               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * ucf_trc/zfn2
670               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * ucf_trc/zfn2
671               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * ucf_trc/zfn2
672               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * ucf_trc/zfn2
673            ENDIF
674         END DO
675         !
676      ENDIF MODULO_NTRD
677
678      ! ======================================================================
679      ! IV. Write trends in the NetCDF file
680      ! ======================================================================
681
682      ! IV.1 Code for dimg mpp output
683      ! -----------------------------
684
685# if defined key_dimgout
686      STOP 'Not implemented'
687# else
688     
689      ! IV.2 Code for IOIPSL/NetCDF output
690      ! ----------------------------------
691
692      IF( lwp .AND. MOD( kt , ntrd_trc ) == 0 ) THEN
693         WRITE(numout,*) ' '
694         WRITE(numout,*) 'trd_mld_trc : write passive tracer trends in the NetCDF file :'
695         WRITE(numout,*) '~~~~~~~~~~~ '
696         WRITE(numout,*) '          ', trim(clhstnam), ' at kt = ', kt
697         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
698         WRITE(numout,*) ' '
699      ENDIF
700         
701      it = kt - nit000 + 1
702
703      NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
704         !
705
706         DO jn = 1, jptra
707            !
708            IF( luttrd(jn) ) THEN
709               !-- Specific treatment for EIV trends
710               !   WARNING : When eiv is switched on but key_diaeiv is not, we do NOT diagnose
711               !   u_eiv, v_eiv, and w_eiv : the exact eiv advective trends thus cannot be computed,
712               !   only their sum makes sense => mask directional contrib. to avoid confusion
713               z2d(:,:) = tmltrd_trc(:,:,jpmld_trc_xei,jn) + tmltrd_trc(:,:,jpmld_trc_yei,jn) &
714                    &   + tmltrd_trc(:,:,jpmld_trc_zei,jn)
715#if ( defined key_trcldf_eiv && defined key_diaeiv )
716               tmltrd_trc(:,:,jpmld_trc_xei,jn) = -999.
717               tmltrd_trc(:,:,jpmld_trc_yei,jn) = -999.
718               tmltrd_trc(:,:,jpmld_trc_zei,jn) = -999.
719#endif   
720               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
721               !-- Output the fields
722               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
723               CALL histwrite( nidtrd(jn), clvar         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
724               CALL histwrite( nidtrd(jn), clvar//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
725               CALL histwrite( nidtrd(jn), clvar//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
726           
727               DO jl = 1, jpltrd_trc - 2
728                  CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)),             &
729                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
730               END DO
731
732               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
733                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
734
735               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
736                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
737                     
738               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)),     &  ! now total EIV : jpltrd_trc + 1
739                    &          it, z2d(:,:), ndimtrd1, ndextrd1 )                     
740            !
741            ENDIF
742         END DO
743
744         IF( kt == nitend ) THEN
745            DO jn = 1, jptra
746               IF( luttrd(jn) )  CALL histclo( nidtrd(jn) )
747            END DO
748         ENDIF
749
750      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
751         
752                 
753         DO jn = 1, jptra
754            !
755            IF( luttrd(jn) ) THEN
756               !-- Specific treatment for EIV trends
757               !   WARNING : see above
758               z2d(:,:) = ztmltrd2(:,:,jpmld_trc_xei,jn) + ztmltrd2(:,:,jpmld_trc_yei,jn) &
759                   &   + ztmltrd2(:,:,jpmld_trc_zei,jn)
760
761#if ( defined key_trcldf_eiv && defined key_diaeiv )
762               ztmltrd2(:,:,jpmld_trc_xei,jn) = -999.
763               ztmltrd2(:,:,jpmld_trc_yei,jn) = -999.
764               ztmltrd2(:,:,jpmld_trc_zei,jn) = -999.
765#endif
766               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 
767               !-- Output the fields
768               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
769
770               CALL histwrite( nidtrd(jn), clvar         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
771               CALL histwrite( nidtrd(jn), clvar//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
772               CALL histwrite( nidtrd(jn), clvar//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
773
774               DO jl = 1, jpltrd_trc - 2
775                  CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)),           &
776                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
777               END DO
778           
779               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
780                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
781
782               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
783                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
784
785               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)),    &  ! now total EIV : jpltrd_trc + 1
786                 &          it, z2d(:,:), ndimtrd1, ndextrd1 )
787
788            ENDIF 
789            !
790         END DO
791         IF( kt == nitend ) THEN
792            DO jn = 1, jptra
793               IF( luttrd(jn) )  CALL histclo( nidtrd(jn) )
794            END DO
795         ENDIF
796
797         !
798      ENDIF NETCDF_OUTPUT
799         
800      ! Compute the control surface (for next time step) : flag = on
801      icount = 1
802
803# endif /* key_dimgout */
804
805      IF( MOD( kt, ntrd_trc ) == 0 ) THEN
806         !
807         ! Reset cumulative arrays to zero
808         ! -------------------------------         
809         nmoymltrd = 0
810         tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
811         tmlradm_trc        (:,:,:)   = 0.e0   ;   tml_sum_trc        (:,:,:)   = 0.e0
812         tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0
813         rmld_sum_trc       (:,:)     = 0.e0
814         !
815      ENDIF
816
817      ! ======================================================================
818      ! V. Write restart file
819      ! ======================================================================
820
821      IF( lrst_trc )   CALL trd_mld_trc_rst_write( kt )  ! this must be after the array swap above (III.3)
822
823   END SUBROUTINE trd_mld_trc
824
825
826   REAL FUNCTION sum2d( ztab )
827      !!----------------------------------------------------------------------
828      !! CD ??? prevoir d'utiliser plutot prtctl
829      !!----------------------------------------------------------------------
830      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab     
831      !!----------------------------------------------------------------------
832      sum2d = SUM(ztab(2:jpi-1,2:jpj-1))
833   END FUNCTION sum2d
834
835   SUBROUTINE trd_mld_trc_init
836      !!----------------------------------------------------------------------
837      !!                  ***  ROUTINE trd_mld_init  ***
838      !!
839      !! ** Purpose :   computation of vertically integrated T and S budgets
840      !!      from ocean surface down to control surface (NetCDF output)
841      !!
842      !! ** Method/usage :
843      !!
844      !!----------------------------------------------------------------------
845      INTEGER :: ilseq, jl, jn
846      REAL(wp) ::   zjulian, zsto, zout
847      CHARACTER (LEN=21) ::   clold ='OLD'         ! open specifier (direct access files)
848      CHARACTER (LEN=21) ::   clunf ='UNFORMATTED' ! open specifier (direct access files)
849      CHARACTER (LEN=21) ::   clseq ='SEQUENTIAL'  ! open specifier (direct access files)
850      CHARACTER (LEN=40) ::   clop, cleiv
851      CHARACTER (LEN=15) ::   csuff
852      CHARACTER (LEN=12) ::   clmxl
853      CHARACTER (LEN=16) ::   cltrcu
854      CHARACTER (LEN= 5) ::   clvar
855      CHARACTER (LEN=80) ::   clname
856
857      NAMELIST/namtoptrd/ ntrd_trc, nctls_trc, ucf_trc, &
858                          ln_trdmld_trc_restart, ln_trdmld_trc_instant, luttrd
859
860      !!----------------------------------------------------------------------
861
862      ! ======================================================================
863      ! I. initialization
864      ! ======================================================================
865
866      IF(lwp) THEN
867         WRITE(numout,*)
868         WRITE(numout,*) ' trd_mld_trc_init : Mixed-layer trends for passive tracers                '
869         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
870         WRITE(numout,*)
871      ENDIF
872
873     
874      ! I.1 Check consistency of user defined preferences
875      ! -------------------------------------------------
876#if defined key_trcldf_eiv
877      IF( lk_trdmld_trc .AND. ln_trcldf_iso ) THEN
878         WRITE(numout,cform_war)
879         WRITE(numout,*) '                You asked for ML diagnostics with iso-neutral diffusion   '
880         WRITE(numout,*) '                and eiv physics.                                          '
881         WRITE(numout,*) '                Yet, key_diaeiv is NOT switched on, so the eddy induced   '
882         WRITE(numout,*) '                velocity is not diagnosed.                                '
883         WRITE(numout,*) '                Therefore, we cannot deduce the eiv advective trends.     '
884         WRITE(numout,*) '                Only THE SUM of the i,j,k directional contributions then  '
885         WRITE(numout,*) '                makes sense => To avoid any confusion, we choosed to mask '
886         WRITE(numout,*) '                these i,j,k directional contributions (with -999.)        '
887         nwarn = nwarn + 1
888      ENDIF
889#  endif
890
891      IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, ntrd_trc ) /= 0 ) ) THEN
892         WRITE(numout,cform_err)
893         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
894         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
895         WRITE(numout,*) '                          you defined, ntrd_trc   = ', ntrd_trc
896         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
897         WRITE(numout,*) '                You should reconsider this choice.                        ' 
898         WRITE(numout,*) 
899         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
900         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
901         nstop = nstop + 1
902      ENDIF
903
904      IF( ( lk_trdmld_trc ) .AND. ( n_cla == 1 ) ) THEN
905         WRITE(numout,cform_war)
906         WRITE(numout,*) '                You set n_cla = 1. Note that the Mixed-Layer diagnostics  '
907         WRITE(numout,*) '                are not exact along the corresponding straits.            '
908         nwarn = nwarn + 1
909      ENDIF
910
911
912      ! * Debugging information *
913      IF( lldebug ) THEN
914         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
915         WRITE(numout,*) '               ln_trcadv_smolar = '     , ln_trcadv_smolar
916         WRITE(numout,*) '               ln_trdmld_trc_instant = ', ln_trdmld_trc_instant
917      ENDIF
918
919      IF( ln_trcadv_smolar .AND. .NOT. ln_trdmld_trc_instant ) THEN
920         WRITE(numout,cform_err)
921         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer Smolark. '
922         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
923         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
924         WRITE(numout,*) '                Smolarkiewicz scheme is Euler forward.                    '
925         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
926         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
927         WRITE(numout,*) 
928         nstop = nstop + 1
929      ENDIF
930
931      IF( ln_trcadv_muscl .AND. .NOT. ln_trdmld_trc_instant ) THEN
932         WRITE(numout,cform_err)
933         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
934         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
935         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
936         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
937         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
938         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
939         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
940         WRITE(numout,*) 
941         nstop = nstop + 1
942      ENDIF
943
944      IF( ln_trcadv_muscl2 .AND. .NOT. ln_trdmld_trc_instant ) THEN
945         WRITE(numout,cform_err)
946         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL2    '
947         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
948         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
949         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
950         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
951         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
952         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
953         WRITE(numout,*) 
954         nstop = nstop + 1
955      ENDIF
956
957      ! I.2 Initialize arrays to zero or read a restart file
958      ! ----------------------------------------------------
959      nmoymltrd   = 0
960
961      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
962      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
963      tmlradm_trc        (:,:,:)   = 0.e0
964
965      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
966      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
967
968      IF( lrsttr .AND. ln_trdmld_trc_restart ) THEN
969         CALL trd_mld_trc_rst_read
970      ELSE
971         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
972         tmlbn_trc          (:,:,:)   = 0.e0
973
974         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
975         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
976      ENDIF
977
978      ilseq  = 1   ;   icount = 1   ;   ionce  = 1  ! open specifier   
979
980      ! I.3 Read control surface from file ctlsurf_idx
981      ! ----------------------------------------------
982      IF( nctls_trc == 1 ) THEN
983         clname = 'ctlsurf_idx'
984         CALL ctlopn( numbol, clname, clold, clunf, clseq, ilseq, numout, lwp, 1 )
985         REWIND( numbol )
986         READ  ( numbol ) nbol_trc
987      ENDIF
988
989      ! ======================================================================
990      ! II. netCDF output initialization
991      ! ======================================================================
992
993#if defined key_dimgout 
994      ???
995#else
996      ! clmxl = legend root for netCDF output
997      IF( nctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
998         clmxl = 'Mixed Layer '
999      ELSE IF( nctls_trc == 1 ) THEN                              ! control surface = read index from file
1000         clmxl = '      Bowl '
1001      ELSE IF( nctls_trc >= 2 ) THEN                              ! control surface = model level
1002         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls_trc
1003      ENDIF
1004
1005      ! II.1 Define frequency of output and means
1006      ! -----------------------------------------
1007#  if defined key_diainstant
1008      IF( .NOT. ln_trdmld_trc_instant ) THEN
1009         STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...'
1010      ENDIF
1011      zsto = ntrd_trc * rdt
1012      clop ="inst(only(x))"
1013#  else
1014      IF( ln_trdmld_trc_instant ) THEN
1015         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
1016      ELSE
1017         zsto = ntrd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
1018      ENDIF
1019      clop ="ave(only(x))"
1020#  endif
1021      zout = ntrd_trc * rdt
1022
1023      IF(lwp) WRITE (numout,*) '                netCDF initialization'
1024
1025      ! II.2 Compute julian date from starting date of the run
1026      ! ------------------------------------------------------
1027      CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
1028      IF(lwp) WRITE(numout,*)' ' 
1029      IF(lwp) WRITE(numout,*)' Date 0 used :', nit000                  &
1030           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
1031           &   ,'Julian day : ', zjulian
1032
1033      ! II.3 Define the T grid trend file (nidtrd)
1034      ! ------------------------------------------
1035
1036      !-- Define long and short names for the NetCDF output variables
1037      !       ==> choose them according to trdmld_trc_oce.F90 <==
1038
1039#if defined key_diaeiv
1040      cleiv = " (*** only total EIV is meaningful ***)"           ! eiv advec. trends require u_eiv, v_eiv
1041#else
1042      cleiv = " "
1043#endif
1044      ctrd_trc(jpmld_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmld_trc_xad    ,2) = "_xad"
1045      ctrd_trc(jpmld_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmld_trc_yad    ,2) = "_yad"
1046      ctrd_trc(jpmld_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmld_trc_zad    ,2) = "_zad"
1047      ctrd_trc(jpmld_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmld_trc_ldf    ,2) = "_ldf"
1048      ctrd_trc(jpmld_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmld_trc_zdf    ,2) = "_zdf"
1049      ctrd_trc(jpmld_trc_xei    ,1) = " Zonal EIV advection"//cleiv      ;   ctrd_trc(jpmld_trc_xei    ,2) = "_xei"
1050      ctrd_trc(jpmld_trc_yei    ,1) = " Merid. EIV advection"//cleiv     ;   ctrd_trc(jpmld_trc_yei    ,2) = "_yei"
1051      ctrd_trc(jpmld_trc_zei    ,1) = " Vertical EIV advection"//cleiv   ;   ctrd_trc(jpmld_trc_zei    ,2) = "_zei"
1052      ctrd_trc(jpmld_trc_bbc    ,1) = " Geothermal flux"                 ;   ctrd_trc(jpmld_trc_bbc    ,2) = "_bbc"
1053      ctrd_trc(jpmld_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmld_trc_bbl    ,2) = "_bbl"
1054      ctrd_trc(jpmld_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmld_trc_dmp    ,2) = "_dmp"
1055      ctrd_trc(jpmld_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmld_trc_sbc    ,2) = "_sbc"
1056#if defined key_lobster
1057      ctrd_trc(jpmld_trc_sms_sed,1) = " Sources minus sinks : sed"       ;   ctrd_trc(jpmld_trc_sms_sed,2) = "_sms_sed"
1058      ctrd_trc(jpmld_trc_sms_bio,1) = " Sources minus sinks : bio"       ;   ctrd_trc(jpmld_trc_sms_bio,2) = "_sms_bio"
1059      ctrd_trc(jpmld_trc_sms_exp,1) = " Sources minus sinks : exp"       ;   ctrd_trc(jpmld_trc_sms_exp,2) = "_sms_exp"
1060#else
1061      ctrd_trc(jpmld_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmld_trc_sms    ,2) = "_sms"
1062#endif                                                                                                 
1063      ctrd_trc(jpmld_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radb   ,2) = "_radb"
1064      ctrd_trc(jpmld_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radn   ,2) = "_radn"
1065      ctrd_trc(jpmld_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmld_trc_atf    ,2) = "_atf"
1066      ctrd_trc(jpltrd_trc+1     ,1) = " Total EIV"//cleiv                ;   ctrd_trc(jpltrd_trc+1     ,2) = "_tei"
1067
1068      DO jn = 1, jptra     
1069      !-- Create a NetCDF file and enter the define mode
1070         IF( luttrd(jn) ) THEN
1071            csuff="TD_"//ctrcnm(jn)
1072            CALL dia_nam( clhstnam, ntrd_trc, csuff )
1073            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
1074               &        1, jpi, 1, jpj, 0, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )
1075     
1076            !-- Define the ML depth variable
1077            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
1078               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1079
1080         ENDIF
1081      END DO
1082
1083      !-- Define physical units
1084      IF( ucf_trc == 1. ) THEN
1085         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit
1086      ELSEIF ( ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
1087         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf
1088      ELSE
1089         cltrcu = "unknown?"
1090      ENDIF
1091
1092      !-- Define miscellaneous passive tracer mixed-layer variables
1093      IF( jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb ) THEN
1094         STOP 'Error : jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb'  ! see below
1095      ENDIF
1096#if defined key_lobster
1097      IF( lldebug ) THEN
1098         DO jn = 1, jptra
1099            WRITE(numout, *) 'TRC jpdet=', jpdet, ' jpnh4=', jpnh4
1100            WRITE(numout, *) 'TRC short title  ctrcnm  jn=", jn, " : ', ctrcnm(jn)
1101            WRITE(numout, *) 'TRC trim(ctrcnm(jn))//"_tot" = ', trim(ctrcnm(jn))//"ml_tot"  ! tml_tot -> detml_tot
1102         END DO
1103         CALL flush(numout)
1104      ENDIF
1105#else
1106!!      Error : this is not ready (PISCES)
1107#endif     
1108
1109      DO jn = 1, jptra
1110         !
1111         IF( luttrd(jn) ) THEN
1112            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
1113            CALL histdef(nidtrd(jn), clvar,           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
1114              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
1115            CALL histdef(nidtrd(jn), clvar//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
1116              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
1117            CALL histdef(nidtrd(jn), clvar//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
1118              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
1119         
1120            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmld_trc_atf
1121               CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
1122                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1123            END DO                                                                         ! if zsto=rdt above
1124         
1125            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_radb,1), & 
1126              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1127         
1128            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_atf,1),   & 
1129              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1130         
1131            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpltrd_trc+1,2)),  clmxl//" "//clvar//ctrd_trc(jpltrd_trc+1 ,1),   & 
1132              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! Total EIV
1133         !
1134         ENDIF
1135      END DO
1136
1137      !-- Leave IOIPSL/NetCDF define mode
1138      DO jn = 1, jptra
1139         IF( luttrd(jn) )  CALL histend( nidtrd(jn) )
1140      END DO
1141
1142#endif        /* key_dimgout */
1143   END SUBROUTINE trd_mld_trc_init
1144
1145#else
1146   !!----------------------------------------------------------------------
1147   !!   Default option :                                       Empty module
1148   !!----------------------------------------------------------------------
1149
1150CONTAINS
1151
1152   SUBROUTINE trd_mld_trc( kt )                                   ! Empty routine
1153      INTEGER, INTENT( in) ::   kt
1154      WRITE(*,*) 'trd_mld_trc: You should not have seen this print! error?', kt
1155   END SUBROUTINE trd_mld_trc
1156
1157   SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
1158      INTEGER, INTENT( in ) ::   ktrd, kjn                        ! ocean trend index and passive tracer rank
1159      CHARACTER(len=2), INTENT( in ) ::  ctype                    ! surface/bottom (2D) or interior (3D) physics
1160      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmld ! passive trc trend
1161      WRITE(*,*) 'trd_mld_trc_zint: You should not have seen this print! error?', ptrc_trdmld(1,1,1)
1162      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
1163      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
1164   END SUBROUTINE trd_mld_trc_zint
1165
1166   SUBROUTINE trd_mld_trc_init                                    ! Empty routine
1167      WRITE(*,*) 'trd_mld_trc_init: You should not have seen this print! error?'
1168   END SUBROUTINE trd_mld_trc_init
1169#endif
1170
1171   !!======================================================================
1172END MODULE trdmld_trc
Note: See TracBrowser for help on using the repository browser.