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

Last change on this file since 1317 was 1317, checked in by smasson, 15 years ago

nwrite = modulo referenced to nit000 in all ouputs, see ticket:339

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