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

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

CT : UPDATE151 : New trends organization

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