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 @ 1685

Last change on this file since 1685 was 1685, checked in by smasson, 14 years ago

cleaning of logical units (use, flush and close), see ticket:570

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