New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trdmld.F90 in branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 3792

Last change on this file since 3792 was 3792, checked in by gm, 11 years ago

dev_MERGE_2012: #960 : correct bug in trdmld & trdmld_trc (time step control) + KPP/TRC compatibility (zdfkpp.F90)

  • Property svn:keywords set to Id
File size: 48.5 KB
Line 
1MODULE trdmld
2   !!======================================================================
3   !!                       ***  MODULE  trdmld  ***
4   !! Ocean diagnostics:  mixed layer T-S trends
5   !!=====================================================================
6   !! History :       !  95-04  (J. Vialard)    Original code
7   !!                 !  97-02  (E. Guilyardi)  Adaptation global + base cmo
8   !!                 !  99-09  (E. Guilyardi)  Re-writing + netCDF output
9   !!            8.5  !  02-06  (G. Madec)      F90: Free form and module
10   !!            9.0  !  04-08  (C. Talandier)  New trends organization
11   !!                 !  05-05  (C. Deltel)     Diagnose trends of time averaged ML T & S
12   !!----------------------------------------------------------------------
13#if   defined key_trdmld   ||   defined key_esopa
14   !!----------------------------------------------------------------------
15   !!   'key_trdmld'                          mixed layer trend diagnostics
16   !!----------------------------------------------------------------------
17   !!   trd_mld          : T and S cumulated trends averaged over the mixed layer
18   !!   trd_mld_zint     : T and S trends vertical integration
19   !!   trd_mld_init     : initialization step
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and tracers variables
22   USE dom_oce         ! ocean space and time domain variables
23   USE trdmod_oce      ! ocean variables trends
24   USE trdmld_oce      ! ocean variables trends
25   USE ldftra_oce      ! ocean active tracers lateral physics
26   USE zdf_oce         ! ocean vertical physics
27   USE in_out_manager  ! I/O manager
28   USE phycst          ! Define parameters for the routines
29   USE dianam          ! build the name of file (routine)
30   USE ldfslp          ! iso-neutral slopes
31   USE zdfmxl          ! mixed layer depth
32   USE zdfddm          ! ocean vertical physics: double diffusion
33   USE ioipsl          ! NetCDF library
34   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
35   USE diadimg         ! dimg direct access file format output
36   USE trdmld_rst      ! restart for diagnosing the ML trends
37   USE prtctl          ! Print control
38   USE lib_mpp         ! MPP library
39   USE wrk_nemo        ! Memory allocation
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   trd_mld        ! routine called by step.F90
45   PUBLIC   trd_mld_init   ! routine called by opa.F90
46   PUBLIC   trd_mld_zint   ! routine called by tracers routines
47
48   CHARACTER (LEN=40) ::  clhstnam         ! name of the trends NetCDF file
49   INTEGER ::   nh_t, nmoymltrd
50   INTEGER ::   nidtrd
51   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::   ndextrd1
52   INTEGER ::   ndimtrd1                       
53   INTEGER ::   ionce, icount                   
54
55   !! * Substitutions
56#  include "domzgr_substitute.h90"
57#  include "ldftra_substitute.h90"
58#  include "zdfddm_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   INTEGER FUNCTION trd_mld_alloc()
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE trd_mld_alloc  ***
69      !!----------------------------------------------------------------------
70      ALLOCATE( ndextrd1(jpi*jpj) , STAT=trd_mld_alloc )
71      !
72      IF( lk_mpp             )   CALL mpp_sum ( trd_mld_alloc )
73      IF( trd_mld_alloc /= 0 )   CALL ctl_warn('trd_mld_alloc: failed to allocate array ndextrd1')
74   END FUNCTION trd_mld_alloc
75
76
77   SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
78      !!----------------------------------------------------------------------
79      !!                  ***  ROUTINE trd_mld_zint  ***
80      !!
81      !! ** Purpose :   Compute the vertical average of the 3D fields given as arguments
82      !!                to the subroutine. This vertical average is performed from ocean
83      !!                surface down to a chosen control surface.
84      !!
85      !! ** Method/usage :
86      !!      The control surface can be either a mixed layer depth (time varying)
87      !!      or a fixed surface (jk level or bowl).
88      !!      Choose control surface with nn_ctls in namelist NAMTRD :
89      !!        nn_ctls = 0  : use mixed layer with density criterion
90      !!        nn_ctls = 1  : read index from file 'ctlsurf_idx'
91      !!        nn_ctls > 1  : use fixed level surface jk = nn_ctls
92      !!      Note: in the remainder of the routine, the volume between the
93      !!            surface and the control surface is called "mixed-layer"
94      !!----------------------------------------------------------------------
95      !
96      INTEGER                         , INTENT( in ) ::   ktrd       ! ocean trend index
97      CHARACTER(len=2)                , INTENT( in ) ::   ctype      ! 2D surface/bottom or 3D interior physics
98      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pttrdmld   ! temperature trend
99      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   pstrdmld   ! salinity trend
100      !
101      INTEGER ::   ji, jj, jk, isum
102      REAL(wp), POINTER, DIMENSION(:,:)  :: zvlmsk 
103      !!----------------------------------------------------------------------
104
105      CALL wrk_alloc( jpi, jpj, zvlmsk ) 
106
107      ! I. Definition of control surface and associated fields
108      ! ------------------------------------------------------
109      !            ==> only once per time step <==
110
111      IF( icount == 1 ) THEN       
112         !
113         tmltrd(:,:,:) = 0.e0    ;    smltrd(:,:,:) = 0.e0    ! <<< reset trend arrays to zero
114         
115         ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer
116         IF( nn_ctls == 0 ) THEN       ! * control surface = mixed-layer with density criterion
117            nmld(:,:) = nmln(:,:)    ! array nmln computed in zdfmxl.F90
118         ELSE IF( nn_ctls == 1 ) THEN  ! * control surface = read index from file
119            nmld(:,:) = nbol(:,:)
120         ELSE IF( nn_ctls >= 2 ) THEN  ! * control surface = model level
121            nn_ctls = MIN( nn_ctls, jpktrd - 1 )
122            nmld(:,:) = nn_ctls + 1
123         ENDIF
124
125         ! ... Compute ndextrd1 and ndimtrd1 only once
126         IF( ionce == 1 ) THEN
127            !
128            ! Check of validity : nmld(ji,jj) <= jpktrd
129            isum        = 0
130            zvlmsk(:,:) = 0.e0
131
132            IF( jpktrd < jpk ) THEN
133               DO jj = 1, jpj
134                  DO ji = 1, jpi
135                     IF( nmld(ji,jj) <= jpktrd ) THEN
136                        zvlmsk(ji,jj) = tmask(ji,jj,1)
137                     ELSE
138                        isum = isum + 1
139                        zvlmsk(ji,jj) = 0.
140                     ENDIF
141                  END DO
142               END DO
143            ENDIF
144
145            ! Index of ocean points (2D only)
146            IF( isum > 0 ) THEN
147               WRITE(numout,*)' Number of invalid points nmld > jpktrd', isum 
148               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )    ! surface
149            ELSE
150               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )    ! surface
151            ENDIF                               
152
153            ionce = 0                ! no more pass here
154            !
155         END IF
156         
157         ! ... Weights for vertical averaging
158         wkx(:,:,:) = 0.e0
159         DO jk = 1, jpktrd             ! initialize wkx with vertical scale factor in mixed-layer
160            DO jj = 1,jpj
161               DO ji = 1,jpi
162                  IF( jk - nmld(ji,jj) < 0.e0 )   wkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
163               END DO
164            END DO
165         END DO
166         
167         rmld(:,:) = 0.e0                ! compute mixed-layer depth : rmld
168         DO jk = 1, jpktrd
169            rmld(:,:) = rmld(:,:) + wkx(:,:,jk)
170         END DO
171         
172         DO jk = 1, jpktrd             ! compute integration weights
173            wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1., rmld(:,:) )
174         END DO
175
176         icount = 0                    ! <<< flag = off : control surface & integr. weights
177         !                             !     computed only once per time step
178      END IF
179
180      ! II. Vertical integration of trends in the mixed-layer
181      ! -----------------------------------------------------
182
183      SELECT CASE (ctype)
184      CASE ( '3D' )   ! mean T/S trends in the mixed-layer
185         DO jk = 1, jpktrd
186            tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,jk) * wkx(:,:,jk)   ! temperature
187            smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,jk) * wkx(:,:,jk)   ! salinity
188         END DO
189      CASE ( '2D' )   ! forcing at upper boundary of the mixed-layer
190         tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,1) * wkx(:,:,1)        ! non penetrative
191         smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,1) * wkx(:,:,1)           
192      END SELECT
193      !
194      CALL wrk_dealloc( jpi, jpj, zvlmsk ) 
195      !
196   END SUBROUTINE trd_mld_zint
197   
198
199   SUBROUTINE trd_mld( kt )
200      !!----------------------------------------------------------------------
201      !!                  ***  ROUTINE trd_mld  ***
202      !!
203      !! ** Purpose :  Compute and cumulate the mixed layer trends over an analysis
204      !!               period, and write NetCDF (or dimg) outputs.
205      !!
206      !! ** Method/usage :
207      !!          The stored trends can be chosen twofold (according to the ln_trdmld_instant
208      !!          logical namelist variable) :
209      !!          1) to explain the difference between initial and final
210      !!             mixed-layer T & S (where initial and final relate to the
211      !!             current analysis window, defined by nn_trd in the namelist)
212      !!          2) to explain the difference between the current and previous
213      !!             TIME-AVERAGED mixed-layer T & S (where time-averaging is
214      !!             performed over each analysis window).
215      !!
216      !! ** Consistency check :
217      !!        If the control surface is fixed ( nn_ctls > 1 ), the residual term (dh/dt
218      !!        entrainment) should be zero, at machine accuracy. Note that in the case
219      !!        of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO
220      !!        over the first two analysis windows (except if restart).
221      !!        N.B. For ORCA2_LIM, use e.g. nn_trd=5, rn_ucf=1., nn_ctls=8
222      !!             for checking residuals.
223      !!             On a NEC-SX5 computer, this typically leads to:
224      !!                   O(1.e-20) temp. residuals (tml_res) when ln_trdmld_instant=.false.
225      !!                   O(1.e-21) temp. residuals (tml_res) when ln_trdmld_instant=.true.
226      !!
227      !! ** Action :
228      !!       At each time step, mixed-layer averaged trends are stored in the
229      !!       tmltrd(:,:,jpmld_xxx) array (see trdmld_oce.F90 for definitions of jpmld_xxx).
230      !!       This array is known when trd_mld is called, at the end of the stp subroutine,
231      !!       except for the purely vertical K_z diffusion term, which is embedded in the
232      !!       lateral diffusion trend.
233      !!
234      !!       In I), this K_z term is diagnosed and stored, thus its contribution is removed
235      !!       from the lateral diffusion trend.
236      !!       In II), the instantaneous mixed-layer T & S are computed, and misc. cumulative
237      !!       arrays are updated.
238      !!       In III), called only once per analysis window, we compute the total trends,
239      !!       along with the residuals and the Asselin correction terms.
240      !!       In IV), the appropriate trends are written in the trends NetCDF file.
241      !!
242      !! References :
243      !!       - Vialard & al.
244      !!       - See NEMO documentation (in preparation)
245      !!----------------------------------------------------------------------
246      !
247      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
248      !
249      INTEGER :: ji, jj, jk, jl, ik, it, itmod
250      LOGICAL :: lldebug = .TRUE.
251      REAL(wp) :: zavt, zfn, zfn2
252      !                                              ! z(ts)mltot : dT/dt over the anlysis window (including Asselin)
253      !                                              ! z(ts)mlres : residual = dh/dt entrainment term
254      REAL(wp), POINTER, DIMENSION(:,:  ) ::  ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
255      REAL(wp), POINTER, DIMENSION(:,:  ) ::  ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 
256      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztmltrd2, zsmltrd2   ! only needed for mean diagnostics
257#if defined key_dimgout
258      INTEGER ::  iyear,imon,iday
259      CHARACTER(LEN=80) :: cltext, clmode
260#endif
261      !!----------------------------------------------------------------------
262 
263      CALL wrk_alloc( jpi, jpj,         ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf                        )
264      CALL wrk_alloc( jpi, jpj,         ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 ) 
265      CALL wrk_alloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2                                                               )
266
267      ! ======================================================================
268      ! I. Diagnose the purely vertical (K_z) diffusion trend
269      ! ======================================================================
270
271      ! ... These terms can be estimated by flux computation at the lower boundary of the ML
272      !     (we compute (-1/h) * K_z * d_z( T ) and (-1/h) * K_z * d_z( S ))
273      IF( ln_traldf_iso ) THEN
274         DO jj = 1,jpj
275            DO ji = 1,jpi
276               ik = nmld(ji,jj)
277               zavt = avt(ji,jj,ik)
278               tmltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
279                  &                      * ( tsn(ji,jj,ik-1,jp_tem) - tsn(ji,jj,ik,jp_tem) )         &
280                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
281               zavt = fsavs(ji,jj,ik)
282               smltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
283                  &                      * ( tsn(ji,jj,ik-1,jp_sal) - tsn(ji,jj,ik,jp_sal) )         &
284                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
285            END DO
286         END DO
287
288         ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
289         tmltrd(:,:,jpmld_ldf) = tmltrd(:,:,jpmld_ldf) - tmltrd(:,:,jpmld_zdf)
290         smltrd(:,:,jpmld_ldf) = smltrd(:,:,jpmld_ldf) - smltrd(:,:,jpmld_zdf)
291      END IF
292
293      ! ... Lateral boundary conditions
294      DO jl = 1, jpltrd
295         CALL lbc_lnk( tmltrd(:,:,jl), 'T', 1. )
296         CALL lbc_lnk( smltrd(:,:,jl), 'T', 1. )
297      END DO
298
299      ! ======================================================================
300      ! II. Cumulate the trends over the analysis window
301      ! ======================================================================
302
303      ztmltrd2(:,:,:) = 0.e0   ;    zsmltrd2(:,:,:) = 0.e0  ! <<< reset arrays to zero
304      ztmltot2(:,:)   = 0.e0   ;    zsmltot2(:,:)   = 0.e0
305      ztmlres2(:,:)   = 0.e0   ;    zsmlres2(:,:)   = 0.e0
306      ztmlatf2(:,:)   = 0.e0   ;    zsmlatf2(:,:)   = 0.e0
307
308      ! II.1 Set before values of vertically average T and S
309      ! ----------------------------------------------------
310      IF( kt > nit000 ) THEN
311         !   ... temperature ...                    ... salinity ...
312         tmlb   (:,:) = tml   (:,:)           ; smlb   (:,:) = sml   (:,:)
313         tmlatfn(:,:) = tmltrd(:,:,jpmld_atf) ; smlatfn(:,:) = smltrd(:,:,jpmld_atf)
314      END IF
315
316      ! II.2 Vertically averaged T and S
317      ! --------------------------------
318      tml(:,:) = 0.e0   ;   sml(:,:) = 0.e0
319      DO jk = 1, jpktrd - 1
320         tml(:,:) = tml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_tem)
321         sml(:,:) = sml(:,:) + wkx(:,:,jk) * tsn(:,:,jk,jp_sal)
322      END DO
323
324      ! II.3 Initialize mixed-layer "before" arrays for the 1rst analysis window   
325      ! ------------------------------------------------------------------------
326      IF( kt == 2 ) THEN  !  i.e. ( .NOT. ln_rstart ).AND.( kt == nit000 + 1)
327         !
328         !   ... temperature ...                ... salinity ...
329         tmlbb  (:,:) = tmlb   (:,:)   ;   smlbb  (:,:) = smlb   (:,:)
330         tmlbn  (:,:) = tml    (:,:)   ;   smlbn  (:,:) = sml    (:,:)
331         tmlatfb(:,:) = tmlatfn(:,:)   ;   smlatfb(:,:) = smlatfn(:,:)
332         
333         tmltrd_csum_ub (:,:,:) = 0.e0  ;   smltrd_csum_ub (:,:,:) = 0.e0
334         tmltrd_atf_sumb(:,:)   = 0.e0  ;   smltrd_atf_sumb(:,:)   = 0.e0
335
336         rmldbn(:,:) = rmld(:,:)
337
338         IF( ln_ctl ) THEN
339            WRITE(numout,*) '             we reach kt == nit000 + 1 = ', nit000+1
340            CALL prt_ctl(tab2d_1=tmlbb   , clinfo1=' tmlbb   -   : ', mask1=tmask, ovlap=1)
341            CALL prt_ctl(tab2d_1=tmlbn   , clinfo1=' tmlbn   -   : ', mask1=tmask, ovlap=1)
342            CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb -   : ', mask1=tmask, ovlap=1)
343         END IF
344         !
345      END IF
346
347      IF( ( ln_rstart ) .AND. ( kt == nit000 ) .AND. ( ln_ctl ) ) THEN
348         IF( ln_trdmld_instant ) THEN
349            WRITE(numout,*) '             restart from kt == nit000 = ', nit000
350            CALL prt_ctl(tab2d_1=tmlbb   , clinfo1=' tmlbb   -   : ', mask1=tmask, ovlap=1)
351            CALL prt_ctl(tab2d_1=tmlbn   , clinfo1=' tmlbn   -   : ', mask1=tmask, ovlap=1)
352            CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb -   : ', mask1=tmask, ovlap=1)
353         ELSE
354            WRITE(numout,*) '             restart from kt == nit000 = ', nit000
355            CALL prt_ctl(tab2d_1=tmlbn          , clinfo1=' tmlbn           -  : ', mask1=tmask, ovlap=1)
356            CALL prt_ctl(tab2d_1=rmldbn         , clinfo1=' rmldbn          -  : ', mask1=tmask, ovlap=1)
357            CALL prt_ctl(tab2d_1=tml_sumb       , clinfo1=' tml_sumb        -  : ', mask1=tmask, ovlap=1)
358            CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb -  : ', mask1=tmask, ovlap=1)
359            CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub  -  : ', mask1=tmask, ovlap=1, kdim=1)
360         END IF
361      END IF
362
363      ! II.4 Cumulated trends over the analysis period
364      ! ----------------------------------------------
365      !
366      !         [  1rst analysis window ] [     2nd analysis window     ]                       
367      !
368      !     o---[--o-----o-----o-----o--]-[--o-----o-----o-----o-----o--]---o-----o--> time steps
369      !                          nn_trd                           2*nn_trd       etc.
370      !     1      2     3     4    =5 e.g.                          =10
371      !
372      IF( ( kt >= 2 ).OR.( ln_rstart ) ) THEN
373         !
374         nmoymltrd = nmoymltrd + 1
375         
376         ! ... Cumulate over BOTH physical contributions AND over time steps
377         DO jl = 1, jpltrd
378            tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl)
379            smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl)
380         END DO
381
382         ! ... Special handling of the Asselin trend
383         tmlatfm(:,:) = tmlatfm(:,:) + tmlatfn(:,:)
384         smlatfm(:,:) = smlatfm(:,:) + smlatfn(:,:)
385
386         ! ... Trends associated with the time mean of the ML T/S
387         tmltrd_sum    (:,:,:) = tmltrd_sum    (:,:,:) + tmltrd    (:,:,:) ! tem
388         tmltrd_csum_ln(:,:,:) = tmltrd_csum_ln(:,:,:) + tmltrd_sum(:,:,:)
389         tml_sum       (:,:)   = tml_sum       (:,:)   + tml       (:,:)
390         smltrd_sum    (:,:,:) = smltrd_sum    (:,:,:) + smltrd    (:,:,:) ! sal
391         smltrd_csum_ln(:,:,:) = smltrd_csum_ln(:,:,:) + smltrd_sum(:,:,:)
392         sml_sum       (:,:)   = sml_sum       (:,:)   + sml       (:,:)
393         rmld_sum      (:,:)   = rmld_sum      (:,:)   + rmld      (:,:)   ! rmld
394         !
395      END IF
396
397      ! ======================================================================
398      ! III. Prepare fields for output (get here ONCE PER ANALYSIS PERIOD)
399      ! ======================================================================
400
401      ! Convert to appropriate physical units
402      ! N.B. It may be useful to check IOIPSL time averaging with :
403      !      tmltrd (:,:,:) = 1. ; smltrd (:,:,:) = 1.
404      tmltrd(:,:,:) = tmltrd(:,:,:) * rn_ucf   ! (actually needed for 1:jpltrd-1, but trdmld(:,:,jpltrd)
405      smltrd(:,:,:) = smltrd(:,:,:) * rn_ucf   !  is no longer used, and is reset to 0. at next time step)
406     
407      ! define time axis
408      it = kt
409      itmod = kt - nit000 + 1
410
411      MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN        ! nitend MUST be multiple of nn_trd
412         !
413         ztmltot (:,:) = 0.e0   ;   zsmltot (:,:) = 0.e0   ! reset arrays to zero
414         ztmlres (:,:) = 0.e0   ;   zsmlres (:,:) = 0.e0
415         ztmltot2(:,:) = 0.e0   ;   zsmltot2(:,:) = 0.e0
416         ztmlres2(:,:) = 0.e0   ;   zsmlres2(:,:) = 0.e0
417     
418         zfn  = float(nmoymltrd)    ;    zfn2 = zfn * zfn
419         
420         ! III.1 Prepare fields for output ("instantaneous" diagnostics)
421         ! -------------------------------------------------------------
422         
423         !-- Compute total trends
424         ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) / ( 2.*rdt )
425         zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) / ( 2.*rdt )
426         
427         !-- Compute residuals
428         ztmlres(:,:) = ztmltot(:,:) - ( tmltrdm(:,:) - tmlatfn(:,:) + tmlatfb(:,:) )
429         zsmlres(:,:) = zsmltot(:,:) - ( smltrdm(:,:) - smlatfn(:,:) + smlatfb(:,:) )
430     
431         !-- Diagnose Asselin trend over the analysis window
432         ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
433         zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
434         
435         !-- Lateral boundary conditions
436         !         ... temperature ...                    ... salinity ...
437         CALL lbc_lnk( ztmltot , 'T', 1. )  ;   CALL lbc_lnk( zsmltot , 'T', 1. )
438         CALL lbc_lnk( ztmlres , 'T', 1. )  ;   CALL lbc_lnk( zsmlres , 'T', 1. )
439         CALL lbc_lnk( ztmlatf , 'T', 1. )  ;   CALL lbc_lnk( zsmlatf , 'T', 1. )
440
441#if defined key_diainstant
442         CALL ctl_stop( 'tml_trd : key_diainstant was never checked within trdmld. Comment this to proceed.')
443#endif
444         ! III.2 Prepare fields for output ("mean" diagnostics)
445         ! ----------------------------------------------------
446         
447         !-- Update the ML depth time sum (to build the Leap-Frog time mean)
448         rmld_sum(:,:) = rmldbn(:,:) + 2 * ( rmld_sum(:,:) - rmld(:,:) ) + rmld(:,:)
449
450         !-- Compute temperature total trends
451         tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
452         ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) /  ( 2.*rdt )    ! now in degC/s
453         
454         !-- Compute salinity total trends
455         sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
456         zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) /  ( 2.*rdt )    ! now in psu/s
457         
458         !-- Compute temperature residuals
459         DO jl = 1, jpltrd
460            ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
461         END DO
462
463         ztmltrdm2(:,:) = 0.e0
464         DO jl = 1, jpltrd
465            ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
466         END DO
467
468         ztmlres2(:,:) =  ztmltot2(:,:)  -       &
469              ( ztmltrdm2(:,:) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:) )
470         
471         !-- Compute salinity residuals
472         DO jl = 1, jpltrd
473            zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
474         END DO
475
476         zsmltrdm2(:,:) = 0.
477         DO jl = 1, jpltrd
478            zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
479         END DO
480
481         zsmlres2(:,:) =  zsmltot2(:,:)  -       &
482              ( zsmltrdm2(:,:) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:) )
483         
484         !-- Diagnose Asselin trend over the analysis window
485         ztmlatf2(:,:) = ztmltrd2(:,:,jpmld_atf) - tmltrd_sum(:,:,jpmld_atf) + tmltrd_atf_sumb(:,:)
486         zsmlatf2(:,:) = zsmltrd2(:,:,jpmld_atf) - smltrd_sum(:,:,jpmld_atf) + smltrd_atf_sumb(:,:)
487
488         !-- Lateral boundary conditions
489         !         ... temperature ...                    ... salinity ...
490         CALL lbc_lnk( ztmltot2, 'T', 1. )  ;   CALL lbc_lnk( zsmltot2, 'T', 1. )
491         CALL lbc_lnk( ztmlres2, 'T', 1. )  ;   CALL lbc_lnk( zsmlres2, 'T', 1. )
492         DO jl = 1, jpltrd
493            CALL lbc_lnk( ztmltrd2(:,:,jl), 'T', 1. ) ! \  these will be output
494            CALL lbc_lnk( zsmltrd2(:,:,jl), 'T', 1. ) ! /  in the NetCDF trends file
495         END DO
496         
497         ! III.3 Time evolution array swap
498         ! -------------------------------
499         
500         ! For T/S instantaneous diagnostics
501         !   ... temperature ...               ... salinity ...
502         tmlbb  (:,:) = tmlb   (:,:)  ;   smlbb  (:,:) = smlb   (:,:)
503         tmlbn  (:,:) = tml    (:,:)  ;   smlbn  (:,:) = sml    (:,:)
504         tmlatfb(:,:) = tmlatfn(:,:)  ;   smlatfb(:,:) = smlatfn(:,:)
505
506         ! For T mean diagnostics
507         tmltrd_csum_ub (:,:,:) = zfn * tmltrd_sum(:,:,:) - tmltrd_csum_ln(:,:,:)
508         tml_sumb       (:,:)   = tml_sum(:,:)
509         tmltrd_atf_sumb(:,:)   = tmltrd_sum(:,:,jpmld_atf)
510         
511         ! For S mean diagnostics
512         smltrd_csum_ub (:,:,:) = zfn * smltrd_sum(:,:,:) - smltrd_csum_ln(:,:,:)
513         sml_sumb       (:,:)   = sml_sum(:,:)
514         smltrd_atf_sumb(:,:)   = smltrd_sum(:,:,jpmld_atf)
515         
516         ! ML depth
517         rmldbn         (:,:)   = rmld    (:,:)
518         
519         IF( ln_ctl ) THEN
520            IF( ln_trdmld_instant ) THEN
521               CALL prt_ctl(tab2d_1=tmlbb   , clinfo1=' tmlbb   -   : ', mask1=tmask, ovlap=1)
522               CALL prt_ctl(tab2d_1=tmlbn   , clinfo1=' tmlbn   -   : ', mask1=tmask, ovlap=1)
523               CALL prt_ctl(tab2d_1=tmlatfb , clinfo1=' tmlatfb -   : ', mask1=tmask, ovlap=1)
524            ELSE
525               CALL prt_ctl(tab2d_1=tmlbn          , clinfo1=' tmlbn           -  : ', mask1=tmask, ovlap=1)
526               CALL prt_ctl(tab2d_1=rmldbn         , clinfo1=' rmldbn          -  : ', mask1=tmask, ovlap=1)
527               CALL prt_ctl(tab2d_1=tml_sumb       , clinfo1=' tml_sumb        -  : ', mask1=tmask, ovlap=1)
528               CALL prt_ctl(tab2d_1=tmltrd_atf_sumb, clinfo1=' tmltrd_atf_sumb -  : ', mask1=tmask, ovlap=1)
529               CALL prt_ctl(tab3d_1=tmltrd_csum_ub , clinfo1=' tmltrd_csum_ub  -  : ', mask1=tmask, ovlap=1, kdim=1)
530            END IF
531         END IF
532
533         ! III.4 Convert to appropriate physical units
534         ! -------------------------------------------
535
536         !    ... temperature ...                         ... salinity ...
537         ztmltot (:,:)   = ztmltot(:,:)   * rn_ucf/zfn  ; zsmltot (:,:)   = zsmltot(:,:)   * rn_ucf/zfn
538         ztmlres (:,:)   = ztmlres(:,:)   * rn_ucf/zfn  ; zsmlres (:,:)   = zsmlres(:,:)   * rn_ucf/zfn
539         ztmlatf (:,:)   = ztmlatf(:,:)   * rn_ucf/zfn  ; zsmlatf (:,:)   = zsmlatf(:,:)   * rn_ucf/zfn
540
541         tml_sum (:,:)   = tml_sum (:,:)  /  (2*zfn) ; sml_sum (:,:)   = sml_sum (:,:)  /  (2*zfn)
542         ztmltot2(:,:)   = ztmltot2(:,:)  * rn_ucf/zfn2 ; zsmltot2(:,:)   = zsmltot2(:,:)  * rn_ucf/zfn2
543         ztmltrd2(:,:,:) = ztmltrd2(:,:,:)* rn_ucf/zfn2 ; zsmltrd2(:,:,:) = zsmltrd2(:,:,:)* rn_ucf/zfn2
544         ztmlatf2(:,:)   = ztmlatf2(:,:)  * rn_ucf/zfn2 ; zsmlatf2(:,:)   = zsmlatf2(:,:)  * rn_ucf/zfn2
545         ztmlres2(:,:)   = ztmlres2(:,:)  * rn_ucf/zfn2 ; zsmlres2(:,:)   = zsmlres2(:,:)  * rn_ucf/zfn2
546
547         rmld_sum(:,:)   = rmld_sum(:,:)  /  (2*zfn)  ! similar to tml_sum and sml_sum
548
549         ! * Debugging information *
550         IF( lldebug ) THEN
551            !
552            WRITE(numout,*)
553            WRITE(numout,*) 'trd_mld : write trends in the Mixed Layer for debugging process:'
554            WRITE(numout,*) '~~~~~~~  '
555            WRITE(numout,*) '          TRA kt = ', kt, 'nmoymltrd = ', nmoymltrd
556            WRITE(numout,*)
557            WRITE(numout,*) '          >>>>>>>>>>>>>>>>>>  TRA TEMPERATURE <<<<<<<<<<<<<<<<<<'
558            WRITE(numout,*) '          TRA ztmlres    : ', SUM(ztmlres(:,:))
559            WRITE(numout,*) '          TRA ztmltot    : ', SUM(ztmltot(:,:))
560            WRITE(numout,*) '          TRA tmltrdm    : ', SUM(tmltrdm(:,:))
561            WRITE(numout,*) '          TRA tmlatfb    : ', SUM(tmlatfb(:,:))
562            WRITE(numout,*) '          TRA tmlatfn    : ', SUM(tmlatfn(:,:))
563            DO jl = 1, jpltrd
564               WRITE(numout,*) '          * TRA TREND INDEX jpmld_xxx = jl = ', jl, &
565                    & ' tmltrd : ', SUM(tmltrd(:,:,jl))
566            END DO
567            WRITE(numout,*) '          TRA ztmlres (jpi/2,jpj/2) : ', ztmlres (jpi/2,jpj/2)
568            WRITE(numout,*) '          TRA ztmlres2(jpi/2,jpj/2) : ', ztmlres2(jpi/2,jpj/2)
569            WRITE(numout,*)
570            WRITE(numout,*) '          >>>>>>>>>>>>>>>>>>  TRA SALINITY <<<<<<<<<<<<<<<<<<'
571            WRITE(numout,*) '          TRA zsmlres    : ', SUM(zsmlres(:,:))
572            WRITE(numout,*) '          TRA zsmltot    : ', SUM(zsmltot(:,:))
573            WRITE(numout,*) '          TRA smltrdm    : ', SUM(smltrdm(:,:))
574            WRITE(numout,*) '          TRA smlatfb    : ', SUM(smlatfb(:,:))
575            WRITE(numout,*) '          TRA smlatfn    : ', SUM(smlatfn(:,:))
576            DO jl = 1, jpltrd
577               WRITE(numout,*) '          * TRA TREND INDEX jpmld_xxx = jl = ', jl, &
578                    & ' smltrd : ', SUM(smltrd(:,:,jl))
579            END DO
580            WRITE(numout,*) '          TRA zsmlres (jpi/2,jpj/2) : ', zsmlres (jpi/2,jpj/2)
581            WRITE(numout,*) '          TRA zsmlres2(jpi/2,jpj/2) : ', zsmlres2(jpi/2,jpj/2)
582            !
583         END IF
584         !
585      END IF MODULO_NTRD
586
587      ! ======================================================================
588      ! IV. Write trends in the NetCDF file
589      ! ======================================================================
590
591      ! IV.1 Code for dimg mpp output
592      ! -----------------------------
593
594#if defined key_dimgout
595
596      IF( MOD( itmod, nn_trd ) == 0 ) THEN
597         iyear =  ndastp/10000
598         imon  = (ndastp-iyear*10000)/100
599         iday  =  ndastp - imon*100 - iyear*10000
600         WRITE(clname,9000) TRIM(cexper),'MLDiags',iyear,imon,iday
601         WRITE(clmode,'(f5.1,a)') nn_trd*rdt/86400.,' days average'
602         cltext = TRIM(cexper)//' mld diags'//TRIM(clmode)
603         CALL dia_wri_dimg (clname, cltext, smltrd, jpltrd, '2')
604      END IF
605
6069000  FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")
607
608#else
609     
610      ! IV.2 Code for IOIPSL/NetCDF output
611      ! ----------------------------------
612
613      IF( lwp .AND. MOD( itmod , nn_trd ) == 0 ) THEN
614         WRITE(numout,*) ' '
615         WRITE(numout,*) 'trd_mld : write trends in the NetCDF file :'
616         WRITE(numout,*) '~~~~~~~  '
617         WRITE(numout,*) '          ', TRIM(clhstnam), ' at kt = ', kt
618         WRITE(numout,*) '          N.B. nmoymltrd = ', nmoymltrd
619         WRITE(numout,*) ' '
620      END IF
621         
622      !-- Write the trends for T/S instantaneous diagnostics
623      IF( ln_trdmld_instant ) THEN           
624
625         CALL histwrite( nidtrd, "mxl_depth", it, rmld(:,:), ndimtrd1, ndextrd1 )
626         
627         !................................. ( ML temperature ) ...................................
628         
629         !-- Output the fields
630         CALL histwrite( nidtrd, "tml"     , it, tml    (:,:), ndimtrd1, ndextrd1 ) 
631         CALL histwrite( nidtrd, "tml_tot" , it, ztmltot(:,:), ndimtrd1, ndextrd1 ) 
632         CALL histwrite( nidtrd, "tml_res" , it, ztmlres(:,:), ndimtrd1, ndextrd1 ) 
633         
634         DO jl = 1, jpltrd - 1
635            CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)),            &
636                 &          it, tmltrd (:,:,jl), ndimtrd1, ndextrd1 )
637         END DO
638         
639         CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)),        &
640              &          it, ztmlatf(:,:), ndimtrd1, ndextrd1 )
641         
642         !.................................. ( ML salinity ) .....................................
643         
644         !-- Output the fields
645         CALL histwrite( nidtrd, "sml"     , it, sml    (:,:), ndimtrd1, ndextrd1 ) 
646         CALL histwrite( nidtrd, "sml_tot" , it, zsmltot(:,:), ndimtrd1, ndextrd1 ) 
647         CALL histwrite( nidtrd, "sml_res" , it, zsmlres(:,:), ndimtrd1, ndextrd1 ) 
648         
649         DO jl = 1, jpltrd - 1
650            CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)),            &
651                 &          it, smltrd(:,:,jl), ndimtrd1, ndextrd1 )
652         END DO
653         
654         CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)),        &
655              &          it, zsmlatf(:,:), ndimtrd1, ndextrd1 )
656         
657         IF( kt == nitend )   CALL histclo( nidtrd )
658
659      !-- Write the trends for T/S mean diagnostics
660      ELSE
661         
662         CALL histwrite( nidtrd, "mxl_depth", it, rmld_sum(:,:), ndimtrd1, ndextrd1 ) 
663         
664         !................................. ( ML temperature ) ...................................
665         
666         !-- Output the fields
667         CALL histwrite( nidtrd, "tml"     , it, tml_sum (:,:), ndimtrd1, ndextrd1 ) 
668         CALL histwrite( nidtrd, "tml_tot" , it, ztmltot2(:,:), ndimtrd1, ndextrd1 ) 
669         CALL histwrite( nidtrd, "tml_res" , it, ztmlres2(:,:), ndimtrd1, ndextrd1 ) 
670         
671         DO jl = 1, jpltrd - 1
672            CALL histwrite( nidtrd, trim("tml"//ctrd(jl,2)),            &
673                 &          it, ztmltrd2(:,:,jl), ndimtrd1, ndextrd1 )
674         END DO
675         
676         CALL histwrite( nidtrd, trim("tml"//ctrd(jpmld_atf,2)),        &
677              &          it, ztmlatf2(:,:), ndimtrd1, ndextrd1 )
678         
679         !.................................. ( ML salinity ) .....................................
680                     
681         !-- Output the fields
682         CALL histwrite( nidtrd, "sml"     , it, sml_sum (:,:), ndimtrd1, ndextrd1 ) 
683         CALL histwrite( nidtrd, "sml_tot" , it, zsmltot2(:,:), ndimtrd1, ndextrd1 ) 
684         CALL histwrite( nidtrd, "sml_res" , it, zsmlres2(:,:), ndimtrd1, ndextrd1 ) 
685         
686         DO jl = 1, jpltrd - 1
687            CALL histwrite( nidtrd, trim("sml"//ctrd(jl,2)),            &
688                 &          it, zsmltrd2(:,:,jl), ndimtrd1, ndextrd1 )
689         END DO
690         
691         CALL histwrite( nidtrd, trim("sml"//ctrd(jpmld_atf,2)),        &
692              &          it, zsmlatf2(:,:), ndimtrd1, ndextrd1 )
693         
694         IF( kt == nitend )   CALL histclo( nidtrd )
695
696      END IF
697     
698      ! Compute the control surface (for next time step) : flag = on
699      icount = 1
700      !
701#endif
702
703      IF( MOD( itmod, nn_trd ) == 0 ) THEN
704         !
705         ! III.5 Reset cumulative arrays to zero
706         ! -------------------------------------
707         nmoymltrd = 0
708         
709         !   ... temperature ...               ... salinity ...
710         tmltrdm        (:,:)   = 0.e0   ;   smltrdm        (:,:)   = 0.e0
711         tmlatfm        (:,:)   = 0.e0   ;   smlatfm        (:,:)   = 0.e0
712         tml_sum        (:,:)   = 0.e0   ;   sml_sum        (:,:)   = 0.e0
713         tmltrd_csum_ln (:,:,:) = 0.e0   ;   smltrd_csum_ln (:,:,:) = 0.e0
714         tmltrd_sum     (:,:,:) = 0.e0   ;   smltrd_sum     (:,:,:) = 0.e0
715
716         rmld_sum       (:,:)   = 0.e0
717         !
718      END IF
719
720      ! ======================================================================
721      ! V. Write restart file
722      ! ======================================================================
723
724      IF( lrst_oce )   CALL trd_mld_rst_write( kt ) 
725
726      CALL wrk_dealloc( jpi, jpj,         ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf                        )
727      CALL wrk_dealloc( jpi, jpj,         ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2 ) 
728      CALL wrk_dealloc( jpi, jpj, jpltrd, ztmltrd2, zsmltrd2                                                               )
729      !
730   END SUBROUTINE trd_mld
731
732
733   SUBROUTINE trd_mld_init
734      !!----------------------------------------------------------------------
735      !!                  ***  ROUTINE trd_mld_init  ***
736      !!
737      !! ** Purpose :   computation of vertically integrated T and S budgets
738      !!      from ocean surface down to control surface (NetCDF output)
739      !!----------------------------------------------------------------------
740      INTEGER :: jl
741      INTEGER :: inum   ! logical unit
742      REAL(wp) ::   zjulian, zsto, zout
743      CHARACTER (LEN=40) ::   clop
744      CHARACTER (LEN=12) ::   clmxl, cltu, clsu
745      !!----------------------------------------------------------------------
746
747      ! ======================================================================
748      ! I. initialization
749      ! ======================================================================
750
751      IF(lwp) THEN
752         WRITE(numout,*)
753         WRITE(numout,*) ' trd_mld_init : Mixed-layer trends'
754         WRITE(numout,*) ' ~~~~~~~~~~~~~'
755         WRITE(numout,*) '                namelist namtrd read in trd_mod_init                        '
756         WRITE(numout,*)
757      END IF
758
759      ! I.1 Check consistency of user defined preferences
760      ! -------------------------------------------------
761
762      IF( ( lk_trdmld ) .AND. ( MOD( nitend-nit000+1, nn_trd ) /= 0 ) ) THEN
763         WRITE(numout,cform_err)
764         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
765         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
766         WRITE(numout,*) '                          you defined, nn_trd   = ', nn_trd
767         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
768         WRITE(numout,*) '                You should reconsider this choice.                        ' 
769         WRITE(numout,*) 
770         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
771         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
772         nstop = nstop + 1
773      END IF
774
775      IF( ( lk_trdmld ) .AND. ( nn_cla == 1 ) ) THEN
776         WRITE(numout,cform_war)
777         WRITE(numout,*) '                You set n_cla = 1. Note that the Mixed-Layer diagnostics  '
778         WRITE(numout,*) '                are not exact along the corresponding straits.            '
779         nwarn = nwarn + 1
780      END IF
781
782      !                                   ! allocate trdmld arrays
783      IF( trd_mld_alloc()    /= 0 )   CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld     arrays' )
784      IF( trdmld_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_mld_init : unable to allocate trdmld_oce arrays' )
785
786      ! I.2 Initialize arrays to zero or read a restart file
787      ! ----------------------------------------------------
788
789      nmoymltrd = 0
790
791      !     ... temperature ...                  ... salinity ...
792      tml            (:,:)   = 0.e0    ;    sml            (:,:)   = 0.e0     ! inst.
793      tmltrdm        (:,:)   = 0.e0    ;    smltrdm        (:,:)   = 0.e0
794      tmlatfm        (:,:)   = 0.e0    ;    smlatfm        (:,:)   = 0.e0
795      tml_sum        (:,:)   = 0.e0    ;    sml_sum        (:,:)   = 0.e0     ! mean
796      tmltrd_sum     (:,:,:) = 0.e0    ;    smltrd_sum     (:,:,:) = 0.e0
797      tmltrd_csum_ln (:,:,:) = 0.e0    ;    smltrd_csum_ln (:,:,:) = 0.e0
798
799      rmld           (:,:)   = 0.e0           
800      rmld_sum       (:,:)   = 0.e0
801
802      IF( ln_rstart .AND. ln_trdmld_restart ) THEN
803         CALL trd_mld_rst_read
804      ELSE
805         !     ... temperature ...                  ... salinity ...
806         tmlb           (:,:)   = 0.e0    ;    smlb           (:,:)   = 0.e0  ! inst.
807         tmlbb          (:,:)   = 0.e0    ;    smlbb          (:,:)   = 0.e0 
808         tmlbn          (:,:)   = 0.e0    ;    smlbn          (:,:)   = 0.e0 
809         tml_sumb       (:,:)   = 0.e0    ;    sml_sumb       (:,:)   = 0.e0  ! mean
810         tmltrd_csum_ub (:,:,:) = 0.e0    ;    smltrd_csum_ub (:,:,:) = 0.e0
811         tmltrd_atf_sumb(:,:)   = 0.e0    ;    smltrd_atf_sumb(:,:)   = 0.e0 
812      END IF
813
814      icount = 1   ;   ionce  = 1                            ! open specifier
815
816      ! I.3 Read control surface from file ctlsurf_idx
817      ! ----------------------------------------------
818 
819      IF( nn_ctls == 1 ) THEN
820         CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp )
821         READ ( inum ) nbol
822         CLOSE( inum )
823      END IF
824
825      ! ======================================================================
826      ! II. netCDF output initialization
827      ! ======================================================================
828
829#if defined key_dimgout 
830      ???
831#else
832      ! clmxl = legend root for netCDF output
833      IF( nn_ctls == 0 ) THEN      ! control surface = mixed-layer with density criterion
834         clmxl = 'Mixed Layer '  !                   (array nmln computed in zdfmxl.F90)
835      ELSE IF( nn_ctls == 1 ) THEN ! control surface = read index from file
836         clmxl = '      Bowl '
837      ELSE IF( nn_ctls >= 2 ) THEN ! control surface = model level
838         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls
839      END IF
840
841      ! II.1 Define frequency of output and means
842      ! -----------------------------------------
843      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
844      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
845      ENDIF
846#  if defined key_diainstant
847      IF( .NOT. ln_trdmld_instant ) THEN
848         CALL ctl_stop( 'trd_mld : this was never checked. Comment this line to proceed...' )
849      END IF
850      zsto = nn_trd * rdt
851      clop = "inst("//TRIM(clop)//")"
852#  else
853      IF( ln_trdmld_instant ) THEN
854         zsto = rdt                 ! inst. diags : we use IOIPSL time averaging
855      ELSE
856         zsto = nn_trd * rdt          ! mean  diags : we DO NOT use any IOIPSL time averaging
857      END IF
858      clop = "ave("//TRIM(clop)//")"
859#  endif
860      zout = nn_trd * rdt
861
862      IF(lwp) WRITE (numout,*) '                netCDF initialization'
863
864      ! II.2 Compute julian date from starting date of the run
865      ! ------------------------------------------------------
866      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
867      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
868      IF(lwp) WRITE(numout,*)' ' 
869      IF(lwp) WRITE(numout,*)'                Date 0 used :',nit000,    &
870         &                   ' YEAR ', nyear,' MONTH '      , nmonth,   &
871         &                   ' DAY ' , nday, 'Julian day : ', zjulian
872
873
874      ! II.3 Define the T grid trend file (nidtrd)
875      ! ------------------------------------------
876      !-- Define long and short names for the NetCDF output variables
877      !       ==> choose them according to trdmld_oce.F90 <==
878
879      ctrd(jpmld_xad,1) = " Zonal advection"                  ;   ctrd(jpmld_xad,2) = "_xad"
880      ctrd(jpmld_yad,1) = " Meridional advection"             ;   ctrd(jpmld_yad,2) = "_yad"
881      ctrd(jpmld_zad,1) = " Vertical advection"               ;   ctrd(jpmld_zad,2) = "_zad"
882      ctrd(jpmld_ldf,1) = " Lateral diffusion"                ;   ctrd(jpmld_ldf,2) = "_ldf"
883      ctrd(jpmld_for,1) = " Forcing"                          ;   ctrd(jpmld_for,2) = "_for"
884      ctrd(jpmld_zdf,1) = " Vertical diff. (Kz)"              ;   ctrd(jpmld_zdf,2) = "_zdf"
885      ctrd(jpmld_bbc,1) = " Geothermal flux"                  ;   ctrd(jpmld_bbc,2) = "_bbc"
886      ctrd(jpmld_bbl,1) = " Adv/diff. Bottom boundary layer"  ;   ctrd(jpmld_bbl,2) = "_bbl"
887      ctrd(jpmld_dmp,1) = " Tracer damping"                   ;   ctrd(jpmld_dmp,2) = "_dmp"
888      ctrd(jpmld_npc,1) = " Non penetrative convec. adjust."  ;   ctrd(jpmld_npc,2) = "_npc"
889      ctrd(jpmld_atf,1) = " Asselin time filter"              ;   ctrd(jpmld_atf,2) = "_atf"
890                                                                 
891      !-- Create a NetCDF file and enter the define mode
892      CALL dia_nam( clhstnam, nn_trd, 'trends' )
893      IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam
894      CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
895      &             1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_t, nidtrd, domain_id=nidom, snc4chunks=snc4set )
896
897      !-- Define the ML depth variable
898      CALL histdef(nidtrd, "mxl_depth", clmxl//"  Mixed Layer Depth"              , "m",         &
899                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
900
901      !-- Define physical units
902      IF     ( rn_ucf == 1.       ) THEN   ;   cltu = "degC/s"     ;   clsu = "p.s.u./s"
903      ELSEIF ( rn_ucf == 3600.*24.) THEN   ;   cltu = "degC/day"   ;   clsu = "p.s.u./day"
904      ELSE                                 ;   cltu = "unknown?"   ;   clsu = "unknown?"
905      END IF
906
907
908      !-- Define miscellaneous T and S mixed-layer variables
909
910      IF( jpltrd /= jpmld_atf ) CALL ctl_stop( 'Error : jpltrd /= jpmld_atf' ) ! see below
911
912      !................................. ( ML temperature ) ...................................
913
914      CALL histdef(nidtrd, "tml"      , clmxl//" T Mixed Layer Temperature"       ,  "C",        &
915                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
916      CALL histdef(nidtrd, "tml_tot",   clmxl//" T Total trend"                   , cltu,        & 
917                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )             
918      CALL histdef(nidtrd, "tml_res",   clmxl//" T dh/dt Entrainment (Resid.)"    , cltu,        & 
919                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
920     
921      DO jl = 1, jpltrd - 1      ! <== only true if jpltrd == jpmld_atf
922         CALL histdef(nidtrd, trim("tml"//ctrd(jl,2)), clmxl//" T"//ctrd(jl,1), cltu,            & 
923                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
924      END DO                                                                 ! if zsto=rdt above
925     
926      CALL histdef(nidtrd, trim("tml"//ctrd(jpmld_atf,2)), clmxl//" T"//ctrd(jpmld_atf,1), cltu, & 
927                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
928     
929      !.................................. ( ML salinity ) .....................................
930     
931      CALL histdef(nidtrd, "sml"      , clmxl//" S Mixed Layer Salinity"          , "p.s.u.",       &
932                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
933      CALL histdef(nidtrd, "sml_tot",   clmxl//" S Total trend"                   , clsu,        & 
934                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )             
935      CALL histdef(nidtrd, "sml_res",   clmxl//" S dh/dt Entrainment (Resid.)"    , clsu,        & 
936                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
937     
938      DO jl = 1, jpltrd - 1      ! <== only true if jpltrd == jpmld_atf
939         CALL histdef(nidtrd, trim("sml"//ctrd(jl,2)), clmxl//" S"//ctrd(jl,1), clsu,            & 
940                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
941      END DO                                                                 ! if zsto=rdt above
942     
943      CALL histdef(nidtrd, trim("sml"//ctrd(jpmld_atf,2)), clmxl//" S"//ctrd(jpmld_atf,1), clsu, & 
944                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
945
946      !-- Leave IOIPSL/NetCDF define mode
947      CALL histend( nidtrd, snc4set )
948
949#endif        /* key_dimgout */
950   END SUBROUTINE trd_mld_init
951
952#else
953   !!----------------------------------------------------------------------
954   !!   Default option :                                       Empty module
955   !!----------------------------------------------------------------------
956CONTAINS
957   SUBROUTINE trd_mld( kt )             ! Empty routine
958      INTEGER, INTENT( in) ::   kt
959      WRITE(*,*) 'trd_mld: You should not have seen this print! error?', kt
960   END SUBROUTINE trd_mld
961   SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
962      REAL, DIMENSION(:,:,:), INTENT( in ) ::   &
963         pttrdmld, pstrdmld                   ! Temperature and Salinity trends
964      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index
965      CHARACTER(len=2), INTENT( in ) ::   & 
966         ctype                                ! surface/bottom (2D arrays) or
967         !                                    ! interior (3D arrays) physics
968      WRITE(*,*) 'trd_mld_zint: You should not have seen this print! error?', pttrdmld(1,1,1)
969      WRITE(*,*) '  "      "  : You should not have seen this print! error?', pstrdmld(1,1,1)
970      WRITE(*,*) '  "      "  : You should not have seen this print! error?', ctype
971      WRITE(*,*) '  "      "  : You should not have seen this print! error?', ktrd
972   END SUBROUTINE trd_mld_zint
973   SUBROUTINE trd_mld_init              ! Empty routine
974      WRITE(*,*) 'trd_mld_init: You should not have seen this print! error?'
975   END SUBROUTINE trd_mld_init
976#endif
977
978   !!======================================================================
979END MODULE trdmld
Note: See TracBrowser for help on using the repository browser.