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/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trdmxl_trc.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496

  • Property svn:keywords set to Id
File size: 51.5 KB
Line 
1MODULE trdmxl_trc
2   !!======================================================================
3   !!                       ***  MODULE  trdmxl_trc  ***
4   !! Ocean diagnostics:  mixed layer passive tracer trends
5   !!======================================================================
6   !! History :  9.0  !  06-08  (C. Deltel)  Original code (from trdmxl.F90)
7   !!                 !  07-04  (C. Deltel)  Bug fix : add trcrad trends
8   !!                 !  07-06  (C. Deltel)  key_gyre : do not call lbc_lnk
9   !!----------------------------------------------------------------------
10#if   defined key_top   &&   defined key_trdmxl_trc
11   !!----------------------------------------------------------------------
12   !!   'key_trdmxl_trc'                      mixed layer trend diagnostics
13   !!----------------------------------------------------------------------
14   !!   trd_mxl_trc      : passive tracer cumulated trends averaged over ML
15   !!   trd_mxl_trc_zint : passive tracer trends vertical integration
16   !!   trd_mxl_trc_init : initialization step
17   !!----------------------------------------------------------------------
18   USE trc               ! tracer definitions (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#if defined key_mpi3
287                  CALL lbc_lnk_nc_multi( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1. )        ! lateral boundary conditions
288#else
289                  CALL lbc_lnk( 'trdmxl_trc', tmltrd_trc(:,:,jl,jn), 'T', 1. )        ! lateral boundary conditions
290#endif
291               END DO
292            ENDIF
293         END DO
294!!gm      ENDIF
295     
296      ! ======================================================================
297      ! II. Cumulate the trends over the analysis window
298      ! ======================================================================
299
300      ztmltrd2(:,:,:,:) = 0.e0   ;   ztmltot2(:,:,:)   = 0.e0     ! <<< reset arrays to zero
301      ztmlres2(:,:,:)   = 0.e0   ;   ztmlatf2(:,:,:)   = 0.e0
302      ztmlrad2(:,:,:)   = 0.e0
303
304      ! II.1 Set before values of vertically averages passive tracers
305      ! -------------------------------------------------------------
306      IF( kt > nittrc000 ) THEN
307         DO jn = 1, jptra
308            IF( ln_trdtrc(jn) ) THEN
309               tmlb_trc   (:,:,jn) = tml_trc   (:,:,jn)
310               tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_atf,jn)
311               tmlradn_trc(:,:,jn) = tmltrd_trc(:,:,jpmxl_trc_radb,jn)
312            ENDIF
313         END DO
314      ENDIF
315
316      ! II.2 Vertically averaged passive tracers
317      ! ----------------------------------------
318      tml_trc(:,:,:) = 0.e0
319      DO jk = 1, jpktrd_trc ! - 1 ???
320         DO jn = 1, jptra
321            IF( ln_trdtrc(jn) ) &
322               tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * tr(:,:,jk,jn,Kmm)
323         END DO
324      END DO
325
326      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window   
327      ! ------------------------------------------------------------------------
328      IF( kt == nittrc000 + 1 ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)    ???
329         !
330         DO jn = 1, jptra
331            IF( ln_trdtrc(jn) ) THEN
332               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
333               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
334               
335               tmltrd_csum_ub_trc (:,:,:,jn) = 0.e0   ;   tmltrd_atf_sumb_trc  (:,:,jn) = 0.e0
336               tmltrd_rad_sumb_trc  (:,:,jn) = 0.e0
337            ENDIF
338         END DO
339         
340         rmldbn_trc(:,:) = rmld_trc(:,:)
341         !
342      ENDIF
343
344      ! II.4 Cumulated trends over the analysis period
345      ! ----------------------------------------------
346      !
347      !         [  1rst analysis window ] [     2nd analysis window     ]                       
348      !
349      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
350      !                            ntrd                             2*ntrd       etc.
351      !     1      2     3     4    =5 e.g.                          =10
352      !
353      IF( ( kt >= 2 ).OR.( ln_rsttr ) ) THEN                        ! ???
354         !
355         nmoymltrd = nmoymltrd + 1
356
357
358         ! ... Cumulate over BOTH physical contributions AND over time steps
359         DO jn = 1, jptra
360            IF( ln_trdtrc(jn) ) THEN
361               DO jl = 1, jpltrd_trc
362                  tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn)
363               END DO
364            ENDIF
365         END DO
366
367         DO jn = 1, jptra
368            IF( ln_trdtrc(jn) ) THEN
369               ! ... Special handling of the Asselin trend
370               tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn)
371               tmlradm_trc(:,:,jn) = tmlradm_trc(:,:,jn) + tmlradn_trc(:,:,jn)
372
373               ! ... Trends associated with the time mean of the ML passive tracers
374               tmltrd_sum_trc    (:,:,:,jn) = tmltrd_sum_trc    (:,:,:,jn) + tmltrd_trc    (:,:,:,jn)
375               tmltrd_csum_ln_trc(:,:,:,jn) = tmltrd_csum_ln_trc(:,:,:,jn) + tmltrd_sum_trc(:,:,:,jn)
376               tml_sum_trc       (:,:,jn)   = tml_sum_trc       (:,:,jn)   + tml_trc       (:,:,jn)
377            ENDIF
378         ENDDO
379
380         rmld_sum_trc      (:,:)     = rmld_sum_trc      (:,:)     + rmld_trc      (:,:)
381         !
382      ENDIF
383
384      ! ======================================================================
385      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
386      ! ======================================================================
387
388      ! Convert to appropriate physical units
389      tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc
390
391      itmod = kt - nittrc000 + 1
392      it    = kt
393
394      MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN           ! nitend MUST be multiple of nn_trd_trc
395         !
396         ztmltot (:,:,:) = 0.e0                                   ! reset arrays to zero
397         ztmlres (:,:,:) = 0.e0
398         ztmltot2(:,:,:) = 0.e0
399         ztmlres2(:,:,:) = 0.e0
400     
401         zfn  = FLOAT( nmoymltrd )    ;    zfn2 = zfn * zfn
402         
403         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
404         ! -------------------------------------------------------------
405
406         DO jn = 1, jptra
407            IF( ln_trdtrc(jn) ) THEN
408               !-- Compute total trends
409               IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN  ! EULER-FORWARD schemes
410                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rn_Dt
411               ELSE                                                                     ! LEAP-FROG schemes
412                  ztmltot(:,:,jn) =  ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) + tmlb_trc(:,:,jn) - tmlbb_trc(:,:,jn))/(2.*rn_Dt)
413               ENDIF
414               
415               !-- Compute residuals
416               ztmlres(:,:,jn) = ztmltot(:,:,jn) - ( tmltrdm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn) &
417                  &                                                 - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)  )
418               
419               !-- Diagnose Asselin trend over the analysis window
420               ztmlatf(:,:,jn) = tmlatfm_trc(:,:,jn) - tmlatfn_trc(:,:,jn) + tmlatfb_trc(:,:,jn)
421               ztmlrad(:,:,jn) = tmlradm_trc(:,:,jn) - tmlradn_trc(:,:,jn) + tmlradb_trc(:,:,jn)
422               
423         !-- Lateral boundary conditions
424               IF ( cn_cfg .NE. 'gyre' ) THEN
425#if defined key_mpi3
426                  CALL lbc_lnk_nc_multi( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1. , ztmlres(:,:,jn) , 'T', 1., &
427                     &                ztmlatf(:,:,jn) , 'T', 1. , ztmlrad(:,:,jn) , 'T', 1. )
428#else
429                  CALL lbc_lnk_multi( 'trdmxl_trc', ztmltot(:,:,jn) , 'T', 1. , ztmlres(:,:,jn) , 'T', 1., &
430                     &                ztmlatf(:,:,jn) , 'T', 1. , ztmlrad(:,:,jn) , 'T', 1. )
431#endif
432               ENDIF
433
434
435#if defined key_diainstant
436               CALL ctl_stop( 'STOP', 'tmltrd_trc : key_diainstant was never checked within trdmxl. Comment this to proceed.' )
437#endif
438            ENDIF
439         END DO
440
441         ! III.2 Prepare fields for output ("mean" diagnostics)
442         ! ----------------------------------------------------
443         
444         !-- Update the ML depth time sum (to build the Leap-Frog time mean)
445         rmld_sum_trc(:,:) = rmldbn_trc(:,:) + 2 * ( rmld_sum_trc(:,:) - rmld_trc(:,:) ) + rmld_trc(:,:)
446
447               !-- Compute passive tracer total trends
448         DO jn = 1, jptra
449            IF( ln_trdtrc(jn) ) THEN
450               tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn)
451               ztmltot2   (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) /  ( 2.*rn_Dt )    ! now tracer unit is /sec
452            ENDIF
453         END DO
454
455         !-- Compute passive tracer residuals
456         DO jn = 1, jptra
457            IF( ln_trdtrc(jn) ) THEN
458               !
459               DO jl = 1, jpltrd_trc
460                  ztmltrd2(:,:,jl,jn) = tmltrd_csum_ub_trc(:,:,jl,jn) + tmltrd_csum_ln_trc(:,:,jl,jn)
461               END DO
462               
463               ztmltrdm2(:,:,jn) = 0.e0
464               DO jl = 1, jpltrd_trc
465                  ztmltrdm2(:,:,jn) = ztmltrdm2(:,:,jn) + ztmltrd2(:,:,jl,jn)
466               END DO
467               
468               ztmlres2(:,:,jn) =  ztmltot2(:,:,jn)  -       &
469                  & ( ztmltrdm2(:,:,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) + tmltrd_atf_sumb_trc(:,:,jn) &
470                  &                     - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) + tmltrd_rad_sumb_trc(:,:,jn) )
471               !
472
473               !-- Diagnose Asselin trend over the analysis window
474               ztmlatf2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_atf ,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn) &
475                  &                                               + tmltrd_atf_sumb_trc(:,:,jn)
476               ztmlrad2(:,:,jn) = ztmltrd2(:,:,jpmxl_trc_radb,jn) - tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn) &
477                  &                                               + tmltrd_rad_sumb_trc(:,:,jn)
478
479         !-- Lateral boundary conditions
480               IF ( cn_cfg .NE. 'gyre' ) THEN            ! other than GYRE configuration   
481#if defined key_mpi3
482                  CALL lbc_lnk_nc_multi( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1., ztmlres2(:,:,jn), 'T', 1. )
483#else
484                  CALL lbc_lnk_multi( 'trdmxl_trc', ztmltot2(:,:,jn), 'T', 1., ztmlres2(:,:,jn), 'T', 1. )
485#endif
486                  DO jl = 1, jpltrd_trc
487#if defined key_mpi3
488                     CALL lbc_lnk_nc_multi( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
489#else
490                     CALL lbc_lnk( 'trdmxl_trc', ztmltrd2(:,:,jl,jn), 'T', 1. )       ! will be output in the NetCDF trends file
491#endif
492                  END DO
493               ENDIF
494
495            ENDIF
496         END DO
497
498         ! * Debugging information *
499         IF( lldebug ) THEN
500            !
501            WRITE(numout,*) 'trd_mxl_trc : write trends in the Mixed Layer for debugging process:'
502            WRITE(numout,*) '~~~~~~~~~~~  '
503            WRITE(numout,*)
504            WRITE(numout,*) 'TRC kt = ', kt, '    nmoymltrd = ', nmoymltrd
505
506            DO jn = 1, jptra
507
508               IF( ln_trdtrc(jn) ) THEN
509                  WRITE(numout, *)
510                  WRITE(numout, *) '>>>>>>>>>>>>>>>>>>  TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<'
511                 
512                  WRITE(numout, *)
513                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres     : ', SUM2D(ztmlres(:,:,jn))
514                  !CD??? PREVOIR: z2d = ztmlres(:,:,jn)   ;   CALL prt_ctl(tab2d_1=z2d, clinfo1=' ztmlres   -   : ')
515                 
516                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres): ', SUM2D(ABS(ztmlres(:,:,jn)))
517                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres is computed from ------------- '
518                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot    : ', SUM2D(+ztmltot    (:,:,jn))
519                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrdm_trc: ', SUM2D(+tmltrdm_trc(:,:,jn))
520                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlatfn_trc: ', SUM2D(-tmlatfn_trc(:,:,jn))
521                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlatfb_trc: ', SUM2D(+tmlatfb_trc(:,:,jn))
522                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlradn_trc: ', SUM2D(-tmlradn_trc(:,:,jn))
523                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlradb_trc: ', SUM2D(+tmlradb_trc(:,:,jn))
524                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
525                 
526                  WRITE(numout, *)
527                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmlres2    : ', SUM2D(ztmlres2(:,:,jn))
528                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ABS(ztmlres2):', SUM2D(ABS(ztmlres2(:,:,jn)))
529                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmlres2 is computed from ------------ '
530                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltot2            : ', SUM2D(+ztmltot2(:,:,jn))
531                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +ztmltrdm2           : ', SUM2D(+ztmltrdm2(:,:,jn)) 
532                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_atf,jn)) 
533                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_atf_sumb_trc : ', SUM2D(+tmltrd_atf_sumb_trc(:,:,jn))
534                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmltrd_sum_trc      : ', SUM2D(-tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn))
535                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmltrd_rad_sumb_trc : ', SUM2D(+tmltrd_rad_sumb_trc(:,:,jn) )
536                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
537                 
538                  WRITE(numout, *)
539                  WRITE(numout,98) 'TRC jn =', jn, ' SUM ztmltot     : ', SUM2D(ztmltot    (:,:,jn))
540                  WRITE(numout, '(3x,a)') ' -->>>------------------- ztmltot is computed from ------------- '
541                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tml_trc    : ', SUM2D(tml_trc    (:,:,jn))
542                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbn_trc  : ', SUM2D(tmlbn_trc  (:,:,jn))
543                  WRITE(numout,98) 'TRC jn =', jn, ' SUM +tmlb_trc   : ', SUM2D(tmlb_trc   (:,:,jn))
544                  WRITE(numout,98) 'TRC jn =', jn, ' SUM -tmlbb_trc  : ', SUM2D(tmlbb_trc  (:,:,jn))
545                  WRITE(numout, '(3x,a)') ' --<<<----------------------------------------------------------- '
546                 
547                  WRITE(numout, *)
548                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmltrdm_trc : ', SUM2D(tmltrdm_trc(:,:,jn))
549                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfb_trc : ', SUM2D(tmlatfb_trc(:,:,jn))
550                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlatfn_trc : ', SUM2D(tmlatfn_trc(:,:,jn))
551                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradb_trc : ', SUM2D(tmlradb_trc(:,:,jn))
552                  WRITE(numout,98) 'TRC jn =', jn, ' SUM tmlradn_trc : ', SUM2D(tmlradn_trc(:,:,jn))
553                 
554                  WRITE(numout, *)
555                  DO jl = 1, jpltrd_trc
556                     WRITE(numout,97) 'TRC jn =', jn, ' TREND INDEX jpmxl_trc_xxx = ', jl, &
557                        & ' SUM tmltrd_trc : ', SUM2D(tmltrd_trc(:,:,jl,jn))
558                  END DO
559               
560                  WRITE(numout,*) 
561                  WRITE(numout,*) ' *********************** ZTMLRES, ZTMLRES2 *********************** '
562                  WRITE(numout,*)
563                  WRITE(numout,*) 'TRC ztmlres (jpi/2,jpj/2,:) : ', ztmlres (jpi/2,jpj/2,jn)
564                  WRITE(numout,*)
565                  WRITE(numout,*) 'TRC ztmlres2(jpi/2,jpj/2,:) : ', ztmlres2(jpi/2,jpj/2,jn)
566                 
567                  WRITE(numout,*) 
568                  WRITE(numout,*) ' *********************** ZTMLRES *********************** '
569                  WRITE(numout,*)
570                 
571                  WRITE(numout,*) '...................................................'
572                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres (1:10,1:5,jn) : '
573                  DO jj = 5, 1, -1
574                     WRITE(numout,99) jj, ( ztmlres (ji,jj,jn), ji=1,10 )
575                  END DO
576                 
577                  WRITE(numout,*) 
578                  WRITE(numout,*) ' *********************** ZTMLRES2 *********************** '
579                  WRITE(numout,*)
580                 
581                  WRITE(numout,*) '...................................................'
582                  WRITE(numout,*) 'TRC jn =', jn, ' ztmlres2 (1:10,1:5,jn) : '
583                  DO jj = 5, 1, -1
584                     WRITE(numout,99) jj, ( ztmlres2 (ji,jj,jn), ji=1,10 )
585                  END DO
586                  !
587               ENDIF
588               !
589            END DO
590
591
59297            FORMAT(a10, i3, 2x, a30, i3, a20, 2x, g20.10)
59398            FORMAT(a10, i3, 2x, a30, 2x, g20.10)
59499            FORMAT('TRC jj =', i3,' : ', 10(g10.3,2x))
595              WRITE(numout,*)
596            !
597         ENDIF
598
599         ! III.3 Time evolution array swap
600         ! -------------------------------
601         ! ML depth
602         rmldbn_trc(:,:)   = rmld_trc(:,:)
603         rmld_sum_trc(:,:)     = rmld_sum_trc(:,:)     /      (2*zfn)  ! similar to tml_sum and sml_sum
604         DO jn = 1, jptra
605            IF( ln_trdtrc(jn) ) THEN       
606               ! For passive tracer instantaneous diagnostics
607               tmlbb_trc  (:,:,jn) = tmlb_trc   (:,:,jn)   ;   tmlbn_trc  (:,:,jn) = tml_trc    (:,:,jn)
608               tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn)   ;   tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn)
609               
610               ! For passive tracer mean diagnostics
611               tmltrd_csum_ub_trc (:,:,:,jn) = zfn * tmltrd_sum_trc(:,:,:,jn) - tmltrd_csum_ln_trc(:,:,:,jn)
612               tml_sumb_trc       (:,:,jn)   = tml_sum_trc(:,:,jn)
613               tmltrd_atf_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_atf ,jn)
614               tmltrd_rad_sumb_trc(:,:,jn)   = tmltrd_sum_trc(:,:,jpmxl_trc_radb,jn)
615               
616               
617               ! III.4 Convert to appropriate physical units
618               ! -------------------------------------------
619               ztmltot     (:,:,jn)   = ztmltot     (:,:,jn)   * rn_ucf_trc/zfn   ! instant diags
620               ztmlres     (:,:,jn)   = ztmlres     (:,:,jn)   * rn_ucf_trc/zfn
621               ztmlatf     (:,:,jn)   = ztmlatf     (:,:,jn)   * rn_ucf_trc/zfn
622               ztmlrad     (:,:,jn)   = ztmlrad     (:,:,jn)   * rn_ucf_trc/zfn
623               tml_sum_trc (:,:,jn)   = tml_sum_trc (:,:,jn)   /      (2*zfn)  ! mean diags
624               ztmltot2    (:,:,jn)   = ztmltot2    (:,:,jn)   * rn_ucf_trc/zfn2
625               ztmltrd2    (:,:,:,jn) = ztmltrd2    (:,:,:,jn) * rn_ucf_trc/zfn2
626               ztmlatf2    (:,:,jn)   = ztmlatf2    (:,:,jn)   * rn_ucf_trc/zfn2
627               ztmlrad2    (:,:,jn)   = ztmlrad2    (:,:,jn)   * rn_ucf_trc/zfn2
628               ztmlres2    (:,:,jn)   = ztmlres2    (:,:,jn)   * rn_ucf_trc/zfn2
629            ENDIF
630         END DO
631         !
632      ENDIF MODULO_NTRD
633
634      ! ======================================================================
635      ! IV. Write trends in the NetCDF file
636      ! ======================================================================
637
638      ! IV.1 Code for IOIPSL/NetCDF output
639      ! ----------------------------------
640
641      IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN
642         WRITE(numout,*) ' '
643         WRITE(numout,*) 'trd_mxl_trc : write passive tracer trends in the NetCDF file :'
644         WRITE(numout,*) '~~~~~~~~~~~ '
645         WRITE(numout,*) '          ', trim(clhstnam), ' at kt = ', kt
646         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
647         WRITE(numout,*) ' '
648      ENDIF
649         
650      NETCDF_OUTPUT : IF( ln_trdmxl_trc_instant ) THEN            ! <<< write the trends for passive tracer instant. diags
651         !
652
653         DO jn = 1, jptra
654            !
655            IF( ln_trdtrc(jn) ) THEN
656               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 )
657               !-- Output the fields
658               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
659               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_trc(:,:,jn), ndimtrd1, ndextrd1 ) 
660               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it, ztmltot(:,:,jn), ndimtrd1, ndextrd1 ) 
661               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it, ztmlres(:,:,jn), ndimtrd1, ndextrd1 ) 
662           
663               DO jl = 1, jpltrd_trc - 2
664                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),             &
665                    &          it, tmltrd_trc(:,:,jl,jn), ndimtrd1, ndextrd1 )
666               END DO
667
668               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),    &  ! now trcrad    : jpltrd_trc - 1
669                    &          it, ztmlrad(:,:,jn), ndimtrd1, ndextrd1 )
670
671               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),     &  ! now Asselin   : jpltrd_trc
672                    &          it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 )
673                     
674            ENDIF
675         END DO
676
677         IF( kt == nitend ) THEN
678            DO jn = 1, jptra
679               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
680            END DO
681         ENDIF
682
683      ELSE                                                        ! <<< write the trends for passive tracer mean diagnostics
684         
685         DO jn = 1, jptra
686            !
687            IF( ln_trdtrc(jn) ) THEN
688               CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 
689               !-- Output the fields
690               clvar = trim(ctrcnm(jn))//"ml"                        ! e.g. detml, zooml, nh4ml, etc.
691
692               CALL histwrite( nidtrd(jn), trim(clvar)         , it, tml_sum_trc(:,:,jn), ndimtrd1, ndextrd1 )
693               CALL histwrite( nidtrd(jn), trim(clvar)//"_tot" , it,    ztmltot2(:,:,jn), ndimtrd1, ndextrd1 ) 
694               CALL histwrite( nidtrd(jn), trim(clvar)//"_res" , it,    ztmlres2(:,:,jn), ndimtrd1, ndextrd1 ) 
695
696               DO jl = 1, jpltrd_trc - 2
697                  CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)),           &
698                    &          it, ztmltrd2(:,:,jl,jn), ndimtrd1, ndextrd1 )
699               END DO
700           
701               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)),   &  ! now trcrad    : jpltrd_trc - 1
702                 &          it, ztmlrad2(:,:,jn), ndimtrd1, ndextrd1 )
703
704               CALL histwrite( nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)),    &  ! now Asselin   : jpltrd_trc
705                 &          it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 )
706
707            ENDIF 
708            !
709         END DO
710         IF( kt == nitend ) THEN
711            DO jn = 1, jptra
712               IF( ln_trdtrc(jn) )  CALL histclo( nidtrd(jn) )
713            END DO
714         ENDIF
715
716         !
717      ENDIF NETCDF_OUTPUT
718         
719      ! Compute the control surface (for next time step) : flag = on
720      icount = 1
721
722      IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN
723         !
724         ! Reset cumulative arrays to zero
725         ! -------------------------------         
726         nmoymltrd = 0
727         tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
728         tmlradm_trc        (:,:,:)   = 0.e0   ;   tml_sum_trc        (:,:,:)   = 0.e0
729         tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0
730         rmld_sum_trc       (:,:)     = 0.e0
731         !
732      ENDIF
733
734      ! ======================================================================
735      ! V. Write restart file
736      ! ======================================================================
737
738      IF( lrst_trc )   CALL trd_mxl_trc_rst_write( kt )  ! this must be after the array swap above (III.3)
739      !
740   END SUBROUTINE trd_mxl_trc
741
742   REAL FUNCTION sum2d( ztab )
743      !!----------------------------------------------------------------------
744      !! CD ??? prevoir d'utiliser plutot prtctl
745      !!----------------------------------------------------------------------
746      REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) ::  ztab     
747      !!----------------------------------------------------------------------
748      sum2d = SUM( ztab(2:jpi-1,2:jpj-1) )
749   END FUNCTION sum2d
750
751
752   SUBROUTINE trd_mxl_trc_init
753      !!----------------------------------------------------------------------
754      !!                  ***  ROUTINE trd_mxl_init  ***
755      !!
756      !! ** Purpose :   computation of vertically integrated T and S budgets
757      !!      from ocean surface down to control surface (NetCDF output)
758      !!
759      !! ** Method/usage :
760      !!
761      !!----------------------------------------------------------------------
762      INTEGER :: inum   ! logical unit
763      INTEGER :: ilseq, jl, jn, iiter
764      REAL(wp) ::   zjulian, zsto, zout
765      CHARACTER (LEN=40) ::   clop
766      CHARACTER (LEN=15) ::   csuff
767      CHARACTER (LEN=12) ::   clmxl
768      CHARACTER (LEN=16) ::   cltrcu
769      CHARACTER (LEN=10) ::   clvar
770
771      !!----------------------------------------------------------------------
772
773      ! ======================================================================
774      ! I. initialization
775      ! ======================================================================
776
777      IF(lwp) THEN
778         WRITE(numout,*)
779         WRITE(numout,*) ' trd_mxl_trc_init : Mixed-layer trends for passive tracers                '
780         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
781         WRITE(numout,*)
782      ENDIF
783
784     
785      ! I.1 Check consistency of user defined preferences
786      ! -------------------------------------------------
787
788      IF( ( lk_trdmxl_trc ) .AND. ( MOD( nitend-nittrc000+1, nn_trd_trc ) /= 0 ) ) THEN
789         WRITE(ctmp1,*) '                Your nitend parameter, nitend = ', nitend
790         WRITE(ctmp2,*) '                is no multiple of the trends diagnostics frequency        '
791         WRITE(ctmp3,*) '                          you defined, nn_trd_trc   = ', nn_trd_trc
792         WRITE(ctmp4,*) '                This will not allow you to restart from this simulation.  '
793         WRITE(ctmp5,*) '                You should reconsider this choice.                        ' 
794         WRITE(ctmp6,*) 
795         WRITE(ctmp7,*) '                N.B. the nitend parameter is also constrained to be a     '
796         WRITE(ctmp8,*) '                multiple of the sea-ice frequency parameter (typically 5) '
797         CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4, ctmp5, ctmp6, ctmp7, ctmp8 )
798      ENDIF
799
800      ! * Debugging information *
801      IF( lldebug ) THEN
802         WRITE(numout,*) '               ln_trcadv_muscl = '      , ln_trcadv_muscl
803         WRITE(numout,*) '               ln_trdmxl_trc_instant = ', ln_trdmxl_trc_instant
804      ENDIF
805
806      IF( ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) .AND. .NOT. ln_trdmxl_trc_instant ) THEN
807         WRITE(ctmp1,*) '                Currently, you can NOT use simultaneously tracer MUSCL    '
808         WRITE(ctmp2,*) '                advection and window averaged diagnostics of ML trends.   '
809         WRITE(ctmp3,*) '                WHY? Everything in trdmxl_trc is coded for leap-frog, and '
810         WRITE(ctmp4,*) '                MUSCL scheme is Euler forward for passive tracers (note   '
811         WRITE(ctmp5,*) '                that MUSCL is leap-frog for active tracers T/S).          '
812         WRITE(ctmp6,*) '                In particuliar, entrainment trend would be FALSE. However '
813         WRITE(ctmp7,*) '                this residual is correct for instantaneous ML diagnostics.'
814         CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4, ctmp5, ctmp6, ctmp7 )
815      ENDIF
816
817      ! I.2 Initialize arrays to zero or read a restart file
818      ! ----------------------------------------------------
819      nmoymltrd   = 0
820
821      rmld_trc           (:,:)     = 0.e0   ;   tml_trc            (:,:,:)   = 0.e0       ! inst.
822      tmltrdm_trc        (:,:,:)   = 0.e0   ;   tmlatfm_trc        (:,:,:)   = 0.e0
823      tmlradm_trc        (:,:,:)   = 0.e0
824
825      tml_sum_trc        (:,:,:)   = 0.e0   ;   tmltrd_sum_trc     (:,:,:,:) = 0.e0       ! mean
826      tmltrd_csum_ln_trc (:,:,:,:) = 0.e0   ;   rmld_sum_trc       (:,:)     = 0.e0
827
828      IF( ln_rsttr .AND. ln_trdmxl_trc_restart ) THEN
829         CALL trd_mxl_trc_rst_read
830      ELSE
831         tmlb_trc           (:,:,:)   = 0.e0   ;   tmlbb_trc          (:,:,:)   = 0.e0     ! inst.
832         tmlbn_trc          (:,:,:)   = 0.e0
833
834         tml_sumb_trc       (:,:,:)   = 0.e0   ;   tmltrd_csum_ub_trc (:,:,:,:) = 0.e0     ! mean
835         tmltrd_atf_sumb_trc(:,:,:)   = 0.e0   ;   tmltrd_rad_sumb_trc(:,:,:)   = 0.e0 
836
837       ENDIF
838
839      icount = 1   ;   ionce  = 1  ! open specifier   
840
841
842      ! I.3 Read control surface from file ctlsurf_idx
843      ! ----------------------------------------------
844      IF( nn_ctls_trc == 1 ) THEN
845         CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
846         READ ( inum ) nbol_trc
847         CLOSE( inum )
848      ENDIF
849
850      ! ======================================================================
851      ! II. netCDF output initialization
852      ! ======================================================================
853
854      ! clmxl = legend root for netCDF output
855      IF( nn_ctls_trc == 0 ) THEN                                   ! control surface = mixed-layer with density criterion
856         clmxl = 'Mixed Layer '
857      ELSE IF( nn_ctls_trc == 1 ) THEN                              ! control surface = read index from file
858         clmxl = '      Bowl '
859      ELSE IF( nn_ctls_trc >= 2 ) THEN                              ! control surface = model level
860         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc
861      ENDIF
862
863      ! II.1 Define frequency of output and means
864      ! -----------------------------------------
865      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
866      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cp time)
867      ENDIF
868#  if defined key_diainstant
869      IF( .NOT. ln_trdmxl_trc_instant ) THEN
870         CALL ctl_stop( 'STOP', 'trd_mxl_trc : this was never checked. Comment this line to proceed...' )
871      ENDIF
872      zsto = nn_trd_trc * rn_Dt
873      clop = "inst("//TRIM(clop)//")"
874#  else
875      IF( ln_trdmxl_trc_instant ) THEN
876         zsto = rn_Dt                                               ! inst. diags : we use IOIPSL time averaging
877      ELSE
878         zsto = nn_trd_trc * rn_Dt                                    ! mean  diags : we DO NOT use any IOIPSL time averaging
879      ENDIF
880      clop = "ave("//TRIM(clop)//")"
881#  endif
882      zout = nn_trd_trc * rn_Dt
883      iiter = nittrc000 - 1
884
885      IF(lwp) WRITE (numout,*) '                netCDF initialization'
886
887      ! II.2 Compute julian date from starting date of the run
888      ! ------------------------------------------------------
889      CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
890      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
891      IF(lwp) WRITE(numout,*)' ' 
892      IF(lwp) WRITE(numout,*)' Date 0 used :', nittrc000               &
893           &   ,' YEAR ', nyear, ' MONTH ', nmonth,' DAY ', nday       &
894           &   ,'Julian day : ', zjulian
895
896      ! II.3 Define the T grid trend file (nidtrd)
897      ! ------------------------------------------
898
899      !-- Define long and short names for the NetCDF output variables
900      !       ==> choose them according to trdmxl_trc_oce.F90 <==
901
902      ctrd_trc(jpmxl_trc_xad    ,1) = " Zonal advection"                 ;   ctrd_trc(jpmxl_trc_xad    ,2) = "_xad"
903      ctrd_trc(jpmxl_trc_yad    ,1) = " Meridional advection"            ;   ctrd_trc(jpmxl_trc_yad    ,2) = "_yad"
904      ctrd_trc(jpmxl_trc_zad    ,1) = " Vertical advection"              ;   ctrd_trc(jpmxl_trc_zad    ,2) = "_zad"
905      ctrd_trc(jpmxl_trc_ldf    ,1) = " Lateral diffusion"               ;   ctrd_trc(jpmxl_trc_ldf    ,2) = "_ldf"
906      ctrd_trc(jpmxl_trc_zdf    ,1) = " Vertical diff. (Kz)"             ;   ctrd_trc(jpmxl_trc_zdf    ,2) = "_zdf"
907      ctrd_trc(jpmxl_trc_bbl    ,1) = " Adv/diff. Bottom boundary layer" ;   ctrd_trc(jpmxl_trc_bbl    ,2) = "_bbl"
908      ctrd_trc(jpmxl_trc_dmp    ,1) = " Tracer damping"                  ;   ctrd_trc(jpmxl_trc_dmp    ,2) = "_dmp"
909      ctrd_trc(jpmxl_trc_sbc    ,1) = " Surface boundary cond."          ;   ctrd_trc(jpmxl_trc_sbc    ,2) = "_sbc"
910      ctrd_trc(jpmxl_trc_sms,    1) = " Sources minus sinks"             ;   ctrd_trc(jpmxl_trc_sms    ,2) = "_sms"
911      ctrd_trc(jpmxl_trc_radb   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radb   ,2) = "_radb"
912      ctrd_trc(jpmxl_trc_radn   ,1) = " Correct negative concentrations" ;   ctrd_trc(jpmxl_trc_radn   ,2) = "_radn"
913      ctrd_trc(jpmxl_trc_atf    ,1) = " Asselin time filter"             ;   ctrd_trc(jpmxl_trc_atf    ,2) = "_atf"
914
915      DO jn = 1, jptra     
916      !-- Create a NetCDF file and enter the define mode
917         IF( ln_trdtrc(jn) ) THEN
918            csuff="ML_"//ctrcnm(jn)
919            CALL dia_nam( clhstnam, nn_trd_trc, csuff )
920            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
921               &        1, jpi, 1, jpj, iiter, zjulian, rn_Dt, nh_t(jn), nidtrd(jn), domain_id=nidom, snc4chunks=snc4set )
922     
923            !-- Define the ML depth variable
924            CALL histdef(nidtrd(jn), "mxl_depth", clmxl//" Mixed Layer Depth", "m",                        &
925               &        jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )
926
927         ENDIF
928      END DO
929
930      !-- Define physical units
931      IF( rn_ucf_trc == 1. ) THEN
932         cltrcu = "(mmole-N/m3)/sec"                              ! all passive tracers have the same unit
933      ELSEIF ( rn_ucf_trc == 3600.*24.) THEN                         ! ??? trop long : seulement (mmole-N/m3)
934         cltrcu = "(mmole-N/m3)/day"                              ! ??? apparait dans les sorties netcdf
935      ELSE
936         cltrcu = "unknown?"
937      ENDIF
938
939      !-- Define miscellaneous passive tracer mixed-layer variables
940      IF( jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb ) THEN
941         CALL ctl_stop( 'STOP', 'Error : jpltrd_trc /= jpmxl_trc_atf .OR.  jpltrd_trc - 1 /= jpmxl_trc_radb' ) ! see below
942      ENDIF
943
944      DO jn = 1, jptra
945         !
946         IF( ln_trdtrc(jn) ) THEN
947            clvar = trim(ctrcnm(jn))//"ml"                           ! e.g. detml, zooml, no3ml, etc.
948            CALL histdef(nidtrd(jn), trim(clvar),           clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ",                         &
949              & "mmole-N/m3", jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
950            CALL histdef(nidtrd(jn), trim(clvar)//"_tot"  , clmxl//" "//trim(ctrcnm(jn))//" Total trend ",                         & 
951              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) 
952            CALL histdef(nidtrd(jn), trim(clvar)//"_res"  , clmxl//" "//trim(ctrcnm(jn))//" dh/dt Entrainment (Resid.)",           & 
953              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
954         
955            DO jl = 1, jpltrd_trc - 2                                ! <== only true if jpltrd_trc == jpmxl_trc_atf
956               CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jl,2)), clmxl//" "//clvar//ctrd_trc(jl,1),                      & 
957                 &    cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
958            END DO                                                                         ! if zsto=rn_Dt above
959         
960            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_radb,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_radb,1), & 
961              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
962         
963            CALL histdef(nidtrd(jn), trim(clvar)//trim(ctrd_trc(jpmxl_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmxl_trc_atf,1),   & 
964              &       cltrcu, jpi, jpj, nh_t(jn), 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
965         !
966         ENDIF
967      END DO
968
969      !-- Leave IOIPSL/NetCDF define mode
970      DO jn = 1, jptra
971         IF( ln_trdtrc(jn) )  CALL histend( nidtrd(jn), snc4set )
972      END DO
973
974      IF(lwp) WRITE(numout,*)
975
976   END SUBROUTINE trd_mxl_trc_init
977
978#else
979   !!----------------------------------------------------------------------
980   !!   Default option :                                       Empty module
981   !!----------------------------------------------------------------------
982CONTAINS
983   SUBROUTINE trd_mxl_trc( kt, Kmm )                                   ! Empty routine
984      INTEGER, INTENT( in) ::   kt
985      INTEGER, INTENT( in) ::   Kmm            ! time level index
986      WRITE(*,*) 'trd_mxl_trc: You should not have seen this print! error?', kt
987   END SUBROUTINE trd_mxl_trc
988   SUBROUTINE trd_mxl_trc_zint( ptrc_trdmxl, ktrd, ctype, kjn, Kmm )
989      INTEGER               , INTENT( in ) ::  ktrd, kjn              ! ocean trend index and passive tracer rank
990      INTEGER               , INTENT( in ) ::  Kmm                    ! time level index
991      CHARACTER(len=2)      , INTENT( in ) ::  ctype                  ! surface/bottom (2D) or interior (3D) physics
992      REAL, DIMENSION(:,:,:), INTENT( in ) ::  ptrc_trdmxl            ! passive trc trend
993      WRITE(*,*) 'trd_mxl_trc_zint: You should not have seen this print! error?', ptrc_trdmxl(1,1,1)
994      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ctype
995      WRITE(*,*) '  "      "      : You should not have seen this print! error?', ktrd
996      WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
997   END SUBROUTINE trd_mxl_trc_zint
998   SUBROUTINE trd_mxl_trc_init                                    ! Empty routine
999      WRITE(*,*) 'trd_mxl_trc_init: You should not have seen this print! error?'
1000   END SUBROUTINE trd_mxl_trc_init
1001#endif
1002
1003   !!======================================================================
1004END MODULE trdmxl_trc
Note: See TracBrowser for help on using the repository browser.