New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trdmxl_trc.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trdmxl_trc.F90 @ 11480

Last change on this file since 11480 was 10966, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert TOP routines in TOP/TRP directory and all knock on effects of these conversions. SETTE tested (GYRE_PISCES only)
  • Property svn:keywords set to Id
File size: 51.0 KB
Line 
1MODULE trdmxl_trc
2   !!======================================================================
3   !!                       ***  MODULE  trdmxl_trc  ***
4   !! Ocean diagnostics:  mixed layer passive tracer trends
5   !!======================================================================
6   !! History :  9.0  !  06-08  (C. Deltel)  Original code (from trdmxl.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_trdmxl_trc
11   !!----------------------------------------------------------------------
12   !!   'key_trdmxl_trc'                      mixed layer trend diagnostics
13   !!----------------------------------------------------------------------
14   !!   trd_mxl_trc      : passive tracer cumulated trends averaged over ML
15   !!   trd_mxl_trc_zint : passive tracer trends vertical integration
16   !!   trd_mxl_trc_init : initialization step
17   !!----------------------------------------------------------------------
18   USE trc               ! tracer definitions (tr etc.)
19   USE trc_oce, ONLY :   nn_dttrc  ! frequency of step on passive tracers
20   USE dom_oce           ! domain definition
21   USE zdfmxl  , ONLY : nmln ! number of level in the mixed layer
22   USE zdf_oce , ONLY : avs  ! vert. diffusivity coef. at w-point for temp 
23   USE trdtrc_oce    ! definition of main arrays used for trends computations
24   USE in_out_manager    ! I/O manager
25   USE dianam            ! build the name of file (routine)
26   USE ldfslp            ! iso-neutral slopes
27   USE ioipsl            ! NetCDF library
28   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
29   USE lib_mpp           ! MPP library
30   USE trdmxl_trc_rst    ! restart for diagnosing the ML trends
31   USE prtctl            ! print control
32   USE sms_pisces        ! PISCES bio-model
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC trd_mxl_trc
38   PUBLIC trd_mxl_trc_alloc
39   PUBLIC trd_mxl_trc_init
40   PUBLIC trd_mxl_trc_zint
41
42   CHARACTER (LEN=40) ::  clhstnam                                ! name of the trends NetCDF file
43   INTEGER ::   nmoymltrd
44   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::   ndextrd1, nidtrd, nh_t
45   INTEGER ::   ndimtrd1                       
46   INTEGER, SAVE ::  ionce, icount
47   LOGICAL :: llwarn  = .TRUE.                                    ! this should always be .TRUE.
48   LOGICAL :: lldebug = .TRUE.
49
50   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztmltrd2   !
51
52   !!----------------------------------------------------------------------
53   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
54   !! $Id$
55   !! Software governed by the CeCILL license (see ./LICENSE)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   INTEGER FUNCTION trd_mxl_trc_alloc()
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE trd_mxl_trc_alloc  ***
62      !!----------------------------------------------------------------------
63      ALLOCATE( ztmltrd2(jpi,jpj,jpltrd_trc,jptra) ,      &
64         &      ndextrd1(jpi*jpj), nidtrd(jptra), nh_t(jptra),  STAT=trd_mxl_trc_alloc)
65         !
66      CALL mpp_sum ( 'trdmxl_trc', trd_mxl_trc_alloc )
67      IF( trd_mxl_trc_alloc /=0 )   CALL ctl_stop( 'STOP', 'trd_mxl_trc_alloc: failed to allocate arrays' )
68      !
69   END FUNCTION trd_mxl_trc_alloc
70
71
72   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn, Kmm )
73      !!----------------------------------------------------------------------
74      !!                  ***  ROUTINE trd_mxl_trc_zint  ***
75      !!
76      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments
77      !!                to the subroutine. This vertical average is performed from ocean
78      !!                surface down to a chosen control surface.
79      !!
80      !! ** Method/usage :
81      !!      The control surface can be either a mixed layer depth (time varying)
82      !!      or a fixed surface (jk level or bowl).
83      !!      Choose control surface with nctls_trc in namelist NAMTRD :
84      !!        nctls_trc = -2  : use isopycnal surface
85      !!        nctls_trc = -1  : use euphotic layer with light criterion
86      !!        nctls_trc =  0  : use mixed layer with density criterion
87      !!        nctls_trc =  1  : read index from file 'ctlsurf_idx'
88      !!        nctls_trc >  1  : use fixed level surface jk = nctls_trc
89      !!      Note: in the remainder of the routine, the volume between the
90      !!            surface and the control surface is called "mixed-layer"
91      !!----------------------------------------------------------------------
92      !!
93      INTEGER, INTENT( in ) ::   ktrd, kjn                        ! ocean trend index and passive tracer rank
94      INTEGER, INTENT( in ) ::   Kmm                              ! time level index
95      CHARACTER(len=2), INTENT( in ) ::  ctype                    ! surface/bottom (2D) or interior (3D) physics
96      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::  ptrc_trdmxl ! passive tracer trend
97      !
98      INTEGER ::   ji, jj, jk, isum
99      REAL(wp), DIMENSION(jpi,jpj) :: zvlmsk
100      !!----------------------------------------------------------------------
101
102      ! I. Definition of control surface and integration weights
103      ! --------------------------------------------------------
104
105      ONCE_PER_TIME_STEP : IF( icount == 1 ) THEN
106         !
107         tmltrd_trc(:,:,:,:) = 0.e0                               ! <<< reset trend arrays to zero
108         
109         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
110         SELECT CASE ( nn_ctls_trc )                                ! choice of the control surface
111            CASE ( -2  )   ;   STOP 'trdmxl_trc : not ready '     !     -> isopycnal surface (see ???)
112            CASE ( -1  )   ;   nmld_trc(:,:) = neln(:,:)          !     -> euphotic layer with light criterion
113            CASE (  0  )   ;   nmld_trc(:,:) = nmln(:,:)          !     -> ML with density criterion (see zdfmxl)
114            CASE (  1  )   ;   nmld_trc(:,:) = nbol_trc(:,:)          !     -> read index from file
115            CASE (  2: )   ;   nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 )
116                               nmld_trc(:,:) = nn_ctls_trc + 1      !     -> model level
117         END SELECT
118
119         ! ... Compute ndextrd1 and ndimtrd1  ??? role de jpktrd_trc
120         IF( ionce == 1 ) THEN
121            !
122            isum  = 0   ;   zvlmsk(:,:) = 0.e0
123
124            IF( jpktrd_trc < jpk ) THEN                           ! description ???
125               DO jj = 1, jpj
126                  DO ji = 1, jpi
127                     IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
128                        zvlmsk(ji,jj) = tmask(ji,jj,1)
129                     ELSE
130                        isum = isum + 1
131                        zvlmsk(ji,jj) = 0.e0
132                     ENDIF
133                  END DO
134               END DO
135            ENDIF
136
137            IF( isum > 0 ) THEN                                   ! index of ocean points (2D only)
138               WRITE(numout,*)' tmltrd_trc : Number of invalid points nmld_trc > jpktrd', isum 
139               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )
140            ELSE
141               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )
142            ENDIF                               
143
144            ionce = 0                                             ! no more pass here
145            !
146         ENDIF   ! ionce == 1
147         
148         ! ... Weights for vertical averaging
149         wkx_trc(:,:,:) = 0.e0
150         DO jk = 1, jpktrd_trc                                    ! initialize wkx_trc with vertical scale factor in mixed-layer
151            DO jj = 1, jpj
152               DO ji = 1, jpi
153                  IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
154               END DO
155            END DO
156         END DO
157         
158         rmld_trc(:,:) = 0.e0
159         DO jk = 1, jpktrd_trc                                    ! compute mixed-layer depth : rmld_trc
160            rmld_trc(:,:) = rmld_trc(:,:) + wkx_trc(:,:,jk)
161         END DO
162         
163         DO jk = 1, jpktrd_trc                                    ! compute integration weights
164            wkx_trc(:,:,jk) = wkx_trc(:,:,jk) / MAX( 1., rmld_trc(:,:) )
165         END DO
166
167         icount = 0                                               ! <<< flag = off : control surface & integr. weights
168         !                                                        !     computed only once per time step
169      ENDIF ONCE_PER_TIME_STEP
170
171      ! II. Vertical integration of trends in the mixed-layer
172      ! -----------------------------------------------------
173
174      SELECT CASE ( ctype )
175         CASE ( '3D' )                                            ! mean passive tracer trends in the mixed-layer
176            DO jk = 1, jpktrd_trc
177               tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,jk) * wkx_trc(:,:,jk)   
178            END DO
179         CASE ( '2D' )                                            ! forcing at upper boundary of the mixed-layer
180            tmltrd_trc(:,:,ktrd,kjn) = tmltrd_trc(:,:,ktrd,kjn) + ptrc_trdmxl(:,:,1) * wkx_trc(:,:,1)  ! non penetrative
181      END SELECT
182      !
183   END SUBROUTINE trd_mxl_trc_zint
184
185
186   SUBROUTINE trd_mxl_trc( kt, Kmm )
187      !!----------------------------------------------------------------------
188      !!                  ***  ROUTINE trd_mxl_trc  ***
189      !!
190      !! ** Purpose :  Compute and cumulate the mixed layer trends over an analysis
191      !!               period, and write NetCDF outputs.
192      !!
193      !! ** Method/usage :
194      !!          The stored trends can be chosen twofold (according to the ln_trdmxl_trc_instant
195      !!          logical namelist variable) :
196      !!          1) to explain the difference between initial and final
197      !!             mixed-layer T & S (where initial and final relate to the
198      !!             current analysis window, defined by ntrc_trc in the namelist)
199      !!          2) to explain the difference between the current and previous
200      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
201      !!             performed over each analysis window).
202      !!
203      !! ** Consistency check :
204      !!        If the control surface is fixed ( nctls_trc > 1 ), the residual term (dh/dt
205      !!        entrainment) should be zero, at machine accuracy. Note that in the case
206      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
207      !!        over the first two analysis windows (except if restart).
208      !!        N.B. For ORCA2_ICE, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8
209      !!             for checking residuals.
210      !!             On a NEC-SX5 computer, this typically leads to:
211      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.false.
212      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmxl_trc_instant=.true.
213      !!
214      !! ** Action :
215      !!       At each time step, mixed-layer averaged trends are stored in the
216      !!       tmltrd(:,:,jpmxl_xxx) array (see trdmxl_oce.F90 for definitions of jpmxl_xxx).
217      !!       This array is known when trd_mld is called, at the end of the stp subroutine,
218      !!       except for the purely vertical K_z diffusion term, which is embedded in the
219      !!       lateral diffusion trend.
220      !!
221      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
222      !!       from the lateral diffusion trend.
223      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
224      !!       arrays are updated.
225      !!       In III), called only once per analysis window, we compute the total trends,
226      !!       along with the residuals and the Asselin correction terms.
227      !!       In IV), the appropriate trends are written in the trends NetCDF file.
228      !!
229      !! References :
230      !!       - Vialard & al.
231      !!       - See NEMO documentation (in preparation)
232      !!----------------------------------------------------------------------
233      !
234      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
235      INTEGER, INTENT(in) ::   Kmm                              ! time level index
236      !
237      INTEGER ::   ji, jj, jk, jl, ik, it, itmod, jn
238      REAL(wp) ::   zavt, zfn, zfn2
239      !
240      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltot             ! d(trc)/dt over the anlysis window (incl. Asselin)
241      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlres             ! residual = dh/dt entrainment term
242      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlatf             ! for storage only
243      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlrad             ! for storage only (for trb<0 corr in trcrad)
244      !
245      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltot2            ! -+
246      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlres2            !  | working arrays to diagnose the trends
247      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmltrdm2           !  | associated with the time meaned ML
248      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlatf2            !  | passive tracers
249      REAL(wp), DIMENSION(jpi,jpj,jptra) ::   ztmlrad2            !  | (-> for trb<0 corr in trcrad)
250      !
251      CHARACTER (LEN=10) ::   clvar
252      !!----------------------------------------------------------------------
253
254
255      IF( nn_dttrc  /= 1  )   CALL ctl_stop( " Be careful, trends diags never validated " )
256
257      ! ======================================================================
258      ! I. Diagnose the purely vertical (K_z) diffusion trend
259      ! ======================================================================
260
261      ! ... These terms can be estimated by flux computation at the lower boundary of the ML
262      !     (we compute (-1/h) * K_z * d_z( tracer ))
263
264      IF( ln_trcldf_iso ) THEN
265         !
266         DO jn = 1, jptra
267            DO jj = 1, jpj
268               DO ji = 1, jpi
269                  ik = nmld_trc(ji,jj)
270                  IF( ln_trdtrc(jn) )    &
271                  tmltrd_trc(ji,jj,jpmxl_trc_zdf,jn) = - avs(ji,jj,ik) / e3w(ji,jj,ik,Kmm) * tmask(ji,jj,ik)  &
272                       &                    * ( tr(ji,jj,ik-1,jn,Kmm) - tr(ji,jj,ik,jn,Kmm) )            &
273                       &                    / MAX( 1., rmld_trc(ji,jj) ) * tmask(ji,jj,1)
274               END DO
275            END DO
276         END DO
277
278         DO jn = 1, jptra
279            ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
280            IF( ln_trdtrc(jn) ) &
281                 tmltrd_trc(:,:,jpmxl_trc_ldf,jn) = tmltrd_trc(:,:,jpmxl_trc_ldf,jn) - tmltrd_trc(:,:,jpmxl_trc_zdf,jn)
282   
283         END DO
284         !     
285      ENDIF
286
287!!gm Test removed, nothing specific to a configuration should survive out of usrdef modules
288!!gm      IF ( cn_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration
289!!gm      ! GYRE : for diagnostic fields, are needed if cyclic B.C. are present, but not for purely MPI comm.
290!!gm      ! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.)
291         DO jn = 1, jptra
292            IF( ln_trdtrc(jn) ) THEN
293               DO jl = 1, jpltrd_trc
294                  CALL lbc_lnk( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1. )        ! lateral boundary conditions
295               END DO
296            ENDIF
297         END DO
298!!gm      ENDIF
299     
300      ! ======================================================================
301      ! II. Cumulate the trends over the analysis window
302      ! ======================================================================
303
304      ztmltrd2(:,:,:,:) = 0.e0   ;   ztmltot2(:,:,:)   = 0.e0     ! <<< reset arrays to zero
305      ztmlres2(:,:,:)   = 0.e0   ;   ztmlatf2(:,:,:)   = 0.e0
306      ztmlrad2(:,:,:)   = 0.e0
307
308      ! II.1 Set before values of vertically averages passive tracers
309      ! -------------------------------------------------------------
310      IF( kt > nittrc000 ) THEN
311         DO jn = 1, jptra
312            IF( ln_trdtrc(jn) ) THEN
313               tmlb_trc   (:,:,jn) = tml_trc   (:,:,jn)
314               tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_atf,jn)
315               tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_radb,jn)
316            ENDIF
317         END DO
318      ENDIF
319
320      ! II.2 Vertically averaged passive tracers
321      ! ----------------------------------------
322      tml_trc(:,:,:) = 0.e0
323      DO jk = 1, jpktrd_trc ! - 1 ???
324         DO jn = 1, jptra
325            IF( ln_trdtrc(jn) ) &
326               tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * tr(:,:,jk,jn,Kmm)
327         END DO
328      END DO
329
330      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window   
331      ! ------------------------------------------------------------------------
332      IF( kt == nittrc000 + nn_dttrc ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)    ???
333         !
334         DO jn = 1, jptra
335            IF( ln_trdtrc(jn) ) THEN
336               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
337               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
338               
339               tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0   ;   tmltrd_atf_sumb_trc  (:,:,jn) = 0.e0
340               tmltrd_rad_sumb_trc  (:,:,jn) = 0.e0
341            ENDIF
342         END DO
343         
344         rmldbn_trc(:,:) = rmld_trc(:,:)
345         !
346      ENDIF
347
348      ! II.4 Cumulated trends over the analysis period
349      ! ----------------------------------------------
350      !
351      !         [  1rst analysis window ] [     2nd analysis window     ]                       
352      !
353      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
354      !                            ntrd                             2*ntrd       etc.
355      !     1      2     3     4    =5 e.g.                          =10
356      !
357      IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN                        ! ???
358         !
359         nmoymltrd = nmoymltrd + 1
360
361
362         ! ... Cumulate over BOTH physical contributions AND over time steps
363         DO jn = 1, jptra
364            IF( ln_trdtrc(jn) ) THEN
365               DO jl = 1, jpltrd_trc
366                  tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
367               END DO
368            ENDIF
369         END DO
370
371         DO jn = 1, jptra
372            IF( ln_trdtrc(jn) ) THEN
373               ! ... Special handling of the Asselin trend
374               tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
375               tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
376
377               ! ... Trends associated with the time mean of the ML passive tracers
378               tmltrd_sum_trc    (:,:,:,jn) = tmltrd_sum_trc    (:,:,:,jn) + tmltrd_trc    (:,:,:,jn)
379               tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
380               tml_sum_trc       (:,:,jn)   = tml_sum_trc       (:,:,jn)   + tml_trc       (:,:,jn)
381            ENDIF
382         ENDDO
383
384         rmld_sum_trc      (:,:)     = rmld_sum_trc      (:,:)     + rmld_trc      (:,:)
385         !
386      ENDIF
387
388      ! ======================================================================
389      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
390      ! ======================================================================
391
392      ! Convert to appropriate physical units
393      tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
394
395      itmod = kt - nittrc000 + 1
396      it    = kt
397
398      MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN           ! nitend MUST be multiple of nn_trd_trc
399         !
400         ztmltot (:,:,:) = 0.e0                                   ! reset arrays to zero
401         ztmlres (:,:,:) = 0.e0
402         ztmltot2(:,:,:) = 0.e0
403         ztmlres2(:,:,:) = 0.e0
404     
405         zfn  = FLOAT( nmoymltrd )    ;    zfn2 = zfn * zfn
406         
407         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
408         ! -------------------------------------------------------------
409
410         DO jn = 1, jptra
411            IF( ln_trdtrc(jn) ) THEN
412               !-- Compute total trends    (use rdttrc instead of rdt ???)
413               IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN  ! EULER-FORWARD schemes
414                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt
415               ELSE                                                                     ! LEAP-FROG schemes
416                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rdt)
417               ENDIF
418               
419               !-- Compute residuals
420               ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
421                  &                                                 - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)  )
422               
423               !-- Diagnose Asselin trend over the analysis window
424               ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
425               ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
426               
427         !-- Lateral boundary conditions
428               IF ( cn_cfg .NE. 'gyre' ) THEN
429                  CALL lbc_lnk_multi( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1. , ztmlres(:,:,jn) , 'T', 1., &
430                     &                ztmlatf(:,:,jn) , 'T', 1. , ztmlrad(:,:,jn) , 'T', 1. )
431               ENDIF
432
433
434#if defined key_diainstant
435               STOP 'tmltrd_trc : key_diainstant was never checked within trdmxl. Comment this to proceed.'
436#endif
437            ENDIF
438         END DO
439
440         ! III.2 Prepare fields for output ("mean" diagnostics)
441         ! ----------------------------------------------------
442         
443         !-- Update the ML depth time sum (to build the Leap-Frog time mean)
444         rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
445
446               !-- Compute passive tracer total trends
447         DO jn = 1, jptra
448            IF( ln_trdtrc(jn) ) THEN
449               tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
450               ztmltot2   (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) /  ( 2.*rdt )    ! now tracer unit is /sec
451            ENDIF
452         END DO
453
454         !-- Compute passive tracer residuals
455         DO jn = 1, jptra
456            IF( ln_trdtrc(jn) ) THEN
457               !
458               DO jl = 1, jpltrd_trc
459                  ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
460               END DO
461               
462               ztmltrdm2(:,:,jn) = 0.e0
463               DO jl = 1, jpltrd_trc
464                  ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
465               END DO
466               
467               ztmlres2(:,:,jn) =  ztmltot2(:,:,jn)  -       &
468                  & ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
469                  &                     - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
470               !
471
472               !-- Diagnose Asselin trend over the analysis window
473               ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) &
474                  &                                               + tmltrd_atf_sumb_trc(:,:,jn)
475               ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) &
476                  &                                               + tmltrd_rad_sumb_trc(:,:,jn)
477
478         !-- Lateral boundary conditions
479               IF ( cn_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration   
480                  CALL lbc_lnk_multi( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1., ztmlres2(:,:,jn), 'T', 1. )
481                  DO jl = 1, jpltrd_trc
482                     CALL lbc_lnk( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
483                  END DO
484               ENDIF
485
486            ENDIF
487         END DO
488
489         ! * Debugging information *
490         IF( lldebug ) THEN
491            !
492            WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
493            WRITE(numout,*) '~~~~~~~~~~~  '
494            WRITE(numout,*)
495            WRITE(numout,*) 'TRC kt = ', kt, '    nmoymltrd = ', nmoymltrd
496
497            DO jn = 1, jptra
498
499               IF( ln_trdtrc(jn) ) THEN
500                  WRITE(numout, *)
501                  WRITE(numout, *) '>>>>>>>>>>>>>>>>>>  TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
502                 
503                  WRITE(numout, *)
504                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres     : ', SUM2D(ztmlres(:,:,jn))
505                  !CD??? PREVOIR: z2d = ztmlres(:,:,jn)   ;   CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres   -   : ')
506                 
507                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
508                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
509                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot    : ', SUM2D(+ztmltot    (:,:,jn))
510                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
511                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
512                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
513                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
514                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
515                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
516                 
517                  WRITE(numout, *)
518                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2    : ', SUM2D(ztmlres2(:,:,jn))
519                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
520                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
521                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2            : ', SUM2D(+ztmltot2(:,:,jn))
522                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2           : ', SUM2D(+ztmltrdm2(:,:,jn)) 
523                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn)) 
524                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
525                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
526                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
527                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
528                 
529                  WRITE(numout, *)
530                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot     : ', SUM2D(ztmltot    (:,:,jn))
531                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
532                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc    : ', SUM2D(tml_trc    (:,:,jn))
533                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc  : ', SUM2D(tmlbn_trc  (:,:,jn))
534                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc   : ', SUM2D(tmlb_trc   (:,:,jn))
535                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc  : ', SUM2D(tmlbb_trc  (:,:,jn))
536                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
537                 
538                  WRITE(numout, *)
539                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
540                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
541                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
542                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
543                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
544                 
545                  WRITE(numout, *)
546                  DO jl = 1, jpltrd_trc
547                     WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
548                        & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
549                  END DO
550               
551                  WRITE(numout,*) 
552                  WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** '
553                  WRITE(numout,*)
554                  WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
555                  WRITE(numout,*)
556                  WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
557                 
558                  WRITE(numout,*) 
559                  WRITE(numout,*) ' *********************** ZTMLRES *********************** '
560                  WRITE(numout,*)
561                 
562                  WRITE(numout,*) '...................................................'
563                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
564                  DO jj = 5, 1, -1
565                     WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
566                  END DO
567                 
568                  WRITE(numout,*) 
569                  WRITE(numout,*) ' *********************** ZTMLRES2 *********************** '
570                  WRITE(numout,*)
571                 
572                  WRITE(numout,*) '...................................................'
573                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
574                  DO jj = 5, 1, -1
575                     WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
576                  END DO
577                  !
578               ENDIF
579               !
580            END DO
581
582
58397            FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
58498            FORMAT(a10, i3, 2x, a30, 2x, g20.10)
58599            FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
586              WRITE(numout,*)
587            !
588         ENDIF
589
590         ! III.3 Time evolution array swap
591         ! -------------------------------
592         ! ML depth
593         rmldbn_trc(:,:)   = rmld_trc(:,:)
594         rmld_sum_trc(:,:)     = rmld_sum_trc(:,:)     /      (2*zfn)  ! similar to tml_sum and sml_sum
595         DO jn = 1, jptra
596            IF( ln_trdtrc(jn) ) THEN       
597               ! For passive tracer instantaneous diagnostics
598               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
599               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
600               
601               ! For passive tracer mean diagnostics
602               tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
603               tml_sumb_trc       (:,:,jn)   = tml_sum_trc(:,:,jn)
604               tmltrd_atf_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
605               tmltrd_rad_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
606               
607               
608               ! III.4 Convert to appropriate physical units
609               ! -------------------------------------------
610               ztmltot     (:,:,jn)   = ztmltot     (:,:,jn)   * rn_ucf_trc/zfn   ! instant diags
611               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * rn_ucf_trc/zfn
612               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * rn_ucf_trc/zfn
613               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * rn_ucf_trc/zfn
614               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
615               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * rn_ucf_trc/zfn2
616               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * rn_ucf_trc/zfn2
617               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * rn_ucf_trc/zfn2
618               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * rn_ucf_trc/zfn2
619               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * rn_ucf_trc/zfn2
620            ENDIF
621         END DO
622         !
623      ENDIF MODULO_NTRD
624
625      ! ======================================================================
626      ! IV. Write trends in the NetCDF file
627      ! ======================================================================
628
629      ! IV.1 Code for IOIPSL/NetCDF output
630      ! ----------------------------------
631
632      IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
633         WRITE(numout,*) ' '
634         WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
635         WRITE(numout,*) '~~~~~~~~~~~ '
636         WRITE(numout,*) '          ', trim(clhstnam), ' at kt = ', kt
637         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
638         WRITE(numout,*) ' '
639      ENDIF
640         
641      NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
642         !
643
644         DO jn = 1, jptra
645            !
646            IF( ln_trdtrc(jn) ) THEN
647               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
648               !-- Output the fields
649               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
650               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
651               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
652               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
653           
654               DO jl = 1, jpltrd_trc - 2
655                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),             &
656                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
657               END DO
658
659               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
660                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
661
662               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
663                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
664                     
665            ENDIF
666         END DO
667
668         IF( kt == nitend ) THEN
669            DO jn = 1, jptra
670               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
671            END DO
672         ENDIF
673
674      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
675         
676         DO jn = 1, jptra
677            !
678            IF( ln_trdtrc(jn) ) THEN
679               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 
680               !-- Output the fields
681               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
682
683               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
684               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
685               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
686
687               DO jl = 1, jpltrd_trc - 2
688                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),           &
689                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
690               END DO
691           
692               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
693                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
694
695               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
696                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
697
698            ENDIF 
699            !
700         END DO
701         IF( kt == nitend ) THEN
702            DO jn = 1, jptra
703               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
704            END DO
705         ENDIF
706
707         !
708      ENDIF NETCDF_OUTPUT
709         
710      ! Compute the control surface (for next time step) : flag = on
711      icount = 1
712
713      IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
714         !
715         ! Reset cumulative arrays to zero
716         ! -------------------------------         
717         nmoymltrd = 0
718         tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
719         tmlradm_trc        (:,:,:)   = 0.e0   ;   tml_sum_trc        (:,:,:)   = 0.e0
720         tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0
721         rmld_sum_trc       (:,:)     = 0.e0
722         !
723      ENDIF
724
725      ! ======================================================================
726      ! V. Write restart file
727      ! ======================================================================
728
729      IF( lrst_trc )   CALL trd_mxl_trc_rst_write( kt )  ! this must be after the array swap above (III.3)
730      !
731   END SUBROUTINE trd_mxl_trc
732
733   REAL FUNCTION sum2d( ztab )
734      !!----------------------------------------------------------------------
735      !! CD ??? prevoir d'utiliser plutot prtctl
736      !!----------------------------------------------------------------------
737      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab     
738      !!----------------------------------------------------------------------
739      sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
740   END FUNCTION sum2d
741
742
743   SUBROUTINE trd_mxl_trc_init
744      !!----------------------------------------------------------------------
745      !!                  ***  ROUTINE trd_mxl_init  ***
746      !!
747      !! ** Purpose :   computation of vertically integrated T and S budgets
748      !!      from ocean surface down to control surface (NetCDF output)
749      !!
750      !! ** Method/usage :
751      !!
752      !!----------------------------------------------------------------------
753      INTEGER :: inum   ! logical unit
754      INTEGER :: ilseq, jl, jn, iiter
755      REAL(wp) ::   zjulian, zsto, zout
756      CHARACTER (LEN=40) ::   clop
757      CHARACTER (LEN=15) ::   csuff
758      CHARACTER (LEN=12) ::   clmxl
759      CHARACTER (LEN=16) ::   cltrcu
760      CHARACTER (LEN=10) ::   clvar
761
762      !!----------------------------------------------------------------------
763
764      ! ======================================================================
765      ! I. initialization
766      ! ======================================================================
767
768      IF(lwp) THEN
769         WRITE(numout,*)
770         WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers                '
771         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
772         WRITE(numout,*)
773      ENDIF
774
775     
776      ! I.1 Check consistency of user defined preferences
777      ! -------------------------------------------------
778
779      IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
780         WRITE(ctmp1,*) '                Your nitend parameter, nitend = ', nitend
781         WRITE(ctmp2,*) '                is no multiple of the trends diagnostics frequency        '
782         WRITE(ctmp3,*) '                          you defined, nn_trd_trc   = ', nn_trd_trc
783         WRITE(ctmp4,*) '                This will not allow you to restart from this simulation.  '
784         WRITE(ctmp5,*) '                You should reconsider this choice.                        ' 
785         WRITE(ctmp6,*) 
786         WRITE(ctmp7,*) '                N.B. the nitend parameter is also constrained to be a     '
787         WRITE(ctmp8,*) '                multiple of the sea-ice frequency parameter (typically 5) '
788         CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4, ctmp5, ctmp6, ctmp7, ctmp8 )
789      ENDIF
790
791      ! * Debugging information *
792      IF( lldebug ) THEN
793         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
794         WRITE(numout,*) '               ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
795      ENDIF
796
797      IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
798         WRITE(ctmp1,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
799         WRITE(ctmp2,*) '                advection and window averaged diagnostics of ML trends.   '
800         WRITE(ctmp3,*) '                WHY? Everything in trdmxl_trc is coded for leap-frog, and '
801         WRITE(ctmp4,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
802         WRITE(ctmp5,*) '                that MUSCL is leap-frog for active tracers T/S).          '
803         WRITE(ctmp6,*) '                In particuliar, entrainment trend would be FALSE. However '
804         WRITE(ctmp7,*) '                this residual is correct for instantaneous ML diagnostics.'
805         CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4, ctmp5, ctmp6, ctmp7 )
806      ENDIF
807
808      ! I.2 Initialize arrays to zero or read a restart file
809      ! ----------------------------------------------------
810      nmoymltrd   = 0
811
812      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
813      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
814      tmlradm_trc        (:,:,:)   = 0.e0
815
816      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
817      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
818
819      IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
820         CALL trd_mxl_trc_rst_read
821      ELSE
822         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
823         tmlbn_trc          (:,:,:)   = 0.e0
824
825         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
826         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
827
828       ENDIF
829
830      icount = 1   ;   ionce  = 1  ! open specifier   
831
832
833      ! I.3 Read control surface from file ctlsurf_idx
834      ! ----------------------------------------------
835      IF( nn_ctls_trc == 1 ) THEN
836         CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
837         READ ( inum ) nbol_trc
838         CLOSE( inum )
839      ENDIF
840
841      ! ======================================================================
842      ! II. netCDF output initialization
843      ! ======================================================================
844
845      ! clmxl = legend root for netCDF output
846      IF( nn_ctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
847         clmxl = 'Mixed Layer '
848      ELSE IF( nn_ctls_trc == 1 ) THEN                              ! control surface = read index from file
849         clmxl = '      Bowl '
850      ELSE IF( nn_ctls_trc >= 2 ) THEN                              ! control surface = model level
851         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
852      ENDIF
853
854      ! II.1 Define frequency of output and means
855      ! -----------------------------------------
856      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
857      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cp time)
858      ENDIF
859#  if defined key_diainstant
860      IF( .NOT. ln_trdmxl_trc_instant ) THEN
861         STOP 'trd_mxl_trc : this was never checked. Comment this line to proceed...'
862      ENDIF
863      zsto = nn_trd_trc * rdt
864      clop = "inst("//TRIM(clop)//")"
865#  else
866      IF( ln_trdmxl_trc_instant ) THEN
867         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
868      ELSE
869         zsto = nn_trd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
870      ENDIF
871      clop = "ave("//TRIM(clop)//")"
872#  endif
873      zout = nn_trd_trc * rdt
874      iiter = ( nittrc000 - 1 ) / nn_dttrc
875
876      IF(lwp) WRITE (numout,*) '                netCDF initialization'
877
878      ! II.2 Compute julian date from starting date of the run
879      ! ------------------------------------------------------
880      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
881      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
882      IF(lwp) WRITE(numout,*)' ' 
883      IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000               &
884           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
885           &   ,'Julian day : ', zjulian
886
887      ! II.3 Define the T grid trend file (nidtrd)
888      ! ------------------------------------------
889
890      !-- Define long and short names for the NetCDF output variables
891      !       ==> choose them according to trdmxl_trc_oce.F90 <==
892
893      ctrd_trc(jpmxl_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmxl_trc_xad    ,2) = "_xad"
894      ctrd_trc(jpmxl_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmxl_trc_yad    ,2) = "_yad"
895      ctrd_trc(jpmxl_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmxl_trc_zad    ,2) = "_zad"
896      ctrd_trc(jpmxl_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmxl_trc_ldf    ,2) = "_ldf"
897      ctrd_trc(jpmxl_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmxl_trc_zdf    ,2) = "_zdf"
898      ctrd_trc(jpmxl_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmxl_trc_bbl    ,2) = "_bbl"
899      ctrd_trc(jpmxl_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmxl_trc_dmp    ,2) = "_dmp"
900      ctrd_trc(jpmxl_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmxl_trc_sbc    ,2) = "_sbc"
901      ctrd_trc(jpmxl_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmxl_trc_sms    ,2) = "_sms"
902      ctrd_trc(jpmxl_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radb   ,2) = "_radb"
903      ctrd_trc(jpmxl_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radn   ,2) = "_radn"
904      ctrd_trc(jpmxl_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmxl_trc_atf    ,2) = "_atf"
905
906      DO jn = 1, jptra     
907      !-- Create a NetCDF file and enter the define mode
908         IF( ln_trdtrc(jn) ) THEN
909            csuff="ML_"//ctrcnm(jn)
910            CALL dia_nam( clhstnam, nn_trd_trc, csuff )
911            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
912               &        1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
913     
914            !-- Define the ML depth variable
915            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
916               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
917
918         ENDIF
919      END DO
920
921      !-- Define physical units
922      IF( rn_ucf_trc == 1. ) THEN
923         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit
924      ELSEIF ( rn_ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
925         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf
926      ELSE
927         cltrcu = "unknown?"
928      ENDIF
929
930      !-- Define miscellaneous passive tracer mixed-layer variables
931      IF( jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
932         STOP 'Error : jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb'  ! see below
933      ENDIF
934
935      DO jn = 1, jptra
936         !
937         IF( ln_trdtrc(jn) ) THEN
938            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
939            CALL histdef(nidtrd(jn), trim(clvar),           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
940              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
941            CALL histdef(nidtrd(jn), trim(clvar)//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
942              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
943            CALL histdef(nidtrd(jn), trim(clvar)//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
944              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
945         
946            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmxl_trc_atf
947               CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
948                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
949            END DO                                                                         ! if zsto=rdt above
950         
951            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), & 
952              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
953         
954            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1),   & 
955              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
956         !
957         ENDIF
958      END DO
959
960      !-- Leave IOIPSL/NetCDF define mode
961      DO jn = 1, jptra
962         IF( ln_trdtrc(jn) )  CALL histend( nidtrd(jn), snc4set )
963      END DO
964
965      IF(lwp) WRITE(numout,*)
966
967   END SUBROUTINE trd_mxl_trc_init
968
969#else
970   !!----------------------------------------------------------------------
971   !!   Default option :                                       Empty module
972   !!----------------------------------------------------------------------
973CONTAINS
974   SUBROUTINE trd_mxl_trc( kt, Kmm )                                   ! Empty routine
975      INTEGER, INTENT( in) ::   kt
976      INTEGER, INTENT( in) ::   Kmm            ! time level index
977      WRITE(*,*) 'trd_mxl_trc: You should not have seen this print! error?', kt
978   END SUBROUTINE trd_mxl_trc
979   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn, Kmm )
980      INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
981      INTEGER               , INTENT( in ) ::  Kmm                    ! time level index
982      CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
983      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmxl            ! passive trc trend
984      WRITE(*,*) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
985      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
986      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
987      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
988   END SUBROUTINE trd_mxl_trc_zint
989   SUBROUTINE trd_mxl_trc_init                                    ! Empty routine
990      WRITE(*,*) 'trd_mxl_trc_init: You should not have seen this print! error?'
991   END SUBROUTINE trd_mxl_trc_init
992#endif
993
994   !!======================================================================
995END MODULE trdmxl_trc
Note: See TracBrowser for help on using the repository browser.