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
RevLine 
[2148]1MODULE diahsb
2   !!======================================================================
3   !!                       ***  MODULE  diahsb  ***
[2334]4   !! Ocean diagnostics: Heat, salt and volume budgets
[2148]5   !!======================================================================
[2334]6   !! History :  3.3  ! 2010-09  (M. Leclair)  Original code
[4161]7   !!                 ! 2012-10  (C. Rousset)  add iom_put
[2148]8   !!----------------------------------------------------------------------
[2334]9
10   !!----------------------------------------------------------------------
[2148]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
[2337]18   USE trabbc          ! bottom boundary condition
[2148]19   USE lib_mpp         ! distributed memory computing library
20   USE trabbc          ! bottom boundary condition
21   USE bdy_par         ! (for lk_bdy)
[3294]22   USE timing          ! preformance summary
[4161]23   USE iom             ! I/O manager
24   USE lib_fortran     ! glob_sum
25   USE restart         ! ocean restart
26   USE wrk_nemo         ! work arrays
[4166]27   USE sbcrnf         ! river runoffd
[2148]28
29   IMPLICIT NONE
30   PRIVATE
31
[2334]32   PUBLIC   dia_hsb        ! routine called by step.F90
[4161]33   PUBLIC   dia_hsb_init   ! routine called by nemogcm.F90
34   PUBLIC   dia_hsb_rst    ! routine called by step.F90
[2148]35
[4147]36   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
[2148]37
[4333]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     !
[2148]42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46
47   !!----------------------------------------------------------------------
[2287]48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[2281]49   !! $Id$
[2334]50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[2148]51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE dia_hsb( kt )
56      !!---------------------------------------------------------------------------
57      !!                  ***  ROUTINE dia_hsb  ***
58      !!     
[2334]59      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
[2148]60      !!
61      !! ** Method : - Compute the deviation of heat content, salt content and volume
[2334]62      !!             at the current time step from their values at nit000
63      !!             - Compute the contribution of forcing and remove it from these deviations
64      !!
[2148]65      !!---------------------------------------------------------------------------
66      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
67      !!
68      INTEGER    ::   jk                          ! dummy loop indice
[4333]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              !
[2148]77      !!---------------------------------------------------------------------------
[4161]78      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')     
[2148]79
[4161]80      CALL wrk_alloc( jpi, jpj, zsurf )
81 
82      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
83     
[2148]84      ! ------------------------- !
85      ! 1 - Trends due to forcing !
86      ! ------------------------- !
[4161]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
[4162]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
[2148]94      ! Add penetrative solar radiation
[4161]95      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * zsurf(:,:) )
[2148]96      ! Add geothermal heat flux
[4161]97      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * zsurf(:,:) )
98      !
[2148]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
[4161]103      ! ------------------------ !
104      ! 2a -  Content variations !
105      ! ------------------------ !
[4333]106      zdiff_v2 = 0._wp
107      zdiff_hc = 0._wp
108      zdiff_sc = 0._wp
[2148]109      ! volume variation (calculated with ssh)
[4161]110      zdiff_v1 = glob_sum( zsurf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
[2148]111      DO jk = 1, jpkm1
[4161]112         ! volume variation (calculated with scale factors)
[4333]113         zdiff_v2 = zdiff_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
[2148]114         ! heat content variation
[4333]115         zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
[2148]116            &                           - hc_loc_ini(:,:,jk) ) )
117         ! salt content variation
[4333]118         zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
[2148]119            &                           - sc_loc_ini(:,:,jk) ) )
120      ENDDO
121
122      ! Substract forcing from heat content, salt content and volume variations
[4161]123      !frc_v = zdiff_v2 - frc_v
124      !frc_t = zdiff_hc - frc_t
125      !frc_s = zdiff_sc - frc_s
[2148]126     
[4161]127      ! add ssh if not vvl
[4162]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
[4161]135      !
[2148]136      ! ----------------------- !
[4161]137      ! 2b -  Content           !
138      ! ----------------------- !
[4333]139      z_v2 = 0._wp
140      z_hc = 0._wp
141      z_sc = 0._wp
[4161]142      ! volume (calculated with ssh)
143      z_v1 = glob_sum( zsurf(:,:) * sshn(:,:) )
144      DO jk = 1, jpkm1
145         ! volume (calculated with scale factors)
[4333]146         z_v2 = z_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
[4161]147         ! heat content
[4333]148         z_hc = z_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) )
[4161]149         ! salt content
[4333]150         z_sc = z_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) )
[4161]151      ENDDO
152      ! add ssh if not vvl
[4162]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
[4161]158
159      ! ----------------------- !
[2148]160      ! 3 - Diagnostics writing !
161      ! ----------------------- !
162      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt )
[4161]163!
164      CALL iom_put( 'bgtemper' , z_hc / z_v2 )                      ! Temperature (C)
165      CALL iom_put( 'bgsaline' , z_sc / z_v2 )                      ! Salinity (psu)
[4333]166      CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J)
[4161]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)
[4333]172      CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_wp ) ! hc  - surface forcing (heat content)
[4161]173      CALL iom_put( 'bgfrcsal' , frc_s * 1.e-9 )                     ! sc  - surface forcing (salt content)
174      !
175      CALL wrk_dealloc( jpi, jpj, zsurf )
176      !
[3294]177      IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
[4161]178!
[2148]179   END SUBROUTINE dia_hsb
180
[2334]181
[2148]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
[2334]196      !!
[2148]197      NAMELIST/namhsb/ ln_diahsb
[4166]198      !
199      INTEGER  ::   ios
[2148]200      !!----------------------------------------------------------------------
[4166]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
[2334]217      !
[2148]218      IF(lwp) THEN                   ! Control print
219         WRITE(numout,*)
[4333]220         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
221         WRITE(numout,*) '~~~~~~~~~~~~'
[2148]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
[4161]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"
[4328]260         IF( lk_bdy ) THEN
[4161]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
[2148]267
[4161]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 )
[2148]290
[4161]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
[4333]306          frc_v = 0._wp                                           
307          frc_t = 0._wp                                           
308          frc_s = 0._wp                                                 
[4161]309       ENDIF
[2148]310
[4161]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
[2148]328
329   !!======================================================================
330END MODULE diahsb
Note: See TracBrowser for help on using the repository browser.