New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trdmld.F90 in trunk/NEMO/OPA_SRC/TRD – NEMO

source: trunk/NEMO/OPA_SRC/TRD/trdmld.F90 @ 1334

Last change on this file since 1334 was 1334, checked in by smasson, 15 years ago

complete work on time origin in outputs (ticket:335) + downward vertical axis (ticket:357)

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