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

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

source: NEMO/trunk/src/OCE/TRD/trdmxl.F90 @ 13286

Last change on this file since 13286 was 13237, checked in by smasson, 4 years ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

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