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 branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by gm, 11 years ago

dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation

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