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 branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 8882

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864

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