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_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 @ 3938

Last change on this file since 3938 was 3938, checked in by flavoni, 11 years ago

dev_r3406_CNRS_LIM3: update LIM3, see ticket #1116
  • Property svn:keywords set to Id
File size: 15.9 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 obc_par         ! (for lk_obc)
22   USE bdy_par         ! (for lk_bdy)
23   USE timing          ! preformance summary
24   USE iom             ! I/O manager
25   USE lib_fortran     ! glob_sum
26   USE restart         ! ocean restart
27   USE wrk_nemo         ! work arrays
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  = .FALSE.   !: 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)   ::   z1_rau0                     ! local scalars
74      REAL(wp)   ::   zdeltat                     !    -     -
75      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
76      REAL(wp)   ::   z_frc_trd_v                 !    -     -
77      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsurf              !
78      !!---------------------------------------------------------------------------
79      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')     
80
81      CALL wrk_alloc( jpi, jpj, zsurf )
82 
83      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
84     
85      ! ------------------------- !
86      ! 1 - Trends due to forcing !
87      ! ------------------------- !
88      z1_rau0 = 1.e0 / rau0
89      z_frc_trd_v = z1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * zsurf(:,:) ) ! volume fluxes
90      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * zsurf(:,:) )       ! heat fluxes
91      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * zsurf(:,:) )       ! salt fluxes
92      ! Add penetrative solar radiation
93      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + ro0cpr * glob_sum( qsr     (:,:) * zsurf(:,:) )
94      ! Add geothermal heat flux
95      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + ro0cpr * glob_sum( qgh_trd0(:,:) * zsurf(:,:) )
96      !
97      frc_v = frc_v + z_frc_trd_v * rdt
98      frc_t = frc_t + z_frc_trd_t * rdt
99      frc_s = frc_s + z_frc_trd_s * rdt
100
101      ! ------------------------ !
102      ! 2a -  Content variations !
103      ! ------------------------ !
104      zdiff_v2 = 0._wp
105      zdiff_hc = 0._wp
106      zdiff_sc = 0._wp
107      ! volume variation (calculated with ssh)
108      zdiff_v1 = glob_sum( zsurf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
109      DO jk = 1, jpkm1
110         ! volume variation (calculated with scale factors)
111         zdiff_v2 = zdiff_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
112         ! heat content variation
113         zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
114            &                           - hc_loc_ini(:,:,jk) ) )
115         ! salt content variation
116         zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
117            &                           - sc_loc_ini(:,:,jk) ) )
118      ENDDO
119
120      ! Substract forcing from heat content, salt content and volume variations
121      !frc_v = zdiff_v2 - frc_v
122      !frc_t = zdiff_hc - frc_t
123      !frc_s = zdiff_sc - frc_s
124     
125      ! add ssh if not vvl
126#ifndef key_vvl
127     zdiff_v2 = zdiff_v2 + zdiff_v1
128     zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_tem)   &
129            &                           - hcssh_loc_ini(:,:) ) )
130     zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_sal)   &
131            &                           - scssh_loc_ini(:,:) ) )
132#endif 
133      !
134      ! ----------------------- !
135      ! 2b -  Content           !
136      ! ----------------------- !
137      z_v2 = 0._wp
138      z_hc = 0._wp
139      z_sc = 0._wp
140      ! volume (calculated with ssh)
141      z_v1 = glob_sum( zsurf(:,:) * sshn(:,:) )
142      DO jk = 1, jpkm1
143         ! volume (calculated with scale factors)
144         z_v2 = z_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
145         ! heat content
146         z_hc = z_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) )
147         ! salt content
148         z_sc = z_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) )
149      ENDDO
150      ! add ssh if not vvl
151#ifndef key_vvl
152     z_v2 = z_v2 + z_v1
153     z_hc = z_hc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_tem) )
154     z_sc = z_sc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_sal) )
155#endif 
156
157      ! ----------------------- !
158      ! 3 - Diagnostics writing !
159      ! ----------------------- !
160      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt )
161!
162      CALL iom_put( 'bgtemper',z_hc / z_v2 )               ! Temperature (C)
163      CALL iom_put( 'bgsaline',z_sc / z_v2 )               ! Salinity (psu)
164      !CALL iom_put( 'bgheatco',zdiff_hc*fact1*zdeltat )      ! Equivalent heat flux (W/m2)
165      !CALL iom_put( 'bgsaltco',zdiff_sc*fact21*zdeltat )     ! Equivalent water flux (mm/s)
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      !
200      REWIND ( numnam )              ! Read Namelist namhsb
201      READ   ( numnam, namhsb )
202      !
203      IF(lwp) THEN                   ! Control print
204         WRITE(numout,*)
205         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
206         WRITE(numout,*) '~~~~~~~~~~~~'
207         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
208         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
209      ENDIF
210
211      IF( .NOT. ln_diahsb )   RETURN
212
213         ! ------------------- !
214         ! 1 - Allocate memory !
215         ! ------------------- !
216         ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror )
217         IF( ierror > 0 ) THEN
218            CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
219         ENDIF
220         ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror )
221         IF( ierror > 0 ) THEN
222            CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' )   ;   RETURN
223         ENDIF
224         ALLOCATE( hcssh_loc_ini(jpi,jpj), STAT=ierror )
225         IF( ierror > 0 ) THEN
226            CALL ctl_stop( 'dia_hsb: unable to allocate hcssh_loc_ini' )   ;   RETURN
227         ENDIF
228         ALLOCATE( scssh_loc_ini(jpi,jpj), STAT=ierror )
229         IF( ierror > 0 ) THEN
230            CALL ctl_stop( 'dia_hsb: unable to allocate scssh_loc_ini' )   ;   RETURN
231         ENDIF
232         ALLOCATE( e3t_ini(jpi,jpj,jpk)   , STAT=ierror )
233         IF( ierror > 0 ) THEN
234            CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' )      ;   RETURN
235         ENDIF
236         ALLOCATE( ssh_ini(jpi,jpj)       , STAT=ierror )
237         IF( ierror > 0 ) THEN
238            CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' )      ;   RETURN
239         ENDIF
240         
241         ! ----------------------------------------------- !
242         ! 2 - Time independant variables and file opening !
243         ! ----------------------------------------------- !
244         WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
245         IF( lk_obc .or. lk_bdy ) THEN
246            CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
247         ENDIF
248         
249         ! ---------------------------------- !
250         ! 4 - initial conservation variables !
251         ! ---------------------------------- !
252         !ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
253         !DO jk = 1, jpk
254         !   e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors
255         !   hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content
256         !   sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
257         !END DO
258         !hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
259         !scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
260         !frc_v = 0._wp                                           ! volume       trend due to forcing
261         !frc_t = 0._wp                                           ! heat content   -    -   -    -   
262         !frc_s = 0._wp                                           ! salt content   -    -   -    -         
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.