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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdmxl.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.

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