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

Last change on this file since 1317 was 1317, checked in by smasson, 12 years ago

nwrite = modulo referenced to nit000 in all ouputs, see ticket:339

File size: 77.9 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, 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.( lrsttr ) ) 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      it = kt - nit000 + 1
578
579      MODULO_NTRD : IF( MOD( it, 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( it , 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      NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
831         !
832
833         DO jn = 1, jptra
834            !
835            IF( luttrd(jn) ) THEN
836               !-- Specific treatment for EIV trends
837               !   WARNING : When eiv is switched on but key_diaeiv is not, we do NOT diagnose
838               !   u_eiv, v_eiv, and w_eiv : the exact eiv advective trends thus cannot be computed,
839               !   only their sum makes sense => mask directional contrib. to avoid confusion
840               z2d(:,:) = tmltrd_trc(:,:,jpmld_trc_xei,jn) + tmltrd_trc(:,:,jpmld_trc_yei,jn) &
841                    &   + tmltrd_trc(:,:,jpmld_trc_zei,jn)
842#if ( defined key_trcldf_eiv && defined key_diaeiv )
843               tmltrd_trc(:,:,jpmld_trc_xei,jn) = -999.
844               tmltrd_trc(:,:,jpmld_trc_yei,jn) = -999.
845               tmltrd_trc(:,:,jpmld_trc_zei,jn) = -999.
846#endif   
847               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
848               !-- Output the fields
849               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
850               CALL histwrite( nidtrd(jn), clvar         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
851               CALL histwrite( nidtrd(jn), clvar//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
852               CALL histwrite( nidtrd(jn), clvar//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
853           
854               DO jl = 1, jpltrd_trc - 2
855                  CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)),             &
856                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
857               END DO
858
859               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
860                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
861
862               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
863                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
864                     
865               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)),     &  ! now total EIV : jpltrd_trc + 1
866                    &          it, z2d(:,:), ndimtrd1, ndextrd1 )                     
867            !
868            ENDIF
869         END DO
870
871         IF( kt == nitend ) THEN
872            DO jn = 1, jptra
873               IF( luttrd(jn) )  CALL histclo( nidtrd(jn) )
874            END DO
875         ENDIF
876
877      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
878         
879                 
880         DO jn = 1, jptra
881            !
882            IF( luttrd(jn) ) THEN
883               !-- Specific treatment for EIV trends
884               !   WARNING : see above
885               z2d(:,:) = ztmltrd2(:,:,jpmld_trc_xei,jn) + ztmltrd2(:,:,jpmld_trc_yei,jn) &
886                   &   + ztmltrd2(:,:,jpmld_trc_zei,jn)
887
888#if ( defined key_trcldf_eiv && defined key_diaeiv )
889               ztmltrd2(:,:,jpmld_trc_xei,jn) = -999.
890               ztmltrd2(:,:,jpmld_trc_yei,jn) = -999.
891               ztmltrd2(:,:,jpmld_trc_zei,jn) = -999.
892#endif
893               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 
894               !-- Output the fields
895               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
896
897               CALL histwrite( nidtrd(jn), clvar         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
898               CALL histwrite( nidtrd(jn), clvar//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
899               CALL histwrite( nidtrd(jn), clvar//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
900
901               DO jl = 1, jpltrd_trc - 2
902                  CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jl,2)),           &
903                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
904               END DO
905           
906               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
907                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
908
909               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
910                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
911
912               CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)),    &  ! now total EIV : jpltrd_trc + 1
913                 &          it, z2d(:,:), ndimtrd1, ndextrd1 )
914
915            ENDIF 
916            !
917         END DO
918         IF( kt == nitend ) THEN
919            DO jn = 1, jptra
920               IF( luttrd(jn) )  CALL histclo( nidtrd(jn) )
921            END DO
922         ENDIF
923
924         !
925      ENDIF NETCDF_OUTPUT
926         
927      ! Compute the control surface (for next time step) : flag = on
928      icount = 1
929
930# endif /* key_dimgout */
931
932      IF( MOD( it, ntrd_trc ) == 0 ) THEN
933         !
934         ! Reset cumulative arrays to zero
935         ! -------------------------------         
936         nmoymltrd = 0
937         tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
938         tmlradm_trc        (:,:,:)   = 0.e0   ;   tml_sum_trc        (:,:,:)   = 0.e0
939         tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0
940         rmld_sum_trc       (:,:)     = 0.e0
941         !
942      ENDIF
943
944      ! ======================================================================
945      ! V. Write restart file
946      ! ======================================================================
947
948      IF( lrst_trc )   CALL trd_mld_trc_rst_write( kt )  ! this must be after the array swap above (III.3)
949
950   END SUBROUTINE trd_mld_trc
951
952    SUBROUTINE trd_mld_bio( kt )
953      !!----------------------------------------------------------------------
954      !!                  ***  ROUTINE trd_mld  ***
955      !!
956      !! ** Purpose :  Compute and cumulate the mixed layer biological trends over an analysis
957      !!               period, and write NetCDF (or dimg) outputs.
958      !!
959      !! ** Method/usage :
960      !!          The stored trends can be chosen twofold (according to the ln_trdmld_trc_instant
961      !!          logical namelist variable) :
962      !!          1) to explain the difference between initial and final
963      !!             mixed-layer T & S (where initial and final relate to the
964      !!             current analysis window, defined by ntrd in the namelist)
965      !!          2) to explain the difference between the current and previous
966      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
967      !!             performed over each analysis window).
968      !!
969      !! ** Consistency check :
970      !!        If the control surface is fixed ( nctls > 1 ), the residual term (dh/dt
971      !!        entrainment) should be zero, at machine accuracy. Note that in the case
972      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
973      !!        over the first two analysis windows (except if restart).
974      !!        N.B. For ORCA2_LIM, use e.g. ntrd=5, ucf=1., nctls=8
975      !!             for checking residuals.
976      !!             On a NEC-SX5 computer, this typically leads to:
977      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmld_trc_instant=.false.
978      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmld_trc_instant=.true.
979      !!
980      !! ** Action :
981      !!       At each time step, mixed-layer averaged trends are stored in the
982      !!       tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
983      !!       This array is known when trd_mld is called, at the end of the stp subroutine,
984      !!       except for the purely vertical K_z diffusion term, which is embedded in the
985      !!       lateral diffusion trend.
986      !!
987      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
988      !!       from the lateral diffusion trend.
989      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
990      !!       arrays are updated.
991      !!       In III), called only once per analysis window, we compute the total trends,
992      !!       along with the residuals and the Asselin correction terms.
993      !!       In IV), the appropriate trends are written in the trends NetCDF file.
994      !!
995      !! References :
996      !!       - Vialard & al.
997      !!       - See NEMO documentation (in preparation)
998      !!----------------------------------------------------------------------
999      INTEGER, INTENT( in ) ::   kt                       ! ocean time-step index
1000#if defined key_lobster
1001      INTEGER  ::  jl, it
1002      LOGICAL  :: llwarn  = .TRUE., lldebug = .TRUE.
1003      REAL(wp), DIMENSION(jpi,jpj,jpdiabio) ::  ztmltrdbio2  ! only needed for mean diagnostics
1004      REAL(wp) :: zfn, zfn2
1005#if defined key_dimgout
1006      INTEGER ::  iyear,imon,iday
1007      CHARACTER(LEN=80) :: cltext, clmode
1008#endif
1009      !!----------------------------------------------------------------------
1010      ! ... Warnings
1011      IF( llwarn ) THEN
1012         IF(      ( nittrc000 /= nit000   ) &
1013              .OR.( ndttrc    /= 1        )    ) THEN
1014
1015            WRITE(numout,*) 'Be careful, trends diags never validated'
1016            STOP 'Uncomment this line to proceed'
1017         END IF
1018      END IF
1019
1020      ! ======================================================================
1021      ! II. Cumulate the trends over the analysis window
1022      ! ======================================================================
1023
1024      ztmltrdbio2(:,:,:) = 0.e0  ! <<< reset arrays to zero
1025
1026      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window
1027      ! ------------------------------------------------------------------------
1028      IF( kt == 2 ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
1029         !
1030         tmltrd_csum_ub_bio (:,:,:) = 0.e0
1031         !
1032      END IF
1033
1034      ! II.4 Cumulated trends over the analysis period
1035      ! ----------------------------------------------
1036      !
1037      !         [  1rst analysis window ] [     2nd analysis window     ]
1038      !
1039      !
1040      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
1041      !                            ntrd                             2*ntrd       etc.
1042      !     1      2     3     4    =5 e.g.                          =10
1043      !
1044      IF( ( kt >= 2 ).OR.( lrsttr ) ) THEN
1045         !
1046         nmoymltrdbio = nmoymltrdbio + 1
1047
1048         ! ... Trends associated with the time mean of the ML passive tracers
1049         tmltrd_sum_bio    (:,:,:) = tmltrd_sum_bio    (:,:,:) + tmltrd_bio    (:,:,:)
1050         tmltrd_csum_ln_bio(:,:,:) = tmltrd_csum_ln_bio(:,:,:) + tmltrd_sum_bio(:,:,:)
1051         !
1052      END IF
1053
1054      ! ======================================================================
1055      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
1056      ! ======================================================================
1057
1058      ! Convert to appropriate physical units
1059      tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * ucf_trc
1060
1061      MODULO_NTRD : IF( MOD( kt, ntrd_trc ) == 0 ) THEN      ! nitend MUST be multiple of ntrd
1062         !
1063         zfn  = float(nmoymltrdbio)    ;    zfn2 = zfn * zfn
1064
1065         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
1066         ! -------------------------------------------------------------
1067
1068#if defined key_diainstant
1069         STOP 'tmltrd_bio : key_diainstant was never checked within trdmld. Comment this to proceed.'
1070#endif
1071         ! III.2 Prepare fields for output ("mean" diagnostics)
1072         ! ----------------------------------------------------
1073
1074         ztmltrdbio2(:,:,:) = tmltrd_csum_ub_bio(:,:,:) + tmltrd_csum_ln_bio(:,:,:)
1075
1076         !-- Lateral boundary conditions
1077#if ! defined key_gyre
1078         ! ES_B27_CD_WARN : lbc inutile GYRE, cf. + haut
1079         DO jn = 1, jpdiabio
1080           CALL lbc_lnk( ztmltrdbio2(:,:,jn), 'T', 1. )
1081         ENDDO
1082#endif
1083         IF( lldebug ) THEN
1084            !
1085            WRITE(numout,*) 'trd_mld_bio : write trends in the Mixed Layer for debugging process:'
1086            WRITE(numout,*) '~~~~~~~~~~~  '
1087            WRITE(numout,*) 'TRC kt = ', kt, 'nmoymltrdbio = ', nmoymltrdbio
1088            WRITE(numout,*)
1089
1090            DO jl = 1, jpdiabio
1091              IF( ln_trdmld_trc_instant ) THEN
1092                  WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX  = ', jl, &
1093                     & ' SUM tmltrd_bio : ', SUM2D(tmltrd_bio(:,:,jl))
1094              ELSE
1095                  WRITE(numout,97) 'TRC jl =', jl, ' bio TREND INDEX  = ', jl, &
1096                     & ' SUM ztmltrdbio2 : ', SUM2D(ztmltrdbio2(:,:,jl))
1097              endif
1098            END DO
1099
110097          FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
110198          FORMAT(a10, i3, 2x, a30, 2x, g20.10)
110299          FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
1103            WRITE(numout,*)
1104            !
1105         ENDIF
1106
1107         ! III.3 Time evolution array swap
1108         ! -------------------------------
1109
1110         ! For passive tracer mean diagnostics
1111         tmltrd_csum_ub_bio (:,:,:) = zfn * tmltrd_sum_bio(:,:,:) - tmltrd_csum_ln_bio(:,:,:)
1112
1113         ! III.4 Convert to appropriate physical units
1114         ! -------------------------------------------
1115         ztmltrdbio2    (:,:,:) = ztmltrdbio2    (:,:,:) * ucf_trc/zfn2
1116
1117      END IF MODULO_NTRD
1118
1119      ! ======================================================================
1120      ! IV. Write trends in the NetCDF file
1121      ! ======================================================================
1122
1123      ! IV.1 Code for dimg mpp output
1124      ! -----------------------------
1125
1126# if defined key_dimgout
1127      STOP 'Not implemented'
1128# else
1129
1130      ! IV.2 Code for IOIPSL/NetCDF output
1131      ! ----------------------------------
1132
1133      ! define time axis
1134      it = kt - nit000 + 1
1135
1136      IF( lwp .AND. MOD( it , ntrd_trc ) == 0 ) THEN
1137         WRITE(numout,*) ' '
1138         WRITE(numout,*) 'trd_mld_bio : write ML bio trends in the NetCDF file :'
1139         WRITE(numout,*) '~~~~~~~~~~~ '
1140         WRITE(numout,*) '          ', TRIM(clhstnam), ' at kt = ', kt
1141         WRITE(numout,*) '          N.B. nmoymltrdbio = ', nmoymltrdbio
1142         WRITE(numout,*) ' '
1143      END IF
1144
1145
1146      ! 2. Start writing data
1147      ! ---------------------
1148
1149      NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN    ! <<< write the trends for passive tracer instant. diags
1150         !
1151            DO jl = 1, jpdiabio
1152               CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) ,            &
1153                    &          it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
1154            END DO
1155
1156
1157         IF( kt == nitend )   CALL histclo( nidtrdbio )
1158
1159      ELSE    ! <<< write the trends for passive tracer mean diagnostics
1160
1161            DO jl = 1, jpdiabio
1162               CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) ,            &
1163                    &          it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
1164            END DO
1165
1166            IF( kt == nitend )   CALL histclo( nidtrdbio )
1167            !
1168      END IF NETCDF_OUTPUT
1169
1170      ! Compute the control surface (for next time step) : flag = on
1171      icountbio = 1
1172
1173
1174# endif /* key_dimgout */
1175
1176      IF( MOD( it, ntrd_trc ) == 0 ) THEN
1177         !
1178         ! III.5 Reset cumulative arrays to zero
1179         ! -------------------------------------
1180         nmoymltrdbio = 0
1181         tmltrd_csum_ln_bio (:,:,:) = 0.e0
1182         tmltrd_sum_bio     (:,:,:) = 0.e0
1183      END IF
1184
1185      ! ======================================================================
1186      ! Write restart file
1187      ! ======================================================================
1188
1189! restart write is done in trd_mld_trc_write which is called by trd_mld_bio (Marina)
1190!
1191#endif
1192   END SUBROUTINE trd_mld_bio
1193
1194   REAL FUNCTION sum2d( ztab )
1195      !!----------------------------------------------------------------------
1196      !! CD ??? prevoir d'utiliser plutot prtctl
1197      !!----------------------------------------------------------------------
1198      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab     
1199      !!----------------------------------------------------------------------
1200      sum2d = SUM(ztab(2:jpi-1,2:jpj-1))
1201   END FUNCTION sum2d
1202
1203   SUBROUTINE trd_mld_trc_init
1204      !!----------------------------------------------------------------------
1205      !!                  ***  ROUTINE trd_mld_init  ***
1206      !!
1207      !! ** Purpose :   computation of vertically integrated T and S budgets
1208      !!      from ocean surface down to control surface (NetCDF output)
1209      !!
1210      !! ** Method/usage :
1211      !!
1212      !!----------------------------------------------------------------------
1213      INTEGER :: ilseq, jl, jn
1214      REAL(wp) ::   zjulian, zsto, zout
1215      CHARACTER (LEN=21) ::   clold ='OLD'         ! open specifier (direct access files)
1216      CHARACTER (LEN=21) ::   clunf ='UNFORMATTED' ! open specifier (direct access files)
1217      CHARACTER (LEN=21) ::   clseq ='SEQUENTIAL'  ! open specifier (direct access files)
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      CHARACTER (LEN=80) ::   clname
1224
1225      NAMELIST/namtoptrd/ ntrd_trc, nctls_trc, ucf_trc, &
1226                          ln_trdmld_trc_restart, ln_trdmld_trc_instant, luttrd
1227
1228      !!----------------------------------------------------------------------
1229
1230      ! ======================================================================
1231      ! I. initialization
1232      ! ======================================================================
1233
1234      IF(lwp) THEN
1235         WRITE(numout,*)
1236         WRITE(numout,*) ' trd_mld_trc_init : Mixed-layer trends for passive tracers                '
1237         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
1238         WRITE(numout,*)
1239      ENDIF
1240
1241     
1242      ! I.1 Check consistency of user defined preferences
1243      ! -------------------------------------------------
1244#if defined key_trcldf_eiv
1245      IF( lk_trdmld_trc .AND. ln_trcldf_iso ) THEN
1246         WRITE(numout,cform_war)
1247         WRITE(numout,*) '                You asked for ML diagnostics with iso-neutral diffusion   '
1248         WRITE(numout,*) '                and eiv physics.                                          '
1249         WRITE(numout,*) '                Yet, key_diaeiv is NOT switched on, so the eddy induced   '
1250         WRITE(numout,*) '                velocity is not diagnosed.                                '
1251         WRITE(numout,*) '                Therefore, we cannot deduce the eiv advective trends.     '
1252         WRITE(numout,*) '                Only THE SUM of the i,j,k directional contributions then  '
1253         WRITE(numout,*) '                makes sense => To avoid any confusion, we choosed to mask '
1254         WRITE(numout,*) '                these i,j,k directional contributions (with -999.)        '
1255         nwarn = nwarn + 1
1256      ENDIF
1257#  endif
1258
1259      IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, ntrd_trc ) /= 0 ) ) THEN
1260         WRITE(numout,cform_err)
1261         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
1262         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
1263         WRITE(numout,*) '                          you defined, ntrd_trc   = ', ntrd_trc
1264         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
1265         WRITE(numout,*) '                You should reconsider this choice.                        ' 
1266         WRITE(numout,*) 
1267         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
1268         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
1269         nstop = nstop + 1
1270      ENDIF
1271
1272      IF( ( lk_trdmld_trc ) .AND. ( n_cla == 1 ) ) THEN
1273         WRITE(numout,cform_war)
1274         WRITE(numout,*) '                You set n_cla = 1. Note that the Mixed-Layer diagnostics  '
1275         WRITE(numout,*) '                are not exact along the corresponding straits.            '
1276         nwarn = nwarn + 1
1277      ENDIF
1278
1279
1280      ! * Debugging information *
1281      IF( lldebug ) THEN
1282         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
1283         WRITE(numout,*) '               ln_trcadv_smolar = '     , ln_trcadv_smolar
1284         WRITE(numout,*) '               ln_trdmld_trc_instant = ', ln_trdmld_trc_instant
1285      ENDIF
1286
1287      IF( ln_trcadv_smolar .AND. .NOT. ln_trdmld_trc_instant ) THEN
1288         WRITE(numout,cform_err)
1289         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer Smolark. '
1290         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
1291         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
1292         WRITE(numout,*) '                Smolarkiewicz scheme is Euler forward.                    '
1293         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
1294         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
1295         WRITE(numout,*) 
1296         nstop = nstop + 1
1297      ENDIF
1298
1299      IF( ln_trcadv_muscl .AND. .NOT. ln_trdmld_trc_instant ) THEN
1300         WRITE(numout,cform_err)
1301         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
1302         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
1303         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
1304         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
1305         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
1306         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
1307         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
1308         WRITE(numout,*) 
1309         nstop = nstop + 1
1310      ENDIF
1311
1312      IF( ln_trcadv_muscl2 .AND. .NOT. ln_trdmld_trc_instant ) THEN
1313         WRITE(numout,cform_err)
1314         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL2    '
1315         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
1316         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
1317         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
1318         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
1319         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
1320         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
1321         WRITE(numout,*) 
1322         nstop = nstop + 1
1323      ENDIF
1324
1325      ! I.2 Initialize arrays to zero or read a restart file
1326      ! ----------------------------------------------------
1327      nmoymltrd   = 0
1328
1329      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
1330      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
1331      tmlradm_trc        (:,:,:)   = 0.e0
1332
1333      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
1334      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
1335
1336#if defined key_lobster
1337      nmoymltrdbio   = 0
1338      tmltrd_sum_bio     (:,:,:) = 0.e0     ;   tmltrd_csum_ln_bio (:,:,:) = 0.e0
1339      DO jl = 1, jp_lobster_trd
1340          ctrd_bio(jl,1) = ctrbil(jl)   ! long name
1341          ctrd_bio(jl,2) = ctrbio(jl)   ! short name
1342       ENDDO
1343#endif
1344
1345      IF( lrsttr .AND. ln_trdmld_trc_restart ) THEN
1346         CALL trd_mld_trc_rst_read
1347      ELSE
1348         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
1349         tmlbn_trc          (:,:,:)   = 0.e0
1350
1351         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
1352         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
1353#if defined key_lobster
1354         tmltrd_csum_ub_bio (:,:,:) = 0.e0
1355#endif
1356
1357       ENDIF
1358
1359      ilseq  = 1   ;   icount = 1   ;   ionce  = 1  ! open specifier   
1360
1361#if defined key_lobster
1362      icountbio = 1   ;   ioncebio  = 1  ! open specifier
1363#endif
1364
1365      ! I.3 Read control surface from file ctlsurf_idx
1366      ! ----------------------------------------------
1367      IF( nctls_trc == 1 ) THEN
1368         clname = 'ctlsurf_idx'
1369         CALL ctlopn( numbol, clname, clold, clunf, clseq, ilseq, numout, lwp, 1 )
1370         REWIND( numbol )
1371         READ  ( numbol ) nbol_trc
1372      ENDIF
1373
1374      ! ======================================================================
1375      ! II. netCDF output initialization
1376      ! ======================================================================
1377
1378#if defined key_dimgout 
1379      ???
1380#else
1381      ! clmxl = legend root for netCDF output
1382      IF( nctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
1383         clmxl = 'Mixed Layer '
1384      ELSE IF( nctls_trc == 1 ) THEN                              ! control surface = read index from file
1385         clmxl = '      Bowl '
1386      ELSE IF( nctls_trc >= 2 ) THEN                              ! control surface = model level
1387         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls_trc
1388      ENDIF
1389
1390      ! II.1 Define frequency of output and means
1391      ! -----------------------------------------
1392      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
1393      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
1394      ENDIF
1395#  if defined key_diainstant
1396      IF( .NOT. ln_trdmld_trc_instant ) THEN
1397         STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...'
1398      ENDIF
1399      zsto = ntrd_trc * rdt
1400      clop = "inst("//TRIM(clop)//")"
1401#  else
1402      IF( ln_trdmld_trc_instant ) THEN
1403         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
1404      ELSE
1405         zsto = ntrd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
1406      ENDIF
1407      clop = "ave("//TRIM(clop)//")"
1408#  endif
1409      zout = ntrd_trc * rdt
1410
1411      IF(lwp) WRITE (numout,*) '                netCDF initialization'
1412
1413      ! II.2 Compute julian date from starting date of the run
1414      ! ------------------------------------------------------
1415      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
1416      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
1417      IF(lwp) WRITE(numout,*)' ' 
1418      IF(lwp) WRITE(numout,*)' Date 0 used :', nit000                  &
1419           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
1420           &   ,'Julian day : ', zjulian
1421
1422      ! II.3 Define the T grid trend file (nidtrd)
1423      ! ------------------------------------------
1424
1425      !-- Define long and short names for the NetCDF output variables
1426      !       ==> choose them according to trdmld_trc_oce.F90 <==
1427
1428#if defined key_diaeiv
1429      cleiv = " (*** only total EIV is meaningful ***)"           ! eiv advec. trends require u_eiv, v_eiv
1430#else
1431      cleiv = " "
1432#endif
1433      ctrd_trc(jpmld_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmld_trc_xad    ,2) = "_xad"
1434      ctrd_trc(jpmld_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmld_trc_yad    ,2) = "_yad"
1435      ctrd_trc(jpmld_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmld_trc_zad    ,2) = "_zad"
1436      ctrd_trc(jpmld_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmld_trc_ldf    ,2) = "_ldf"
1437      ctrd_trc(jpmld_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmld_trc_zdf    ,2) = "_zdf"
1438      ctrd_trc(jpmld_trc_xei    ,1) = " Zonal EIV advection"//cleiv      ;   ctrd_trc(jpmld_trc_xei    ,2) = "_xei"
1439      ctrd_trc(jpmld_trc_yei    ,1) = " Merid. EIV advection"//cleiv     ;   ctrd_trc(jpmld_trc_yei    ,2) = "_yei"
1440      ctrd_trc(jpmld_trc_zei    ,1) = " Vertical EIV advection"//cleiv   ;   ctrd_trc(jpmld_trc_zei    ,2) = "_zei"
1441      ctrd_trc(jpmld_trc_bbc    ,1) = " Geothermal flux"                 ;   ctrd_trc(jpmld_trc_bbc    ,2) = "_bbc"
1442      ctrd_trc(jpmld_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmld_trc_bbl    ,2) = "_bbl"
1443      ctrd_trc(jpmld_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmld_trc_dmp    ,2) = "_dmp"
1444      ctrd_trc(jpmld_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmld_trc_sbc    ,2) = "_sbc"
1445      ctrd_trc(jpmld_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmld_trc_sms    ,2) = "_sms"
1446      ctrd_trc(jpmld_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radb   ,2) = "_radb"
1447      ctrd_trc(jpmld_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radn   ,2) = "_radn"
1448      ctrd_trc(jpmld_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmld_trc_atf    ,2) = "_atf"
1449      ctrd_trc(jpltrd_trc+1     ,1) = " Total EIV"//cleiv                ;   ctrd_trc(jpltrd_trc+1     ,2) = "_tei"
1450
1451      DO jn = 1, jptra     
1452      !-- Create a NetCDF file and enter the define mode
1453         IF( luttrd(jn) ) THEN
1454            csuff="ML_"//ctrcnm(jn)
1455            CALL dia_nam( clhstnam, ntrd_trc, csuff )
1456            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
1457               &        1, jpi, 1, jpj, 0, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )
1458     
1459            !-- Define the ML depth variable
1460            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
1461               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1462
1463         ENDIF
1464      END DO
1465
1466#if defined key_lobster
1467          !-- Create a NetCDF file and enter the define mode
1468          CALL dia_nam( clhstnam, ntrd_trc, 'trdbio' )
1469          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
1470             &             1, jpi, 1, jpj, 0, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom )
1471#endif
1472
1473      !-- Define physical units
1474      IF( ucf_trc == 1. ) THEN
1475         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit
1476      ELSEIF ( ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
1477         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf
1478      ELSE
1479         cltrcu = "unknown?"
1480      ENDIF
1481
1482      !-- Define miscellaneous passive tracer mixed-layer variables
1483      IF( jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb ) THEN
1484         STOP 'Error : jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb'  ! see below
1485      ENDIF
1486
1487      DO jn = 1, jptra
1488         !
1489         IF( luttrd(jn) ) THEN
1490            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
1491            CALL histdef(nidtrd(jn), clvar,           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
1492              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
1493            CALL histdef(nidtrd(jn), clvar//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
1494              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
1495            CALL histdef(nidtrd(jn), clvar//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
1496              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
1497         
1498            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmld_trc_atf
1499               CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
1500                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1501            END DO                                                                         ! if zsto=rdt above
1502         
1503            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_radb,1), & 
1504              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1505         
1506            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_atf,1),   & 
1507              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1508         
1509            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpltrd_trc+1,2)),  clmxl//" "//clvar//ctrd_trc(jpltrd_trc+1 ,1),   & 
1510              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! Total EIV
1511         !
1512         ENDIF
1513      END DO
1514
1515#if defined key_lobster
1516      DO jl = 1, jp_lobster_trd
1517         CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1))   ,            &
1518             &    cltrcu, jpi, jpj, nh_tb, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1519      END DO                                                                         ! if zsto=rdt above
1520#endif
1521
1522      !-- Leave IOIPSL/NetCDF define mode
1523      DO jn = 1, jptra
1524         IF( luttrd(jn) )  CALL histend( nidtrd(jn) )
1525      END DO
1526
1527#if defined key_lobster
1528      !-- Leave IOIPSL/NetCDF define mode
1529      CALL histend( nidtrdbio )
1530
1531      IF(lwp) WRITE(numout,*)
1532       IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
1533#endif
1534
1535#endif        /* key_dimgout */
1536   END SUBROUTINE trd_mld_trc_init
1537
1538#else
1539   !!----------------------------------------------------------------------
1540   !!   Default option :                                       Empty module
1541   !!----------------------------------------------------------------------
1542
1543   INTERFACE trd_mod_trc
1544      MODULE PROCEDURE trd_mod_trc_trp, trd_mod_trc_bio
1545   END INTERFACE
1546
1547CONTAINS
1548
1549   SUBROUTINE trd_mld_trc( kt )                                   ! Empty routine
1550      INTEGER, INTENT( in) ::   kt
1551      WRITE(*,*) 'trd_mld_trc: You should not have seen this print! error?', kt
1552   END SUBROUTINE trd_mld_trc
1553
1554   SUBROUTINE trd_mld_bio( kt )
1555      INTEGER, INTENT( in) ::   kt
1556      WRITE(*,*) 'trd_mld_bio: You should not have seen this print! error?', kt
1557   END SUBROUTINE trd_mld_bio
1558
1559   SUBROUTINE trd_mod_trc_bio( ptrbio, ktrd, kt )
1560      INTEGER               , INTENT( in )     ::   kt      ! time step
1561      INTEGER               , INTENT( in )     ::   ktrd    ! bio trend index
1562      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrbio  ! Bio trend
1563      WRITE(*,*) 'trd_mod_trc_bio : You should not have seen this print! error?', ptrbio(1,1,1)
1564      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
1565      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
1566   END SUBROUTINE trd_mod_trc_bio
1567
1568   SUBROUTINE trd_mod_trc_trp( ptrtrd, kjn, ktrd, kt )
1569      INTEGER               , INTENT( in )     ::   kt      ! time step
1570      INTEGER               , INTENT( in )     ::   kjn     ! tracer index
1571      INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
1572      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
1573      WRITE(*,*) 'trd_mod_trc_trp : You should not have seen this print! error?', ptrtrd(1,1,1)
1574      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
1575      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
1576      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
1577   END SUBROUTINE trd_mod_trc_trp
1578
1579   SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
1580      INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
1581      CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
1582      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmld            ! passive trc trend
1583      WRITE(*,*) 'trd_mld_trc_zint: You should not have seen this print! error?', ptrc_trdmld(1,1,1)
1584      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
1585      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
1586      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
1587   END SUBROUTINE trd_mld_trc_zint
1588
1589   SUBROUTINE trd_mld_trc_init                                    ! Empty routine
1590      WRITE(*,*) 'trd_mld_trc_init: You should not have seen this print! error?'
1591   END SUBROUTINE trd_mld_trc_init
1592#endif
1593
1594   !!======================================================================
1595END MODULE trdmld_trc
Note: See TracBrowser for help on using the repository browser.