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/2020/r12377_ticket2386/src/TOP/TRP – NEMO

source: NEMO/branches/2020/r12377_ticket2386/src/TOP/TRP/trdmxl_trc.F90 @ 12511

Last change on this file since 12511 was 12511, checked in by andmirek, 4 years ago

ticket #2386: update trunk@12493 to have AGRIF sette working

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