source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 18 months ago

The Dr Hook changes from my perl code.

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