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_merge_2017/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 9050

Last change on this file since 9050 was 9050, checked in by cetlod, 6 years ago

First use of lbc_lnk_multi in TOP routines

  • Property svn:keywords set to Id
File size: 51.2 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_multi( ztmltot(:,:,jn) , 'T', 1. , ztmlres(:,:,jn) , 'T', 1., &
436                     &                ztmlatf(:,:,jn) , 'T', 1. , 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_multi( ztmltot2(:,:,jn), 'T', 1., ztmlres2(:,:,jn), 'T', 1. )
487                  DO jl = 1, jpltrd_trc
488                     CALL lbc_lnk( ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
489                  END DO
490               ENDIF
491
492            ENDIF
493         END DO
494
495         ! * Debugging information *
496         IF( lldebug ) THEN
497            !
498            WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
499            WRITE(numout,*) '~~~~~~~~~~~  '
500            WRITE(numout,*)
501            WRITE(numout,*) 'TRC kt = ', kt, '    nmoymltrd = ', nmoymltrd
502
503            DO jn = 1, jptra
504
505               IF( ln_trdtrc(jn) ) THEN
506                  WRITE(numout, *)
507                  WRITE(numout, *) '>>>>>>>>>>>>>>>>>>  TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
508                 
509                  WRITE(numout, *)
510                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres     : ', SUM2D(ztmlres(:,:,jn))
511                  !CD??? PREVOIR: z2d = ztmlres(:,:,jn)   ;   CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres   -   : ')
512                 
513                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
514                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
515                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot    : ', SUM2D(+ztmltot    (:,:,jn))
516                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
517                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
518                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
519                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
520                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
521                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
522                 
523                  WRITE(numout, *)
524                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2    : ', SUM2D(ztmlres2(:,:,jn))
525                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
526                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
527                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2            : ', SUM2D(+ztmltot2(:,:,jn))
528                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2           : ', SUM2D(+ztmltrdm2(:,:,jn)) 
529                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn)) 
530                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
531                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
532                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
533                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
534                 
535                  WRITE(numout, *)
536                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot     : ', SUM2D(ztmltot    (:,:,jn))
537                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
538                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc    : ', SUM2D(tml_trc    (:,:,jn))
539                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc  : ', SUM2D(tmlbn_trc  (:,:,jn))
540                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc   : ', SUM2D(tmlb_trc   (:,:,jn))
541                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc  : ', SUM2D(tmlbb_trc  (:,:,jn))
542                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
543                 
544                  WRITE(numout, *)
545                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
546                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
547                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
548                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
549                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
550                 
551                  WRITE(numout, *)
552                  DO jl = 1, jpltrd_trc
553                     WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
554                        & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
555                  END DO
556               
557                  WRITE(numout,*) 
558                  WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** '
559                  WRITE(numout,*)
560                  WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
561                  WRITE(numout,*)
562                  WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
563                 
564                  WRITE(numout,*) 
565                  WRITE(numout,*) ' *********************** ZTMLRES *********************** '
566                  WRITE(numout,*)
567                 
568                  WRITE(numout,*) '...................................................'
569                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
570                  DO jj = 5, 1, -1
571                     WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
572                  END DO
573                 
574                  WRITE(numout,*) 
575                  WRITE(numout,*) ' *********************** ZTMLRES2 *********************** '
576                  WRITE(numout,*)
577                 
578                  WRITE(numout,*) '...................................................'
579                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
580                  DO jj = 5, 1, -1
581                     WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
582                  END DO
583                  !
584               ENDIF
585               !
586            END DO
587
588
58997            FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
59098            FORMAT(a10, i3, 2x, a30, 2x, g20.10)
59199            FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
592              WRITE(numout,*)
593            !
594         ENDIF
595
596         ! III.3 Time evolution array swap
597         ! -------------------------------
598         ! ML depth
599         rmldbn_trc(:,:)   = rmld_trc(:,:)
600         rmld_sum_trc(:,:)     = rmld_sum_trc(:,:)     /      (2*zfn)  ! similar to tml_sum and sml_sum
601         DO jn = 1, jptra
602            IF( ln_trdtrc(jn) ) THEN       
603               ! For passive tracer instantaneous diagnostics
604               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
605               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
606               
607               ! For passive tracer mean diagnostics
608               tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
609               tml_sumb_trc       (:,:,jn)   = tml_sum_trc(:,:,jn)
610               tmltrd_atf_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
611               tmltrd_rad_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
612               
613               
614               ! III.4 Convert to appropriate physical units
615               ! -------------------------------------------
616               ztmltot     (:,:,jn)   = ztmltot     (:,:,jn)   * rn_ucf_trc/zfn   ! instant diags
617               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * rn_ucf_trc/zfn
618               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * rn_ucf_trc/zfn
619               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * rn_ucf_trc/zfn
620               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
621               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * rn_ucf_trc/zfn2
622               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * rn_ucf_trc/zfn2
623               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * rn_ucf_trc/zfn2
624               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * rn_ucf_trc/zfn2
625               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * rn_ucf_trc/zfn2
626            ENDIF
627         END DO
628         !
629      ENDIF MODULO_NTRD
630
631      ! ======================================================================
632      ! IV. Write trends in the NetCDF file
633      ! ======================================================================
634
635      ! IV.1 Code for IOIPSL/NetCDF output
636      ! ----------------------------------
637
638      IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
639         WRITE(numout,*) ' '
640         WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
641         WRITE(numout,*) '~~~~~~~~~~~ '
642         WRITE(numout,*) '          ', trim(clhstnam), ' at kt = ', kt
643         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
644         WRITE(numout,*) ' '
645      ENDIF
646         
647      NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
648         !
649
650         DO jn = 1, jptra
651            !
652            IF( ln_trdtrc(jn) ) THEN
653               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
654               !-- Output the fields
655               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
656               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
657               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
658               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
659           
660               DO jl = 1, jpltrd_trc - 2
661                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),             &
662                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
663               END DO
664
665               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
666                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
667
668               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
669                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
670                     
671            ENDIF
672         END DO
673
674         IF( kt == nitend ) THEN
675            DO jn = 1, jptra
676               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
677            END DO
678         ENDIF
679
680      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
681         
682         DO jn = 1, jptra
683            !
684            IF( ln_trdtrc(jn) ) THEN
685               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 
686               !-- Output the fields
687               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
688
689               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
690               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
691               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
692
693               DO jl = 1, jpltrd_trc - 2
694                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),           &
695                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
696               END DO
697           
698               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
699                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
700
701               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
702                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
703
704            ENDIF 
705            !
706         END DO
707         IF( kt == nitend ) THEN
708            DO jn = 1, jptra
709               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
710            END DO
711         ENDIF
712
713         !
714      ENDIF NETCDF_OUTPUT
715         
716      ! Compute the control surface (for next time step) : flag = on
717      icount = 1
718
719      IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
720         !
721         ! Reset cumulative arrays to zero
722         ! -------------------------------         
723         nmoymltrd = 0
724         tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
725         tmlradm_trc        (:,:,:)   = 0.e0   ;   tml_sum_trc        (:,:,:)   = 0.e0
726         tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0
727         rmld_sum_trc       (:,:)     = 0.e0
728         !
729      ENDIF
730
731      ! ======================================================================
732      ! V. Write restart file
733      ! ======================================================================
734
735      IF( lrst_trc )   CALL trd_mxl_trc_rst_write( kt )  ! this must be after the array swap above (III.3)
736
737      CALL wrk_dealloc( jpi, jpj, jptra, ztmltot , ztmlres , ztmlatf , ztmlrad             )
738      CALL wrk_dealloc( jpi, jpj, jptra, ztmltot2, ztmlres2, ztmlatf2, ztmlrad2, ztmltrdm2 )
739      !
740   END SUBROUTINE trd_mxl_trc
741
742   REAL FUNCTION sum2d( ztab )
743      !!----------------------------------------------------------------------
744      !! CD ??? prevoir d'utiliser plutot prtctl
745      !!----------------------------------------------------------------------
746      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab     
747      !!----------------------------------------------------------------------
748      sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
749   END FUNCTION sum2d
750
751
752   SUBROUTINE trd_mxl_trc_init
753      !!----------------------------------------------------------------------
754      !!                  ***  ROUTINE trd_mxl_init  ***
755      !!
756      !! ** Purpose :   computation of vertically integrated T and S budgets
757      !!      from ocean surface down to control surface (NetCDF output)
758      !!
759      !! ** Method/usage :
760      !!
761      !!----------------------------------------------------------------------
762      INTEGER :: inum   ! logical unit
763      INTEGER :: ilseq, jl, jn, iiter
764      REAL(wp) ::   zjulian, zsto, zout
765      CHARACTER (LEN=40) ::   clop
766      CHARACTER (LEN=15) ::   csuff
767      CHARACTER (LEN=12) ::   clmxl
768      CHARACTER (LEN=16) ::   cltrcu
769      CHARACTER (LEN=10) ::   clvar
770
771      !!----------------------------------------------------------------------
772
773      ! ======================================================================
774      ! I. initialization
775      ! ======================================================================
776
777      IF(lwp) THEN
778         WRITE(numout,*)
779         WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers                '
780         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
781         WRITE(numout,*)
782      ENDIF
783
784     
785      ! I.1 Check consistency of user defined preferences
786      ! -------------------------------------------------
787
788      IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
789         WRITE(numout,cform_err)
790         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
791         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
792         WRITE(numout,*) '                          you defined, nn_trd_trc   = ', nn_trd_trc
793         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
794         WRITE(numout,*) '                You should reconsider this choice.                        ' 
795         WRITE(numout,*) 
796         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
797         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
798         nstop = nstop + 1
799      ENDIF
800
801      ! * Debugging information *
802      IF( lldebug ) THEN
803         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
804         WRITE(numout,*) '               ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
805      ENDIF
806
807      IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
808         WRITE(numout,cform_err)
809         WRITE(numout,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
810         WRITE(numout,*) '                advection and window averaged diagnostics of ML trends.   '
811         WRITE(numout,*) '                WHY? Everything in trdmxl_trc is coded for leap-frog, and '
812         WRITE(numout,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
813         WRITE(numout,*) '                that MUSCL is leap-frog for active tracers T/S).          '
814         WRITE(numout,*) '                In particuliar, entrainment trend would be FALSE. However '
815         WRITE(numout,*) '                this residual is correct for instantaneous ML diagnostics.'
816         WRITE(numout,*) 
817         nstop = nstop + 1
818      ENDIF
819
820      ! I.2 Initialize arrays to zero or read a restart file
821      ! ----------------------------------------------------
822      nmoymltrd   = 0
823
824      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
825      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
826      tmlradm_trc        (:,:,:)   = 0.e0
827
828      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
829      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
830
831      IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
832         CALL trd_mxl_trc_rst_read
833      ELSE
834         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
835         tmlbn_trc          (:,:,:)   = 0.e0
836
837         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
838         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
839
840       ENDIF
841
842      icount = 1   ;   ionce  = 1  ! open specifier   
843
844
845      ! I.3 Read control surface from file ctlsurf_idx
846      ! ----------------------------------------------
847      IF( nn_ctls_trc == 1 ) THEN
848         CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
849         READ ( inum ) nbol_trc
850         CLOSE( inum )
851      ENDIF
852
853      ! ======================================================================
854      ! II. netCDF output initialization
855      ! ======================================================================
856
857      ! clmxl = legend root for netCDF output
858      IF( nn_ctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
859         clmxl = 'Mixed Layer '
860      ELSE IF( nn_ctls_trc == 1 ) THEN                              ! control surface = read index from file
861         clmxl = '      Bowl '
862      ELSE IF( nn_ctls_trc >= 2 ) THEN                              ! control surface = model level
863         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
864      ENDIF
865
866      ! II.1 Define frequency of output and means
867      ! -----------------------------------------
868      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
869      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cp time)
870      ENDIF
871#  if defined key_diainstant
872      IF( .NOT. ln_trdmxl_trc_instant ) THEN
873         STOP 'trd_mxl_trc : this was never checked. Comment this line to proceed...'
874      ENDIF
875      zsto = nn_trd_trc * rdt
876      clop = "inst("//TRIM(clop)//")"
877#  else
878      IF( ln_trdmxl_trc_instant ) THEN
879         zsto = rdt                                               ! inst. diags : we use IOIPSL time averaging
880      ELSE
881         zsto = nn_trd_trc * rdt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
882      ENDIF
883      clop = "ave("//TRIM(clop)//")"
884#  endif
885      zout = nn_trd_trc * rdt
886      iiter = ( nittrc000 - 1 ) / nn_dttrc
887
888      IF(lwp) WRITE (numout,*) '                netCDF initialization'
889
890      ! II.2 Compute julian date from starting date of the run
891      ! ------------------------------------------------------
892      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
893      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
894      IF(lwp) WRITE(numout,*)' ' 
895      IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000               &
896           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
897           &   ,'Julian day : ', zjulian
898
899      ! II.3 Define the T grid trend file (nidtrd)
900      ! ------------------------------------------
901
902      !-- Define long and short names for the NetCDF output variables
903      !       ==> choose them according to trdmxl_trc_oce.F90 <==
904
905      ctrd_trc(jpmxl_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmxl_trc_xad    ,2) = "_xad"
906      ctrd_trc(jpmxl_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmxl_trc_yad    ,2) = "_yad"
907      ctrd_trc(jpmxl_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmxl_trc_zad    ,2) = "_zad"
908      ctrd_trc(jpmxl_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmxl_trc_ldf    ,2) = "_ldf"
909      ctrd_trc(jpmxl_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmxl_trc_zdf    ,2) = "_zdf"
910      ctrd_trc(jpmxl_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmxl_trc_bbl    ,2) = "_bbl"
911      ctrd_trc(jpmxl_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmxl_trc_dmp    ,2) = "_dmp"
912      ctrd_trc(jpmxl_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmxl_trc_sbc    ,2) = "_sbc"
913      ctrd_trc(jpmxl_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmxl_trc_sms    ,2) = "_sms"
914      ctrd_trc(jpmxl_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radb   ,2) = "_radb"
915      ctrd_trc(jpmxl_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radn   ,2) = "_radn"
916      ctrd_trc(jpmxl_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmxl_trc_atf    ,2) = "_atf"
917
918      DO jn = 1, jptra     
919      !-- Create a NetCDF file and enter the define mode
920         IF( ln_trdtrc(jn) ) THEN
921            csuff="ML_"//ctrcnm(jn)
922            CALL dia_nam( clhstnam, nn_trd_trc, csuff )
923            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
924               &        1, jpi, 1, jpj, iiter, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
925     
926            !-- Define the ML depth variable
927            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
928               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
929
930         ENDIF
931      END DO
932
933      !-- Define physical units
934      IF( rn_ucf_trc == 1. ) THEN
935         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit
936      ELSEIF ( rn_ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
937         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf
938      ELSE
939         cltrcu = "unknown?"
940      ENDIF
941
942      !-- Define miscellaneous passive tracer mixed-layer variables
943      IF( jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
944         STOP 'Error : jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb'  ! see below
945      ENDIF
946
947      DO jn = 1, jptra
948         !
949         IF( ln_trdtrc(jn) ) THEN
950            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
951            CALL histdef(nidtrd(jn), trim(clvar),           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
952              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
953            CALL histdef(nidtrd(jn), trim(clvar)//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
954              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
955            CALL histdef(nidtrd(jn), trim(clvar)//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
956              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
957         
958            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmxl_trc_atf
959               CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
960                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
961            END DO                                                                         ! if zsto=rdt above
962         
963            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), & 
964              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
965         
966            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1),   & 
967              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
968         !
969         ENDIF
970      END DO
971
972      !-- Leave IOIPSL/NetCDF define mode
973      DO jn = 1, jptra
974         IF( ln_trdtrc(jn) )  CALL histend( nidtrd(jn), snc4set )
975      END DO
976
977      IF(lwp) WRITE(numout,*)
978
979   END SUBROUTINE trd_mxl_trc_init
980
981#else
982   !!----------------------------------------------------------------------
983   !!   Default option :                                       Empty module
984   !!----------------------------------------------------------------------
985CONTAINS
986   SUBROUTINE trd_mxl_trc( kt )                                   ! Empty routine
987      INTEGER, INTENT( in) ::   kt
988      WRITE(*,*) 'trd_mxl_trc: You should not have seen this print! error?', kt
989   END SUBROUTINE trd_mxl_trc
990   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn )
991      INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
992      CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
993      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmxl            ! passive trc trend
994      WRITE(*,*) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
995      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
996      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
997      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
998   END SUBROUTINE trd_mxl_trc_zint
999   SUBROUTINE trd_mxl_trc_init                                    ! Empty routine
1000      WRITE(*,*) 'trd_mxl_trc_init: You should not have seen this print! error?'
1001   END SUBROUTINE trd_mxl_trc_init
1002#endif
1003
1004   !!======================================================================
1005END MODULE trdmxl_trc
Note: See TracBrowser for help on using the repository browser.