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

Last change on this file since 1204 was 1204, checked in by ctlod, 16 years ago

add subroutines for dummy module, see ticket: #259

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