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/UKMO/r8395_mix-lyr_diag/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/UKMO/r8395_mix-lyr_diag/NEMOGCM/NEMO/TOP_SRC/TRP/trdmxl_trc.F90 @ 11290

Last change on this file since 11290 was 11290, checked in by jcastill, 5 years ago

Remove svn keywords

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