New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trdmld.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 13 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

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