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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 @ 4731

Last change on this file since 4731 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
  • Property svn:keywords set to Id
File size: 17.6 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
[4558]38   REAL(dp)                                ::   surf_tot                !
39   REAL(dp)                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
40   REAL(dp)                                ::   frc_wn_t      , frc_wn_s ! global forcing trends
41   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf      , ssh_ini              !
42   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
43   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini
[2148]44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48
49   !!----------------------------------------------------------------------
[2287]50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[2281]51   !! $Id$
[2334]52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[2148]53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE dia_hsb( kt )
58      !!---------------------------------------------------------------------------
59      !!                  ***  ROUTINE dia_hsb  ***
60      !!     
[2334]61      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
[2148]62      !!
63      !! ** Method : - Compute the deviation of heat content, salt content and volume
[2334]64      !!             at the current time step from their values at nit000
65      !!             - Compute the contribution of forcing and remove it from these deviations
66      !!
[2148]67      !!---------------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
69      !!
70      INTEGER    ::   jk                          ! dummy loop indice
[4558]71      REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
72      REAL(dp)   ::   zdiff_hc1   , zdiff_sc1     ! -   -   -   -   -   -   -   -
73      REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
74      REAL(dp)   ::   zerr_hc1    , zerr_sc1       ! heat and salt content misfit
75      REAL(dp)   ::   zvol_tot                    ! volume
76      REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
77      REAL(dp)   ::   z_frc_trd_v                 !    -     -
78      REAL(dp)   ::   z_wn_trd_t , z_wn_trd_s   !    -     -
79      REAL(dp)   ::   z_ssh_hc , z_ssh_sc   !    -     -
[2148]80      !!---------------------------------------------------------------------------
[4161]81      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')     
[2148]82
83      ! ------------------------- !
84      ! 1 - Trends due to forcing !
85      ! ------------------------- !
[4558]86      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) ) ! volume fluxes
87      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )       ! heat fluxes
88      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )       ! salt fluxes
89      ! Add runoff heat & salt input
90      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) )
91      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) )
[4162]92
[2148]93      ! Add penetrative solar radiation
[4558]94      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) )
[2148]95      ! Add geothermal heat flux
[4558]96      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * surf(:,:) )
[4161]97      !
[4558]98      IF( .NOT. lk_vvl ) THEN
99         z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) )
100         z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) )
101      ENDIF
102
[2148]103      frc_v = frc_v + z_frc_trd_v * rdt
104      frc_t = frc_t + z_frc_trd_t * rdt
105      frc_s = frc_s + z_frc_trd_s * rdt
[4558]106      !                                          ! Advection flux through fixed surface (z=0)
107      IF( .NOT. lk_vvl ) THEN
108         frc_wn_t = frc_wn_t + z_wn_trd_t * rdt
109         frc_wn_s = frc_wn_s + z_wn_trd_s * rdt
110      ENDIF
[2148]111
[4161]112      ! ------------------------ !
[4558]113      ! 2 -  Content variations !
[4161]114      ! ------------------------ !
[4558]115      zdiff_v2 = 0.d0
116      zdiff_hc = 0.d0
117      zdiff_sc = 0.d0
118
[2148]119      ! volume variation (calculated with ssh)
[4558]120      zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
121
122      ! heat & salt content variation (associated with ssh)
123      IF( .NOT. lk_vvl ) THEN
124         z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) )
125         z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) )
126      ENDIF
127
[2148]128      DO jk = 1, jpkm1
[4161]129         ! volume variation (calculated with scale factors)
[4558]130         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) &
131            &                           * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
[2148]132         ! heat content variation
[4558]133         zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * tmask(:,:,jk) & 
134            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
[2148]135         ! salt content variation
[4558]136         zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * tmask(:,:,jk)   &
137            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
[2148]138      ENDDO
139
140      ! Substract forcing from heat content, salt content and volume variations
[4558]141      zdiff_v1 = zdiff_v1 - frc_v
142      IF( lk_vvl )   zdiff_v2 = zdiff_v2 - frc_v
143      zdiff_hc = zdiff_hc - frc_t
144      zdiff_sc = zdiff_sc - frc_s
[4162]145      IF( .NOT. lk_vvl ) THEN
[4558]146         zdiff_hc1 = zdiff_hc + z_ssh_hc 
147         zdiff_sc1 = zdiff_sc + z_ssh_sc
148         zerr_hc1  = z_ssh_hc - frc_wn_t
149         zerr_sc1  = z_ssh_sc - frc_wn_s
[4162]150      ENDIF
[4558]151
[2148]152      ! ----------------------- !
[4558]153      ! 3 - Diagnostics writing !
[4161]154      ! ----------------------- !
[4558]155      zvol_tot   = 0.d0                                                   ! total ocean volume
[4161]156      DO jk = 1, jpkm1
[4558]157         zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
158      END DO
159
160      IF( lk_vvl ) THEN
161        CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature variation (C)
162        CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    variation (psu)
163        CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content variation (1.e20 J)
164        CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content variation (psu*km3)
165        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3) 
166        CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t variation (km3) 
167        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)
168        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C)
169        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu)
170      ELSE
171        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content variation (C)
172        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu)
173        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J)
174        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3)
175        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3) 
176        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)
177        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C)
178        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu)
179        CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
180        CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
[4162]181      ENDIF
[4558]182      !
183      IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' )
[4161]184
[3294]185      IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
[4161]186!
[2148]187   END SUBROUTINE dia_hsb
188
[2334]189
[2148]190   SUBROUTINE dia_hsb_init
191      !!---------------------------------------------------------------------------
192      !!                  ***  ROUTINE dia_hsb  ***
193      !!     
194      !! ** Purpose: Initialization for the heat salt volume budgets
195      !!
196      !! ** Method : Compute initial heat content, salt content and volume
197      !!
198      !! ** Action : - Compute initial heat content, salt content and volume
199      !!             - Initialize forcing trends
200      !!             - Compute coefficients for conversion
201      !!---------------------------------------------------------------------------
202      INTEGER            ::   jk       ! dummy loop indice
203      INTEGER            ::   ierror   ! local integer
[2334]204      !!
[2148]205      NAMELIST/namhsb/ ln_diahsb
[4166]206      !
207      INTEGER  ::   ios
[2148]208      !!----------------------------------------------------------------------
[4166]209
210      IF(lwp) THEN
211         WRITE(numout,*)
212         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
213         WRITE(numout,*) '~~~~~~~~ '
214      ENDIF
215
216      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
217      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
218901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )
219
220      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
221      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
222902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )
[4624]223      IF(lwm) WRITE ( numond, namhsb )
[4166]224
[2334]225      !
[2148]226      IF(lwp) THEN                   ! Control print
227         WRITE(numout,*)
[4333]228         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
229         WRITE(numout,*) '~~~~~~~~~~~~'
[2148]230         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
231         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
[4558]232         WRITE(numout,*)
[2148]233      ENDIF
234
235      IF( .NOT. ln_diahsb )   RETURN
[4558]236!      IF( .NOT. lk_mpp_rep ) &
237!        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', &
238!             &         ' whereas the global sum to be precise must be done in double precision ',&
239!             &         ' please add key_mpp_rep')
[2148]240
[4558]241      ! ------------------- !
242      ! 1 - Allocate memory !
243      ! ------------------- !
244      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &
245         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
246      IF( ierror > 0 ) THEN
247         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
248      ENDIF
249
250      IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
251      IF( ierror > 0 ) THEN
252         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
253      ENDIF
254
255      ! ----------------------------------------------- !
256      ! 2 - Time independant variables and file opening !
257      ! ----------------------------------------------- !
258      IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
259      IF(lwp) WRITE(numout,*) '~~~~~~~'
260      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
261      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area
262
263      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
[4161]264      !
[4558]265      ! ---------------------------------- !
266      ! 4 - initial conservation variables !
267      ! ---------------------------------- !
268      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
269      !
[4161]270   END SUBROUTINE dia_hsb_init
[2148]271
[4161]272   SUBROUTINE dia_hsb_rst( kt, cdrw )
273     !!---------------------------------------------------------------------
274     !!                   ***  ROUTINE limdia_rst  ***
275     !!                     
276     !! ** Purpose :   Read or write DIA file in restart file
277     !!
278     !! ** Method  :   use of IOM library
279     !!----------------------------------------------------------------------
280     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
281     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
282     !
283     INTEGER ::   jk   !
284     INTEGER ::   id1   ! local integers
285     !!----------------------------------------------------------------------
286     !
287     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
288        IF( ln_rstart ) THEN                   !* Read the restart file
289           !id1 = iom_varid( numror, 'frc_vol'  , ldstop = .FALSE. )
290           !
[4558]291           IF(lwp) WRITE(numout,*) '~~~~~~~'
292           IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
293           IF(lwp) WRITE(numout,*) '~~~~~~~'
[4161]294           CALL iom_get( numror, 'frc_v', frc_v )
295           CALL iom_get( numror, 'frc_t', frc_t )
296           CALL iom_get( numror, 'frc_s', frc_s )
[4558]297           IF( .NOT. lk_vvl ) THEN
298              CALL iom_get( numror, 'frc_wn_t', frc_wn_t )
299              CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
300           ENDIF
[4161]301           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
302           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
303           CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini )
304           CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini )
[4558]305           IF( .NOT. lk_vvl ) THEN
306              CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
307              CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
308           ENDIF
[4161]309       ELSE
[4558]310          IF(lwp) WRITE(numout,*) '~~~~~~~'
311          IF(lwp) WRITE(numout,*) ' dia_hsb at initial state '
312          IF(lwp) WRITE(numout,*) '~~~~~~~'
[4161]313          ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
314          DO jk = 1, jpk
315             e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors
316             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content
317             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
318          END DO
[4558]319          frc_v = 0.d0                                           ! volume       trend due to forcing
320          frc_t = 0.d0                                           ! heat content   -    -   -    -   
321          frc_s = 0.d0                                           ! salt content   -    -   -    -       
322          IF( .NOT. lk_vvl ) THEN
323             ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
324             ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
325             frc_wn_t = 0.d0                                       ! initial heat content misfit due to free surface
326             frc_wn_s = 0.d0                                       ! initial salt content misfit due to free surface
327          ENDIF
[4161]328       ENDIF
[2148]329
[4161]330     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
331        !                                   ! -------------------
[4558]332        IF(lwp) WRITE(numout,*) '~~~~~~~'
333        IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
334        IF(lwp) WRITE(numout,*) '~~~~~~~'
335
[4161]336        CALL iom_rstput( kt, nitrst, numrow, 'frc_v'   , frc_v     )
337        CALL iom_rstput( kt, nitrst, numrow, 'frc_t'   , frc_t     )
338        CALL iom_rstput( kt, nitrst, numrow, 'frc_s'   , frc_s     )
[4558]339        IF( .NOT. lk_vvl ) THEN
340           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t )
341           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
342        ENDIF
[4161]343        CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
344        CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
345        CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
346        CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
[4558]347        IF( .NOT. lk_vvl ) THEN
348           CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
349           CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
350        ENDIF
[4161]351        !
352     ENDIF
353     !
354   END SUBROUTINE dia_hsb_rst
[2148]355
356   !!======================================================================
357END MODULE diahsb
Note: See TracBrowser for help on using the repository browser.