source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trdmxl_trc.F90 @ 10420

Last change on this file since 10420 was 10420, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: force STOP when fail to allocate array, see #2133

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