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

source: branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago

FASTNEt code modifications

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