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 tags/nemo_v1_13_dev7/NEMO/OPA_SRC/TRD – NEMO

source: tags/nemo_v1_13_dev7/NEMO/OPA_SRC/TRD/trdmld.F90 @ 1596

Last change on this file since 1596 was 503, checked in by opalod, 18 years ago

nemo_v1_update_064 : CT : general trends update including the addition of mean windows analysis possibility in the mixed layer

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