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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

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