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

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

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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.3 KB
Line 
1MODULE trdmld
2   !!======================================================================
3   !!                       ***  MODULE  trdmld  ***
4   !! Ocean diagnostics:  mixed layer T-S trends
5   !!=====================================================================
6#if   defined key_trdmld   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_trdmld'                          mixed layer trend diagnostics
9   !!----------------------------------------------------------------------
10   !!   trd_mld          : T and S cumulated trends averaged over the mixed layer
11   !!   trd_mld_zint     : T and S trends vertical integration
12   !!   trd_mld_init     : initialization step
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers variables
16   USE dom_oce         ! ocean space and time domain variables
17   USE trdmod_oce      ! ocean variables trends
18   USE ldftra_oce      ! ocean active tracers lateral physics
19   USE zdf_oce         ! ocean vertical physics
20   USE in_out_manager  ! I/O manager
21   USE phycst          ! Define parameters for the routines
22   USE daymod          ! calendar
23   USE dianam          ! build the name of file (routine)
24   USE ldfslp          ! iso-neutral slopes
25   USE zdfmxl          ! mixed layer depth
26   USE zdfddm          ! ocean vertical physics: double diffusion
27   USE ioipsl          ! NetCDF library
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
29   USE diadimg         ! dimg direct access file format output
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Accessibility
35   PUBLIC trd_mld        ! routine called by step.F90
36   PUBLIC trd_mld_init   ! routine called by opa.F90
37   PUBLIC trd_mld_zint   ! routine called by tracers routines
38
39   !! * Shared module variables
40   LOGICAL, PUBLIC, PARAMETER ::   lk_trdmld = .TRUE.    !: momentum trend flag
41
42   !! * Module variables
43   INTEGER ::   &
44      nh_t, nmoymltrd,             &  ! ???
45      nidtrd,                      &
46      ndextrd1(jpi*jpj),           &
47      ndimtrd1
48   INTEGER, SAVE ::   &
49      ionce, icount,               &
50      idebug                          ! (0/1) set it to 1 in case of problem to have more print
51
52   INTEGER, DIMENSION(jpi,jpj) ::   &
53      nmld,                         & ! mixed layer depth
54      nbol               
55
56   REAL(wp), DIMENSION(jpi,jpj) ::   &
57      rmld   ,          &  ! mld depth (m) corresponding to nmld
58      tml    , sml  ,   &  ! average T and S over mixed layer
59      tmlb   , smlb ,   &  ! before tml and sml (kt-1)
60      tmlbb  , smlbb,   &  ! tml and sml at begining of the nwrite-1
61      !                    ! timestep averaging period
62      tmlbn  , smlbn,   &  ! after tml and sml at time step after the
63      !                    ! begining of the NWRITE-1 timesteps
64      tmltrdm, smltrdm     !
65
66   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
67      tmltrd ,          &  ! total cumulative trends of temperature and
68      smltrd ,          &  ! salinity over nwrite-1 time steps
69      wkx
70
71   CHARACTER(LEN=80) :: clname
72
73   !! * Substitutions
74#  include "domzgr_substitute.h90"
75#  include "ldftra_substitute.h90"
76#  include "zdfddm_substitute.h90"
77   !!----------------------------------------------------------------------
78   !!   OPA 9.0 , LOCEAN-IPSL (2005)
79   !! $Header$
80   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
81   !!----------------------------------------------------------------------
82
83CONTAINS
84
85SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
86      !!----------------------------------------------------------------------
87      !!                  ***  ROUTINE trd_mld_zint  ***
88      !!
89      !! ** Purpose :   computation of vertically integrated T and S budgets
90      !!                from ocean surface down to control surface
91      !!
92      !! ** Method/usage :
93      !!      integration done over nwrite-1 time steps
94      !!      Control surface can be either a mixed layer depth (time varying)
95      !!      or a fixed surface (jk level or bowl).
96      !!      Choose control surface with nctls in namelist NAMDIA.
97      !!      nctls = 0  : use mixed layer with density criterion
98      !!      nctls = 1  : read index from file 'ctlsurf_idx'
99      !!      nctls > 1  : use fixed level surface jk = nctls
100      !!      Note: in the remainder of the routine, the volume between the
101      !!            surface and the control surface is called "mixed-layer"
102      !!      Method check : if the control surface is fixed, the residual dh/dt
103      !!                     entrainment should be zero
104      !!
105      !! ** Action :
106      !!            /commld/   : rmld         mld depth corresponding to nmld
107      !!                         tml          average T over mixed layer
108      !!                         tmlb         tml at kt-1
109      !!                         tmlbb        tml at begining of the NWRITE-1
110      !!                                      time steps averaging period
111      !!                         tmlbn        tml at time step after the
112      !!                                      begining of the NWRITE-1 time
113      !!                                      steps averaging period
114      !!
115      !!                  mixed layer trends :
116      !!
117      !!                  tmltrd (,,1) = zonal advection
118      !!                  tmltrd (,,2) = meridional advection
119      !!                  tmltrd (,,3) = vertical advection
120      !!                  tmltrd (,,4) = lateral diffusion (horiz. component+Beckman)
121      !!                  tmltrd (,,5) = forcing
122      !!                  tmltrd (,,6) = entrainment due to vertical diffusion (TKE)
123      !!          if iso  tmltrd (,,7) = lateral diffusion (vertical component)
124      !!                  tmltrd (,,8) = eddy induced zonal advection
125      !!                  tmltrd (,,9) = eddy induced meridional advection
126      !!                  tmltrd (,,10) = eddy induced vertical advection
127      !!
128      !!                  tmltrdm(,) : total cumulative trends over nwrite-1 time steps
129      !!                  ztmltot(,) : dT/dt over the NWRITE-1 time steps
130      !!                               averaging period (including Asselin
131      !!                               terms)
132      !!                  ztmlres(,) : residual = dh/dt entrainment
133      !!
134      !!      trends output in netCDF format using ioipsl
135      !!
136      !! History :
137      !!        !  95-04  (J. Vialard)  Original code
138      !!        !  97-02  (E. Guilyardi)  Adaptation global + base cmo
139      !!        !  99-09  (E. Guilyardi)  Re-writing + netCDF output
140      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
141      !!   9.0  !  04-08  (C. Talandier) New trends organization
142      !!----------------------------------------------------------------------
143      !! * Arguments
144      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index
145
146      CHARACTER(len=2), INTENT( in ) ::   &
147         ctype                                ! surface/bottom (2D arrays) or
148                                              ! interior (3D arrays) physics
149
150      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
151         pttrdmld,                         &  ! Temperature trend
152         pstrdmld                             ! Salinity    trend
153
154      !! * Local declarations
155      INTEGER ::   ji, jj, jk, isum
156# if defined key_trabbl_dif
157      INTEGER ::   ikb
158# endif
159
160      REAL(wp), DIMENSION(jpi,jpj) ::   &
161         zvlmsk
162      !!----------------------------------------------------------------------
163
164      IF( icount == 1 ) THEN       
165
166         zvlmsk(:,:)   = 0.e0
167         tmltrd(:,:,:) = 0.e0
168         smltrd(:,:,:) = 0.e0
169         
170         ! This computation should be done only once per time step
171
172         !  ========================================================
173         !   I. definition of control surface and associated fields
174         !  ========================================================
175
176         !    I.1 set nmld(ji,jj) = index of first T point below control surface
177         !    -------------------                       or outside mixed-layer
178
179         IF( nctls == 0 ) THEN
180            ! control surface = mixed-layer with density criterion
181            ! (array nmln computed in zdfmxl.F90)
182            nmld(:,:) = nmln(:,:)
183         ELSE IF( nctls == 1 ) THEN
184            ! control surface = read index from file
185            nmld(:,:) = nbol(:,:)
186         ELSE IF( nctls >= 2 ) THEN
187            ! control surface = model level
188            nctls = MIN( nctls, jpktrd - 1 )
189            nmld(:,:) = nctls + 1
190         ENDIF
191
192         IF( ionce == 1 ) THEN  ! compute ndextrd1 and ndimtrd1 only once
193            ! Check of validity : nmld(ji,jj) =< jpktrd
194            isum = 0
195
196            IF( jpktrd < jpk ) THEN
197               DO jj = 1, jpj
198                  DO ji = 1, jpi
199                     IF( nmld(ji,jj) <= jpktrd ) THEN
200                        zvlmsk(ji,jj) = tmask(ji,jj,1)
201                     ELSE
202                        isum = isum + 1
203                        zvlmsk(ji,jj) = 0.
204                     ENDIF
205                  END DO
206               END DO
207            ENDIF
208
209            ! Index of ocean points (2D only)
210            IF( isum > 0 ) THEN
211               WRITE(numout,*)' Number of invalid points nmld > jpktrd', isum 
212               CALL wheneq( jpi*jpj, zvlmsk(:,:) , 1, 1., ndextrd1, ndimtrd1 )    ! surface
213            ELSE
214               CALL wheneq( jpi*jpj, tmask(:,:,1), 1, 1., ndextrd1, ndimtrd1 )    ! surface
215            ENDIF                               
216
217            ! no more pass here
218            ionce = 0
219
220         ENDIF
221         
222         IF( idebug /= 0 ) THEN
223            ! CALL prihre (zvlmsk,jpi,jpj,1,jpi,2,1,jpj,2,3,numout)
224            WRITE(numout,*) ' debuging trd_mld_zint: I.1 done ' 
225            CALL FLUSH(numout)
226         ENDIF
227
228
229         ! I.2 probability density function of presence in mixed-layer
230         ! --------------------------------
231         ! (i.e. weight of each grid point in vertical integration : wkx(ji,jj,jk)
232
233
234         ! initialize wkx with vertical scale factor in mixed-layer
235
236         wkx(:,:,:) = 0.e0
237         DO jk = 1, jpktrd
238            DO jj = 1,jpj
239               DO ji = 1,jpi
240                  IF( jk - nmld(ji,jj) < 0. )   wkx(ji,jj,jk) = fse3t(ji,jj,jk) * tmask(ji,jj,jk)
241               END DO
242            END DO
243         END DO
244         
245         ! compute mixed-layer depth : rmld
246         
247         rmld(:,:) = 0.
248         DO jk = 1, jpktrd
249            rmld(:,:) = rmld(:,:) + wkx(:,:,jk)
250         END DO
251         
252         ! compute PDF
253
254         DO jk = 1, jpktrd
255            wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1., rmld(:,:) )
256         END DO
257
258         IF( idebug /= 0 ) THEN
259            WRITE(numout,*) ' debuging trd_mld_zint: I.2 done ' 
260            CALL FLUSH(numout)
261         ENDIF
262
263         ! Set counter icount to 0 to avoid this part at each time step
264         icount = 0
265
266      ENDIF
267
268
269      !  ====================================================
270      !   II. vertical integration of trends in mixed-layer
271      !  ====================================================
272
273      ! II.1 vertical integration of 3D and 2D trends
274      ! ---------------------------------------------
275
276      SELECT CASE (ctype)
277
278      CASE ('3D')       ! 3D treatment
279
280         ! trends terms in the mixed-layer
281         DO jk = 1, jpktrd
282            ! Temperature
283            tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,jk) * wkx(:,:,jk)   
284
285            ! Salinity
286            smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,jk) * wkx(:,:,jk)   
287         ENDDO
288
289      CASE ('2D')       ! 2D treatment
290
291         SELECT CASE (ktrd) 
292
293         CASE (jpmldldf)
294
295# if defined key_trabbl_dif
296               ! trends terms from Beckman over-flow parameterization
297               DO jj = 1,jpj
298                  DO ji = 1,jpi
299                     ikb = MAX( mbathy(ji,jj)-1, 1 )
300                     ! beckmann component -> horiz. part of lateral diffusion
301                     tmltrd(ji,jj,ktrd) = tmltrd(ji,jj,ktrd) + pttrdmld(ji,jj,1) * wkx(ji,jj,ikb)
302                     smltrd(ji,jj,ktrd) = smltrd(ji,jj,ktrd) + pstrdmld(ji,jj,1) * wkx(ji,jj,ikb)
303                  END DO
304               END DO
305# endif
306
307         CASE DEFAULT
308
309            ! trends terms at upper boundary of mixed-layer
310
311            ! forcing term (non penetrative)
312            ! Temperature
313            tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + pttrdmld(:,:,1) * wkx(:,:,1)   
314
315            ! forcing term
316            ! Salinity
317            smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + pstrdmld(:,:,1) * wkx(:,:,1)   
318
319         END SELECT
320
321      END SELECT
322
323      IF( idebug /= 0 ) THEN
324         IF(lwp) WRITE(numout,*) ' debuging trd_mld_zint: II.1 done' 
325         CALL FLUSH(numout)
326      ENDIF
327
328   END SUBROUTINE trd_mld_zint
329
330
331
332   SUBROUTINE trd_mld( kt )
333      !!----------------------------------------------------------------------
334      !!                  ***  ROUTINE trd_mld  ***
335      !!
336      !! ** Purpose :  computation of cumulated trends over analysis period
337      !!               and make outputs (NetCDF or DIMG format)
338      !!
339      !! ** Method/usage :
340      !!
341      !! History :
342      !!   9.0  !  04-08  (C. Talandier) New trends organization
343      !!----------------------------------------------------------------------
344      !! * Arguments
345      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
346
347      !! * Local declarations
348      INTEGER :: ji, jj, jk, jl, ik, it
349
350      REAL(wp) :: zmean, zavt
351
352      REAL(wp) ,DIMENSION(jpi,jpj) ::   &
353         ztmltot, ztmlres,              &
354         zsmltot, zsmlres,              & 
355         z2d
356
357#if defined key_dimgout
358      INTEGER ::  iyear,imon,iday
359      CHARACTER(LEN=80) :: cltext, clmode
360#endif
361      !!----------------------------------------------------------------------
362
363      ! I. trends terms at lower boundary of mixed-layer
364      ! ------------------------------------------------
365
366      DO jj = 1,jpj
367         DO ji = 1,jpi
368           
369            ik = nmld(ji,jj)
370           
371            ! Temperature
372            ! entrainment due to vertical diffusion
373            !       - due to vertical mixing scheme (TKE)
374            zavt = avt(ji,jj,ik)
375            tmltrd(ji,jj,jpmldevd) = - 1. * zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)   &
376               &                   * ( tn(ji,jj,ik-1) - tn(ji,jj,ik) )   &
377               &                   / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
378            ! Salinity
379            ! entrainment due to vertical diffusion
380            !       - due to vertical mixing scheme (TKE)
381            zavt = fsavs(ji,jj,ik)
382            smltrd(ji,jj,jpmldevd) = -1. * zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik)   &
383               &                  * ( sn(ji,jj,ik-1) - sn(ji,jj,ik) )   &
384               &                  / MAX( 1., rmld(ji,jj) ) * tmask(ji,jj,1)
385         END DO
386      END DO
387
388      IF( l_traldf_iso ) THEN
389         ! We substract to the TOTAL vertical diffusion tmltrd(:,:,jpmldzdf)
390         ! computed in subroutines trazdf_iso.F90 or trazdf_imp.F90
391         ! the vertical part du to the Kz in order to keep only the vertical
392         ! isopycnal diffusion (i.e the isopycnal diffusion componant on the vertical):
393         tmltrd(:,:,jpmldzdf) = tmltrd(:,:,jpmldzdf) - tmltrd(:,:,jpmldevd)   ! - due to isopycnal mixing scheme (implicit part)
394         smltrd(:,:,jpmldzdf) = smltrd(:,:,jpmldzdf) - smltrd(:,:,jpmldevd)   ! - due to isopycnal mixing scheme (implicit part)
395      ENDIF
396
397      ! Boundary conditions
398      DO jk = 1, jpltrd
399         CALL lbc_lnk( tmltrd, 'T', 1. )
400         CALL lbc_lnk( smltrd, 'T', 1. )
401      END DO
402
403      IF( idebug /= 0 ) THEN
404         WRITE(numout,*) ' debuging trd_mld: I. done' 
405         CALL FLUSH(numout)
406      ENDIF
407
408      !  =================================
409      !   II. Cumulated trends
410      !  =================================
411
412      ! II.1 set before values of vertically average T and S
413      ! ---------------------------------------------------
414
415      IF( kt > nit000 ) THEN
416         tmlb(:,:) = tml(:,:)
417         smlb(:,:) = sml(:,:)
418      ENDIF
419
420      ! II.2 vertically integrated T and S
421      ! ---------------------------------
422
423      tml(:,:) = 0.
424      sml(:,:) = 0.
425
426      DO jk = 1, jpktrd - 1
427         tml(:,:) = tml(:,:) + wkx(:,:,jk) * tn(:,:,jk)
428         sml(:,:) = sml(:,:) + wkx(:,:,jk) * sn(:,:,jk) 
429      END DO
430
431      IF(idebug /= 0) THEN
432         WRITE(numout,*) ' debuging trd_mld: II.2 done' 
433         CALL FLUSH(numout)
434      ENDIF
435
436      ! II.3 set `before' mixed layer values for kt = nit000+1
437      ! --------------------------------------------------------
438
439      IF( kt == nit000+1 ) THEN
440         tmlbb(:,:) = tmlb(:,:)
441         tmlbn(:,:) = tml (:,:)
442         smlbb(:,:) = smlb(:,:)
443         smlbn(:,:) = sml (:,:)
444      ENDIF
445
446      IF( idebug /= 0 ) THEN
447         WRITE(numout,*) ' debuging trd_mld: II.3 done' 
448         CALL FLUSH(numout)
449      ENDIF
450
451      ! II.4 cumulated trends over analysis period (kt=2 to nwrite)
452      ! -----------------------------------------------------------
453
454      ! trends cumulated over nwrite-2 time steps
455
456      IF( kt >= nit000+2 ) THEN
457         nmoymltrd = nmoymltrd + 1
458         DO jl = 1, jpltrd
459            tmltrdm(:,:) = tmltrdm(:,:) + tmltrd(:,:,jl)
460            smltrdm(:,:) = smltrdm(:,:) + smltrd(:,:,jl)
461         END DO
462      ENDIF
463
464      IF( idebug /= 0 ) THEN
465         WRITE(numout,*) ' debuging trd_mld: II.4 done' 
466         CALL FLUSH(numout)
467      ENDIF
468
469      !  =============================================
470      !   III. Output in netCDF + residual computation
471      !  =============================================
472
473      ztmltot(:,:) = 0.
474      zsmltot(:,:) = 0.
475      ztmlres(:,:) = 0.
476      zsmlres(:,:) = 0.
477
478      IF( MOD( kt - nit000+1, nwrite ) == 0 ) THEN
479
480         ! III.1 compute total trend
481         ! ------------------------
482
483         zmean = float(nmoymltrd)
484         
485         ztmltot(:,:) = ( tml(:,:) - tmlbn(:,:) + tmlb(:,:) - tmlbb(:,:) ) /  (zmean * 2. * rdt)
486         zsmltot(:,:) = ( sml(:,:) - smlbn(:,:) + smlb(:,:) - smlbb(:,:) ) /  (zmean * 2. * rdt)
487
488         IF(idebug /= 0) THEN
489            WRITE(numout,*) ' zmean = ',zmean 
490            WRITE(numout,*) ' debuging trd_mld: III.1 done' 
491            CALL FLUSH(numout)
492         ENDIF
493         
494
495         ! III.2 compute residual
496         ! ---------------------
497
498         ztmlres(:,:) = ztmltot(:,:) - tmltrdm(:,:) / zmean
499         zsmlres(:,:) = zsmltot(:,:) - smltrdm(:,:) / zmean
500
501
502         ! Boundary conditions
503
504         CALL lbc_lnk( ztmltot, 'T', 1. )
505         CALL lbc_lnk( ztmlres, 'T', 1. )
506         CALL lbc_lnk( zsmltot, 'T', 1. )
507         CALL lbc_lnk( zsmlres, 'T', 1. )
508
509         IF( idebug /= 0 ) THEN
510            WRITE(numout,*) ' debuging trd_mld: III.2 done' 
511            CALL FLUSH(numout)
512         ENDIF
513
514
515         ! III.3 time evolution array swap
516         ! ------------------------------
517
518         tmlbb(:,:) = tmlb(:,:)
519         tmlbn(:,:) = tml (:,:)
520         smlbb(:,:) = smlb(:,:)
521         smlbn(:,:) = sml (:,:)
522
523         IF( idebug /= 0 ) THEN
524            WRITE(numout,*) ' debuging trd_mld: III.3 done' 
525            CALL FLUSH(numout)
526         ENDIF
527
528
529         ! III.4 zero cumulative array
530         ! ---------------------------
531
532          nmoymltrd = 0
533
534          tmltrdm(:,:) = 0.
535          smltrdm(:,:) = 0.
536
537          IF(idebug /= 0) THEN
538              WRITE(numout,*) ' debuging trd_mld: III.4 done' 
539              CALL FLUSH(numout)
540          ENDIF
541         
542      ENDIF
543
544      ! III.5 write trends to output
545      ! ---------------------------
546
547#if defined key_dimgout
548! code for dimg mpp output
549      IF ( MOD(kt,nwrite) == 0 ) THEN
550         WRITE(clmode,'(f5.1,a)' ) nwrite*rdt/86400.,' days average'
551         iyear = ndastp/10000
552         imon = (ndastp-iyear*10000)/100
553         iday = ndastp - imon*100 - iyear*10000
554         WRITE(clname,9000) TRIM(cexper),'MLDiags',iyear,imon,iday
555         cltext=TRIM(cexper)//' mld diags'//TRIM(clmode)
556         CALL dia_wri_dimg (clname, cltext, smltrd, jpltrd, '2')
557  9000   FORMAT(a,"_",a,"_y",i4.4,"m",i2.2,"d",i2.2,".dimgproc")
558       END IF
559
560#else
561      IF( kt >=  nit000+1 ) THEN
562
563         ! define time axis
564         it= kt-nit000+1
565         IF( lwp .AND. MOD( kt, nwrite ) == 0 ) THEN
566            WRITE(numout,*) '     trd_mld : write NetCDF fields'
567         ENDIF
568         
569         CALL histwrite( nidtrd,"somlttml",it,rmld          ,ndimtrd1,ndextrd1) ! Mixed-layer depth
570         
571         ! Temperature trends
572         ! ------------------
573         CALL histwrite( nidtrd,"somltemp",it,tml           ,ndimtrd1,ndextrd1) ! Mixed-layer temperature
574         CALL histwrite( nidtrd,"somlttto",it,ztmltot       ,ndimtrd1,ndextrd1) ! total
575         CALL histwrite( nidtrd,"somlttax",it,tmltrd(:,:, 1),ndimtrd1,ndextrd1) ! i- adv.
576         CALL histwrite( nidtrd,"somlttay",it,tmltrd(:,:, 2),ndimtrd1,ndextrd1) ! j- adv.
577         CALL histwrite( nidtrd,"somlttaz",it,tmltrd(:,:, 3),ndimtrd1,ndextrd1) ! vertical adv.
578         CALL histwrite( nidtrd,"somlttdh",it,tmltrd(:,:, 4),ndimtrd1,ndextrd1) ! hor. lateral diff.
579         CALL histwrite( nidtrd,"somlttfo",it,tmltrd(:,:, 5),ndimtrd1,ndextrd1) ! forcing
580
581         CALL histwrite( nidtrd,"somlbtdz",it,tmltrd(:,:, 6),ndimtrd1,ndextrd1) ! vert. diffusion
582         CALL histwrite( nidtrd,"somlbtdt",it,ztmlres       ,ndimtrd1,ndextrd1) ! dh/dt entrainment (residual)
583         IF( l_traldf_iso ) THEN
584            CALL histwrite( nidtrd,"somlbtdv",it,tmltrd(:,:, 7),ndimtrd1,ndextrd1) ! vert. lateral diff.
585         ENDIF
586#if defined key_traldf_eiv
587         CALL histwrite( nidtrd,"somlgtax",it,tmltrd(:,:, 8),ndimtrd1,ndextrd1) ! i- adv. (eiv)
588         CALL histwrite( nidtrd,"somlgtay",it,tmltrd(:,:, 9),ndimtrd1,ndextrd1) ! j- adv. (eiv)
589         CALL histwrite( nidtrd,"somlgtaz",it,tmltrd(:,:,10),ndimtrd1,ndextrd1) ! vert. adv. (eiv)
590         z2d(:,:) = tmltrd(:,:,8) + tmltrd(:,:,9) + tmltrd(:,:,10)
591         CALL histwrite( nidtrd,"somlgtat",it,z2d           ,ndimtrd1,ndextrd1) ! total adv. (eiv)
592#endif   
593
594         ! Salinity trends
595         ! ---------------
596         CALL histwrite( nidtrd,"somlsalt",it,sml           ,ndimtrd1,ndextrd1) ! Mixed-layer salinity
597         CALL histwrite( nidtrd,"somltsto",it,zsmltot       ,ndimtrd1,ndextrd1) ! total
598         CALL histwrite( nidtrd,"somltsax",it,smltrd(:,:, 1),ndimtrd1,ndextrd1) ! i- adv.
599         CALL histwrite( nidtrd,"somltsay",it,smltrd(:,:, 2),ndimtrd1,ndextrd1) ! j- adv.
600         CALL histwrite( nidtrd,"somltsaz",it,smltrd(:,:, 3),ndimtrd1,ndextrd1) ! vert. adv.
601         CALL histwrite( nidtrd,"somltsdh",it,smltrd(:,:, 4),ndimtrd1,ndextrd1) ! hor. lateral diff.
602         CALL histwrite( nidtrd,"somltsfo",it,smltrd(:,:, 5),ndimtrd1,ndextrd1) ! forcing
603         CALL histwrite( nidtrd,"somlbsdz",it,smltrd(:,:, 6),ndimtrd1,ndextrd1) ! vert. diff.
604         CALL histwrite( nidtrd,"somlbsdt",it,zsmlres       ,ndimtrd1,ndextrd1) ! dh/dt entrainment (residual)
605         IF( l_traldf_iso ) THEN
606            CALL histwrite( nidtrd,"somlbsdv",it,smltrd(:,:, 7),ndimtrd1,ndextrd1) ! vert. lateral diff;
607         ENDIF
608#if defined key_traldf_eiv
609         CALL histwrite( nidtrd,"somlgsax",it,smltrd(:,:, 8),ndimtrd1,ndextrd1) ! i-adv. (eiv)
610         CALL histwrite( nidtrd,"somlgsay",it,smltrd(:,:, 9),ndimtrd1,ndextrd1) ! j-adv. (eiv)
611         CALL histwrite( nidtrd,"somlgsaz",it,smltrd(:,:,10),ndimtrd1,ndextrd1) ! vert. adv. (eiv)
612         z2d(:,:) = smltrd(:,:,8) + smltrd(:,:,9) + smltrd(:,:,10)
613         CALL histwrite( nidtrd,"somlgsat",it,z2d           ,ndimtrd1,ndextrd1) ! total adv. (eiv)
614#endif
615
616         IF( idebug /= 0 ) THEN
617            WRITE(numout,*) ' debuging trd_mld: III.5 done' 
618            CALL FLUSH(numout)
619         ENDIF
620
621         ! set counter icount to one to allow the calculation
622         ! of the surface control in the next time step in the trd_mld_zint subroutine
623         icount = 1
624
625      ENDIF
626
627      ! At the end of the 1st time step, set icount to 1 to be
628      ! able to compute the surface control at the beginning of
629      ! the second time step
630      IF( kt == nit000 )   icount = 1
631
632      IF( kt == nitend )   CALL histclo( nidtrd )
633#endif
634
635   END SUBROUTINE trd_mld
636
637
638
639   SUBROUTINE trd_mld_init
640      !!----------------------------------------------------------------------
641      !!                  ***  ROUTINE trd_mld_init  ***
642      !!
643      !! ** Purpose :   computation of vertically integrated T and S budgets
644      !!      from ocean surface down to control surface (NetCDF output)
645      !!
646      !! ** Method/usage :
647      !!
648      !! History :
649      !!        !  95-04  (J. Vialard)  Original code
650      !!        !  97-02  (E. Guilyardi)  Adaptation global + base cmo
651      !!        !  99-09  (E. Guilyardi)  Re-writing + netCDF output
652      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
653      !!   9.0  !  04-08  (C. Talandier) New trends organization
654      !!----------------------------------------------------------------------
655      !! * Local declarations
656      INTEGER :: ilseq
657
658      REAL(wp) ::   zjulian, zsto, zout
659
660      CHARACTER (LEN=21) ::   &
661         clold ='OLD'        , & ! open specifier (direct access files)
662         clunf ='UNFORMATTED', & ! open specifier (direct access files)
663         clseq ='SEQUENTIAL'     ! open specifier (direct access files)
664      CHARACTER (LEN=40) ::   clhstnam
665      CHARACTER (LEN=40) ::   clop
666      CHARACTER (LEN=12) ::   clmxl
667
668      NAMELIST/namtrd/ ntrd, nctls
669      !!----------------------------------------------------------------------
670
671      !  ===================
672      !   I. initialization
673      !  ===================
674
675      ! Open specifier
676      ilseq  = 1
677      idebug = 0      ! set it to 1 in case of problem to have more print
678      icount = 1     
679      ionce  = 1
680
681      ! namelist namtrd : trend diagnostic
682      REWIND( numnam )
683      READ  ( numnam, namtrd )
684
685      IF(lwp) THEN
686         WRITE(numout,*) ' '
687         WRITE(numout,*) 'trd_mld_init: mixed layer heat & freshwater budget trends'
688         WRITE(numout,*) '~~~~~~~~~~~~~'
689         WRITE(numout,*) ' '
690         WRITE(numout,*) '          Namelist namtrd : '
691         WRITE(numout,*) '             control surface for trends      nctls = ',nctls
692         WRITE(numout,*) ' '
693      ENDIF
694
695      ! cumulated trends array init
696      nmoymltrd = 0
697      tmltrdm(:,:) = 0.e0
698      smltrdm(:,:) = 0.e0
699
700      !  read control surface from file ctlsurf_idx
701
702      IF( nctls == 1 ) THEN
703         clname ='ctlsurf_idx'
704         CALL ctlopn(numbol,clname,clold,clunf,clseq,   &
705              ilseq,numout,lwp,1)
706         REWIND (numbol)
707         READ(numbol) nbol
708      ENDIF
709
710
711      IF( idebug /= 0 ) THEN
712         WRITE(numout,*) ' debuging trd_mld_init: 0. done ' 
713         CALL FLUSH(numout)
714      ENDIF
715
716      !  ===================================
717      !   II. netCDF output initialization
718      !  ===================================
719
720#if defined key_dimgout 
721
722#else
723      !     clmxl = legend root for netCDF output
724      IF( nctls == 0 ) THEN
725         ! control surface = mixed-layer with density criterion
726         ! (array nmln computed in zdfmxl.F90)
727         clmxl = 'Mixed Layer '
728      ELSE IF( nctls == 1 ) THEN
729         ! control surface = read index from file
730         clmxl = '      Bowl '
731      ELSE IF( nctls >= 2 ) THEN
732         ! control surface = model level
733         WRITE(clmxl,'(A9,I2,1X)') 'Levels 1-', nctls
734      ENDIF
735
736      !-----------------------------------------
737      ! II.1 Define frequency of output and means
738      ! -----------------------------------------
739
740#if defined key_diainstant
741      zsto = nwrite*rdt
742      clop ="inst(x)"
743#else
744      zsto = rdt
745      clop ="ave(x)"
746#endif
747      zout = nwrite*rdt
748
749      IF(lwp) WRITE (numout,*) ' trdmld_ncinit: netCDF initialization'
750
751      ! II.2 Compute julian date from starting date of the run
752      ! ------------------------
753
754      CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
755      IF (lwp) WRITE(numout,*)' ' 
756      IF (lwp) WRITE(numout,*)' Date 0 used :',nit000   &
757           ,' YEAR ', nyear,' MONTH ', nmonth,' DAY ', nday   &
758           ,'Julian day : ', zjulian
759
760
761      ! II.3 Define the T grid trend file (nidtrd)
762      ! ---------------------------------
763
764      CALL dia_nam( clhstnam, nwrite, 'trends' )                  ! filename
765      IF(lwp) WRITE(numout,*) ' Name of NETCDF file ', clhstnam
766      CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,1, jpi,   &  ! Horizontal grid : glamt and gphit
767           1, jpj, 0, zjulian, rdt, nh_t, nidtrd)
768
769      ! Declare output fields as netCDF variables
770
771      ! Mixed layer Depth
772      CALL histdef( nidtrd, "somlttml", clmxl//"Depth"              , "m"   ,   &  ! hmlp
773         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
774
775      ! Temperature
776      CALL histdef( nidtrd, "somltemp", clmxl//"Temperature"        , "C"   ,   &  ! ???
777         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
778      ! Temperature trends
779      CALL histdef( nidtrd, "somlttto", clmxl//"T Total"             , "C/s",   &  ! total
780         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )
781      CALL histdef( nidtrd, "somlttax", clmxl//"T Zonal Advection", "C/s",       & ! i-adv.
782         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
783      CALL histdef( nidtrd, "somlttay", clmxl//"T Meridional Advection", "C/s",   & ! j-adv.
784         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
785      CALL histdef( nidtrd, "somlttaz", clmxl//"T Vertical Advection", "C/s",   & ! vert. adv.
786         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
787      CALL histdef( nidtrd, "somlttdh", clmxl//"T Horizontal Diffusion ", "C/s",   & ! hor. lateral diff.
788         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
789      CALL histdef( nidtrd, "somlttfo", clmxl//"T Forcing", "C/s",   & ! forcing
790         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
791      CALL histdef( nidtrd, "somlbtdz", clmxl//"T Vertical Diffusion", "C/s",   & ! vert. diff.
792         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
793      CALL histdef( nidtrd, "somlbtdt", clmxl//"T dh/dt Entrainment (Residual)", "C/s",   & ! T * dh/dt
794         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zout, zout )
795      IF( l_traldf_iso ) THEN
796      CALL histdef( nidtrd, "somlbtdv", clmxl//"T Vert. lateral Diffusion","C/s",   & ! vertical diffusion entrainment (ISO)
797         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
798      ENDIF
799#if defined key_traldf_eiv
800      CALL histdef( nidtrd, "somlgtax", clmxl//"T Zonal EIV Advection", "C/s",   & ! i-adv. (eiv)
801         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
802      CALL histdef( nidtrd, "somlgtay", clmxl//"T Meridional EIV Advection", "C/s",   & ! j-adv. (eiv)
803         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
804      CALL histdef( nidtrd, "somlgtaz", clmxl//"T Vertical EIV Advection", "C/s",   & ! vert. adv. (eiv)
805         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
806      CALL histdef( nidtrd, "somlgtat", clmxl//"T Total EIV Advection", "C/s",   & ! total advection (eiv)
807         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
808#endif
809      ! Salinity
810      CALL histdef( nidtrd, "somlsalt", clmxl//"Salinity", "PSU",   & ! Mixed-layer salinity
811         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
812      ! Salinity trends
813      CALL histdef( nidtrd, "somltsto", clmxl//"S Total", "PSU/s",   & ! total
814         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
815      CALL histdef( nidtrd, "somltsax", clmxl//"S Zonal Advection", "PSU/s",   & ! i-advection
816         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
817      CALL histdef( nidtrd, "somltsay", clmxl//"S Meridional Advection", "PSU/s",   & ! j-advection
818         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
819      CALL histdef( nidtrd, "somltsaz", clmxl//"S Vertical Advection", "PSU/s",   & ! vertical advection
820         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
821      CALL histdef( nidtrd, "somltsdh", clmxl//"S Horizontal Diffusion ", "PSU/s",   & ! hor. lat. diff.
822         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
823      CALL histdef( nidtrd, "somltsfo", clmxl//"S Forcing", "PSU/s",   & ! forcing
824         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
825
826      CALL histdef( nidtrd, "somlbsdz", clmxl//"S Vertical Diffusion", "PSU/s",   & ! vert. diff.
827         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
828      CALL histdef( nidtrd, "somlbsdt", clmxl//"S dh/dt Entrainment (Residual)", "PSU/s",   & ! S * dh/dt
829         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
830      IF( l_traldf_iso ) THEN
831      ! vertical diffusion entrainment (ISO)
832      CALL histdef( nidtrd, "somlbsdv", clmxl//"S Vertical lateral Diffusion", "PSU/s",   & ! vert. lat. diff.
833         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
834      ENDIF
835#if defined key_traldf_eiv
836      CALL histdef( nidtrd, "somlgsax", clmxl//"S Zonal EIV Advection", "PSU/s",   & ! i-advection (eiv)
837         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
838      CALL histdef( nidtrd, "somlgsay", clmxl//"S Meridional EIV Advection", "PSU/s",   & ! j-advection (eiv)
839         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
840      CALL histdef( nidtrd, "somlgsaz", clmxl//"S Vertical EIV Advection", "PSU/s",   & ! vert. adv. (eiv)
841         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
842      CALL histdef( nidtrd, "somlgsat", clmxl//"S Total EIV Advection", "PSU/s",   & ! total adv. (eiv)
843         &          jpi, jpj, nh_t, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
844#endif
845      CALL histend( nidtrd )
846#endif
847
848      IF( idebug /= 0 ) THEN
849         WRITE(numout,*) ' debuging trd_mld_init: II. done' 
850         CALL FLUSH(numout)
851      ENDIF
852
853
854      END SUBROUTINE trd_mld_init
855
856#else
857   !!----------------------------------------------------------------------
858   !!   Default option :                                       Empty module
859   !!----------------------------------------------------------------------
860   LOGICAL, PUBLIC, PARAMETER ::   lk_trdmld = .FALSE.   !: momentum trend flag
861CONTAINS
862   SUBROUTINE trd_mld( kt )             ! Empty routine
863      INTEGER, INTENT( in) ::   kt
864      WRITE(*,*) 'trd_mld: You should not have seen this print! error?', kt
865   END SUBROUTINE trd_mld
866   SUBROUTINE trd_mld_zint( pttrdmld, pstrdmld, ktrd, ctype )
867      REAL, DIMENSION(:,:,:), INTENT( in ) ::   &
868         pttrdmld, pstrdmld                   ! Temperature and Salinity trends
869      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index
870      CHARACTER(len=2), INTENT( in ) ::   & 
871         ctype                                ! surface/bottom (2D arrays) or
872         !                                    ! interior (3D arrays) physics
873      WRITE(*,*) 'trd_mld_zint: You should not have seen this print! error?', pttrdmld(1,1,1)
874      WRITE(*,*) '  "      "  : You should not have seen this print! error?', pstrdmld(1,1,1)
875      WRITE(*,*) '  "      "  : You should not have seen this print! error?', ctype
876      WRITE(*,*) '  "      "  : You should not have seen this print! error?', ktrd
877   END SUBROUTINE trd_mld_zint
878   SUBROUTINE trd_mld_init              ! Empty routine
879      WRITE(*,*) 'trd_mld_init: You should not have seen this print! error?'
880   END SUBROUTINE trd_mld_init
881#endif
882
883   !!======================================================================
884END MODULE trdmld
Note: See TracBrowser for help on using the repository browser.