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 branches/CMIP5_IPSL/NEMO/OPA_SRC/DIA – NEMO

source: branches/CMIP5_IPSL/NEMO/OPA_SRC/DIA/diagap.F90 @ 1846

Last change on this file since 1846 was 1846, checked in by mafoipsl, 14 years ago

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