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

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

add a namelist logical to mask land points in NetCDF outputs, see ticket:322

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      MODULO_NTRD : IF( MOD( kt, ntrd_trc ) == 0 ) THEN           ! nitend MUST be multiple of ntrd_trc
578         !
579         ztmltot (:,:,:) = 0.e0                                   ! reset arrays to zero
580         ztmlres (:,:,:) = 0.e0
581         ztmltot2(:,:,:) = 0.e0
582         ztmlres2(:,:,:) = 0.e0
583     
584         zfn  = FLOAT( nmoymltrd )    ;    zfn2 = zfn * zfn
585         
586         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
587         ! -------------------------------------------------------------
588
589         DO jn = 1, jptra
590            IF( luttrd(jn) ) THEN
591               !-- Compute total trends    (use rdttrc instead of rdt ???)
592               IF ( ln_trcadv_smolar .OR. ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN  ! EULER-FORWARD schemes
593                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
594               ELSE                                                                     ! LEAP-FROG schemes
595                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
596               ENDIF
597               
598               !-- Compute residuals
599               ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
600                  &                                                 - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)  )
601               
602               !-- Diagnose Asselin trend over the analysis window
603               ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
604               ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
605               
606         !-- Lateral boundary conditions
607#if ! defined key_gyre
608
609               CALL lbc_lnk( ztmltot(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlres(:,:,jn) , 'T', 1. )
610               CALL lbc_lnk( ztmlatf(:,:,jn) , 'T', 1. )   ;   CALL lbc_lnk( ztmlrad(:,:,jn) , 'T', 1. )
611
612#endif
613
614#if defined key_diainstant
615               STOP 'tmltrd_trc : key_diainstant was never checked within trdmld. Comment this to proceed.'
616#endif
617            ENDIF
618         END DO
619
620         ! III.2 Prepare fields for output ("mean" diagnostics)
621         ! ----------------------------------------------------
622         
623         !-- Update the ML depth time sum (to build the Leap-Frog time mean)
624         rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
625
626               !-- Compute passive tracer total trends
627         DO jn = 1, jptra
628            IF( luttrd(jn) ) THEN
629               tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
630               ztmltot2   (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) /  ( 2.*rdt )    ! now tracer unit is /sec
631            ENDIF
632         END DO
633
634         !-- Compute passive tracer residuals
635         DO jn = 1, jptra
636            IF( luttrd(jn) ) THEN
637               !
638               DO jl = 1, jpltrd_trc
639                  ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
640               END DO
641               
642               ztmltrdm2(:,:,jn) = 0.e0
643               DO jl = 1, jpltrd_trc
644                  ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
645               END DO
646               
647               ztmlres2(:,:,jn) =  ztmltot2(:,:,jn)  -       &
648                  & ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
649                  &                     - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
650               !
651
652               !-- Diagnose Asselin trend over the analysis window
653               ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn) &
654                  &                                               + tmltrd_atf_sumb_trc(:,:,jn)
655               ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmld_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmld_trc_radb,jn) &
656                  &                                               + tmltrd_rad_sumb_trc(:,:,jn)
657
658         !-- Lateral boundary conditions
659#if ! defined key_gyre         
660               CALL lbc_lnk( ztmltot2(:,:,jn), 'T', 1. )
661               CALL lbc_lnk( ztmlres2(:,:,jn), 'T', 1. )
662               DO jl = 1, jpltrd_trc
663                  CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
664               END DO
665#endif
666            ENDIF
667         END DO
668
669         ! * Debugging information *
670         IF( lldebug ) THEN
671            !
672            WRITE(numout,*) 'trd_mld_trc : write trends in the Mixed Layer for debugging process:'
673            WRITE(numout,*) '~~~~~~~~~~~  '
674            WRITE(numout,*)
675            WRITE(numout,*) 'TRC kt = ', kt, '    nmoymltrd = ', nmoymltrd
676
677            DO jn = 1, jptra
678
679               IF( luttrd(jn) ) THEN
680                  WRITE(numout, *)
681                  WRITE(numout, *) '>>>>>>>>>>>>>>>>>>  TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
682                 
683                  WRITE(numout, *)
684                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres     : ', SUM2D(ztmlres(:,:,jn))
685                  !CD??? PREVOIR: z2d = ztmlres(:,:,jn)   ;   CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres   -   : ')
686                 
687                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
688                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
689                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot    : ', SUM2D(+ztmltot    (:,:,jn))
690                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
691                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
692                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
693                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
694                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
695                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
696                 
697                  WRITE(numout, *)
698                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2    : ', SUM2D(ztmlres2(:,:,jn))
699                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
700                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
701                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2            : ', SUM2D(+ztmltot2(:,:,jn))
702                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2           : ', SUM2D(+ztmltrdm2(:,:,jn)) 
703                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_atf,jn)) 
704                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
705                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmld_trc_radb,jn))
706                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
707                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
708                 
709                  WRITE(numout, *)
710                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot     : ', SUM2D(ztmltot    (:,:,jn))
711                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
712                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc    : ', SUM2D(tml_trc    (:,:,jn))
713                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc  : ', SUM2D(tmlbn_trc  (:,:,jn))
714                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc   : ', SUM2D(tmlb_trc   (:,:,jn))
715                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc  : ', SUM2D(tmlbb_trc  (:,:,jn))
716                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
717                 
718                  WRITE(numout, *)
719                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
720                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
721                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
722                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
723                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
724                 
725                  WRITE(numout, *)
726                  DO jl = 1, jpltrd_trc
727                     WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmld_trc_xxx = ', jl, &
728                        & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
729                  END DO
730               
731                  WRITE(numout,*) 
732                  WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** '
733                  WRITE(numout,*)
734                  WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
735                  WRITE(numout,*)
736                  WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
737                 
738                  WRITE(numout,*) 
739                  WRITE(numout,*) ' *********************** ZTMLRES *********************** '
740                  WRITE(numout,*)
741                 
742                  WRITE(numout,*) '...................................................'
743                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
744                  DO jj = 5, 1, -1
745                     WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
746                  END DO
747                 
748                  WRITE(numout,*) 
749                  WRITE(numout,*) ' *********************** ZTMLRES2 *********************** '
750                  WRITE(numout,*)
751                 
752                  WRITE(numout,*) '...................................................'
753                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
754                  DO jj = 5, 1, -1
755                     WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
756                  END DO
757                  !
758               ENDIF
759               !
760            END DO
761
762
76397            FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
76498            FORMAT(a10, i3, 2x, a30, 2x, g20.10)
76599            FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
766              WRITE(numout,*)
767            !
768         ENDIF
769
770         ! III.3 Time evolution array swap
771         ! -------------------------------
772         ! ML depth
773         rmldbn_trc(:,:)   = rmld_trc(:,:)
774         rmld_sum_trc(:,:)     = rmld_sum_trc(:,:)     /      (2*zfn)  ! similar to tml_sum and sml_sum
775         DO jn = 1, jptra
776            IF( luttrd(jn) ) THEN       
777               ! For passive tracer instantaneous diagnostics
778               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
779               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
780               
781               ! For passive tracer mean diagnostics
782               tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
783               tml_sumb_trc       (:,:,jn)   = tml_sum_trc(:,:,jn)
784               tmltrd_atf_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmld_trc_atf ,jn)
785               tmltrd_rad_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmld_trc_radb,jn)
786               
787               
788               ! III.4 Convert to appropriate physical units
789               ! -------------------------------------------
790               ztmltot     (:,:,jn)   = ztmltot     (:,:,jn)   * ucf_trc/zfn   ! instant diags
791               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * ucf_trc/zfn
792               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * ucf_trc/zfn
793               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * ucf_trc/zfn
794               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
795               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * ucf_trc/zfn2
796               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * ucf_trc/zfn2
797               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * ucf_trc/zfn2
798               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * ucf_trc/zfn2
799               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * ucf_trc/zfn2
800            ENDIF
801         END DO
802         !
803      ENDIF MODULO_NTRD
804
805      ! ======================================================================
806      ! IV. Write trends in the NetCDF file
807      ! ======================================================================
808
809      ! IV.1 Code for dimg mpp output
810      ! -----------------------------
811
812# if defined key_dimgout
813      STOP 'Not implemented'
814# else
815     
816      ! IV.2 Code for IOIPSL/NetCDF output
817      ! ----------------------------------
818
819      IF( lwp .AND. MOD( kt , ntrd_trc ) == 0 ) THEN
820         WRITE(numout,*) ' '
821         WRITE(numout,*) 'trd_mld_trc : write passive tracer trends in the NetCDF file :'
822         WRITE(numout,*) '~~~~~~~~~~~ '
823         WRITE(numout,*) '          ', trim(clhstnam), ' at kt = ', kt
824         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
825         WRITE(numout,*) ' '
826      ENDIF
827         
828      it = kt - nit000 + 1
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( kt, 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      IF( lwp .AND. MOD( kt , ntrd_trc ) == 0 ) THEN
1134         WRITE(numout,*) ' '
1135         WRITE(numout,*) 'trd_mld_bio : write ML bio trends in the NetCDF file :'
1136         WRITE(numout,*) '~~~~~~~~~~~ '
1137         WRITE(numout,*) '          ', TRIM(clhstnam), ' at kt = ', kt
1138         WRITE(numout,*) '          N.B. nmoymltrdbio = ', nmoymltrdbio
1139         WRITE(numout,*) ' '
1140      END IF
1141
1142
1143      ! define time axis
1144      it = kt - nit000 + 1
1145
1146
1147      ! 2. Start writing data
1148      ! ---------------------
1149
1150      NETCDF_OUTPUT : IF( ln_trdmld_trc_instant ) THEN    ! <<< write the trends for passive tracer instant. diags
1151         !
1152            DO jl = 1, jpdiabio
1153               CALL histwrite( nidtrdbio,TRIM("ML_"//ctrd_bio(jl,2)) ,            &
1154                    &          it, tmltrd_bio(:,:,jl), ndimtrd1, ndextrd1 )
1155            END DO
1156
1157
1158         IF( kt == nitend )   CALL histclo( nidtrdbio )
1159
1160      ELSE    ! <<< write the trends for passive tracer mean diagnostics
1161
1162            DO jl = 1, jpdiabio
1163               CALL histwrite( nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)) ,            &
1164                    &          it, ztmltrdbio2(:,:,jl), ndimtrd1, ndextrd1 )
1165            END DO
1166
1167            IF( kt == nitend )   CALL histclo( nidtrdbio )
1168            !
1169      END IF NETCDF_OUTPUT
1170
1171      ! Compute the control surface (for next time step) : flag = on
1172      icountbio = 1
1173
1174
1175# endif /* key_dimgout */
1176
1177      IF( MOD( kt, ntrd_trc ) == 0 ) THEN
1178         !
1179         ! III.5 Reset cumulative arrays to zero
1180         ! -------------------------------------
1181         nmoymltrdbio = 0
1182         tmltrd_csum_ln_bio (:,:,:) = 0.e0
1183         tmltrd_sum_bio     (:,:,:) = 0.e0
1184      END IF
1185
1186      ! ======================================================================
1187      ! Write restart file
1188      ! ======================================================================
1189
1190! restart write is done in trd_mld_trc_write which is called by trd_mld_bio (Marina)
1191!
1192#endif
1193   END SUBROUTINE trd_mld_bio
1194
1195   REAL FUNCTION sum2d( ztab )
1196      !!----------------------------------------------------------------------
1197      !! CD ??? prevoir d'utiliser plutot prtctl
1198      !!----------------------------------------------------------------------
1199      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab     
1200      !!----------------------------------------------------------------------
1201      sum2d = SUM(ztab(2:jpi-1,2:jpj-1))
1202   END FUNCTION sum2d
1203
1204   SUBROUTINE trd_mld_trc_init
1205      !!----------------------------------------------------------------------
1206      !!                  ***  ROUTINE trd_mld_init  ***
1207      !!
1208      !! ** Purpose :   computation of vertically integrated T and S budgets
1209      !!      from ocean surface down to control surface (NetCDF output)
1210      !!
1211      !! ** Method/usage :
1212      !!
1213      !!----------------------------------------------------------------------
1214      INTEGER :: ilseq, jl, jn
1215      REAL(wp) ::   zjulian, zsto, zout
1216      CHARACTER (LEN=21) ::   clold ='OLD'         ! open specifier (direct access files)
1217      CHARACTER (LEN=21) ::   clunf ='UNFORMATTED' ! open specifier (direct access files)
1218      CHARACTER (LEN=21) ::   clseq ='SEQUENTIAL'  ! open specifier (direct access files)
1219      CHARACTER (LEN=40) ::   clop, cleiv
1220      CHARACTER (LEN=15) ::   csuff
1221      CHARACTER (LEN=12) ::   clmxl
1222      CHARACTER (LEN=16) ::   cltrcu
1223      CHARACTER (LEN= 5) ::   clvar
1224      CHARACTER (LEN=80) ::   clname
1225
1226      NAMELIST/namtoptrd/ ntrd_trc, nctls_trc, ucf_trc, &
1227                          ln_trdmld_trc_restart, ln_trdmld_trc_instant, luttrd
1228
1229      !!----------------------------------------------------------------------
1230
1231      ! ======================================================================
1232      ! I. initialization
1233      ! ======================================================================
1234
1235      IF(lwp) THEN
1236         WRITE(numout,*)
1237         WRITE(numout,*) ' trd_mld_trc_init : Mixed-layer trends for passive tracers                '
1238         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
1239         WRITE(numout,*)
1240      ENDIF
1241
1242     
1243      ! I.1 Check consistency of user defined preferences
1244      ! -------------------------------------------------
1245#if defined key_trcldf_eiv
1246      IF( lk_trdmld_trc .AND. ln_trcldf_iso ) THEN
1247         WRITE(numout,cform_war)
1248         WRITE(numout,*) '                You asked for ML diagnostics with iso-neutral diffusion   '
1249         WRITE(numout,*) '                and eiv physics.                                          '
1250         WRITE(numout,*) '                Yet, key_diaeiv is NOT switched on, so the eddy induced   '
1251         WRITE(numout,*) '                velocity is not diagnosed.                                '
1252         WRITE(numout,*) '                Therefore, we cannot deduce the eiv advective trends.     '
1253         WRITE(numout,*) '                Only THE SUM of the i,j,k directional contributions then  '
1254         WRITE(numout,*) '                makes sense => To avoid any confusion, we choosed to mask '
1255         WRITE(numout,*) '                these i,j,k directional contributions (with -999.)        '
1256         nwarn = nwarn + 1
1257      ENDIF
1258#  endif
1259
1260      IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, ntrd_trc ) /= 0 ) ) THEN
1261         WRITE(numout,cform_err)
1262         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
1263         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
1264         WRITE(numout,*) '                          you defined, ntrd_trc   = ', ntrd_trc
1265         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
1266         WRITE(numout,*) '                You should reconsider this choice.                        ' 
1267         WRITE(numout,*) 
1268         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
1269         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
1270         nstop = nstop + 1
1271      ENDIF
1272
1273      IF( ( lk_trdmld_trc ) .AND. ( n_cla == 1 ) ) THEN
1274         WRITE(numout,cform_war)
1275         WRITE(numout,*) '                You set n_cla = 1. Note that the Mixed-Layer diagnostics  '
1276         WRITE(numout,*) '                are not exact along the corresponding straits.            '
1277         nwarn = nwarn + 1
1278      ENDIF
1279
1280
1281      ! * Debugging information *
1282      IF( lldebug ) THEN
1283         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
1284         WRITE(numout,*) '               ln_trcadv_smolar = '     , ln_trcadv_smolar
1285         WRITE(numout,*) '               ln_trdmld_trc_instant = ', ln_trdmld_trc_instant
1286      ENDIF
1287
1288      IF( ln_trcadv_smolar .AND. .NOT. ln_trdmld_trc_instant ) THEN
1289         WRITE(numout,cform_err)
1290         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer Smolark. '
1291         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
1292         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
1293         WRITE(numout,*) '                Smolarkiewicz scheme is Euler forward.                    '
1294         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
1295         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
1296         WRITE(numout,*) 
1297         nstop = nstop + 1
1298      ENDIF
1299
1300      IF( ln_trcadv_muscl .AND. .NOT. ln_trdmld_trc_instant ) THEN
1301         WRITE(numout,cform_err)
1302         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
1303         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
1304         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
1305         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
1306         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
1307         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
1308         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
1309         WRITE(numout,*) 
1310         nstop = nstop + 1
1311      ENDIF
1312
1313      IF( ln_trcadv_muscl2 .AND. .NOT. ln_trdmld_trc_instant ) THEN
1314         WRITE(numout,cform_err)
1315         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL2    '
1316         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
1317         WRITE(numout,*) '                WHY? Everything in trdmld_trc is coded for leap-frog, and '
1318         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
1319         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
1320         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
1321         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
1322         WRITE(numout,*) 
1323         nstop = nstop + 1
1324      ENDIF
1325
1326      ! I.2 Initialize arrays to zero or read a restart file
1327      ! ----------------------------------------------------
1328      nmoymltrd   = 0
1329
1330      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
1331      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
1332      tmlradm_trc        (:,:,:)   = 0.e0
1333
1334      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
1335      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
1336
1337#if defined key_lobster
1338      nmoymltrdbio   = 0
1339      tmltrd_sum_bio     (:,:,:) = 0.e0     ;   tmltrd_csum_ln_bio (:,:,:) = 0.e0
1340      DO jl = 1, jp_lobster_trd
1341          ctrd_bio(jl,1) = ctrbil(jl)   ! long name
1342          ctrd_bio(jl,2) = ctrbio(jl)   ! short name
1343       ENDDO
1344#endif
1345
1346      IF( lrsttr .AND. ln_trdmld_trc_restart ) THEN
1347         CALL trd_mld_trc_rst_read
1348      ELSE
1349         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
1350         tmlbn_trc          (:,:,:)   = 0.e0
1351
1352         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
1353         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
1354#if defined key_lobster
1355         tmltrd_csum_ub_bio (:,:,:) = 0.e0
1356#endif
1357
1358       ENDIF
1359
1360      ilseq  = 1   ;   icount = 1   ;   ionce  = 1  ! open specifier   
1361
1362#if defined key_lobster
1363      icountbio = 1   ;   ioncebio  = 1  ! open specifier
1364#endif
1365
1366      ! I.3 Read control surface from file ctlsurf_idx
1367      ! ----------------------------------------------
1368      IF( nctls_trc == 1 ) THEN
1369         clname = 'ctlsurf_idx'
1370         CALL ctlopn( numbol, clname, clold, clunf, clseq, ilseq, numout, lwp, 1 )
1371         REWIND( numbol )
1372         READ  ( numbol ) nbol_trc
1373      ENDIF
1374
1375      ! ======================================================================
1376      ! II. netCDF output initialization
1377      ! ======================================================================
1378
1379#if defined key_dimgout 
1380      ???
1381#else
1382      ! clmxl = legend root for netCDF output
1383      IF( nctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
1384         clmxl = 'Mixed Layer '
1385      ELSE IF( nctls_trc == 1 ) THEN                              ! control surface = read index from file
1386         clmxl = '      Bowl '
1387      ELSE IF( nctls_trc >= 2 ) THEN                              ! control surface = model level
1388         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls_trc
1389      ENDIF
1390
1391      ! II.1 Define frequency of output and means
1392      ! -----------------------------------------
1393      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
1394      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
1395      ENDIF
1396#  if defined key_diainstant
1397      IF( .NOT. ln_trdmld_trc_instant ) THEN
1398         STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...'
1399      ENDIF
1400      zsto = ntrd_trc * rdt
1401      clop = "inst("//TRIM(clop)//")"
1402#  else
1403      IF( ln_trdmld_trc_instant ) THEN
1404         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
1405      ELSE
1406         zsto = ntrd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
1407      ENDIF
1408      clop = "ave("//TRIM(clop)//")"
1409#  endif
1410      zout = ntrd_trc * rdt
1411
1412      IF(lwp) WRITE (numout,*) '                netCDF initialization'
1413
1414      ! II.2 Compute julian date from starting date of the run
1415      ! ------------------------------------------------------
1416      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
1417      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
1418      IF(lwp) WRITE(numout,*)' ' 
1419      IF(lwp) WRITE(numout,*)' Date 0 used :', nit000                  &
1420           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
1421           &   ,'Julian day : ', zjulian
1422
1423      ! II.3 Define the T grid trend file (nidtrd)
1424      ! ------------------------------------------
1425
1426      !-- Define long and short names for the NetCDF output variables
1427      !       ==> choose them according to trdmld_trc_oce.F90 <==
1428
1429#if defined key_diaeiv
1430      cleiv = " (*** only total EIV is meaningful ***)"           ! eiv advec. trends require u_eiv, v_eiv
1431#else
1432      cleiv = " "
1433#endif
1434      ctrd_trc(jpmld_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmld_trc_xad    ,2) = "_xad"
1435      ctrd_trc(jpmld_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmld_trc_yad    ,2) = "_yad"
1436      ctrd_trc(jpmld_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmld_trc_zad    ,2) = "_zad"
1437      ctrd_trc(jpmld_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmld_trc_ldf    ,2) = "_ldf"
1438      ctrd_trc(jpmld_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmld_trc_zdf    ,2) = "_zdf"
1439      ctrd_trc(jpmld_trc_xei    ,1) = " Zonal EIV advection"//cleiv      ;   ctrd_trc(jpmld_trc_xei    ,2) = "_xei"
1440      ctrd_trc(jpmld_trc_yei    ,1) = " Merid. EIV advection"//cleiv     ;   ctrd_trc(jpmld_trc_yei    ,2) = "_yei"
1441      ctrd_trc(jpmld_trc_zei    ,1) = " Vertical EIV advection"//cleiv   ;   ctrd_trc(jpmld_trc_zei    ,2) = "_zei"
1442      ctrd_trc(jpmld_trc_bbc    ,1) = " Geothermal flux"                 ;   ctrd_trc(jpmld_trc_bbc    ,2) = "_bbc"
1443      ctrd_trc(jpmld_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmld_trc_bbl    ,2) = "_bbl"
1444      ctrd_trc(jpmld_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmld_trc_dmp    ,2) = "_dmp"
1445      ctrd_trc(jpmld_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmld_trc_sbc    ,2) = "_sbc"
1446      ctrd_trc(jpmld_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmld_trc_sms    ,2) = "_sms"
1447      ctrd_trc(jpmld_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radb   ,2) = "_radb"
1448      ctrd_trc(jpmld_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmld_trc_radn   ,2) = "_radn"
1449      ctrd_trc(jpmld_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmld_trc_atf    ,2) = "_atf"
1450      ctrd_trc(jpltrd_trc+1     ,1) = " Total EIV"//cleiv                ;   ctrd_trc(jpltrd_trc+1     ,2) = "_tei"
1451
1452      DO jn = 1, jptra     
1453      !-- Create a NetCDF file and enter the define mode
1454         IF( luttrd(jn) ) THEN
1455            csuff="ML_"//ctrcnm(jn)
1456            CALL dia_nam( clhstnam, ntrd_trc, csuff )
1457            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
1458               &        1, jpi, 1, jpj, 0, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )
1459     
1460            !-- Define the ML depth variable
1461            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
1462               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1463
1464         ENDIF
1465      END DO
1466
1467#if defined key_lobster
1468          !-- Create a NetCDF file and enter the define mode
1469          CALL dia_nam( clhstnam, ntrd_trc, 'trdbio' )
1470          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
1471             &             1, jpi, 1, jpj, 0, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom )
1472#endif
1473
1474      !-- Define physical units
1475      IF( ucf_trc == 1. ) THEN
1476         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit
1477      ELSEIF ( ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
1478         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf
1479      ELSE
1480         cltrcu = "unknown?"
1481      ENDIF
1482
1483      !-- Define miscellaneous passive tracer mixed-layer variables
1484      IF( jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb ) THEN
1485         STOP 'Error : jpltrd_trc /= jpmld_trc_atf .OR.  jpltrd_trc - 1 /= jpmld_trc_radb'  ! see below
1486      ENDIF
1487
1488      DO jn = 1, jptra
1489         !
1490         IF( luttrd(jn) ) THEN
1491            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
1492            CALL histdef(nidtrd(jn), clvar,           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
1493              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
1494            CALL histdef(nidtrd(jn), clvar//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
1495              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
1496            CALL histdef(nidtrd(jn), clvar//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
1497              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
1498         
1499            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmld_trc_atf
1500               CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
1501                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1502            END DO                                                                         ! if zsto=rdt above
1503         
1504            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_radb,1), & 
1505              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1506         
1507            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_atf,1),   & 
1508              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
1509         
1510            CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpltrd_trc+1,2)),  clmxl//" "//clvar//ctrd_trc(jpltrd_trc+1 ,1),   & 
1511              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! Total EIV
1512         !
1513         ENDIF
1514      END DO
1515
1516#if defined key_lobster
1517      DO jl = 1, jp_lobster_trd
1518         CALL histdef(nidtrdbio, TRIM("ML_"//ctrd_bio(jl,2)), TRIM(clmxl//" ML_"//ctrd_bio(jl,1))   ,            &
1519             &    cltrcu, jpi, jpj, nh_tb, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
1520      END DO                                                                         ! if zsto=rdt above
1521#endif
1522
1523      !-- Leave IOIPSL/NetCDF define mode
1524      DO jn = 1, jptra
1525         IF( luttrd(jn) )  CALL histend( nidtrd(jn) )
1526      END DO
1527
1528#if defined key_lobster
1529      !-- Leave IOIPSL/NetCDF define mode
1530      CALL histend( nidtrdbio )
1531
1532      IF(lwp) WRITE(numout,*)
1533       IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization for ML bio trends'
1534#endif
1535
1536#endif        /* key_dimgout */
1537   END SUBROUTINE trd_mld_trc_init
1538
1539#else
1540   !!----------------------------------------------------------------------
1541   !!   Default option :                                       Empty module
1542   !!----------------------------------------------------------------------
1543
1544   INTERFACE trd_mod_trc
1545      MODULE PROCEDURE trd_mod_trc_trp, trd_mod_trc_bio
1546   END INTERFACE
1547
1548CONTAINS
1549
1550   SUBROUTINE trd_mld_trc( kt )                                   ! Empty routine
1551      INTEGER, INTENT( in) ::   kt
1552      WRITE(*,*) 'trd_mld_trc: You should not have seen this print! error?', kt
1553   END SUBROUTINE trd_mld_trc
1554
1555   SUBROUTINE trd_mld_bio( kt )
1556      INTEGER, INTENT( in) ::   kt
1557      WRITE(*,*) 'trd_mld_bio: You should not have seen this print! error?', kt
1558   END SUBROUTINE trd_mld_bio
1559
1560   SUBROUTINE trd_mod_trc_bio( ptrbio, ktrd, kt )
1561      INTEGER               , INTENT( in )     ::   kt      ! time step
1562      INTEGER               , INTENT( in )     ::   ktrd    ! bio trend index
1563      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrbio  ! Bio trend
1564      WRITE(*,*) 'trd_mod_trc_bio : You should not have seen this print! error?', ptrbio(1,1,1)
1565      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
1566      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
1567   END SUBROUTINE trd_mod_trc_bio
1568
1569   SUBROUTINE trd_mod_trc_trp( ptrtrd, kjn, ktrd, kt )
1570      INTEGER               , INTENT( in )     ::   kt      ! time step
1571      INTEGER               , INTENT( in )     ::   kjn     ! tracer index
1572      INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
1573      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
1574      WRITE(*,*) 'trd_mod_trc_trp : You should not have seen this print! error?', ptrtrd(1,1,1)
1575      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
1576      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
1577      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kt
1578   END SUBROUTINE trd_mod_trc_trp
1579
1580   SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn )
1581      INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
1582      CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
1583      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmld            ! passive trc trend
1584      WRITE(*,*) 'trd_mld_trc_zint: You should not have seen this print! error?', ptrc_trdmld(1,1,1)
1585      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
1586      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
1587      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
1588   END SUBROUTINE trd_mld_trc_zint
1589
1590   SUBROUTINE trd_mld_trc_init                                    ! Empty routine
1591      WRITE(*,*) 'trd_mld_trc_init: You should not have seen this print! error?'
1592   END SUBROUTINE trd_mld_trc_init
1593#endif
1594
1595   !!======================================================================
1596END MODULE trdmld_trc
Note: See TracBrowser for help on using the repository browser.