Ignore:
Timestamp:
2015-03-25T17:13:51+01:00 (6 years ago)
Author:
clem
Message:

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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r5167 r5176  
    88   !!            3.5  ! 2011-02  (G. Madec)  add mpp considerations 
    99   !!             -   ! 2014-05  (C. Rousset) add lim_cons_hsm 
     10   !!             -   ! 2015-03  (C. Rousset) add lim_cons_final 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_lim3 
     
    157158 
    158159   SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
    159       !!------------------------------------------------------------------- 
    160       !!               ***  ROUTINE lim_cons_hsm *** 
    161       !! 
    162       !! ** Purpose : Test the conservation of heat, salt and mass for each routine 
    163       !! 
    164       !! ** Method  : 
    165       !!--------------------------------------------------------------------- 
    166       INTEGER         , INTENT(in)    :: icount      ! determine wether this is the beggining of the routine (0) or the end (1) 
    167       CHARACTER(len=*), INTENT(in)    :: cd_routine  ! name of the routine 
     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 
    168175      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    169176      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft 
    170177      REAL(wp)                        :: zvmin, zamin, zamax  
    171178      REAL(wp)                        :: zvtrp, zetrp 
    172       REAL(wp), PARAMETER             :: zconv = 1.e-9 
     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 
    173181 
    174182      IF( icount == 0 ) THEN 
    175183 
     184         ! salt flux 
    176185         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    177186            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
    178             &                ) *  e12t(:,:) * tmask(:,:,1) ) 
    179  
     187            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
     188 
     189         ! water flux 
    180190         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    181191            &                  wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    182             &                ) *  e12t(:,:) * tmask(:,:,1) ) 
    183  
     192            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
     193 
     194         ! heat flux 
    184195         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    185196            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    186197            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv ) 
    187198 
    188          zvi_b  = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) ) 
    189  
    190          zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic ) 
     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 ) 
    191202 
    192203         zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
    193204            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    & 
    194                             ) * e12t(:,:) * tmask(:,:,1) * zconv ) 
     205                            ) * e12t * tmask(:,:,1) * zconv ) 
    195206 
    196207      ELSEIF( icount == 1 ) THEN 
    197208 
     209         ! salt flux 
    198210         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    199211            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &  
    200             &              ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b 
    201  
     212            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b 
     213 
     214         ! water flux 
    202215         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
    203216            &                wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    & 
    204             &              ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b 
    205  
     217            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfw_b 
     218 
     219         ! heat flux 
    206220         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    207221            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    208222            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b 
    209223  
    210          zvi  = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 )  & 
    211             &                    * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw  
    212  
    213          zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * rhoic ) - zsmv_b ) * r1_rdtice + zfs 
     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 
    214230 
    215231         zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  & 
    216232            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    & 
    217             &                ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft 
    218  
    219          zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e12t(:,:) * tmask(:,:,1) )  
    220          zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e12t(:,:) * tmask(:,:,1) * zconv )  
     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 
    221239         zvmin = glob_min( v_i ) 
    222240         zamax = glob_max( SUM( a_i, dim=3 ) ) 
    223241         zamin = glob_min( a_i ) 
    224242 
    225         
     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 
    226249         IF(lwp) THEN 
    227             IF ( ABS( zvi  * rday ) >  0.5 * 1.e9 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday) 
    228             IF ( ABS( zsmv * rday ) >  5.  * 1.e9 ) WRITE(numout,*) 'violation saline [psu*kg/day] (',cd_routine,') = ',(zsmv * rday) 
    229             IF ( ABS( zei         ) >  2.  * 1.e9 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei) 
    230             IF ( zvmin <  -epsi10          ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
    231             IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN 
    232                                              WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     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 
    233256            ENDIF 
    234             IF ( zamin <  -epsi10          ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    235             IF( cd_routine == 'limtrp' .AND. ABS( zvtrp * rday ) > 0.5*1.e9 ) THEN 
    236                                              WRITE(numout,*) 'violation vtrp [kg/day]        (',cd_routine,') = ',(zvtrp * rday) 
    237                                              WRITE(numout,*) 'violation etrp [GW]            (',cd_routine,') = ',(zetrp ) 
     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 
    238260            ENDIF 
     261            IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    239262         ENDIF 
    240263 
     
    244267 
    245268   SUBROUTINE lim_cons_final( cd_routine ) 
    246       CHARACTER(len=*), INTENT(in)    :: cd_routine  ! name of the 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 
    247282      REAL(wp)                        :: zhfx, zsfx, zvfx 
    248       REAL(wp), PARAMETER             :: zconv = 1.e-9 
    249  
    250       zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t(:,:) * tmask(:,:,1) * zconv )  
    251       zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t(:,:) * tmask(:,:,1) ) * rday 
    252       zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e12t(:,:) * tmask(:,:,1) ) * rday  
    253  
    254       ! if error > 1 mm / 100 years over the Arctic Basin 
    255       IF( ABS( zvfx ) > 0.5 * 1.e9    ) WRITE(numout,*) 'violation vfx [kg/day]       (',cd_routine,') = ',(zvfx) 
    256       ! if error > 1 mm / 100 years over the Arctic Basin (ice with latent heat = 3e6 J/kg)  
    257       IF( ABS( zhfx ) > 2.  * 1.e9   ) WRITE(numout,*) 'violation hfx [GW]           (',cd_routine,') = ',(zhfx) 
    258       ! if error > 1 mm / 100 years over the Arctic Basin (ice of salinity = 10 pss) 
    259       IF( ABS( zsfx ) > 5.  * 1.e9   ) WRITE(numout,*) 'violation sfx [psu*kg/day]   (',cd_routine,') = ',(zsfx) 
     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) 
    260302 
    261303   END SUBROUTINE lim_cons_final 
Note: See TracChangeset for help on using the changeset viewer.