source: NEMO/trunk/src/TOP/TRP/trdmxl_trc.F90

Last change on this file was 13497, checked in by techene, 2 months ago

re-introduce comments that have been erased by loops transformation see #2525

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