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

source: trunk/NEMO/OPA_SRC/TRD/trdmld.F90 @ 526

Last change on this file since 526 was 521, checked in by opalod, 18 years ago

nemo_v1_update_73 : CT : build Mixed Layer restart files using iom

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 46.7 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      ! 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 = .FALSE.
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      IF( ln_traldf_iso ) THEN
256         DO jj = 1,jpj
257            DO ji = 1,jpi
258               ik = nmld(ji,jj)
259               zavt = avt(ji,jj,ik)
260               tmltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
261                  &                      * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) )         &
262                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
263               zavt = fsavs(ji,jj,ik)
264               smltrd(ji,jj,jpmld_zdf) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)  &
265                  &                      * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) )         &
266                  &                      / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
267            END DO
268         END DO
269
270         ! ... Remove this K_z trend from the iso-neutral diffusion term (if any)
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      ! ======================================================================
701      ! V. Write restart file
702      ! ======================================================================
703
704      CALL trd_mld_rst_write( kt ) 
705
706   END SUBROUTINE trd_mld
707
708
709
710   SUBROUTINE trd_mld_init
711      !!----------------------------------------------------------------------
712      !!                  ***  ROUTINE trd_mld_init  ***
713      !!
714      !! ** Purpose :   computation of vertically integrated T and S budgets
715      !!      from ocean surface down to control surface (NetCDF output)
716      !!
717      !!----------------------------------------------------------------------
718      !! * Local declarations
719      INTEGER :: ilseq, jl
720
721      REAL(wp) ::   zjulian, zsto, zout
722
723      CHARACTER (LEN=21) ::    &
724         clold ='OLD'        , & ! open specifier (direct access files)
725         clunf ='UNFORMATTED', & ! open specifier (direct access files)
726         clseq ='SEQUENTIAL'     ! open specifier (direct access files)
727      CHARACTER (LEN=40) ::   clop
728      CHARACTER (LEN=12) ::   clmxl, cltu, clsu
729
730      !!----------------------------------------------------------------------
731
732      ! ======================================================================
733      ! I. initialization
734      ! ======================================================================
735
736      IF(lwp) THEN
737         WRITE(numout,*)
738         WRITE(numout,*) ' trd_mld_init : Mixed-layer trends'
739         WRITE(numout,*) ' ~~~~~~~~~~~~~'
740         WRITE(numout,*) '                namelist namtrd read in trd_mod_init                        '
741         WRITE(numout,*)
742      END IF
743
744      ! I.1 Check consistency of user defined preferences
745      ! -------------------------------------------------
746
747      IF( ( lk_trdmld ) .AND. ( MOD( nitend, ntrd ) /= 0 ) ) THEN
748         WRITE(numout,cform_err)
749         WRITE(numout,*) '                Your nitend parameter, nitend = ', nitend
750         WRITE(numout,*) '                is no multiple of the trends diagnostics frequency        '
751         WRITE(numout,*) '                          you defined, ntrd   = ', ntrd
752         WRITE(numout,*) '                This will not allow you to restart from this simulation.  '
753         WRITE(numout,*) '                You should reconsider this choice.                        ' 
754         WRITE(numout,*) 
755         WRITE(numout,*) '                N.B. the nitend parameter is also constrained to be a     '
756         WRITE(numout,*) '                multiple of the sea-ice frequency parameter (typically 5) '
757         nstop = nstop + 1
758      END IF
759
760      IF( ( lk_trdmld ) .AND. ( n_cla == 1 ) ) THEN
761         WRITE(numout,cform_war)
762         WRITE(numout,*) '                You set n_cla = 1. Note that the Mixed-Layer diagnostics  '
763         WRITE(numout,*) '                are not exact along the corresponding straits.            '
764         nwarn = nwarn + 1
765      END IF
766
767      ! I.2 Initialize arrays to zero or read a restart file
768      ! ----------------------------------------------------
769
770      nmoymltrd = 0
771
772      !     ... temperature ...                  ... salinity ...
773      tml            (:,:)   = 0.e0    ;    sml            (:,:)   = 0.e0     ! inst.
774      tmltrdm        (:,:)   = 0.e0    ;    smltrdm        (:,:)   = 0.e0
775      tmlatfm        (:,:)   = 0.e0    ;    smlatfm        (:,:)   = 0.e0
776      tml_sum        (:,:)   = 0.e0    ;    sml_sum        (:,:)   = 0.e0     ! mean
777      tmltrd_sum     (:,:,:) = 0.e0    ;    smltrd_sum     (:,:,:) = 0.e0
778      tmltrd_csum_ln (:,:,:) = 0.e0    ;    smltrd_csum_ln (:,:,:) = 0.e0
779
780      rmld           (:,:)   = 0.e0           
781      rmld_sum       (:,:)   = 0.e0
782
783      IF( ln_rstart .AND. ln_trdmld_restart ) THEN
784         CALL trd_mld_rst_read
785      ELSE
786         !     ... temperature ...                  ... salinity ...
787         tmlb           (:,:)   = 0.e0    ;    smlb           (:,:)   = 0.e0  ! inst.
788         tmlbb          (:,:)   = 0.e0    ;    smlbb          (:,:)   = 0.e0 
789         tmlbn          (:,:)   = 0.e0    ;    smlbn          (:,:)   = 0.e0 
790         tml_sumb       (:,:)   = 0.e0    ;    sml_sumb       (:,:)   = 0.e0  ! mean
791         tmltrd_csum_ub (:,:,:) = 0.e0    ;    smltrd_csum_ub (:,:,:) = 0.e0
792         tmltrd_atf_sumb(:,:)   = 0.e0    ;    smltrd_atf_sumb(:,:)   = 0.e0 
793      END IF
794
795      ilseq  = 1   ;   icount = 1   ;   ionce  = 1                            ! open specifier
796
797      ! I.3 Read control surface from file ctlsurf_idx
798      ! ----------------------------------------------
799 
800      IF( nctls == 1 ) THEN
801         clname = 'ctlsurf_idx'
802         CALL ctlopn( numbol, clname, clold, clunf, clseq, ilseq, numout, lwp, 1 )
803         REWIND( numbol )
804         READ  ( numbol ) nbol
805      END IF
806
807      ! ======================================================================
808      ! II. netCDF output initialization
809      ! ======================================================================
810
811#if defined key_dimgout 
812      ???
813#else
814      ! clmxl = legend root for netCDF output
815      IF( nctls == 0 ) THEN      ! control surface = mixed-layer with density criterion
816         clmxl = 'Mixed Layer '  !                   (array nmln computed in zdfmxl.F90)
817      ELSE IF( nctls == 1 ) THEN ! control surface = read index from file
818         clmxl = '      Bowl '
819      ELSE IF( nctls >= 2 ) THEN ! control surface = model level
820         WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nctls
821      END IF
822
823      ! II.1 Define frequency of output and means
824      ! -----------------------------------------
825#  if defined key_diainstant
826      IF( .NOT. ln_trdmld_instant ) THEN
827         CALL ctl_stop( 'trd_mld : this was never checked. Comment this line to proceed...' )
828      END IF
829      zsto = ntrd * rdt
830      clop ="inst(only(x))"
831#  else
832      IF( ln_trdmld_instant ) THEN
833         zsto = rdt                 ! inst. diags : we use IOIPSL time averaging
834      ELSE
835         zsto = ntrd * rdt          ! mean  diags : we DO NOT use any IOIPSL time averaging
836      END IF
837      clop ="ave(only(x))"
838#  endif
839      zout = ntrd * rdt
840
841      IF(lwp) WRITE (numout,*) '                netCDF initialization'
842
843      ! II.2 Compute julian date from starting date of the run
844      ! ------------------------------------------------------
845      CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
846      IF(lwp) WRITE(numout,*)' ' 
847      IF(lwp) WRITE(numout,*)'                Date 0 used :',nit000,    &
848         &                   ' YEAR ', nyear,' MONTH '      , nmonth,   &
849         &                   ' DAY ' , nday, 'Julian day : ', zjulian
850
851
852      ! II.3 Define the T grid trend file (nidtrd)
853      ! ------------------------------------------
854      !-- Define long and short names for the NetCDF output variables
855      !       ==> choose them according to trdmld_oce.F90 <==
856
857      ctrd(jpmld_xad,1) = " Zonal advection"                  ;   ctrd(jpmld_xad,2) = "_xad"
858      ctrd(jpmld_yad,1) = " Meridional advection"             ;   ctrd(jpmld_yad,2) = "_yad"
859      ctrd(jpmld_zad,1) = " Vertical advection"               ;   ctrd(jpmld_zad,2) = "_zad"
860      ctrd(jpmld_ldf,1) = " Lateral diffusion"                ;   ctrd(jpmld_ldf,2) = "_ldf"
861      ctrd(jpmld_for,1) = " Forcing"                          ;   ctrd(jpmld_for,2) = "_for"
862      ctrd(jpmld_zdf,1) = " Vertical diff. (Kz)"              ;   ctrd(jpmld_zdf,2) = "_zdf"
863      ctrd(jpmld_bbc,1) = " Geothermal flux"                  ;   ctrd(jpmld_bbc,2) = "_bbc"
864      ctrd(jpmld_bbl,1) = " Adv/diff. Bottom boundary layer"  ;   ctrd(jpmld_bbl,2) = "_bbl"
865      ctrd(jpmld_dmp,1) = " Tracer damping"                   ;   ctrd(jpmld_dmp,2) = "_dmp"
866      ctrd(jpmld_npc,1) = " Non penetrative convec. adjust."  ;   ctrd(jpmld_npc,2) = "_npc"
867      ctrd(jpmld_atf,1) = " Asselin time filter"              ;   ctrd(jpmld_atf,2) = "_atf"
868                                                                 
869      !-- Create a NetCDF file and enter the define mode
870      CALL dia_nam( clhstnam, ntrd, 'trends' )
871      IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam
872      CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,                                            &
873      &             1, jpi, 1, jpj, 0, zjulian, rdt, nh_t, nidtrd, domain_id=nidom )
874
875      !-- Define the ML depth variable
876      CALL histdef(nidtrd, "mxl_depth", clmxl//"  Mixed Layer Depth"              , "m",         &
877                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
878
879      !-- Define physical units
880      IF( ucf == 1. ) THEN
881         cltu = "degC/s"     ;   clsu = "p.s.u./s"
882      ELSEIF ( ucf == 3600.*24.) THEN
883         cltu = "degC/day"   ;   clsu = "p.s.u./day"
884      ELSE
885         cltu = "unknown?"   ;   clsu = "unknown?"
886      END IF
887
888      !-- Define miscellaneous T and S mixed-layer variables
889
890      IF( jpltrd /= jpmld_atf ) CALL ctl_stop( 'Error : jpltrd /= jpmld_atf' ) ! see below
891
892      !................................. ( ML temperature ) ...................................
893
894      CALL histdef(nidtrd, "tml"      , clmxl//" T Mixed Layer Temperature"       ,  "C",        &
895                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
896      CALL histdef(nidtrd, "tml_tot",   clmxl//" T Total trend"                   , cltu,        & 
897                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )             
898      CALL histdef(nidtrd, "tml_res",   clmxl//" T dh/dt Entrainment (Resid.)"    , cltu,        & 
899                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
900     
901      DO jl = 1, jpltrd - 1      ! <== only true if jpltrd == jpmld_atf
902         CALL histdef(nidtrd, trim("tml"//ctrd(jl,2)), clmxl//" T"//ctrd(jl,1), cltu,            & 
903                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
904      END DO                                                                 ! if zsto=rdt above
905     
906      CALL histdef(nidtrd, trim("tml"//ctrd(jpmld_atf,2)), clmxl//" T"//ctrd(jpmld_atf,1), cltu, & 
907                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
908     
909      !.................................. ( ML salinity ) .....................................
910     
911      CALL histdef(nidtrd, "sml"      , clmxl//" S Mixed Layer Salinity"          , "p.s.u.",       &
912                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )           
913      CALL histdef(nidtrd, "sml_tot",   clmxl//" S Total trend"                   , clsu,        & 
914                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )             
915      CALL histdef(nidtrd, "sml_res",   clmxl//" S dh/dt Entrainment (Resid.)"    , clsu,        & 
916                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )                   
917     
918      DO jl = 1, jpltrd - 1      ! <== only true if jpltrd == jpmld_atf
919         CALL histdef(nidtrd, trim("sml"//ctrd(jl,2)), clmxl//" S"//ctrd(jl,1), clsu,            & 
920                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) ! IOIPSL: time mean
921      END DO                                                                 ! if zsto=rdt above
922     
923      CALL histdef(nidtrd, trim("sml"//ctrd(jpmld_atf,2)), clmxl//" S"//ctrd(jpmld_atf,1), clsu, & 
924                   jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean
925
926      !-- Leave IOIPSL/NetCDF define mode
927      CALL histend( nidtrd )
928
929#endif        /* key_dimgout */
930   END SUBROUTINE trd_mld_init
931
932#else
933   !!----------------------------------------------------------------------
934   !!   Default option :                                       Empty module
935   !!----------------------------------------------------------------------
936CONTAINS
937   SUBROUTINE trd_mld( kt )             ! Empty routine
938      INTEGER, INTENT( in) ::   kt
939      WRITE(*,*) 'trd_mld: You should not have seen this print! error?', kt
940   END SUBROUTINE trd_mld
941   SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
942      REAL, DIMENSION(:,:,:), INTENT( in ) ::   &
943         pttrdmld, pstrdmld                   ! Temperature and Salinity trends
944      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index
945      CHARACTER(len=2), INTENT( in ) ::   & 
946         ctype                                ! surface/bottom (2D arrays) or
947         !                                    ! interior (3D arrays) physics
948      WRITE(*,*) 'trd_mld_zint: You should not have seen this print! error?', pttrdmld(1,1,1)
949      WRITE(*,*) '  "      "  : You should not have seen this print! error?', pstrdmld(1,1,1)
950      WRITE(*,*) '  "      "  : You should not have seen this print! error?', ctype
951      WRITE(*,*) '  "      "  : You should not have seen this print! error?', ktrd
952   END SUBROUTINE trd_mld_zint
953   SUBROUTINE trd_mld_init              ! Empty routine
954      WRITE(*,*) 'trd_mld_init: You should not have seen this print! error?'
955   END SUBROUTINE trd_mld_init
956#endif
957
958   !!======================================================================
959END MODULE trdmld
Note: See TracBrowser for help on using the repository browser.