source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

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