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/2018/dev_r9838_ENHANCE04_RK3/src/TOP/TRP – NEMO

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/TOP/TRP/trdmxl_trc.F90 @ 9939

Last change on this file since 9939 was 9939, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

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