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.
diagap.F90 in trunk/NEMO/OPA_SRC/DIA – NEMO

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

Last change on this file since 135 was 32, checked in by opalod, 20 years ago

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