Changeset 11501


Ignore:
Timestamp:
2019-09-04T19:41:33+02:00 (13 months ago)
Author:
clem
Message:

introduce a point-to-point conservation check, stop the model if it fails and write the issue in a file

Location:
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ice_ref

    r11401 r11501  
    233233!------------------------------------------------------------------------------ 
    234234   ln_icediachk     = .false.         !  check online the heat, mass & salt budgets (T) or not (F) 
     235      rn_icechk     =  1.             !     1. <=> 1mm of ice spuriously gained/lost during 1   year  (at any gridcell) 
     236                                      !                           --                        100 years (over the entire ice cover) 
    235237   ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F) 
    236238   ln_icectl        = .false.         !  ice points output for debug (T or F) 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/ice.F90

    r11397 r11501  
    196196   !                                     !!** ice-diagnostics namelist (namdia) ** 
    197197   LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F) 
     198   REAL(wp), PUBLIC ::   rn_icechk        !: 1. <=> 1mm of ice spuriously gained/lost during 1   year  (at any gridcell) 
     199                                          !                           --                     100 years (globally) 
    198200   LOGICAL , PUBLIC ::   ln_icediahsb     !: flag for ice diag (T) or not (F) 
    199201   LOGICAL , PUBLIC ::   ln_icectl        !: flag for sea-ice points output (T) or not (F) 
     
    351353   !! * Ice diagnostics 
    352354   !!---------------------------------------------------------------------- 
    353    ! thd refers to changes induced by thermodynamics 
    354    ! trp   ''         ''     ''       advection (transport of ice) 
    355    ! 
    356355   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi   !: transport of ice volume 
    357356   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs   !: transport of snw volume 
     
    365364   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw     !: snw volume variation   [m/s]  
    366365 
     366   !!---------------------------------------------------------------------- 
     367   !! * Ice conservation 
     368   !!---------------------------------------------------------------------- 
     369   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_v        !: conservation of ice volume 
     370   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_s        !: conservation of ice salt 
     371   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_t        !: conservation of ice heat 
     372   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fv       !: conservation of ice volume 
     373   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fs       !: conservation of ice salt 
     374   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_ft       !: conservation of ice heat 
    367375   ! 
    368376   !!---------------------------------------------------------------------- 
     
    389397      INTEGER :: ice_alloc 
    390398      ! 
    391       INTEGER :: ierr(15), ii 
     399      INTEGER :: ierr(16), ii 
    392400      !!----------------------------------------------------------------- 
    393401      ierr(:) = 0 
     
    461469         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) ) 
    462470 
     471      ! * Ice conservation 
     472      ii = ii + 1 
     473      ALLOCATE( diag_v (jpi,jpj) , diag_s (jpi,jpj) , diag_t (jpi,jpj),   &  
     474         &      diag_fv(jpi,jpj) , diag_fs(jpi,jpj) , diag_ft(jpi,jpj), STAT=ierr(ii) ) 
     475       
    463476      ! * SIMIP diagnostics 
    464477      ii = ii + 1 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icecor.F90

    r11397 r11501  
    6060      IF( ln_timing    )   CALL timing_start('icecor')                                                             ! timing 
    6161      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     62      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icecor',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    6263      ! 
    6364      IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN 
     
    164165      ! 
    165166      ! controls 
    166       IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    167       IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints 
     167      IF( ln_ctl       )   CALL ice_prt3D   ('icecor')                                                             ! prints 
    168168      IF( ln_icectl .AND. kn == 2 ) & 
    169          &                   CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints 
    170       IF( ln_timing      )   CALL timing_stop ('icecor')                                                             ! timing 
     169         &                 CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints 
     170      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     171      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icecor',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
     172      IF( ln_timing    )   CALL timing_stop ('icecor')                                                             ! timing 
    171173      ! 
    172174   END SUBROUTINE ice_cor 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icectl.F90

    r11382 r11501  
    1212   !!   'key_si3'                                       SI3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!    ice_cons_hsm     : conservation tests on heat, salt and mass 
    15    !!    ice_cons_final   : conservation tests on heat, salt and mass at end of time step 
     14   !!    ice_cons_hsm     : conservation tests on heat, salt and mass during a  time step (global)  
     15   !!    ice_cons_final   : conservation tests on heat, salt and mass at end of time step (global) 
     16   !!    ice_cons2D       : conservation tests on heat, salt and mass at each gridcell 
    1617   !!    ice_ctl          : control prints in case of crash 
    1718   !!    ice_prt          : control prints at a given grid point 
     
    2728   ! 
    2829   USE in_out_manager ! I/O manager 
     30   USE iom            ! I/O manager library 
    2931   USE lib_mpp        ! MPP library 
    3032   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
     
    3739   PUBLIC   ice_cons_hsm 
    3840   PUBLIC   ice_cons_final 
     41   PUBLIC   ice_cons2D 
    3942   PUBLIC   ice_ctl 
    4043   PUBLIC   ice_prt 
    4144   PUBLIC   ice_prt3D 
    4245 
     46   ! thresold values for conservation 
     47   !    these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk 
     48   REAL(wp), PARAMETER ::   zchk_m   = 1.e-5   ! kg/m2/s <=> 1mm of ice per 1   year  spuriously gained/lost 
     49   REAL(wp), PARAMETER ::   zchk2D_m = 1.e-7   !                   --       100 years          -- 
     50   REAL(wp), PARAMETER ::   zchk_s   = 1.e-4   ! g/m2/s  <=> 1mm of ice per 1   year  spuriously gained/lost (considering s=10g/kg) 
     51   REAL(wp), PARAMETER ::   zchk2D_s = 1.e-6   !                   --       100 years          -- 
     52   REAL(wp), PARAMETER ::   zchk_t   = 3.      ! W/m2    <=> 1mm of ice per 1   year  spuriously gained/lost (considering Lf=3e5J/kg) 
     53   REAL(wp), PARAMETER ::   zchk2D_t = 0.03    !                   --       100 years          -- 
     54    
    4355   !! * Substitutions 
    4456#  include "vectopt_loop_substitute.h90" 
     
    5971      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    6072      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    61       !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to 
     73      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine violations are set to 
    6274      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    6375      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
     
    6880      REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
    6981      !! 
    70       REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
    71       REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin 
     82      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
     83         &          zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 
    7284      REAL(wp) ::   zvtrp, zetrp 
    73       REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill 
    74       REAL(wp), PARAMETER ::   zconv = 1.e-9 ! convert W to GW and kg to Mt 
     85      REAL(wp) ::   zarea 
    7586      !!------------------------------------------------------------------- 
    7687      ! 
    7788      IF( icount == 0 ) THEN 
    78          !                          ! water flux 
    79          pdiag_fv = glob_sum( 'icectl',                                                                       & 
    80             &                 -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
    81             &                    wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
    82             &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
    83             &                    wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    84             &                  ) * e1e2t(:,:) ) * zconv 
     89 
     90         pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) 
     91         pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) 
     92         pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 
     93 
     94         ! mass flux 
     95         pdiag_fv = glob_sum( 'icectl',  & 
     96            &                         ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
     97            &                           wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t ) 
     98         ! salt flux 
     99         pdiag_fs = glob_sum( 'icectl',  & 
     100            &                         ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + & 
     101            &                           sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) 
     102         ! heat flux 
     103         pdiag_ft = glob_sum( 'icectl',  & 
     104            &                         (   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  & 
     105            &                           - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) 
     106 
     107      ELSEIF( icount == 1 ) THEN 
     108 
     109         ! -- mass diag -- ! 
     110         zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_rdtice       & 
     111            &         + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn +       & 
     112            &                                 wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & 
     113            &                                 wfx_ice_sub + wfx_spr ) * e1e2t )                                           & 
     114            &         - pdiag_fv 
    85115         ! 
    86          !                          ! salt flux 
    87          pdiag_fs = glob_sum( 'icectl',                                                                     & 
    88             &                  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    89             &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
    90             &                  ) *  e1e2t(:,:) ) * zconv  
     116         ! -- salt diag -- ! 
     117         zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice  & 
     118            &         + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni +           & 
     119            &                                 sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) & 
     120            &         - pdiag_fs 
    91121         ! 
    92          !                          ! heat flux 
    93          pdiag_ft = glob_sum( 'icectl',                                                                    & 
    94             &                  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    95             &                  - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    96             &                  ) *  e1e2t(:,:) ) * zconv 
    97  
    98          pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv ) 
    99  
    100          pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi            , dim=3 ) * e1e2t * zconv ) 
    101  
    102          pdiag_t = glob_sum( 'icectl', (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
    103             &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 
    104  
    105       ELSEIF( icount == 1 ) THEN 
    106  
    107          ! water flux 
    108          zfv = glob_sum( 'icectl',                                                                        & 
    109             &             -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
    110             &                wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
    111             &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
    112             &                wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    113             &              ) * e1e2t(:,:) ) * zconv - pdiag_fv 
    114  
    115          ! salt flux 
    116          zfs = glob_sum( 'icectl',                                                                       & 
    117             &              ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    118             &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &  
    119             &              ) * e1e2t(:,:) ) * zconv - pdiag_fs 
    120  
    121          ! heat flux 
    122          zft = glob_sum( 'icectl',                                                                      & 
    123             &              ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
    124             &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   & 
    125             &              ) * e1e2t(:,:) ) * zconv - pdiag_ft 
    126   
    127          ! outputs 
    128          zv = ( ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv  & 
    129             &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
    130  
    131          zs = ( ( glob_sum( 'icectl', SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) * zconv  & 
    132             &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
    133  
    134          zt = ( glob_sum( 'icectl',                                                                & 
    135             &             (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )                       & 
    136             &              + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   & 
    137             &   - pdiag_t ) * r1_rdtice + zft 
    138  
    139          ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
    140          zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t  ) * zconv * rday  
    141          zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv 
    142  
    143          zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    144          zvmin  = glob_min( 'icectl', v_i ) 
    145          zamin  = glob_min( 'icectl', a_i ) 
    146          zsmin  = glob_min( 'icectl', sv_i ) 
    147          zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
    148          zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
    149  
    150          ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    151          zarea   = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 
    152          zv_sill = zarea * 2.5e-5 
    153          zs_sill = zarea * 25.e-5 
    154          zt_sill = zarea * 10.e-5 
    155  
    156          IF(lwp) THEN 
     122         ! -- heat diag -- ! 
     123         zdiag_heat = ( glob_sum( 'icectl', ( SUM(SUM(e_i, dim=4), dim=3) + SUM(SUM(e_s, dim=4), dim=3) ) * e1e2t ) - pdiag_t & 
     124            &         ) * r1_rdtice                                                                                           & 
     125            &         + glob_sum( 'icectl', (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                      & 
     126            &                                - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )                    & 
     127            &         - pdiag_ft 
     128 
     129         ! -- min/max diag -- ! 
     130         zdiag_amax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     131         zdiag_vmin  = glob_min( 'icectl', v_i ) 
     132         zdiag_amin  = glob_min( 'icectl', a_i ) 
     133         zdiag_smin  = glob_min( 'icectl', sv_i ) 
     134         zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
     135         zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
     136 
     137         ! -- advection scheme is conservative? -- ! 
     138         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 
     139         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0 
     140 
     141         ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     142         zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
     143 
     144         IF( lwp ) THEN 
    157145            ! check conservation issues 
    158             IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv 
    159             IF ( ABS( zs ) > zs_sill )   WRITE(numout,*) 'violation saline [Mkg/day]    (',cd_routine,') = ',zs 
    160             IF ( ABS( zt ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt 
     146            IF( ABS(zdiag_mass) > zchk_m * rn_icechk * zarea ) & 
     147               &                   WRITE(numout,*)   cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 
     148            IF( ABS(zdiag_salt) > zchk_s * rn_icechk * zarea ) & 
     149               &                   WRITE(numout,*)   cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice 
     150            IF( ABS(zdiag_heat) > zchk_t * rn_icechk * zarea ) & 
     151               &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice 
     152            ! check negative values 
     153            IF( zdiag_vmin  < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i < 0         = ',zdiag_vmin 
     154            IF( zdiag_amin  < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i < 0         = ',zdiag_amin 
     155            IF( zdiag_smin  < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i < 0         = ',zdiag_smin 
     156            IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i < 0         = ',zdiag_eimin 
     157            IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s < 0         = ',zdiag_esmin 
    161158            ! check maximum ice concentration 
    162             IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  & 
    163                &                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
    164             ! check negative values 
    165             IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin 
    166             IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
    167             IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin 
    168             IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin 
    169             IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin 
    170 !clem: the following check fails (I think...) 
    171 !            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN 
    172 !                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp 
    173 !                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp 
    174 !            ENDIF 
     159            IF( zdiag_amax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
     160               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
     161            ! check if advection scheme is conservative 
     162            IF( ABS(zvtrp) > zchk_m*rn_icechk*zarea .AND. cd_routine == 'icedyn_adv' ) & 
     163               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
    175164         ENDIF 
    176165         ! 
     
    179168   END SUBROUTINE ice_cons_hsm 
    180169 
    181  
    182170   SUBROUTINE ice_cons_final( cd_routine ) 
    183171      !!------------------------------------------------------------------- 
     
    188176      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    189177      !!              It prints in ocean.output if there is a violation of conservation at each time-step 
    190       !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to 
     178      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine the violation are set to 
    191179      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    192180      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
    193181      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)  
    194182      !!------------------------------------------------------------------- 
    195       CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine 
    196       REAL(wp)                        :: zqmass, zhfx, zsfx, zvfx 
    197       REAL(wp)                        :: zarea, zv_sill, zs_sill, zt_sill 
    198       REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
     183      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine 
     184      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat 
     185      REAL(wp) ::   zarea 
    199186      !!------------------------------------------------------------------- 
    200187 
    201188      ! water flux 
    202       zvfx  = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday 
    203  
    204       ! salt flux 
    205       zsfx  = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) * zconv * rday 
    206  
    207       ! heat flux 
     189      ! -- mass diag -- ! 
     190      zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 
     191 
     192      ! -- salt diag -- ! 
     193      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) 
     194 
     195      ! -- heat diag -- ! 
    208196      ! clem: not the good formulation 
    209       !!zhfx  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  & 
    210       !!   &                        ) * e1e2t ) * zconv 
    211  
    212       ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)  
    213       zarea   = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2 
    214       zv_sill = zarea * 2.5e-5 
    215       zs_sill = zarea * 25.e-5 
    216       zt_sill = zarea * 10.e-5 
    217  
    218       IF(lwp) THEN 
    219          IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx  [Mt/day]       (',cd_routine,') = ',zvfx 
    220          IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx  [Mkg/day]      (',cd_routine,') = ',zsfx 
    221          !!IF( ABS( zhfx ) > zt_sill )   WRITE(numout,*) 'violation hfx  [GW]           (',cd_routine,') = ',zhfx 
     197      !!zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  & 
     198      !!   &                              ) * e1e2t ) 
     199 
     200      ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     201      zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
     202 
     203      IF( lwp ) THEN 
     204         IF( ABS(zdiag_mass) > zchk_m * rn_icechk * zarea ) & 
     205            &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 
     206         IF( ABS(zdiag_salt) > zchk_s * rn_icechk * zarea ) & 
     207            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice 
     208         !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice 
    222209      ENDIF 
    223210      ! 
    224211   END SUBROUTINE ice_cons_final 
    225212 
     213   SUBROUTINE ice_cons2D( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft ) 
     214      !!------------------------------------------------------------------- 
     215      !!                       ***  ROUTINE ice_cons2D *** 
     216      !! 
     217      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine 
     218      !!                     + test if ice concentration and volume are > 0 
     219      !! 
     220      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
     221      !!              It stops the code if there is a violation of conservation at any gridcell 
     222      !!------------------------------------------------------------------- 
     223      INTEGER         , INTENT(in) ::   icount        ! called at: =0 the begining of the routine, =1  the end 
     224      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine 
     225      REAL(wp)        , DIMENSION(jpi,jpj), INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
     226      !! 
     227      REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
     228         &                              zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax   
     229      INTEGER ::   jl, jk 
     230      LOGICAL ::   ll_stop_m = .FALSE. 
     231      LOGICAL ::   ll_stop_s = .FALSE. 
     232      LOGICAL ::   ll_stop_t = .FALSE. 
     233      CHARACTER(len=120) ::   clnam   ! filename for the output 
     234      !!------------------------------------------------------------------- 
     235      ! 
     236      IF( icount == 0 ) THEN 
     237 
     238         pdiag_v = SUM( v_i  * rhoi + v_s * rhos, dim=3 ) 
     239         pdiag_s = SUM( sv_i * rhoi             , dim=3 ) 
     240         pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 
     241 
     242         ! mass flux 
     243         pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd  +  & 
     244            &       wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr 
     245         ! salt flux 
     246         pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam  
     247         ! heat flux 
     248         pdiag_ft =   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  &  
     249            &       - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr 
     250 
     251      ELSEIF( icount == 1 ) THEN 
     252 
     253         ! -- mass diag -- ! 
     254         zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos, dim=3 ) - pdiag_v ) * r1_rdtice                             & 
     255            &         + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
     256            &             wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr )           & 
     257            &         - pdiag_fv 
     258         IF( MAXVAL( ABS(zdiag_mass) ) > zchk2D_m * rn_icechk )   ll_stop_m = .TRUE. 
     259         ! 
     260         ! -- salt diag -- ! 
     261         zdiag_salt =   ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice                                                  & 
     262            &         + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) & 
     263            &         - pdiag_fs 
     264         IF( MAXVAL( ABS(zdiag_salt) ) > zchk2D_s * rn_icechk )   ll_stop_s = .TRUE. 
     265         ! 
     266         ! -- heat diag -- ! 
     267         zdiag_heat =   ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_rdtice & 
     268            &         + (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                                &  
     269            &            - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr )                                        & 
     270            &         - pdiag_ft 
     271         IF( MAXVAL( ABS(zdiag_heat) ) > zchk2D_t * rn_icechk )   ll_stop_t = .TRUE. 
     272         ! 
     273         ! -- other diags -- ! 
     274         ! a_i < 0 
     275         zdiag_amin(:,:) = 0._wp 
     276         DO jl = 1, jpl 
     277            WHERE( a_i(:,:,jl) < 0._wp )   zdiag_amin(:,:) = 1._wp 
     278         ENDDO 
     279         ! v_i < 0 
     280         zdiag_vmin(:,:) = 0._wp 
     281         DO jl = 1, jpl 
     282            WHERE( v_i(:,:,jl) < 0._wp )   zdiag_vmin(:,:) = 1._wp 
     283         ENDDO 
     284         ! s_i < 0 
     285         zdiag_smin(:,:) = 0._wp 
     286         DO jl = 1, jpl 
     287            WHERE( s_i(:,:,jl) < 0._wp )   zdiag_smin(:,:) = 1._wp 
     288         ENDDO 
     289         ! e_i < 0 
     290         zdiag_emin(:,:) = 0._wp 
     291         DO jl = 1, jpl 
     292            DO jk = 1, nlay_i 
     293               WHERE( e_i(:,:,jk,jl) < 0._wp )   zdiag_emin(:,:) = 1._wp 
     294            ENDDO 
     295         ENDDO 
     296         ! a_i > amax 
     297         !WHERE( SUM( a_i, dim=3 ) > ( MAX( rn_amax_n, rn_amax_s ) + epsi10 )   ;   zdiag_amax(:,:) = 1._wp 
     298         !ELSEWHERE                                                             ;   zdiag_amax(:,:) = 0._wp 
     299         !END WHERE 
     300 
     301         IF( ll_stop_m .OR. ll_stop_s .OR. ll_stop_t ) THEN 
     302            clnam = 'diag_ice_conservation_'//cd_routine 
     303            CALL ice_cons_wri( clnam, zdiag_mass, zdiag_salt, zdiag_heat, zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin ) 
     304         ENDIF 
     305 
     306         IF( ll_stop_m )   CALL ctl_stop( 'STOP', cd_routine//': ice mass conservation issue' ) 
     307         IF( ll_stop_s )   CALL ctl_stop( 'STOP', cd_routine//': ice salt conservation issue' ) 
     308         IF( ll_stop_t )   CALL ctl_stop( 'STOP', cd_routine//': ice heat conservation issue' ) 
     309          
     310      ENDIF 
     311 
     312   END SUBROUTINE ice_cons2D 
     313 
     314   SUBROUTINE ice_cons_wri( cdfile_name, pdiag_mass, pdiag_salt, pdiag_heat, pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin ) 
     315      !!--------------------------------------------------------------------- 
     316      !!                 ***  ROUTINE ice_cons_wri  *** 
     317      !!         
     318      !! ** Purpose :   create a NetCDF file named cdfile_name which contains  
     319      !!                the instantaneous fields when conservation issue occurs 
     320      !! 
     321      !! ** Method  :   NetCDF files using ioipsl 
     322      !!---------------------------------------------------------------------- 
     323      CHARACTER(len=*), INTENT( in ) ::   cdfile_name      ! name of the file created 
     324      REAL(wp), DIMENSION(:,:), INTENT( in ) ::   pdiag_mass, pdiag_salt, pdiag_heat, & 
     325         &                                        pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax   
     326      !! 
     327      INTEGER ::   inum 
     328      !!---------------------------------------------------------------------- 
     329      !  
     330      IF(lwp) WRITE(numout,*) 
     331      IF(lwp) WRITE(numout,*) 'ice_cons_wri : single instantaneous ice state' 
     332      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  named :', cdfile_name, '...nc' 
     333      IF(lwp) WRITE(numout,*)                 
     334 
     335      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
     336       
     337      CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 )    ! ice mass spurious lost/gain 
     338      CALL iom_rstput( 0, 0, inum, 'cons_salt', pdiag_salt(:,:) , ktype = jp_r8 )    ! ice salt spurious lost/gain 
     339      CALL iom_rstput( 0, 0, inum, 'cons_heat', pdiag_heat(:,:) , ktype = jp_r8 )    ! ice heat spurious lost/gain 
     340      ! other diags 
     341      CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 )    !  
     342      CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 )    !  
     343      CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 )    !  
     344      CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 )    !  
     345       
     346      CALL iom_close( inum ) 
     347 
     348   END SUBROUTINE ice_cons_wri 
    226349    
    227350   SUBROUTINE ice_ctl( kt ) 
     
    246369      ialert_id = 2 ! reference number of this alert 
    247370      cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert 
    248  
    249371      DO jl = 1, jpl 
    250372         DO jj = 1, jpj 
    251373            DO ji = 1, jpi 
    252374               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    253                   !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    254                   !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    255                   !WRITE(numout,*) ' Point - category', ji, jj, jl 
    256                   !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl) 
    257                   !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl) 
     375                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    258376                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    259377               ENDIF 
     
    269387         DO ji = 1, jpi 
    270388            IF(   h_i(ji,jj,jl)  >  50._wp   ) THEN 
     389               WRITE(numout,*) ' ALERTE 3 :   Very thick ice' 
    271390               !CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    272391               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     
    280399      DO jj = 1, jpj 
    281400         DO ji = 1, jpi 
    282             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
     401            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2.  .AND.  & 
    283402               &  at_i(ji,jj) > 0._wp   ) THEN 
     403               WRITE(numout,*) ' ALERTE 4 :   Very fast ice' 
    284404               !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    285                !WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    286                !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
    287                !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
    288                !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    289                !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    290                !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    291                !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
    292                !WRITE(numout,*)  
     405               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     406            ENDIF 
     407         END DO 
     408      END DO 
     409 
     410      ! Alert on salt flux 
     411      ialert_id = 5 ! reference number of this alert 
     412      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
     413      DO jj = 1, jpj 
     414         DO ji = 1, jpi 
     415            IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
     416               WRITE(numout,*) ' ALERTE 5 :   High salt flux' 
     417               !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    293418               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    294419            ENDIF 
     
    302427         DO ji = 1, jpi 
    303428            IF(   tmask(ji,jj,1) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
     429               WRITE(numout,*) ' ALERTE 6 :   Ice on continents' 
    304430               !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    305                !WRITE(numout,*) ' masks s, u, v        : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1)  
    306                !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    307                !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    308                !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
    309                !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
    310                !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
    311                !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
    312                !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
    313                ! 
    314431               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    315432            ENDIF 
     
    325442            DO ji = 1, jpi 
    326443               IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     444                  WRITE(numout,*) ' ALERTE 7 :   Very fresh ice' 
    327445!                 CALL ice_prt(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    328 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    329 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    330 !                 WRITE(numout,*)  
    331446                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    332447               ENDIF 
     
    335450      END DO 
    336451! 
     452      ! Alert if qns very big 
     453      ialert_id = 8 ! reference number of this alert 
     454      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
     455      DO jj = 1, jpj 
     456         DO ji = 1, jpi 
     457            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     458               ! 
     459               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
     460               !CALL ice_prt( kt, ji, jj, 2, '   ') 
     461               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     462               ! 
     463            ENDIF 
     464         END DO 
     465      END DO 
     466      !+++++ 
    337467 
    338468!     ! Alert if too old ice 
     
    345475                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    346476                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
     477                  WRITE(numout,*) ' ALERTE 9 :   Wrong ice age' 
    347478                  !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    348479                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     
    351482         END DO 
    352483      END DO 
    353   
    354       ! Alert on salt flux 
    355       ialert_id = 5 ! reference number of this alert 
    356       cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
    357       DO jj = 1, jpj 
    358          DO ji = 1, jpi 
    359             IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
    360                !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    361                !DO jl = 1, jpl 
    362                   !WRITE(numout,*) ' Category no: ', jl 
    363                   !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)    
    364                   !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)    
    365                   !WRITE(numout,*) ' ' 
    366                !END DO 
    367                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    368             ENDIF 
    369          END DO 
    370       END DO 
    371  
    372       ! Alert if qns very big 
    373       ialert_id = 8 ! reference number of this alert 
    374       cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
    375       DO jj = 1, jpj 
    376          DO ji = 1, jpi 
    377             IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    378                ! 
    379                !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    380                !WRITE(numout,*) ' ji, jj    : ', ji, jj 
    381                !WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    382                !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    383                !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    384                ! 
    385                !CALL ice_prt( kt, ji, jj, 2, '   ') 
    386                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    387                ! 
    388             ENDIF 
    389          END DO 
    390       END DO 
    391       !+++++ 
    392   
     484   
    393485      ! Alert if very warm ice 
    394486      ialert_id = 10 ! reference number of this alert 
     
    400492               DO ji = 1, jpi 
    401493                  ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0 
    402                   IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    403                      &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    404                      !WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
    405                      !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
    406                      !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
    407                      !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
    408                      !WRITE(numout,*) ' sz_i: ', sz_i(ji,jj,jk,jl) 
    409                      !WRITE(numout,*) ' ztmelts : ', ztmelts 
    410                      inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     494                  IF( t_i(ji,jj,jk,jl) > ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
     495                     &                            .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
     496                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
     497                    inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    411498                  ENDIF 
    412499               END DO 
     
    435522   END SUBROUTINE ice_ctl 
    436523  
    437     
    438524   SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 ) 
    439525      !!------------------------------------------------------------------- 
     
    443529      !!                in ocean.ouput  
    444530      !!                3 possibilities exist  
    445       !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1) 
     531      !!                n = 1/-1 -> simple ice state 
    446532      !!                n = 2    -> exhaustive state 
    447533      !!                n = 3    -> ice/ocean salt fluxes 
     
    482568               WRITE(numout,*) ' - Cell values ' 
    483569               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    484                WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj) 
    485570               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
     571               WRITE(numout,*) ' ato_i         : ', ato_i(ji,jj)        
    486572               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    487573               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
     
    503589               END DO 
    504590            ENDIF 
    505             IF( kn == -1 ) THEN 
    506                WRITE(numout,*) ' Mechanical Check ************** ' 
    507                WRITE(numout,*) ' Check what means ice divergence ' 
    508                WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 
    509                WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj) 
    510                WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj) 
    511                WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 
    512             ENDIF 
    513              
    514591 
    515592            !-------------------- 
     
    525602               WRITE(numout,*) ' - Cell values ' 
    526603               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    527                WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj) 
    528604               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    529605               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     
    624700      !! 
    625701      !!------------------------------------------------------------------- 
    626       CHARACTER(len=*), INTENT(in)  :: cd_routine    ! name of the routine 
    627       INTEGER                       :: jk, jl        ! dummy loop indices 
     702      CHARACTER(len=*), INTENT(in) ::  cd_routine    ! name of the routine 
     703      INTEGER                      ::  jk, jl        ! dummy loop indices 
    628704       
    629705      CALL prt_ctl_info(' ========== ') 
     
    684760       
    685761   END SUBROUTINE ice_prt3D 
    686  
     762       
    687763#else 
    688764   !!---------------------------------------------------------------------- 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedia.F90

    r11371 r11501  
    175175      INTEGER            ::   ios, ierror   ! local integer 
    176176      !! 
    177       NAMELIST/namdia/ ln_icediachk, ln_icediahsb, ln_icectl, iiceprt, jiceprt   
     177      NAMELIST/namdia/ ln_icediachk, rn_icechk, ln_icediahsb, ln_icectl, iiceprt, jiceprt   
    178178      !!---------------------------------------------------------------------- 
    179179      ! 
     
    191191         WRITE(numout,*) ' ~~~~~~~~~~~' 
    192192         WRITE(numout,*) '   Namelist namdia:' 
    193          WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_icediachk = ', ln_icediachk 
    194          WRITE(numout,*) '      Output          heat/mass/salt budget      ln_icediahsb = ', ln_icediahsb 
    195          WRITE(numout,*) '      control prints for a given grid point      ln_icectl    = ', ln_icectl 
    196          WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')' 
     193         WRITE(numout,*) '      Diagnose online heat/mass/salt conservation ln_icediachk = ', ln_icediachk 
     194         WRITE(numout,*) '         threshold for conservation               rn_icechk    = ', rn_icechk 
     195         WRITE(numout,*) '      Output          heat/mass/salt budget       ln_icediahsb = ', ln_icediahsb 
     196         WRITE(numout,*) '      control prints for a given grid point       ln_icectl    = ', ln_icectl 
     197         WRITE(numout,*) '         chosen grid point position          (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')' 
    197198      ENDIF 
    198199      !       
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rdgrft.F90

    r11317 r11501  
    145145      IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing 
    146146      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     147      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    147148 
    148149      IF( kt == nit000 ) THEN 
     
    276277 
    277278      ! controls 
     279      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints 
    278280      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    279       IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints 
     281      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    280282      IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing 
    281283      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg.F90

    r11377 r11501  
    6464      IF( ln_timing    )   CALL timing_start('icedyn_rhg')                                                             ! timing 
    6565      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     66      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    6667      ! 
    6768      IF( kt == nit000 .AND. lwp ) THEN 
     
    8788      ! 
    8889      ! controls 
     90      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints 
    8991      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    90       IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints 
     92      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    9193      IF( ln_timing    )   CALL timing_stop ('icedyn_rhg')                                                             ! timing 
    9294      ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/iceitd.F90

    r11317 r11501  
    8888 
    8989      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     90      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    9091 
    9192      !----------------------------------------------------------------------------------------------- 
     
    316317      ! 
    317318      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     319      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    318320      ! 
    319321   END SUBROUTINE ice_itd_rem 
     
    586588      ! 
    587589      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     590      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    588591      ! 
    589592      jdonor(:,:) = 0 
     
    664667      ! 
    665668      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     669      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    666670      ! 
    667671   END SUBROUTINE ice_itd_reb 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd.F90

    r11317 r11501  
    9595      IF( ln_timing    )   CALL timing_start('icethd')                                                             ! timing 
    9696      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     97      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    9798 
    9899      IF( kt == nit000 .AND. lwp ) THEN 
     
    243244      ! 
    244245      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     246      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    245247      !                    
    246248      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd_do.F90

    r11317 r11501  
    113113 
    114114      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft ) 
     115      IF( ln_icediachk )   CALL ice_cons2D  ( 0, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft ) 
    115116 
    116117      at_i(:,:) = SUM( a_i, dim=3 ) 
     
    420421      ! 
    421422      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
     423      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    422424      ! 
    423425   END SUBROUTINE ice_thd_do 
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyice.F90

    r11210 r11501  
    6363      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing 
    6464      IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     65      IF( ln_icediachk )   CALL ice_cons2D  (0,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    6566      ! 
    6667      CALL ice_var_glo2eqv 
     
    109110      ! 
    110111      ! controls 
     112      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
    111113      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    112       IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints 
     114      IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
    113115      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing 
    114116      ! 
Note: See TracChangeset for help on using the changeset viewer.