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.
limcons.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90 @ 5176

Last change on this file since 5176 was 5176, checked in by clem, 9 years ago

LIM3 online diagnostics (for debugging). Improvements of the online check of heat, salt and mass. The thresholds used to print errors on conservation (in ocean.output) now depend on the ice area (and therefore are independent on the simulation.
  • Property svn:keywords set to Id
File size: 16.2 KB
Line 
1MODULE limcons
2   !!======================================================================
3   !!                   ***  MODULE  limcons  ***
4   !! LIM-3 Sea Ice :   conservation check
5   !!======================================================================
6   !! History :   -   ! Original code from William H. Lipscomb, LANL
7   !!            3.0  ! 2004-06  (M. Vancoppenolle)   Energy Conservation
8   !!            3.5  ! 2011-02  (G. Madec)  add mpp considerations
9   !!             -   ! 2014-05  (C. Rousset) add lim_cons_hsm
10   !!             -   ! 2015-03  (C. Rousset) add lim_cons_final
11   !!----------------------------------------------------------------------
12#if defined key_lim3
13   !!----------------------------------------------------------------------
14   !!   'key_lim3'                                      LIM-3 sea-ice model
15   !!----------------------------------------------------------------------
16   !!    lim_cons     :   checks whether energy, mass and salt are conserved
17   !!----------------------------------------------------------------------
18   USE phycst         ! physical constants
19   USE ice            ! LIM-3 variables
20   USE dom_ice        ! LIM-3 domain
21   USE dom_oce        ! ocean domain
22   USE in_out_manager ! I/O manager
23   USE lib_mpp        ! MPP library
24   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
25   USE sbc_oce , ONLY : sfx  ! Surface boundary condition: ocean fields
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   lim_column_sum
31   PUBLIC   lim_column_sum_energy
32   PUBLIC   lim_cons_check
33   PUBLIC   lim_cons_hsm
34   PUBLIC   lim_cons_final
35
36   !!----------------------------------------------------------------------
37   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE lim_column_sum( ksum, pin, pout )
44      !!-------------------------------------------------------------------
45      !!               ***  ROUTINE lim_column_sum ***
46      !!
47      !! ** Purpose : Compute the sum of xin over nsum categories
48      !!
49      !! ** Method  : Arithmetics
50      !!
51      !! ** Action  : Gets xin(ji,jj,jl) and computes xout(ji,jj)
52      !!---------------------------------------------------------------------
53      INTEGER                   , INTENT(in   ) ::   ksum   ! number of categories/layers
54      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pin    ! input field
55      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   pout   ! output field
56      !
57      INTEGER ::   jl   ! dummy loop indices
58      !!---------------------------------------------------------------------
59      !
60      pout(:,:) = pin(:,:,1)
61      DO jl = 2, ksum
62         pout(:,:) = pout(:,:) + pin(:,:,jl)
63      END DO
64      !
65   END SUBROUTINE lim_column_sum
66
67
68   SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout)
69      !!-------------------------------------------------------------------
70      !!               ***  ROUTINE lim_column_sum_energy ***
71      !!
72      !! ** Purpose : Compute the sum of xin over nsum categories
73      !!              and nlay layers
74      !!
75      !! ** Method  : Arithmetics
76      !!---------------------------------------------------------------------
77      INTEGER                                  , INTENT(in   ) ::   ksum   !: number of categories
78      INTEGER                                  , INTENT(in   ) ::   klay   !: number of vertical layers
79      REAL(wp), DIMENSION(jpi,jpj,nlay_i+1,jpl), INTENT(in   ) ::   pin   !: input field
80      REAL(wp), DIMENSION(jpi,jpj)             , INTENT(  out) ::   pout   !: output field
81      !
82      INTEGER ::   jk, jl   ! dummy loop indices
83      !!---------------------------------------------------------------------
84      !
85      pout(:,:) = 0._wp
86      DO jl = 1, ksum
87         DO jk = 2, klay 
88            pout(:,:) = pout(:,:) + pin(:,:,jk,jl)
89         END DO
90      END DO
91      !
92   END SUBROUTINE lim_column_sum_energy
93
94
95   SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid )
96      !!-------------------------------------------------------------------
97      !!               ***  ROUTINE lim_cons_check ***
98      !!
99      !! ** Purpose : Test the conservation of a certain variable
100      !!              For each physical grid cell, check that initial
101      !!              and final values
102      !!              of a conserved field are equal to within a small value.
103      !!
104      !! ** Method  :
105      !!---------------------------------------------------------------------
106      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px1          !: initial field
107      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px2          !: final field
108      REAL(wp)                , INTENT(in   ) ::   pmax_err     !: max allowed error
109      CHARACTER(len=15)       , INTENT(in   ) ::   cd_fieldid   !: field identifyer
110      !
111      INTEGER  ::   ji, jj          ! dummy loop indices
112      INTEGER  ::   inb_error       ! number of g.c where there is a cons. error
113      LOGICAL  ::   llconserv_err   ! = .true. if conservation check failed
114      REAL(wp) ::   zmean_error     ! mean error on error points
115      !!---------------------------------------------------------------------
116      !
117      IF(lwp) WRITE(numout,*) ' lim_cons_check '
118      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
119
120      llconserv_err = .FALSE.
121      inb_error     = 0
122      zmean_error   = 0._wp
123      IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err )   llconserv_err = .TRUE.
124
125      IF( llconserv_err ) THEN
126         DO jj = 1, jpj 
127            DO ji = 1, jpi
128               IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN
129                  inb_error   = inb_error + 1
130                  zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) )
131                  !
132                  IF(lwp) THEN
133                     WRITE (numout,*) ' ALERTE 99 '
134                     WRITE (numout,*) ' Conservation error: ', cd_fieldid
135                     WRITE (numout,*) ' Point             : ', ji, jj 
136                     WRITE (numout,*) ' lat, lon          : ', gphit(ji,jj), glamt(ji,jj)
137                     WRITE (numout,*) ' Initial value     : ', px1(ji,jj)
138                     WRITE (numout,*) ' Final value       : ', px2(ji,jj)
139                     WRITE (numout,*) ' Difference        : ', px2(ji,jj) - px1(ji,jj)
140                  ENDIF
141               ENDIF
142            END DO
143         END DO
144         !
145      ENDIF
146      IF(lk_mpp)   CALL mpp_sum( inb_error   )
147      IF(lk_mpp)   CALL mpp_sum( zmean_error )
148      !
149      IF( inb_error > 0 .AND. lwp ) THEN
150         zmean_error = zmean_error / REAL( inb_error, wp )
151         WRITE(numout,*) ' Conservation check for : ', cd_fieldid
152         WRITE(numout,*) ' Number of error points : ', inb_error
153         WRITE(numout,*) ' Mean error on these pts: ', zmean_error
154      ENDIF
155      !
156   END SUBROUTINE lim_cons_check
157
158
159   SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
160      !!--------------------------------------------------------------------------------------------------------
161      !!                                        ***  ROUTINE lim_cons_hsm ***
162      !!
163      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
164      !!                     + test if ice concentration and volume are > 0
165      !!
166      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiahsb=true
167      !!              It prints in ocean.output if there is a violation of conservation at each time-step
168      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to
169      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
170      !!              For salt and heat thresholds, ice is considered to have a salinity of 10
171      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)
172      !!--------------------------------------------------------------------------------------------------------
173      INTEGER         , INTENT(in)    :: icount        ! determine wether this is the beggining of the routine (0) or the end (1)
174      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
175      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
176      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft
177      REAL(wp)                        :: zvmin, zamin, zamax 
178      REAL(wp)                        :: zvtrp, zetrp
179      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
180      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
181
182      IF( icount == 0 ) THEN
183
184         ! salt flux
185         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
186            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &
187            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv )
188
189         ! water flux
190         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
191            &                  wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
192            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv )
193
194         ! heat flux
195         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
196            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
197            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv )
198
199         zvi_b  = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e12t * tmask(:,:,1) * zconv )
200
201         zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e12t * tmask(:,:,1) * zconv )
202
203         zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
204            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
205                            ) * e12t * tmask(:,:,1) * zconv )
206
207      ELSEIF( icount == 1 ) THEN
208
209         ! salt flux
210         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
211            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
212            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b
213
214         ! water flux
215         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
216            &                wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
217            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfw_b
218
219         ! heat flux
220         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
221            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
222            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b
223 
224         ! outputs
225         zvi  = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 )  &
226            &                    * e12t * tmask(:,:,1) * zconv ) - zvi_b ) * r1_rdtice - zfw ) * rday
227
228         zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 )  &
229            &                    * e12t * tmask(:,:,1) * zconv ) - zsmv_b ) * r1_rdtice + zfs ) * rday
230
231         zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
232            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
233            &                ) * e12t * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft
234
235         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
236         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t * tmask(:,:,1) * zconv ) * rday 
237         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e12t * tmask(:,:,1) * zconv )
238
239         zvmin = glob_min( v_i )
240         zamax = glob_max( SUM( a_i, dim=3 ) )
241         zamin = glob_min( a_i )
242
243         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)
244         zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2
245         zv_sill = zarea * 2.5e-5
246         zs_sill = zarea * 25.e-5
247         zh_sill = zarea * 10.e-5
248
249         IF(lwp) THEN
250            IF ( ABS( zvi  ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zvi
251            IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv
252            IF ( ABS( zei  ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zei
253            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'limtrp' ) THEN
254                                         WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
255                                         WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
256            ENDIF
257            IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
258            IF (     zamax   > rn_amax+epsi10 .AND. cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN
259                                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
260            ENDIF
261            IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
262         ENDIF
263
264      ENDIF
265
266   END SUBROUTINE lim_cons_hsm
267
268   SUBROUTINE lim_cons_final( cd_routine )
269      !!---------------------------------------------------------------------------------------------------------
270      !!                                   ***  ROUTINE lim_cons_final ***
271      !!
272      !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step
273      !!
274      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiahsb=true
275      !!              It prints in ocean.output if there is a violation of conservation at each time-step
276      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to
277      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
278      !!              For salt and heat thresholds, ice is considered to have a salinity of 10
279      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)
280      !!--------------------------------------------------------------------------------------------------------
281      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
282      REAL(wp)                        :: zhfx, zsfx, zvfx
283      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
284      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
285
286      ! heat flux
287      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t * tmask(:,:,1) * zconv ) 
288      ! salt flux
289      zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday
290      ! water flux
291      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t * tmask(:,:,1) * zconv ) * rday
292
293      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)
294      zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e12t * zconv ) ! in 1.e9 m2
295      zv_sill = zarea * 2.5e-5
296      zs_sill = zarea * 25.e-5
297      zh_sill = zarea * 10.e-5
298
299      IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx)
300      IF( ABS( zsfx ) > zs_sill ) WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx)
301      IF( ABS( zhfx ) > zh_sill ) WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx)
302
303   END SUBROUTINE lim_cons_final
304
305#else
306   !!----------------------------------------------------------------------
307   !!   Default option         Empty module            NO LIM sea-ice model
308   !!----------------------------------------------------------------------
309#endif
310   !!======================================================================
311END MODULE limcons
Note: See TracBrowser for help on using the repository browser.