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 @ 558

Last change on this file since 558 was 557, checked in by opalod, 18 years ago

nemo_v1_bugfix_069: SM+CT+CE: bugfix of mld restart + OFF line compatibiblity

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