source: trunk/NEMO/OPA_SRC/DIA/diagap.F90 @ 1312

Last change on this file since 1312 was 1312, checked in by smasson, 12 years ago

add a namelist logical to mask land points in NetCDF outputs, see ticket:322

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 10.7 KB
Line 
1MODULE diagap
2   !!======================================================================
3   !!                       ***  MODULE  diagap  ***
4   !! Ocean diagnostics : computation of model-data tracer gap
5   !!======================================================================
6#if defined key_diagap
7   !!----------------------------------------------------------------------
8   !!   'key_diagap' :                           model-data gap diagnostics
9   !!----------------------------------------------------------------------
10   !!   dia_gap      : model and data level mean temperature and salinity
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and tracers
14   USE dom_oce         ! ocean space and time domain
15   USE in_out_manager  ! I/O manager
16   USE dtatem          ! ???
17   USE dtasal          ! ???
18   USE daymod          ! calendar
19   USE dianam          ! build name of file (routine)
20   USE lib_mpp         ! distributed memory computing library
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Routine accessibility
26   PUBLIC dia_gap     ! called in step.F90 module
27
28   !! * Shared module variables
29   LOGICAL, PUBLIC, PARAMETER ::   lk_diagap = .TRUE.   !: model-data diagnostics flag
30
31   !! * Module variables
32   INTEGER ::                 &
33      ngap  ,                 &  ! time step frequency
34      nprg                       ! switch for control print
35   ! netcdf files and index common
36   INTEGER ::   &
37      nhoridg, ndepidg,             &
38      ndex(1)
39   REAL(wp) ::     &
40      vol                        ! total ocean volume
41   REAL(wp), DIMENSION(jpk) ::   &
42      volk , volkr,           &  ! level ocean volume and its inverse
43      tdtag, sdtag,           &  ! level mean data temperature & salinity
44      tmodg, smodg               ! level mean model temperature & salinity
45
46   !! * Substitutions
47#  include "domzgr_substitute.h90"
48   !!----------------------------------------------------------------------
49   !!   OPA 9.0 , LOCEAN-IPSL (2005)
50   !! $Id$
51   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE dia_gap( kt )
57      !!----------------------------------------------------------------------
58      !!                ***  ROUTINE dia_gap  ***
59      !!
60      !! ** Purpose :   Compute model and data level mean T and S profiles
61      !!      and output it in numgap (NetCDF or direct access file)
62      !!
63      !! ** Method :
64      !!         tracers (model) : tn, sn
65      !!         tracers (data)  : t_dta, s_dta
66      !!         difference between model and data tracers
67      !!         variance between model and data tracers
68      !!
69      !! ** Action  :   output in file numgap
70      !!
71      !! History :
72      !!   6.0  !  91-10  (G. Madec)  Original code
73      !!   7.0  !  92-07  (M. Imbard)  Add variance and mpp staff
74      !!   8.5  !  02-07  (G. Madec)  Free form, F90
75      !!----------------------------------------------------------------------
76      !! * Modules used
77      USE ioipsl
78
79      !! * Arguments
80      INTEGER, INTENT(in) ::   kt           ! ocean time-step index
81
82      !! * local declarations
83      INTEGER ::   ji, jj, jk   ! dummy loop indices
84      INTEGER ::   it
85      CHARACTER (len=40) ::   clhstnam
86      CHARACTER (len=40) ::   clop
87      REAL(wp) ::   zbt, zdt    ! temporary scalar
88      REAL(wp) ::   zsto, zout, zmax, zjulian
89      REAL(wp), DIMENSION(jpi) ::   zfoo
90      REAL(wp), DIMENSION(jpj) ::   zloo
91
92      NAMELIST/namgap/ ngap, nprg
93      !!----------------------------------------------------------------------
94
95
96      ! 0. initialization
97      ! -----------------
98
99      zdt = rdt
100      IF( nacc == 1 )   zdt = rdtmin
101
102      IF( kt == nit000 ) THEN
103
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) 'dia_gap : level mean model-data gap'
106         IF(lwp) WRITE(numout,*) '~~~~~~~'
107
108         ! Read diagap parameters in namelist namgap
109         ngap = 15
110         nprg = 15
111
112         REWIND( numnam )
113         READ( numnam, namgap )
114
115         IF(lwp) WRITE(numout,*) '          time step frequency for gap     ngap  = ',ngap
116         IF(lwp) WRITE(numout,*) '          switch for control print gap    nprg  = ',nprg
117
118         ! horizontal slab volume (tmask_i to take into account only interior ocean domain)
119         volk(:) = 0.e0
120         DO jk = 1, jpkm1
121            DO jj = 1, jpj
122               DO ji = 1, jpi
123                  zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask_i(ji,jj)
124                  volk(jk) = volk(jk) + tmask(ji,jj,jk) * zbt
125               END DO
126            END DO
127         END DO
128         IF( lk_mpp )   CALL mpp_sum( volk, jpk )   ! sum over the global domain
129
130         volkr(:) = 0.e0
131         DO jk = 1, jpk
132            IF( volk(jk) /= 0.e0 )   volkr(jk) = 1.e0 / volk(jk)
133         END DO
134         IF(lwp) THEN
135            WRITE(numout,*)
136            WRITE(numout,*) '     volume of each horizontal slab (km3) :'
137            DO jk = 1, jpkm1
138               WRITE(numout,9010) jk, volk(jk) * 1.e-9
139            END DO
140         ENDIF
1419010     FORMAT( 9X, 'level ', i3, f15.3 )
142
143         ! basin volume
144         vol = 0.e0
145         DO jk = 1, jpkm1
146            vol = vol + volk(jk)
147         END DO
148         IF(lwp) WRITE(numout,*)
149         IF(lwp) WRITE(numout,*) '     basin volume : ',vol*1.e-9, '  km3'
150
151         ! OPEN netcdf file
152
153         ! Define frequency of output and means
154         zsto = ngap * zdt
155         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!)
156         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time)
157         ENDIF
158         zout = ngap * zdt
159         zmax = FLOAT( nitend - nit000 + 1 ) * zdt
160         zfoo(1:jpi) = 0.e0
161         zloo(1:jpj) = 0.e0
162
163         ! Compute julian date from starting date of the run
164
165         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
166         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
167
168         CALL dia_nam( clhstnam, ngap, 'diagap' )
169         IF(lwp) WRITE(numout,*) 'Name of diagap NETCDF file ', clhstnam
170         ! Horizontal grid : zphi()
171         CALL histbeg( clhstnam, 1, zfoo, 1, zloo,   &
172                       1, 1, 1, 1, 0, zjulian, zdt, nhoridg, numgap, domain_id=nidom )
173         ! Vertical grids : gdept, gdepw
174         CALL histvert( numgap, "deptht", "Vertical T levels",   &
175                         "m", jpk, gdept_0, ndepidg )
176
177         ! define fields to be stored
178         ! mo(temp/sali)(da/mo/va): mean ocean temperature/salinity data/model/differences
179
180          CALL histdef( numgap, "motempda", "Mean Data Temperature","C",   &
181                        1, 1, nhoridg, jpk, 1, jpk, ndepidg, 32, clop,   &
182                        zsto,zout )
183          CALL histdef(numgap, "mosalida", "Mean Data Salinity","PSU",   &
184                        1, 1, nhoridg, jpk, 1, jpk, ndepidg, 32, clop,   &
185                        zsto,zout )
186          CALL histdef(numgap, "motempmo", "Mean Model Temperature","C",   &
187                        1, 1, nhoridg, jpk, 1, jpk, ndepidg, 32, clop,   &
188                        zsto,zout )
189          CALL histdef(numgap, "mosalimo", "Mean Model Salinity","PSU",   &
190                        1, 1, nhoridg, jpk, 1, jpk, ndepidg, 32, clop,   &
191                        zsto,zout )
192          CALL histend( numgap )
193          ndex(1) = 0
194      ENDIF
195
196
197      ! 1. horizontal averaged
198      ! ----------------------
199
200      IF( MOD( kt, ngap ) == 0 ) THEN
201
202         ! initialization
203         tdtag(:) = 0.e0
204         sdtag(:) = 0.e0
205         tmodg(:) = 0.e0
206         smodg(:) = 0.e0
207         
208         !  c a u t i o n  : use of tmask_i : average over the interior domain only
209         
210         ! sum
211         DO jk = 1, jpkm1
212            DO jj = 1, jpj
213              DO ji = 1, jpi
214                zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) * tmask_i(ji,jj)
215                tdtag(jk) = tdtag(jk) + zbt * t_dta(ji,jj,jk) * volkr(jk)
216                sdtag(jk) = sdtag(jk) + zbt * s_dta(ji,jj,jk) * volkr(jk)
217                tmodg(jk) = tmodg(jk) + zbt * tn  (ji,jj,jk) * volkr(jk)
218                smodg(jk) = smodg(jk) + zbt * sn  (ji,jj,jk) * volkr(jk)
219              END DO 
220            END DO 
221         END DO
222
223
224         ! 2. Basin averaged
225         ! -----------------
226         
227         DO jk = 1, jpkm1
228            tdtag(jpk) = tdtag(jpk) + tdtag(jk) * volk(jk) / vol
229            sdtag(jpk) = sdtag(jpk) + sdtag(jk) * volk(jk) / vol
230            tmodg(jpk) = tmodg(jpk) + tmodg(jk) * volk(jk) / vol
231            smodg(jpk) = smodg(jpk) + smodg(jk) * volk(jk) / vol
232         END DO 
233          IF( lk_mpp)   CALL mpp_sum( tdtag, jpk )   ! sum over the global domain
234          IF( lk_mpp)   CALL mpp_sum( sdtag, jpk )
235          IF( lk_mpp)   CALL mpp_sum( tmodg, jpk )
236          IF( lk_mpp)   CALL mpp_sum( smodg, jpk )
237
238          ! 3.  Averaged output in file numgap
239          ! ----------------------------======
240
241          IF( MOD( kt, nprg ) == 0 ) THEN
242              IF(lwp) THEN
243                  WRITE(numout,*) 'dia_gap: time step = ', kt, 'model - data'
244                  WRITE(numout,9300)
245
246                  DO jk = 1, jpk
247                    WRITE(numout,9310) tdtag(jk), tmodg(jk), tdtag(jk) - tmodg(jk), jk, fsdept(1,1,jk),   &
248                                       sdtag(jk), smodg(jk), sdtag(jk) - smodg(jk)
249                  END DO 
250              ENDIF
251          ENDIF
252
253 9300     FORMAT('  T-data  T-model    diff.   level  depth    S-data  S-model    diff.' )
254 9310     FORMAT( 3(f8.3,1x),i5,f9.0,2x,3(f8.3,1x) )
255
256          ! Output in file numgap
257
258          ! Netcdf write
259
260          IF (lwp ) THEN   ! Ok even in mpp after the call to mpp_sum
261             it = kt - nit000 + 1       ! define time axis
262             CALL histwrite( numgap, "motempda", it, tdtag, jpk, ndex )
263             CALL histwrite( numgap, "mosalida", it, sdtag, jpk, ndex )
264             CALL histwrite( numgap, "motempmo", it, tmodg, jpk, ndex )
265             CALL histwrite( numgap, "mosalimo", it, smodg, jpk, ndex )
266          END IF
267      ENDIF
268
269      ! Closing numgap file
270
271      IF( kt == nitend ) THEN
272         CALL histclo( numgap )      !   Netcdf file
273      ENDIF
274
275   END SUBROUTINE dia_gap
276
277#else
278   !!----------------------------------------------------------------------
279   !!   Default option :                                       Dummy module
280   !!----------------------------------------------------------------------
281   LOGICAL, PUBLIC, PARAMETER ::   lk_diagap = .FALSE.    !: diagap flag
282CONTAINS
283   SUBROUTINE dia_gap( kt )           ! Dummy routine
284      WRITE(*,*) 'dia_gap: You should not have seen this print! error?', kt
285   END SUBROUTINE dia_gap
286#endif
287
288   !!======================================================================
289END MODULE diagap
Note: See TracBrowser for help on using the repository browser.