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/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 @ 4810

Last change on this file since 4810 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
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(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
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE dia_hsb( kt )
58      !!---------------------------------------------------------------------------
59      !!                  ***  ROUTINE dia_hsb  ***
60      !!     
61      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
62      !!
63      !! ** Method : - Compute the deviation of heat content, salt content and volume
64      !!             at the current time step from their values at nit000
65      !!             - Compute the contribution of forcing and remove it from these deviations
66      !!
67      !!---------------------------------------------------------------------------
68      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
69      !!
70      INTEGER    ::   jk                          ! dummy loop indice
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   !    -     -
80      !!---------------------------------------------------------------------------
81      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')     
82
83      ! ------------------------- !
84      ! 1 - Trends due to forcing !
85      ! ------------------------- !
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(:,:) )
92
93      ! Add penetrative solar radiation
94      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) )
95      ! Add geothermal heat flux
96      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * surf(:,:) )
97      !
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
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
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
111
112      ! ------------------------ !
113      ! 2 -  Content variations !
114      ! ------------------------ !
115      zdiff_v2 = 0.d0
116      zdiff_hc = 0.d0
117      zdiff_sc = 0.d0
118
119      ! volume variation (calculated with ssh)
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
128      DO jk = 1, jpkm1
129         ! volume variation (calculated with scale factors)
130         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) &
131            &                           * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
132         ! heat content variation
133         zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * tmask(:,:,jk) & 
134            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
135         ! salt content variation
136         zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * tmask(:,:,jk)   &
137            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
138      ENDDO
139
140      ! Substract forcing from heat content, salt content and volume variations
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
145      IF( .NOT. lk_vvl ) THEN
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
150      ENDIF
151
152      ! ----------------------- !
153      ! 3 - Diagnostics writing !
154      ! ----------------------- !
155      zvol_tot   = 0.d0                                                   ! total ocean volume
156      DO jk = 1, jpkm1
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)
181      ENDIF
182      !
183      IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' )
184
185      IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
186!
187   END SUBROUTINE dia_hsb
188
189
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
204      !!
205      NAMELIST/namhsb/ ln_diahsb
206      !
207      INTEGER  ::   ios
208      !!----------------------------------------------------------------------
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 )
223      IF(lwm) WRITE ( numond, namhsb )
224
225      !
226      IF(lwp) THEN                   ! Control print
227         WRITE(numout,*)
228         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
229         WRITE(numout,*) '~~~~~~~~~~~~'
230         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
231         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
232         WRITE(numout,*)
233      ENDIF
234
235      IF( .NOT. ln_diahsb )   RETURN
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')
240
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' )         
264      !
265      ! ---------------------------------- !
266      ! 4 - initial conservation variables !
267      ! ---------------------------------- !
268      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
269      !
270   END SUBROUTINE dia_hsb_init
271
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           !
291           IF(lwp) WRITE(numout,*) '~~~~~~~'
292           IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
293           IF(lwp) WRITE(numout,*) '~~~~~~~'
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 )
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
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 )
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
309       ELSE
310          IF(lwp) WRITE(numout,*) '~~~~~~~'
311          IF(lwp) WRITE(numout,*) ' dia_hsb at initial state '
312          IF(lwp) WRITE(numout,*) '~~~~~~~'
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
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
328       ENDIF
329
330     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
331        !                                   ! -------------------
332        IF(lwp) WRITE(numout,*) '~~~~~~~'
333        IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
334        IF(lwp) WRITE(numout,*) '~~~~~~~'
335
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     )
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
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 )
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
351        !
352     ENDIF
353     !
354   END SUBROUTINE dia_hsb_rst
355
356   !!======================================================================
357END MODULE diahsb
Note: See TracBrowser for help on using the repository browser.