source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 7 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

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