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

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

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

Last change on this file since 11504 was 11504, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Strip out all references to nn_dttrc
and the trcsub.F90 module. Notes:

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