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 trunk/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 7881

Last change on this file since 7881 was 6140, checked in by timgraham, 8 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

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