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

Last change on this file since 13226 was 13226, checked in by orioltp, 3 months ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

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