source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 19 months ago

GMED 450 add flush after prints

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