source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 3325

Last change on this file since 3325 was 3318, checked in by gm, 9 years ago

Ediag branche: #927 split TRA/DYN trd computation

  • Property svn:keywords set to Id
File size: 46.8 KB
Line 
1MODULE trdmld
2   !!======================================================================
3   !!                       ***  MODULE  trdmld  ***
4   !! Ocean diagnostics:  mixed layer T-S trends
5   !!=====================================================================
6   !! History :       !  95-04  (J. Vialard)    Original code
7   !!                 !  97-02  (E. Guilyardi)  Adaptation global + base cmo
8   !!                 !  99-09  (E. Guilyardi)  Re-writing + netCDF output
9   !!            8.5  !  02-06  (G. Madec)      F90: Free form and module
10   !!            9.0  !  04-08  (C. Talandier)  New trends organization
11   !!                 !  05-05  (C. Deltel)     Diagnose trends of time averaged ML T & S
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   trd_mld          : T and S cumulated trends averaged over the mixed layer
16   !!   trd_mld_zint     : T and S trends vertical integration
17   !!   trd_mld_init     : initialization step
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers variables
20   USE dom_oce         ! ocean space and time domain variables
21   USE trd_oce         ! trends: ocean variables
22   USE trdmld_oce      ! ocean variables trends
23   USE ldftra_oce      ! ocean active tracers lateral physics
24   USE zdf_oce         ! ocean vertical physics
25   USE in_out_manager  ! I/O manager
26   USE phycst          ! Define parameters for the routines
27   USE dianam          ! build the name of file (routine)
28   USE ldfslp          ! iso-neutral slopes
29   USE zdfmxl          ! mixed layer depth
30   USE zdfddm          ! ocean vertical physics: double diffusion
31   USE ioipsl          ! NetCDF library
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33   USE diadimg         ! dimg direct access file format output
34   USE trdmld_rst      ! restart for diagnosing the ML trends
35   USE prtctl          ! Print control
36   USE restart         ! for lrst_oce
37   USE lib_mpp         ! MPP library
38   USE wrk_nemo        ! Memory allocation
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   trd_mld        ! routine called by step.F90
44   PUBLIC   trd_mld_init   ! routine called by opa.F90
45   PUBLIC   trd_mld_zint   ! routine called by tracers routines
46
47   CHARACTER (LEN=40) ::  clhstnam         ! name of the trends NetCDF file
48   INTEGER ::   nh_t, nmoymltrd
49   INTEGER ::   nidtrd
50   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::   ndextrd1
51   INTEGER ::   ndimtrd1                       
52   INTEGER ::   ionce, icount                   
53
54   !! * Substitutions
55#  include "domzgr_substitute.h90"
56#  include "ldftra_substitute.h90"
57#  include "zdfddm_substitute.h90"
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
60   !! $Id$
61   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   INTEGER FUNCTION trd_mld_alloc()
66      !!----------------------------------------------------------------------
67      !!                  ***  ROUTINE trd_mld_alloc  ***
68      !!----------------------------------------------------------------------
69      ALLOCATE( ndextrd1(jpi*jpj) , STAT=trd_mld_alloc )
70      !
71      IF( lk_mpp             )   CALL mpp_sum ( trd_mld_alloc )
72      IF( trd_mld_alloc /= 0 )   CALL ctl_warn('trd_mld_alloc: failed to allocate array ndextrd1')
73   END FUNCTION trd_mld_alloc
74
75
76   SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
77      !!----------------------------------------------------------------------
78      !!                  ***  ROUTINE trd_mld_zint  ***
79      !!
80      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments
81      !!                to the subroutine. This vertical average is performed from ocean
82      !!                surface down to a chosen control surface.
83      !!
84      !! ** Method/usage :
85      !!      The control surface can be either a mixed layer depth (time varying)
86      !!      or a fixed surface (jk level or bowl).
87      !!      Choose control surface with nn_ctls in namelist NAMTRD :
88      !!        nn_ctls = 0  : use mixed layer with density criterion
89      !!        nn_ctls = 1  : read index from file 'ctlsurf_idx'
90      !!        nn_ctls > 1  : use fixed level surface jk = nn_ctls
91      !!      Note: in the remainder of the routine, the volume between the
92      !!            surface and the control surface is called "mixed-layer"
93      !!----------------------------------------------------------------------
94      INTEGER                         , INTENT( in ) ::   ktrd       ! ocean trend index
95      CHARACTER(len=2)                , INTENT( in ) ::   ctype      ! 2D surface/bottom or 3D interior physics
96      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pttrdmld   ! temperature trend
97      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pstrdmld   ! salinity trend
98      !
99      INTEGER ::   ji, jj, jk, isum
100      REAL(wp), POINTER, DIMENSION(:,:)  :: zvlmsk 
101      !!----------------------------------------------------------------------
102
103      CALL wrk_alloc( jpi, jpj, zvlmsk ) 
104
105      ! I. Definition of control surface and associated fields
106      ! ------------------------------------------------------
107      !            ==> only once per time step <==
108
109      IF( icount == 1 ) THEN       
110         !
111         tmltrd(:,:,:) = 0._wp    ;    smltrd(:,:,:) = 0._wp    ! <<< reset trend arrays to zero
112         
113         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
114         IF( nn_ctls == 0 ) THEN       ! * control surface = mixed-layer with density criterion
115            nmld(:,:) = nmln(:,:)    ! array nmln computed in zdfmxl.F90
116         ELSEIF( nn_ctls == 1 ) THEN   ! * control surface = read index from file
117            nmld(:,:) = nbol(:,:)
118         ELSEIF( nn_ctls >= 2 ) THEN   ! * control surface = model level
119            nn_ctls = MIN( nn_ctls, jpktrd - 1 )
120            nmld(:,:) = nn_ctls + 1
121         ENDIF
122
123         ! ... Compute ndextrd1 and ndimtrd1 only once
124         IF( ionce == 1 ) THEN
125            !
126            ! Check of validity : nmld(ji,jj) <= jpktrd
127            isum        = 0
128            zvlmsk(:,:) = 0._wp
129
130            IF( jpktrd < jpk ) THEN
131               DO jj = 1, jpj
132                  DO ji = 1, jpi
133                     IF( nmld(ji,jj) <= jpktrd ) THEN
134                        zvlmsk(ji,jj) = tmask(ji,jj,1)
135                     ELSE
136                        isum = isum + 1
137                        zvlmsk(ji,jj) = 0._wp
138                     ENDIF
139                  END DO
140               END DO
141            ENDIF
142
143            ! Index of ocean points (2D only)
144            IF( isum > 0 ) THEN
145               WRITE(numout,*)' Number of invalid points nmld > jpktrd', isum 
146               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )    ! surface
147            ELSE
148               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )    ! surface
149            ENDIF                               
150
151            ionce = 0                ! no more pass here
152            !
153         END IF
154         
155         ! ... Weights for vertical averaging
156         wkx(:,:,:) = 0._wp
157         DO jk = 1, jpktrd             ! initialize wkx with vertical scale factor in mixed-layer
158            DO jj = 1,jpj
159               DO ji = 1,jpi
160                  IF( jk - nmld(ji,jj) < 0.e0 )   wkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
161               END DO
162            END DO
163         END DO
164         
165         rmld(:,:) = 0._wp               ! compute mixed-layer depth : rmld
166         DO jk = 1, jpktrd
167            rmld(:,:) = rmld(:,:) + wkx(:,:,jk)
168         END DO
169         
170         DO jk = 1, jpktrd             ! compute integration weights
171            wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1._wp, rmld(:,:) )
172         END DO
173
174         icount = 0                    ! <<< flag = off : control surface & integr. weights
175         !                             !     computed only once per time step
176      END IF
177
178      ! II. Vertical integration of trends in the mixed-layer
179      ! -----------------------------------------------------
180
181      SELECT CASE (ctype)
182      CASE ( '3D' )   ! mean T/S trends in the mixed-layer
183         DO jk = 1, jpktrd
184            tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,jk) * wkx(:,:,jk)   ! temperature
185            smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,jk) * wkx(:,:,jk)   ! salinity
186         END DO
187      CASE ( '2D' )   ! forcing at upper boundary of the mixed-layer
188         tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,1) * wkx(:,:,1)        ! non penetrative
189         smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,1) * wkx(:,:,1)           
190      END SELECT
191      !
192      CALL wrk_dealloc( jpi, jpj, zvlmsk ) 
193      !
194   END SUBROUTINE trd_mld_zint
195   
196
197   SUBROUTINE trd_mld( kt )
198      !!----------------------------------------------------------------------
199      !!                  ***  ROUTINE trd_mld  ***
200      !!
201      !! ** Purpose :  Compute and cumulate the mixed layer trends over an analysis
202      !!               period, and write NetCDF (or dimg) outputs.
203      !!
204      !! ** Method/usage :
205      !!          The stored trends can be chosen twofold (according to the ln_trdmld_instant
206      !!          logical namelist variable) :
207      !!          1) to explain the difference between initial and final
208      !!             mixed-layer T & S (where initial and final relate to the
209      !!             current analysis window, defined by nn_trd in the namelist)
210      !!          2) to explain the difference between the current and previous
211      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
212      !!             performed over each analysis window).
213      !!
214      !! ** Consistency check :
215      !!        If the control surface is fixed ( nn_ctls > 1 ), the residual term (dh/dt
216      !!        entrainment) should be zero, at machine accuracy. Note that in the case
217      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
218      !!        over the first two analysis windows (except if restart).
219      !!        N.B. For ORCA2_LIM, use e.g. nn_trd=5, rn_ucf=1., nn_ctls=8
220      !!             for checking residuals.
221      !!             On a NEC-SX5 computer, this typically leads to:
222      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmld_instant=.false.
223      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmld_instant=.true.
224      !!
225      !! ** Action :
226      !!       At each time step, mixed-layer averaged trends are stored in the
227      !!       tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
228      !!       This array is known when trd_mld is called, at the end of the stp subroutine,
229      !!       except for the purely vertical K_z diffusion term, which is embedded in the
230      !!       lateral diffusion trend.
231      !!
232      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
233      !!       from the lateral diffusion trend.
234      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
235      !!       arrays are updated.
236      !!       In III), called only once per analysis window, we compute the total trends,
237      !!       along with the residuals and the Asselin correction terms.
238      !!       In IV), the appropriate trends are written in the trends NetCDF file.
239      !!
240      !! References :  Vialard et al.,2001, JPO.
241      !!----------------------------------------------------------------------
242      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
243      !
244      INTEGER :: ji, jj, jk, jl, ik, it, itmod
245      LOGICAL :: lldebug = .TRUE.
246      REAL(wp) :: zavt, zfn, zfn2
247      !                                              ! z(ts)mltot : dT/dt over the anlysis window (including Asselin)
248      !                                              ! z(ts)mlres : residual = dh/dt entrainment term
249      REAL(wp), POINTER, DIMENSION(:,:  ) ::  ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
250      REAL(wp), POINTER, DIMENSION(:,:  ) ::  ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 
251      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztmltrd2, zsmltrd2   ! only needed for mean diagnostics
252#if defined key_dimgout
253      INTEGER ::  iyear,imon,iday
254      CHARACTER(LEN=80) :: cltext, clmode
255#endif
256      !!----------------------------------------------------------------------
257 
258      CALL wrk_alloc( jpi, jpj,         ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf                        )
259      CALL wrk_alloc( jpi, jpj,         ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 ) 
260      CALL wrk_alloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2                                                               )
261
262      ! ======================================================================
263      ! I. Diagnose the purely vertical (K_z) diffusion trend
264      ! ======================================================================
265
266      ! ... These terms can be estimated by flux computation at the lower boundary of the ML
267      !     (we compute (-1/h) * K_z * d_z( T ) and (-1/h) * K_z * d_z( S ))
268      IF( ln_traldf_iso ) THEN
269         DO jj = 1,jpj
270            DO ji = 1,jpi
271               ik = nmld(ji,jj)
272               zavt = avt(ji,jj,ik)
273               tmltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
274                  &                      * ( tsn(ji,jj,ik-1,jp_tem) - tsn(ji,jj,ik,jp_tem) )         &
275                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
276               zavt = fsavs(ji,jj,ik)
277               smltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
278                  &                      * ( tsn(ji,jj,ik-1,jp_sal) - tsn(ji,jj,ik,jp_sal) )         &
279                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
280            END DO
281         END DO
282
283         ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
284         tmltrd(:,:,jpmld_ldf) = tmltrd(:,:,jpmld_ldf) - tmltrd(:,:,jpmld_zdf)
285         smltrd(:,:,jpmld_ldf) = smltrd(:,:,jpmld_ldf) - smltrd(:,:,jpmld_zdf)
286      END IF
287
288      ! ... Lateral boundary conditions
289      DO jl = 1, jpltrd
290         CALL lbc_lnk( tmltrd(:,:,jl), 'T', 1. )
291         CALL lbc_lnk( smltrd(:,:,jl), 'T', 1. )
292      END DO
293
294      ! ======================================================================
295      ! II. Cumulate the trends over the analysis window
296      ! ======================================================================
297
298      ztmltrd2(:,:,:) = 0.e0   ;    zsmltrd2(:,:,:) = 0.e0  ! <<< reset arrays to zero
299      ztmltot2(:,:)   = 0.e0   ;    zsmltot2(:,:)   = 0.e0
300      ztmlres2(:,:)   = 0.e0   ;    zsmlres2(:,:)   = 0.e0
301      ztmlatf2(:,:)   = 0.e0   ;    zsmlatf2(:,:)   = 0.e0
302
303      ! II.1 Set before values of vertically average T and S
304      ! ----------------------------------------------------
305      IF( kt > nit000 ) THEN
306         !   ... temperature ...                    ... salinity ...
307         tmlb   (:,:) = tml   (:,:)           ; smlb   (:,:) = sml   (:,:)
308         tmlatfn(:,:) = tmltrd(:,:,jpmld_atf) ; smlatfn(:,:) = smltrd(:,:,jpmld_atf)
309      END IF
310
311      ! II.2 Vertically averaged T and S
312      ! --------------------------------
313      tml(:,:) = 0.e0   ;   sml(:,:) = 0.e0
314      DO jk = 1, jpktrd - 1
315         tml(:,:) = tml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_tem)
316         sml(:,:) = sml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_sal)
317      END DO
318
319      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window   
320      ! ------------------------------------------------------------------------
321      IF( kt == 2 ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
322         !
323         !   ... temperature ...                ... salinity ...
324         tmlbb  (:,:) = tmlb   (:,:)   ;   smlbb  (:,:) = smlb   (:,:)
325         tmlbn  (:,:) = tml    (:,:)   ;   smlbn  (:,:) = sml    (:,:)
326         tmlatfb(:,:) = tmlatfn(:,:)   ;   smlatfb(:,:) = smlatfn(:,:)
327         
328         tmltrd_csum_ub (:,:,:) = 0.e0  ;   smltrd_csum_ub (:,:,:) = 0.e0
329         tmltrd_atf_sumb(:,:)   = 0.e0  ;   smltrd_atf_sumb(:,:)   = 0.e0
330
331         rmldbn(:,:) = rmld(:,:)
332
333         IF( ln_ctl ) THEN
334            WRITE(numout,*) '             we reach kt == nit000 + 1 = ', nit000+1
335            CALL prt_ctl(tab2d_1=tmlbb   , clinfo1=' tmlbb   -   : ', mask1=tmask, ovlap=1)
336            CALL prt_ctl(tab2d_1=tmlbn   , clinfo1=' tmlbn   -   : ', mask1=tmask, ovlap=1)
337            CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb -   : ', mask1=tmask, ovlap=1)
338         END IF
339         !
340      END IF
341
342      IF( ( ln_rstart ) .AND. ( kt == nit000 ) .AND. ( ln_ctl ) ) THEN
343         IF( ln_trdmld_instant ) THEN
344            WRITE(numout,*) '             restart from kt == nit000 = ', nit000
345            CALL prt_ctl(tab2d_1=tmlbb   , clinfo1=' tmlbb   -   : ', mask1=tmask, ovlap=1)
346            CALL prt_ctl(tab2d_1=tmlbn   , clinfo1=' tmlbn   -   : ', mask1=tmask, ovlap=1)
347            CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb -   : ', mask1=tmask, ovlap=1)
348         ELSE
349            WRITE(numout,*) '             restart from kt == nit000 = ', nit000
350            CALL prt_ctl(tab2d_1=tmlbn          , clinfo1=' tmlbn           -  : ', mask1=tmask, ovlap=1)
351            CALL prt_ctl(tab2d_1=rmldbn         , clinfo1=' rmldbn          -  : ', mask1=tmask, ovlap=1)
352            CALL prt_ctl(tab2d_1=tml_sumb       , clinfo1=' tml_sumb        -  : ', mask1=tmask, ovlap=1)
353            CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb -  : ', mask1=tmask, ovlap=1)
354            CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub  -  : ', mask1=tmask, ovlap=1, kdim=1)
355         END IF
356      END IF
357
358      ! II.4 Cumulated trends over the analysis period
359      ! ----------------------------------------------
360      !
361      !         [  1rst analysis window ] [     2nd analysis window     ]                       
362      !
363      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
364      !                          nn_trd                           2*nn_trd       etc.
365      !     1      2     3     4    =5 e.g.                          =10
366      !
367      IF( ( kt >= 2 ).OR.( ln_rstart ) ) THEN
368         !
369         nmoymltrd = nmoymltrd + 1
370         
371         ! ... Cumulate over BOTH physical contributions AND over time steps
372         DO jl = 1, jpltrd
373            tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl)
374            smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl)
375         END DO
376
377         ! ... Special handling of the Asselin trend
378         tmlatfm(:,:) = tmlatfm(:,:) + tmlatfn(:,:)
379         smlatfm(:,:) = smlatfm(:,:) + smlatfn(:,:)
380
381         ! ... Trends associated with the time mean of the ML T/S
382         tmltrd_sum    (:,:,:) = tmltrd_sum    (:,:,:) + tmltrd    (:,:,:) ! tem
383         tmltrd_csum_ln(:,:,:) = tmltrd_csum_ln(:,:,:) + tmltrd_sum(:,:,:)
384         tml_sum       (:,:)   = tml_sum       (:,:)   + tml       (:,:)
385         smltrd_sum    (:,:,:) = smltrd_sum    (:,:,:) + smltrd    (:,:,:) ! sal
386         smltrd_csum_ln(:,:,:) = smltrd_csum_ln(:,:,:) + smltrd_sum(:,:,:)
387         sml_sum       (:,:)   = sml_sum       (:,:)   + sml       (:,:)
388         rmld_sum      (:,:)   = rmld_sum      (:,:)   + rmld      (:,:)   ! rmld
389         !
390      END IF
391
392      ! ======================================================================
393      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
394      ! ======================================================================
395
396      ! Convert to appropriate physical units
397      ! N.B. It may be useful to check IOIPSL time averaging with :
398      !      tmltrd (:,:,:) = 1. ; smltrd (:,:,:) = 1.
399      tmltrd(:,:,:) = tmltrd(:,:,:) * rn_ucf   ! (actually needed for 1:jpltrd-1, but trdmld(:,:,jpltrd)
400      smltrd(:,:,:) = smltrd(:,:,:) * rn_ucf   !  is no longer used, and is reset to 0. at next time step)
401     
402      ! define time axis
403      it = kt
404      itmod = kt - nit000 + 1
405
406      MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN        ! nitend MUST be multiple of nn_trd
407         !
408         ztmltot (:,:) = 0.e0   ;   zsmltot (:,:) = 0.e0   ! reset arrays to zero
409         ztmlres (:,:) = 0.e0   ;   zsmlres (:,:) = 0.e0
410         ztmltot2(:,:) = 0.e0   ;   zsmltot2(:,:) = 0.e0
411         ztmlres2(:,:) = 0.e0   ;   zsmlres2(:,:) = 0.e0
412     
413         zfn  = float(nmoymltrd)    ;    zfn2 = zfn * zfn
414         
415         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
416         ! -------------------------------------------------------------
417         
418         !-- Compute total trends
419         ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / ( 2.*rdt )
420         zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / ( 2.*rdt )
421         
422         !-- Compute residuals
423         ztmlres(:,:) = ztmltot(:,:) - ( tmltrdm(:,:) - tmlatfn(:,:) + tmlatfb(:,:) )
424         zsmlres(:,:) = zsmltot(:,:) - ( smltrdm(:,:) - smlatfn(:,:) + smlatfb(:,:) )
425     
426         !-- Diagnose Asselin trend over the analysis window
427         ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
428         zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
429         
430         !-- Lateral boundary conditions
431         !         ... temperature ...                    ... salinity ...
432         CALL lbc_lnk( ztmltot , 'T', 1. )  ;   CALL lbc_lnk( zsmltot , 'T', 1. )
433         CALL lbc_lnk( ztmlres , 'T', 1. )  ;   CALL lbc_lnk( zsmlres , 'T', 1. )
434         CALL lbc_lnk( ztmlatf , 'T', 1. )  ;   CALL lbc_lnk( zsmlatf , 'T', 1. )
435
436#if defined key_diainstant
437         CALL ctl_stop( 'tml_trd : key_diainstant was never checked within trdmld. Comment this to proceed.')
438#endif
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(:,:) = rmldbn(:,:) + 2 * ( rmld_sum(:,:) - rmld(:,:) ) + rmld(:,:)
444
445         !-- Compute temperature total trends
446         tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
447         ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) /  ( 2.*rdt )    ! now in degC/s
448         
449         !-- Compute salinity total trends
450         sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
451         zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) /  ( 2.*rdt )    ! now in psu/s
452         
453         !-- Compute temperature residuals
454         DO jl = 1, jpltrd
455            ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
456         END DO
457
458         ztmltrdm2(:,:) = 0.e0
459         DO jl = 1, jpltrd
460            ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
461         END DO
462
463         ztmlres2(:,:) =  ztmltot2(:,:)  -       &
464              ( ztmltrdm2(:,:) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:) )
465         
466         !-- Compute salinity residuals
467         DO jl = 1, jpltrd
468            zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
469         END DO
470
471         zsmltrdm2(:,:) = 0.
472         DO jl = 1, jpltrd
473            zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
474         END DO
475
476         zsmlres2(:,:) =  zsmltot2(:,:)  -       &
477              ( zsmltrdm2(:,:) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:) )
478         
479         !-- Diagnose Asselin trend over the analysis window
480         ztmlatf2(:,:) = ztmltrd2(:,:,jpmld_atf) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:)
481         zsmlatf2(:,:) = zsmltrd2(:,:,jpmld_atf) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:)
482
483         !-- Lateral boundary conditions
484         !         ... temperature ...                    ... salinity ...
485         CALL lbc_lnk( ztmltot2, 'T', 1. )  ;   CALL lbc_lnk( zsmltot2, 'T', 1. )
486         CALL lbc_lnk( ztmlres2, 'T', 1. )  ;   CALL lbc_lnk( zsmlres2, 'T', 1. )
487         DO jl = 1, jpltrd
488            CALL lbc_lnk( ztmltrd2(:,:,jl), 'T', 1. ) ! \  these will be output
489            CALL lbc_lnk( zsmltrd2(:,:,jl), 'T', 1. ) ! /  in the NetCDF trends file
490         END DO
491         
492         ! III.3 Time evolution array swap
493         ! -------------------------------
494         
495         ! For T/S instantaneous diagnostics
496         !   ... temperature ...               ... salinity ...
497         tmlbb  (:,:) = tmlb   (:,:)  ;   smlbb  (:,:) = smlb   (:,:)
498         tmlbn  (:,:) = tml    (:,:)  ;   smlbn  (:,:) = sml    (:,:)
499         tmlatfb(:,:) = tmlatfn(:,:)  ;   smlatfb(:,:) = smlatfn(:,:)
500
501         ! For T mean diagnostics
502         tmltrd_csum_ub (:,:,:) = zfn * tmltrd_sum(:,:,:) - tmltrd_csum_ln(:,:,:)
503         tml_sumb       (:,:)   = tml_sum(:,:)
504         tmltrd_atf_sumb(:,:)   = tmltrd_sum(:,:,jpmld_atf)
505         
506         ! For S mean diagnostics
507         smltrd_csum_ub (:,:,:) = zfn * smltrd_sum(:,:,:) - smltrd_csum_ln(:,:,:)
508         sml_sumb       (:,:)   = sml_sum(:,:)
509         smltrd_atf_sumb(:,:)   = smltrd_sum(:,:,jpmld_atf)
510         
511         ! ML depth
512         rmldbn         (:,:)   = rmld    (:,:)
513         
514         IF( ln_ctl ) THEN
515            IF( ln_trdmld_instant ) THEN
516               CALL prt_ctl(tab2d_1=tmlbb   , clinfo1=' tmlbb   -   : ', mask1=tmask, ovlap=1)
517               CALL prt_ctl(tab2d_1=tmlbn   , clinfo1=' tmlbn   -   : ', mask1=tmask, ovlap=1)
518               CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb -   : ', mask1=tmask, ovlap=1)
519            ELSE
520               CALL prt_ctl(tab2d_1=tmlbn          , clinfo1=' tmlbn           -  : ', mask1=tmask, ovlap=1)
521               CALL prt_ctl(tab2d_1=rmldbn         , clinfo1=' rmldbn          -  : ', mask1=tmask, ovlap=1)
522               CALL prt_ctl(tab2d_1=tml_sumb       , clinfo1=' tml_sumb        -  : ', mask1=tmask, ovlap=1)
523               CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb -  : ', mask1=tmask, ovlap=1)
524               CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub  -  : ', mask1=tmask, ovlap=1, kdim=1)
525            END IF
526         END IF
527
528         ! III.4 Convert to appropriate physical units
529         ! -------------------------------------------
530
531         !    ... temperature ...                         ... salinity ...
532         ztmltot (:,:)   = ztmltot(:,:)   * rn_ucf/zfn  ; zsmltot (:,:)   = zsmltot(:,:)   * rn_ucf/zfn
533         ztmlres (:,:)   = ztmlres(:,:)   * rn_ucf/zfn  ; zsmlres (:,:)   = zsmlres(:,:)   * rn_ucf/zfn
534         ztmlatf (:,:)   = ztmlatf(:,:)   * rn_ucf/zfn  ; zsmlatf (:,:)   = zsmlatf(:,:)   * rn_ucf/zfn
535
536         tml_sum (:,:)   = tml_sum (:,:)  /  (2*zfn) ; sml_sum (:,:)   = sml_sum (:,:)  /  (2*zfn)
537         ztmltot2(:,:)   = ztmltot2(:,:)  * rn_ucf/zfn2 ; zsmltot2(:,:)   = zsmltot2(:,:)  * rn_ucf/zfn2
538         ztmltrd2(:,:,:) = ztmltrd2(:,:,:)* rn_ucf/zfn2 ; zsmltrd2(:,:,:) = zsmltrd2(:,:,:)* rn_ucf/zfn2
539         ztmlatf2(:,:)   = ztmlatf2(:,:)  * rn_ucf/zfn2 ; zsmlatf2(:,:)   = zsmlatf2(:,:)  * rn_ucf/zfn2
540         ztmlres2(:,:)   = ztmlres2(:,:)  * rn_ucf/zfn2 ; zsmlres2(:,:)   = zsmlres2(:,:)  * rn_ucf/zfn2
541
542         rmld_sum(:,:)   = rmld_sum(:,:)  /  (2*zfn)  ! similar to tml_sum and sml_sum
543
544         ! * Debugging information *
545         IF( lldebug ) THEN
546            !
547            WRITE(numout,*)
548            WRITE(numout,*) 'trd_mld : write trends in the Mixed Layer for debugging process:'
549            WRITE(numout,*) '~~~~~~~  '
550            WRITE(numout,*) '          TRA kt = ', kt, 'nmoymltrd = ', nmoymltrd
551            WRITE(numout,*)
552            WRITE(numout,*) '          >>>>>>>>>>>>>>>>>>  TRA TEMPERATURE <<<<<<<<<<<<<<<<<<'
553            WRITE(numout,*) '          TRA ztmlres    : ', SUM(ztmlres(:,:))
554            WRITE(numout,*) '          TRA ztmltot    : ', SUM(ztmltot(:,:))
555            WRITE(numout,*) '          TRA tmltrdm    : ', SUM(tmltrdm(:,:))
556            WRITE(numout,*) '          TRA tmlatfb    : ', SUM(tmlatfb(:,:))
557            WRITE(numout,*) '          TRA tmlatfn    : ', SUM(tmlatfn(:,:))
558            DO jl = 1, jpltrd
559               WRITE(numout,*) '          * TRA TREND INDEX jpmld_xxx = jl = ', jl, &
560                    & ' tmltrd : ', SUM(tmltrd(:,:,jl))
561            END DO
562            WRITE(numout,*) '          TRA ztmlres (jpi/2,jpj/2) : ', ztmlres (jpi/2,jpj/2)
563            WRITE(numout,*) '          TRA ztmlres2(jpi/2,jpj/2) : ', ztmlres2(jpi/2,jpj/2)
564            WRITE(numout,*)
565            WRITE(numout,*) '          >>>>>>>>>>>>>>>>>>  TRA SALINITY <<<<<<<<<<<<<<<<<<'
566            WRITE(numout,*) '          TRA zsmlres    : ', SUM(zsmlres(:,:))
567            WRITE(numout,*) '          TRA zsmltot    : ', SUM(zsmltot(:,:))
568            WRITE(numout,*) '          TRA smltrdm    : ', SUM(smltrdm(:,:))
569            WRITE(numout,*) '          TRA smlatfb    : ', SUM(smlatfb(:,:))
570            WRITE(numout,*) '          TRA smlatfn    : ', SUM(smlatfn(:,:))
571            DO jl = 1, jpltrd
572               WRITE(numout,*) '          * TRA TREND INDEX jpmld_xxx = jl = ', jl, &
573                    & ' smltrd : ', SUM(smltrd(:,:,jl))
574            END DO
575            WRITE(numout,*) '          TRA zsmlres (jpi/2,jpj/2) : ', zsmlres (jpi/2,jpj/2)
576            WRITE(numout,*) '          TRA zsmlres2(jpi/2,jpj/2) : ', zsmlres2(jpi/2,jpj/2)
577            !
578         END IF
579         !
580      END IF MODULO_NTRD
581
582      ! ======================================================================
583      ! IV. Write trends in the NetCDF file
584      ! ======================================================================
585
586      ! IV.1 Code for dimg mpp output
587      ! -----------------------------
588
589#if defined key_dimgout
590
591      IF( MOD( itmod, nn_trd ) == 0 ) THEN
592         iyear =  ndastp/10000
593         imon  = (ndastp-iyear*10000)/100
594         iday  =  ndastp - imon*100 - iyear*10000
595         WRITE(clname,9000) TRIM(cexper),'MLDiags',iyear,imon,iday
596         WRITE(clmode,'(f5.1,a)') nn_trd*rdt/86400.,' days average'
597         cltext = TRIM(cexper)//' mld diags'//TRIM(clmode)
598         CALL dia_wri_dimg (clname, cltext, smltrd, jpltrd, '2')
599      END IF
600
6019000  FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")
602
603#else
604     
605      ! IV.2 Code for IOIPSL/NetCDF output
606      ! ----------------------------------
607
608      IF( lwp .AND. MOD( itmod , nn_trd ) == 0 ) THEN
609         WRITE(numout,*) ' '
610         WRITE(numout,*) 'trd_mld : write trends in the NetCDF file :'
611         WRITE(numout,*) '~~~~~~~  '
612         WRITE(numout,*) '          ', TRIM(clhstnam), ' at kt = ', kt
613         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
614         WRITE(numout,*) ' '
615      END IF
616         
617      !-- Write the trends for T/S instantaneous diagnostics
618      IF( ln_trdmld_instant ) THEN           
619
620         CALL histwrite( nidtrd, "mxl_depth", it, rmld(:,:), ndimtrd1, ndextrd1 )
621         
622         !................................. ( ML temperature ) ...................................
623         
624         !-- Output the fields
625         CALL histwrite( nidtrd, "tml"     , it, tml    (:,:), ndimtrd1, ndextrd1 ) 
626         CALL histwrite( nidtrd, "tml_tot" , it, ztmltot(:,:), ndimtrd1, ndextrd1 ) 
627         CALL histwrite( nidtrd, "tml_res" , it, ztmlres(:,:), ndimtrd1, ndextrd1 ) 
628         
629         DO jl = 1, jpltrd - 1
630            CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)),            &
631                 &          it, tmltrd (:,:,jl), ndimtrd1, ndextrd1 )
632         END DO
633         
634         CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)),        &
635              &          it, ztmlatf(:,:), ndimtrd1, ndextrd1 )
636         
637         !.................................. ( ML salinity ) .....................................
638         
639         !-- Output the fields
640         CALL histwrite( nidtrd, "sml"     , it, sml    (:,:), ndimtrd1, ndextrd1 ) 
641         CALL histwrite( nidtrd, "sml_tot" , it, zsmltot(:,:), ndimtrd1, ndextrd1 ) 
642         CALL histwrite( nidtrd, "sml_res" , it, zsmlres(:,:), ndimtrd1, ndextrd1 ) 
643         
644         DO jl = 1, jpltrd - 1
645            CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)),            &
646                 &          it, smltrd(:,:,jl), ndimtrd1, ndextrd1 )
647         END DO
648         
649         CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)),        &
650              &          it, zsmlatf(:,:), ndimtrd1, ndextrd1 )
651         
652         IF( kt == nitend )   CALL histclo( nidtrd )
653
654      !-- Write the trends for T/S mean diagnostics
655      ELSE
656         
657         CALL histwrite( nidtrd, "mxl_depth", it, rmld_sum(:,:), ndimtrd1, ndextrd1 ) 
658         
659         !................................. ( ML temperature ) ...................................
660         
661         !-- Output the fields
662         CALL histwrite( nidtrd, "tml"     , it, tml_sum (:,:), ndimtrd1, ndextrd1 ) 
663         CALL histwrite( nidtrd, "tml_tot" , it, ztmltot2(:,:), ndimtrd1, ndextrd1 ) 
664         CALL histwrite( nidtrd, "tml_res" , it, ztmlres2(:,:), ndimtrd1, ndextrd1 ) 
665         
666         DO jl = 1, jpltrd - 1
667            CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)),            &
668                 &          it, ztmltrd2(:,:,jl), ndimtrd1, ndextrd1 )
669         END DO
670         
671         CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)),        &
672              &          it, ztmlatf2(:,:), ndimtrd1, ndextrd1 )
673         
674         !.................................. ( ML salinity ) .....................................
675                     
676         !-- Output the fields
677         CALL histwrite( nidtrd, "sml"     , it, sml_sum (:,:), ndimtrd1, ndextrd1 ) 
678         CALL histwrite( nidtrd, "sml_tot" , it, zsmltot2(:,:), ndimtrd1, ndextrd1 ) 
679         CALL histwrite( nidtrd, "sml_res" , it, zsmlres2(:,:), ndimtrd1, ndextrd1 ) 
680         
681         DO jl = 1, jpltrd - 1
682            CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)),            &
683                 &          it, zsmltrd2(:,:,jl), ndimtrd1, ndextrd1 )
684         END DO
685         
686         CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)),        &
687              &          it, zsmlatf2(:,:), ndimtrd1, ndextrd1 )
688         
689         IF( kt == nitend )   CALL histclo( nidtrd )
690
691      END IF
692     
693      ! Compute the control surface (for next time step) : flag = on
694      icount = 1
695      !
696#endif
697
698      IF( MOD( itmod, nn_trd ) == 0 ) THEN
699         !
700         ! III.5 Reset cumulative arrays to zero
701         ! -------------------------------------
702         nmoymltrd = 0
703         
704         !   ... temperature ...               ... salinity ...
705         tmltrdm        (:,:)   = 0.e0   ;   smltrdm        (:,:)   = 0.e0
706         tmlatfm        (:,:)   = 0.e0   ;   smlatfm        (:,:)   = 0.e0
707         tml_sum        (:,:)   = 0.e0   ;   sml_sum        (:,:)   = 0.e0
708         tmltrd_csum_ln (:,:,:) = 0.e0   ;   smltrd_csum_ln (:,:,:) = 0.e0
709         tmltrd_sum     (:,:,:) = 0.e0   ;   smltrd_sum     (:,:,:) = 0.e0
710
711         rmld_sum       (:,:)   = 0.e0
712         !
713      END IF
714
715      ! ======================================================================
716      ! V. Write restart file
717      ! ======================================================================
718
719      IF( lrst_oce )   CALL trd_mld_rst_write( kt ) 
720
721      CALL wrk_dealloc( jpi, jpj,         ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf                        )
722      CALL wrk_dealloc( jpi, jpj,         ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 ) 
723      CALL wrk_dealloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2                                                               )
724      !
725   END SUBROUTINE trd_mld
726
727
728   SUBROUTINE trd_mld_init
729      !!----------------------------------------------------------------------
730      !!                  ***  ROUTINE trd_mld_init  ***
731      !!
732      !! ** Purpose :   computation of vertically integrated T and S budgets
733      !!      from ocean surface down to control surface (NetCDF output)
734      !!----------------------------------------------------------------------
735      INTEGER :: jl
736      INTEGER :: inum   ! logical unit
737      REAL(wp) ::   zjulian, zsto, zout
738      CHARACTER (LEN=40) ::   clop
739      CHARACTER (LEN=12) ::   clmxl, cltu, clsu
740      !!----------------------------------------------------------------------
741
742      ! ======================================================================
743      ! I. initialization
744      ! ======================================================================
745
746      IF(lwp) THEN
747         WRITE(numout,*)
748         WRITE(numout,*) ' trd_mld_init : Mixed-layer trends'
749         WRITE(numout,*) ' ~~~~~~~~~~~~~'
750         WRITE(numout,*) '                namelist namtrd read in trd_mod_init                        '
751         WRITE(numout,*)
752      END IF
753
754      ! I.1 Check consistency of user defined preferences
755      ! -------------------------------------------------
756
757      IF( MOD( nitend, nn_trd ) /= 0 ) THEN
758         WRITE(numout,cform_err)
759         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
760         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
761         WRITE(numout,*) '                          you defined, nn_trd   = ', nn_trd
762         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
763         WRITE(numout,*) '                You should reconsider this choice.                        ' 
764         WRITE(numout,*) 
765         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
766         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
767         nstop = nstop + 1
768      END IF
769
770      IF( nn_cla == 1 )   CALL ctl_warn( '      You set n_cla = 1. Note that the Mixed-Layer diagnostics  ',   &
771         &                               '      are not exact along the corresponding straits.            ')
772
773      !                                   ! allocate trdmld arrays
774      IF( trd_mld_alloc()    /= 0 )   CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld     arrays' )
775      IF( trdmld_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld_oce arrays' )
776
777      ! I.2 Initialize arrays to zero or read a restart file
778      ! ----------------------------------------------------
779
780      nmoymltrd = 0
781
782      !     ... temperature ...                  ... salinity ...
783      tml            (:,:)   = 0.e0    ;    sml            (:,:)   = 0.e0     ! inst.
784      tmltrdm        (:,:)   = 0.e0    ;    smltrdm        (:,:)   = 0.e0
785      tmlatfm        (:,:)   = 0.e0    ;    smlatfm        (:,:)   = 0.e0
786      tml_sum        (:,:)   = 0.e0    ;    sml_sum        (:,:)   = 0.e0     ! mean
787      tmltrd_sum     (:,:,:) = 0.e0    ;    smltrd_sum     (:,:,:) = 0.e0
788      tmltrd_csum_ln (:,:,:) = 0.e0    ;    smltrd_csum_ln (:,:,:) = 0.e0
789
790      rmld           (:,:)   = 0.e0           
791      rmld_sum       (:,:)   = 0.e0
792
793      IF( ln_rstart .AND. ln_trdmld_restart ) THEN
794         CALL trd_mld_rst_read
795      ELSE
796         !     ... temperature ...                  ... salinity ...
797         tmlb           (:,:)   = 0.e0    ;    smlb           (:,:)   = 0.e0  ! inst.
798         tmlbb          (:,:)   = 0.e0    ;    smlbb          (:,:)   = 0.e0 
799         tmlbn          (:,:)   = 0.e0    ;    smlbn          (:,:)   = 0.e0 
800         tml_sumb       (:,:)   = 0.e0    ;    sml_sumb       (:,:)   = 0.e0  ! mean
801         tmltrd_csum_ub (:,:,:) = 0.e0    ;    smltrd_csum_ub (:,:,:) = 0.e0
802         tmltrd_atf_sumb(:,:)   = 0.e0    ;    smltrd_atf_sumb(:,:)   = 0.e0 
803      END IF
804
805      icount = 1   ;   ionce  = 1                            ! open specifier
806
807      ! I.3 Read control surface from file ctlsurf_idx
808      ! ----------------------------------------------
809 
810      IF( nn_ctls == 1 ) THEN
811         CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
812         READ ( inum ) nbol
813         CLOSE( inum )
814      END IF
815
816      ! ======================================================================
817      ! II. netCDF output initialization
818      ! ======================================================================
819
820#if defined key_dimgout 
821      ???
822#else
823      ! clmxl = legend root for netCDF output
824      IF( nn_ctls == 0 ) THEN      ! control surface = mixed-layer with density criterion
825         clmxl = 'Mixed Layer '  !                   (array nmln computed in zdfmxl.F90)
826      ELSE IF( nn_ctls == 1 ) THEN ! control surface = read index from file
827         clmxl = '      Bowl '
828      ELSE IF( nn_ctls >= 2 ) THEN ! control surface = model level
829         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls
830      END IF
831
832      ! II.1 Define frequency of output and means
833      ! -----------------------------------------
834      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
835      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
836      ENDIF
837#  if defined key_diainstant
838      IF( .NOT. ln_trdmld_instant ) THEN
839         CALL ctl_stop( 'trd_mld : this was never checked. Comment this line to proceed...' )
840      END IF
841      zsto = nn_trd * rdt
842      clop = "inst("//TRIM(clop)//")"
843#  else
844      IF( ln_trdmld_instant ) THEN
845         zsto = rdt                 ! inst. diags : we use IOIPSL time averaging
846      ELSE
847         zsto = nn_trd * rdt          ! mean  diags : we DO NOT use any IOIPSL time averaging
848      END IF
849      clop = "ave("//TRIM(clop)//")"
850#  endif
851      zout = nn_trd * rdt
852
853      IF(lwp) WRITE (numout,*) '                netCDF initialization'
854
855      ! II.2 Compute julian date from starting date of the run
856      ! ------------------------------------------------------
857      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
858      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
859      IF(lwp) WRITE(numout,*)' ' 
860      IF(lwp) WRITE(numout,*)'                Date 0 used :',nit000,    &
861         &                   ' YEAR ', nyear,' MONTH '      , nmonth,   &
862         &                   ' DAY ' , nday, 'Julian day : ', zjulian
863
864
865      ! II.3 Define the T grid trend file (nidtrd)
866      ! ------------------------------------------
867      !-- Define long and short names for the NetCDF output variables
868      !       ==> choose them according to trdmld_oce.F90 <==
869
870      ctrd(jpmld_xad,1) = " Zonal advection"                  ;   ctrd(jpmld_xad,2) = "_xad"
871      ctrd(jpmld_yad,1) = " Meridional advection"             ;   ctrd(jpmld_yad,2) = "_yad"
872      ctrd(jpmld_zad,1) = " Vertical advection"               ;   ctrd(jpmld_zad,2) = "_zad"
873      ctrd(jpmld_ldf,1) = " Lateral diffusion"                ;   ctrd(jpmld_ldf,2) = "_ldf"
874      ctrd(jpmld_for,1) = " Forcing"                          ;   ctrd(jpmld_for,2) = "_for"
875      ctrd(jpmld_zdf,1) = " Vertical diff. (Kz)"              ;   ctrd(jpmld_zdf,2) = "_zdf"
876      ctrd(jpmld_bbc,1) = " Geothermal flux"                  ;   ctrd(jpmld_bbc,2) = "_bbc"
877      ctrd(jpmld_bbl,1) = " Adv/diff. Bottom boundary layer"  ;   ctrd(jpmld_bbl,2) = "_bbl"
878      ctrd(jpmld_dmp,1) = " Tracer damping"                   ;   ctrd(jpmld_dmp,2) = "_dmp"
879      ctrd(jpmld_npc,1) = " Non penetrative convec. adjust."  ;   ctrd(jpmld_npc,2) = "_npc"
880      ctrd(jpmld_atf,1) = " Asselin time filter"              ;   ctrd(jpmld_atf,2) = "_atf"
881                                                                 
882      !-- Create a NetCDF file and enter the define mode
883      CALL dia_nam( clhstnam, nn_trd, 'trends' )
884      IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam
885      CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
886      &             1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_t, nidtrd, domain_id=nidom, snc4chunks=snc4set )
887
888      !-- Define the ML depth variable
889      CALL histdef(nidtrd, "mxl_depth", clmxl//"  Mixed Layer Depth"              , "m",         &
890                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
891
892      !-- Define physical units
893      IF     ( rn_ucf == 1.       ) THEN   ;   cltu = "degC/s"     ;   clsu = "p.s.u./s"
894      ELSEIF ( rn_ucf == 3600.*24.) THEN   ;   cltu = "degC/day"   ;   clsu = "p.s.u./day"
895      ELSE                                 ;   cltu = "unknown?"   ;   clsu = "unknown?"
896      END IF
897
898
899      !-- Define miscellaneous T and S mixed-layer variables
900
901      IF( jpltrd /= jpmld_atf ) CALL ctl_stop( 'Error : jpltrd /= jpmld_atf' ) ! see below
902
903      !................................. ( ML temperature ) ...................................
904
905      CALL histdef(nidtrd, "tml"      , clmxl//" T Mixed Layer Temperature"       ,  "C",        &
906                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
907      CALL histdef(nidtrd, "tml_tot",   clmxl//" T Total trend"                   , cltu,        & 
908                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )             
909      CALL histdef(nidtrd, "tml_res",   clmxl//" T dh/dt Entrainment (Resid.)"    , cltu,        & 
910                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
911     
912      DO jl = 1, jpltrd - 1      ! <== only true if jpltrd == jpmld_atf
913         CALL histdef(nidtrd, trim("tml"//ctrd(jl,2)), clmxl//" T"//ctrd(jl,1), cltu,            & 
914                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
915      END DO                                                                 ! if zsto=rdt above
916     
917      CALL histdef(nidtrd, trim("tml"//ctrd(jpmld_atf,2)), clmxl//" T"//ctrd(jpmld_atf,1), cltu, & 
918                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
919     
920      !.................................. ( ML salinity ) .....................................
921     
922      CALL histdef(nidtrd, "sml"      , clmxl//" S Mixed Layer Salinity"          , "p.s.u.",       &
923                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
924      CALL histdef(nidtrd, "sml_tot",   clmxl//" S Total trend"                   , clsu,        & 
925                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )             
926      CALL histdef(nidtrd, "sml_res",   clmxl//" S dh/dt Entrainment (Resid.)"    , clsu,        & 
927                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
928     
929      DO jl = 1, jpltrd - 1      ! <== only true if jpltrd == jpmld_atf
930         CALL histdef(nidtrd, trim("sml"//ctrd(jl,2)), clmxl//" S"//ctrd(jl,1), clsu,            & 
931                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
932      END DO                                                                 ! if zsto=rdt above
933     
934      CALL histdef(nidtrd, trim("sml"//ctrd(jpmld_atf,2)), clmxl//" S"//ctrd(jpmld_atf,1), clsu, & 
935                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
936
937      !-- Leave IOIPSL/NetCDF define mode
938      CALL histend( nidtrd, snc4set )
939
940#endif        /* key_dimgout */
941   END SUBROUTINE trd_mld_init
942
943   !!======================================================================
944END MODULE trdmld
Note: See TracBrowser for help on using the repository browser.