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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90 @ 4460

Last change on this file since 4460 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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