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 @ 11048

Last change on this file since 11048 was 10425, checked in by smasson, 5 years ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

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