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

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

  1. Initial restructuring.
  2. Use fldread to read open boundary data.
  • Property svn:keywords set to Id
File size: 13.0 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   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   USE oce             ! ocean dynamics and tracers
11   USE dom_oce         ! ocean space and time domain
12   USE phycst          ! physical constants
13   USE sbc_oce         ! surface thermohaline fluxes
14   USE in_out_manager  ! I/O manager
15   USE domvvl          ! vertical scale factors
16   USE traqsr          ! penetrative solar radiation
17   USE trabbc          ! bottom boundary condition
18   USE lib_mpp         ! distributed memory computing library
19   USE trabbc          ! bottom boundary condition
20   USE obc_par         ! (for lk_obc)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   dia_hsb        ! routine called by step.F90
26   PUBLIC   dia_hsb_init   ! routine called by opa.F90
27
28   LOGICAL, PUBLIC ::   ln_diahsb  = .FALSE.   !: check the heat and salt budgets
29
30   INTEGER                                 ::   numhsb                           !
31   REAL(dp)                                ::   surf_tot   , vol_tot             !
32   REAL(dp)                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
33   REAL(dp)                                ::   fact1                            ! conversion factors
34   REAL(dp)                                ::   fact21    , fact22               !     -         -
35   REAL(dp)                                ::   fact31    , fact32               !     -         -
36   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf      , ssh_ini              !
37   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48
49CONTAINS
50
51   SUBROUTINE dia_hsb( kt )
52      !!---------------------------------------------------------------------------
53      !!                  ***  ROUTINE dia_hsb  ***
54      !!     
55      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
56      !!
57      !! ** Method : - Compute the deviation of heat content, salt content and volume
58      !!             at the current time step from their values at nit000
59      !!             - Compute the contribution of forcing and remove it from these deviations
60      !!
61      !! ** Action : Write the results in the 'heat_salt_volume_budgets.txt' ASCII file
62      !!---------------------------------------------------------------------------
63      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
64      !!
65      INTEGER    ::   jk                          ! dummy loop indice
66      REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
67      REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
68      REAL(dp)   ::   z1_rau0                     ! local scalars
69      REAL(dp)   ::   zdeltat                     !    -     -
70      REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
71      REAL(dp)   ::   z_frc_trd_v                 !    -     -
72      !!---------------------------------------------------------------------------
73
74      ! ------------------------- !
75      ! 1 - Trends due to forcing !
76      ! ------------------------- !
77      z1_rau0 = 1.e0 / rau0
78      z_frc_trd_v = z1_rau0 * SUM( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) )     ! volume fluxes
79      z_frc_trd_t =           SUM( sbc_tsc(:,:,jp_tem) * surf(:,:) )     ! heat fluxes
80      z_frc_trd_s =           SUM( sbc_tsc(:,:,jp_sal) * surf(:,:) )     ! salt fluxes
81      ! Add penetrative solar radiation
82      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + ro0cpr * SUM( qsr     (:,:) * surf(:,:) )
83      ! Add geothermal heat flux
84      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + ro0cpr * SUM( qgh_trd0(:,:) * surf(:,:) )
85      IF( lk_mpp ) THEN
86         CALL mpp_sum( z_frc_trd_v )
87         CALL mpp_sum( z_frc_trd_t )
88      ENDIF
89      frc_v = frc_v + z_frc_trd_v * rdt
90      frc_t = frc_t + z_frc_trd_t * rdt
91      frc_s = frc_s + z_frc_trd_s * rdt
92
93      ! ----------------------- !
94      ! 2 -  Content variations !
95      ! ----------------------- !
96      zdiff_v2 = 0.d0
97      zdiff_hc = 0.d0
98      zdiff_sc = 0.d0
99      ! volume variation (calculated with ssh)
100      zdiff_v1 = SUM( surf(:,:) * tmask(:,:,1) * ( sshn(:,:) - ssh_ini(:,:) ) )
101      DO jk = 1, jpkm1
102         ! volume variation (calculated with scale factors)
103         zdiff_v2 = zdiff_v2 + SUM( surf(:,:) * tmask(:,:,jk)   &
104            &                       * ( fse3t_n(:,:,jk)         &
105            &                           - e3t_ini(:,:,jk) ) )
106         ! heat content variation
107         zdiff_hc = zdiff_hc + SUM( surf(:,:) * tmask(:,:,jk)          &
108            &                       * ( fse3t_n(:,:,jk) * tn(:,:,jk)   &
109            &                           - hc_loc_ini(:,:,jk) ) )
110         ! salt content variation
111         zdiff_sc = zdiff_sc + SUM( surf(:,:) * tmask(:,:,jk)          &
112            &                       * ( fse3t_n(:,:,jk) * sn(:,:,jk)   &
113            &                           - sc_loc_ini(:,:,jk) ) )
114      ENDDO
115
116      IF( lk_mpp ) THEN
117         CALL mpp_sum( zdiff_hc )
118         CALL mpp_sum( zdiff_sc )
119         CALL mpp_sum( zdiff_v1 )
120         CALL mpp_sum( zdiff_v2 )
121      ENDIF
122
123      ! Substract forcing from heat content, salt content and volume variations
124      zdiff_v1 = zdiff_v1 - frc_v
125      zdiff_v2 = zdiff_v2 - frc_v
126      zdiff_hc = zdiff_hc - frc_t
127      zdiff_sc = zdiff_sc - frc_s
128     
129      ! ----------------------- !
130      ! 3 - Diagnostics writing !
131      ! ----------------------- !
132      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt )
133      WRITE(numhsb , 9020) kt , zdiff_hc / vol_tot , zdiff_hc * fact1  * zdeltat,                                &
134         &                      zdiff_sc / vol_tot , zdiff_sc * fact21 * zdeltat, zdiff_sc * fact22 * zdeltat,   &
135         &                      zdiff_v1           , zdiff_v1 * fact31 * zdeltat, zdiff_v1 * fact32 * zdeltat,   &
136         &                      zdiff_v2           , zdiff_v2 * fact31 * zdeltat, zdiff_v2 * fact32 * zdeltat
137
138      IF ( kt == nitend ) CLOSE( numhsb )
139
1409020  FORMAT(I5,11D15.7)
141      !
142   END SUBROUTINE dia_hsb
143
144
145   SUBROUTINE dia_hsb_init
146      !!---------------------------------------------------------------------------
147      !!                  ***  ROUTINE dia_hsb  ***
148      !!     
149      !! ** Purpose: Initialization for the heat salt volume budgets
150      !!
151      !! ** Method : Compute initial heat content, salt content and volume
152      !!
153      !! ** Action : - Compute initial heat content, salt content and volume
154      !!             - Initialize forcing trends
155      !!             - Compute coefficients for conversion
156      !!---------------------------------------------------------------------------
157      CHARACTER (len=32) ::   cl_name  ! output file name
158      INTEGER            ::   jk       ! dummy loop indice
159      INTEGER            ::   ierror   ! local integer
160      !!
161      NAMELIST/namhsb/ ln_diahsb
162      !!----------------------------------------------------------------------
163      !
164      REWIND ( numnam )              ! Read Namelist namhsb
165      READ   ( numnam, namhsb )
166      !
167      IF(lwp) THEN                   ! Control print
168         WRITE(numout,*)
169         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
170         WRITE(numout,*) '~~~~~~~~~~~~'
171         WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
172         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
173      ENDIF
174
175      IF( .NOT. ln_diahsb )   RETURN
176
177      ! ------------------- !
178      ! 1 - Allocate memory !
179      ! ------------------- !
180      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror )
181      IF( ierror > 0 ) THEN
182         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
183      ENDIF
184      ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror )
185      IF( ierror > 0 ) THEN
186         CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' )   ;   RETURN
187      ENDIF
188      ALLOCATE( e3t_ini(jpi,jpj,jpk)   , STAT=ierror )
189      IF( ierror > 0 ) THEN
190         CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' )      ;   RETURN
191      ENDIF
192      ALLOCATE( surf(jpi,jpj)          , STAT=ierror )
193      IF( ierror > 0 ) THEN
194         CALL ctl_stop( 'dia_hsb: unable to allocate surf' )         ;   RETURN
195      ENDIF
196      ALLOCATE( ssh_ini(jpi,jpj)       , STAT=ierror )
197      IF( ierror > 0 ) THEN
198         CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' )      ;   RETURN
199      ENDIF
200
201      ! ----------------------------------------------- !
202      ! 2 - Time independant variables and file opening !
203      ! ----------------------------------------------- !
204      WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
205      WRITE(numout,*) "~~~~~~~  output written in the 'heat_salt_volume_budgets.txt' ASCII file"
206      IF( lk_obc ) THEN
207         CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
208      ENDIF
209      cl_name    = 'heat_salt_volume_budgets.txt'                         ! name of output file
210      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
211      surf_tot  = SUM( surf(:,:) )                                       ! total ocean surface area
212      vol_tot   = 0.d0                                                   ! total ocean volume
213      DO jk = 1, jpkm1
214         vol_tot  = vol_tot + SUM( surf(:,:) * tmask(:,:,jk)     &
215            &                      * fse3t_n(:,:,jk)         )
216      END DO
217      IF( lk_mpp ) THEN
218         CALL mpp_sum( vol_tot )
219         CALL mpp_sum( surf_tot )
220      ENDIF
221
222      CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 )
223      !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80
224      WRITE( numhsb, 9010 ) "kt   |     heat content budget     |            salt content budget             ",   &
225         !                                                   123456789012345678901234567890123456789012345 -> 45
226         &                                                  "|            volume budget (ssh)             ",   &
227         !                                                   678901234567890123456789012345678901234567890 -> 45
228         &                                                  "|            volume budget (e3t)             "
229      WRITE( numhsb, 9010 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   &
230         &                                                  "|     [m3]         [mmm/s]          [SV]     ",   &
231         &                                                  "|     [m3]         [mmm/s]          [SV]     "
232
233      ! --------------- !
234      ! 3 - Conversions ! (factors will be multiplied by duration afterwards)
235      ! --------------- !
236
237      ! heat content variation   =>   equivalent heat flux:
238      fact1  = rau0 * rcp / surf_tot                                         ! [C*m3]   ->  [W/m2]
239      ! salt content variation   =>   equivalent EMP and equivalent "flow":
240      fact21 = 1.e3  / ( soce * surf_tot )                                   ! [psu*m3] ->  [mm/s]
241      fact22 = 1.e-6 / soce                                                  ! [psu*m3] ->  [Sv]
242      ! volume variation         =>   equivalent EMP and equivalent "flow":
243      fact31 = 1.e3  / surf_tot                                              ! [m3]     ->  [mm/s]
244      fact32 = 1.e-6                                                         ! [m3]     ->  [SV]
245
246      ! ---------------------------------- !
247      ! 4 - initial conservation variables !
248      ! ---------------------------------- !
249      ssh_ini(:,:) = sshn(:,:)                               ! initial ssh
250      DO jk = 1, jpk
251         e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                ! initial vertical scale factors
252         hc_loc_ini(:,:,jk) = tn(:,:,jk) * fse3t_n(:,:,jk)   ! initial heat content
253         sc_loc_ini(:,:,jk) = sn(:,:,jk) * fse3t_n(:,:,jk)   ! initial salt content
254      END DO
255      frc_v = 0.d0                                           ! volume       trend due to forcing
256      frc_t = 0.d0                                           ! heat content   -    -   -    -   
257      frc_s = 0.d0                                           ! salt content   -    -   -    -         
258      !
2599010  FORMAT(A80,A45,A45)
260      !
261   END SUBROUTINE dia_hsb_init
262
263   !!======================================================================
264END MODULE diahsb
Note: See TracBrowser for help on using the repository browser.