source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 10 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

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