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.
diahsb.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago

remove remaining bugs in LIM3, so that it can run in both regional and global config
  • Property svn:keywords set to Id
File size: 15.4 KB
Line 
1MODULE diahsb
2   !!======================================================================
3   !!                       ***  MODULE  diahsb  ***
4   !! Ocean diagnostics: Heat, salt and volume budgets
5   !!======================================================================
6   !! History :  3.3  ! 2010-09  (M. Leclair)  Original code
7   !!                 ! 2012-10  (C. Rousset)  add iom_put
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   USE oce             ! ocean dynamics and tracers
12   USE dom_oce         ! ocean space and time domain
13   USE phycst          ! physical constants
14   USE sbc_oce         ! surface thermohaline fluxes
15   USE in_out_manager  ! I/O manager
16   USE domvvl          ! vertical scale factors
17   USE traqsr          ! penetrative solar radiation
18   USE trabbc          ! bottom boundary condition
19   USE lib_mpp         ! distributed memory computing library
20   USE trabbc          ! bottom boundary condition
21   USE bdy_par         ! (for lk_bdy)
22   USE timing          ! preformance summary
23   USE iom             ! I/O manager
24   USE lib_fortran     ! glob_sum
25   USE restart         ! ocean restart
26   USE wrk_nemo         ! work arrays
27   USE sbcrnf         ! river runoffd
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   dia_hsb        ! routine called by step.F90
33   PUBLIC   dia_hsb_init   ! routine called by nemogcm.F90
34   PUBLIC   dia_hsb_rst    ! routine called by step.F90
35
36   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
37
38   REAL(wp), SAVE                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ssh_ini              !
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcssh_loc_ini, scssh_loc_ini     !
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE dia_hsb( kt )
56      !!---------------------------------------------------------------------------
57      !!                  ***  ROUTINE dia_hsb  ***
58      !!     
59      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
60      !!
61      !! ** Method : - Compute the deviation of heat content, salt content and volume
62      !!             at the current time step from their values at nit000
63      !!             - Compute the contribution of forcing and remove it from these deviations
64      !!
65      !!---------------------------------------------------------------------------
66      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
67      !!
68      INTEGER    ::   jk                          ! dummy loop indice
69      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
70      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
71      REAL(wp)   ::   z_hc        , z_sc          ! heat and salt content
72      REAL(wp)   ::   z_v1        , z_v2          ! volume
73      REAL(wp)   ::   zdeltat                     !    -     -
74      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
75      REAL(wp)   ::   z_frc_trd_v                 !    -     -
76      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsurf              !
77      !!---------------------------------------------------------------------------
78      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')     
79
80      CALL wrk_alloc( jpi, jpj, zsurf )
81 
82      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
83     
84      ! ------------------------- !
85      ! 1 - Trends due to forcing !
86      ! ------------------------- !
87      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * zsurf(:,:) ) ! volume fluxes
88      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * zsurf(:,:) )       ! heat fluxes
89      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * zsurf(:,:) )       ! salt fluxes
90      !
91      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * zsurf(:,:) )
92      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * zsurf(:,:) )
93
94      ! Add penetrative solar radiation
95      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * zsurf(:,:) )
96      ! Add geothermal heat flux
97      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * zsurf(:,:) )
98      !
99      frc_v = frc_v + z_frc_trd_v * rdt
100      frc_t = frc_t + z_frc_trd_t * rdt
101      frc_s = frc_s + z_frc_trd_s * rdt
102
103      ! ------------------------ !
104      ! 2a -  Content variations !
105      ! ------------------------ !
106      zdiff_v2 = 0._wp
107      zdiff_hc = 0._wp
108      zdiff_sc = 0._wp
109      ! volume variation (calculated with ssh)
110      zdiff_v1 = glob_sum( zsurf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
111      DO jk = 1, jpkm1
112         ! volume variation (calculated with scale factors)
113         zdiff_v2 = zdiff_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
114         ! heat content variation
115         zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
116            &                           - hc_loc_ini(:,:,jk) ) )
117         ! salt content variation
118         zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
119            &                           - sc_loc_ini(:,:,jk) ) )
120      ENDDO
121
122      ! Substract forcing from heat content, salt content and volume variations
123      !frc_v = zdiff_v2 - frc_v
124      !frc_t = zdiff_hc - frc_t
125      !frc_s = zdiff_sc - frc_s
126     
127      ! add ssh if not vvl
128      IF( .NOT. lk_vvl ) THEN
129        zdiff_v2 = zdiff_v2 + zdiff_v1
130        zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_tem)   &
131               &                           - hcssh_loc_ini(:,:) ) )
132        zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_sal)   &
133               &                           - scssh_loc_ini(:,:) ) )
134      ENDIF
135      !
136      ! ----------------------- !
137      ! 2b -  Content           !
138      ! ----------------------- !
139      z_v2 = 0._wp
140      z_hc = 0._wp
141      z_sc = 0._wp
142      ! volume (calculated with ssh)
143      z_v1 = glob_sum( zsurf(:,:) * sshn(:,:) )
144      DO jk = 1, jpkm1
145         ! volume (calculated with scale factors)
146         z_v2 = z_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
147         ! heat content
148         z_hc = z_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) )
149         ! salt content
150         z_sc = z_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) )
151      ENDDO
152      ! add ssh if not vvl
153      IF( .NOT. lk_vvl ) THEN
154        z_v2 = z_v2 + z_v1
155        z_hc = z_hc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_tem) )
156        z_sc = z_sc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_sal) )
157      ENDIF
158
159      ! ----------------------- !
160      ! 3 - Diagnostics writing !
161      ! ----------------------- !
162      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt )
163!
164      CALL iom_put( 'bgtemper' , z_hc / z_v2 )                      ! Temperature (C)
165      CALL iom_put( 'bgsaline' , z_sc / z_v2 )                      ! Salinity (psu)
166      CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J)
167      CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9 )                 ! Salt content variation (psu*km3)
168      CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 )                    ! volume ssh (km3) 
169      CALL iom_put( 'bgsshtot' , zdiff_v1 / glob_sum(zsurf) )          ! ssh (m) 
170      CALL iom_put( 'bgvoltot' , zdiff_v2 * 1.e-9 )                 ! volume total (km3)
171      CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 )                     ! vol - surface forcing (volume)
172      CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_wp ) ! hc  - surface forcing (heat content)
173      CALL iom_put( 'bgfrcsal' , frc_s * 1.e-9 )                     ! sc  - surface forcing (salt content)
174      !
175      CALL wrk_dealloc( jpi, jpj, zsurf )
176      !
177      IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
178!
179   END SUBROUTINE dia_hsb
180
181
182   SUBROUTINE dia_hsb_init
183      !!---------------------------------------------------------------------------
184      !!                  ***  ROUTINE dia_hsb  ***
185      !!     
186      !! ** Purpose: Initialization for the heat salt volume budgets
187      !!
188      !! ** Method : Compute initial heat content, salt content and volume
189      !!
190      !! ** Action : - Compute initial heat content, salt content and volume
191      !!             - Initialize forcing trends
192      !!             - Compute coefficients for conversion
193      !!---------------------------------------------------------------------------
194      INTEGER            ::   jk       ! dummy loop indice
195      INTEGER            ::   ierror   ! local integer
196      !!
197      NAMELIST/namhsb/ ln_diahsb
198      !
199      INTEGER  ::   ios
200      !!----------------------------------------------------------------------
201
202      IF(lwp) THEN
203         WRITE(numout,*)
204         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
205         WRITE(numout,*) '~~~~~~~~ '
206      ENDIF
207
208      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
209      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
210901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )
211
212      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
213      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
214902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )
215      WRITE ( numond, namhsb )
216
217      !
218      IF(lwp) THEN                   ! Control print
219         WRITE(numout,*)
220         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
221         WRITE(numout,*) '~~~~~~~~~~~~'
222         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
223         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
224      ENDIF
225
226      IF( .NOT. ln_diahsb )   RETURN
227
228         ! ------------------- !
229         ! 1 - Allocate memory !
230         ! ------------------- !
231         ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror )
232         IF( ierror > 0 ) THEN
233            CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
234         ENDIF
235         ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror )
236         IF( ierror > 0 ) THEN
237            CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' )   ;   RETURN
238         ENDIF
239         ALLOCATE( hcssh_loc_ini(jpi,jpj), STAT=ierror )
240         IF( ierror > 0 ) THEN
241            CALL ctl_stop( 'dia_hsb: unable to allocate hcssh_loc_ini' )   ;   RETURN
242         ENDIF
243         ALLOCATE( scssh_loc_ini(jpi,jpj), STAT=ierror )
244         IF( ierror > 0 ) THEN
245            CALL ctl_stop( 'dia_hsb: unable to allocate scssh_loc_ini' )   ;   RETURN
246         ENDIF
247         ALLOCATE( e3t_ini(jpi,jpj,jpk)   , STAT=ierror )
248         IF( ierror > 0 ) THEN
249            CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' )      ;   RETURN
250         ENDIF
251         ALLOCATE( ssh_ini(jpi,jpj)       , STAT=ierror )
252         IF( ierror > 0 ) THEN
253            CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' )      ;   RETURN
254         ENDIF
255         
256         ! ----------------------------------------------- !
257         ! 2 - Time independant variables and file opening !
258         ! ----------------------------------------------- !
259         IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
260         IF( lk_bdy ) THEN
261            CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
262         ENDIF
263         !
264         CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
265      !
266   END SUBROUTINE dia_hsb_init
267
268   SUBROUTINE dia_hsb_rst( kt, cdrw )
269     !!---------------------------------------------------------------------
270     !!                   ***  ROUTINE limdia_rst  ***
271     !!                     
272     !! ** Purpose :   Read or write DIA file in restart file
273     !!
274     !! ** Method  :   use of IOM library
275     !!----------------------------------------------------------------------
276     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
277     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
278     !
279     INTEGER ::   jk   !
280     INTEGER ::   id1   ! local integers
281     !!----------------------------------------------------------------------
282     !
283     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
284        IF( ln_rstart ) THEN                   !* Read the restart file
285           !id1 = iom_varid( numror, 'frc_vol'  , ldstop = .FALSE. )
286           !
287           CALL iom_get( numror, 'frc_v', frc_v )
288           CALL iom_get( numror, 'frc_t', frc_t )
289           CALL iom_get( numror, 'frc_s', frc_s )
290
291           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
292           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
293           CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini )
294           CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini )
295           CALL iom_get( numror, jpdom_autoglo, 'hcssh_loc_ini', hcssh_loc_ini )
296           CALL iom_get( numror, jpdom_autoglo, 'scssh_loc_ini', scssh_loc_ini )
297       ELSE
298          ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
299          DO jk = 1, jpk
300             e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors
301             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content
302             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
303          END DO
304          hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
305          scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
306          frc_v = 0._wp                                           
307          frc_t = 0._wp                                           
308          frc_s = 0._wp                                                 
309       ENDIF
310
311     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
312        !                                   ! -------------------
313        IF(lwp) WRITE(numout,*) '---- dia-rst ----'
314        CALL iom_rstput( kt, nitrst, numrow, 'frc_v'   , frc_v     )
315        CALL iom_rstput( kt, nitrst, numrow, 'frc_t'   , frc_t     )
316        CALL iom_rstput( kt, nitrst, numrow, 'frc_s'   , frc_s     )
317       
318        CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
319        CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
320        CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
321        CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
322        CALL iom_rstput( kt, nitrst, numrow, 'hcssh_loc_ini', hcssh_loc_ini )
323        CALL iom_rstput( kt, nitrst, numrow, 'scssh_loc_ini', scssh_loc_ini )
324        !
325     ENDIF
326     !
327   END SUBROUTINE dia_hsb_rst
328
329   !!======================================================================
330END MODULE diahsb
Note: See TracBrowser for help on using the repository browser.