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.
Changeset 11586 for NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE – NEMO

Ignore:
Timestamp:
2019-09-20T17:28:02+02:00 (5 years ago)
Author:
gsamson
Message:

dev_r11265_ABL : see #2131

  • merge HPC-13_IRRMANN_BDY_optimization branch @ r11535 (last commit) with dev_r11265_ABL branch @ r11414 (except doc directory)
  • change ORCA2 results due to changes in HPC-13_IRRMANN_BDY_optimization branch
Location:
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/ice.F90

    r11413 r11586  
    110110   !! bv_i        |      -      |    relative brine volume        | ???   |  
    111111   !! at_ip       |      -      |    Total ice pond concentration |       | 
     112   !! hm_ip       |      -      |    Mean ice pond depth          | m     | 
    112113   !! vt_ip       |      -      |    Total ice pond vol. per unit area| m | 
    113114   !!===================================================================== 
    … …  
    188189 
    189190   !                                     !!** ice-ponds namelist (namthd_pnd) 
     191   LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F) 
    190192   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012 
    191193   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
    … …  
    196198   !                                     !!** ice-diagnostics namelist (namdia) ** 
    197199   LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F) 
     200   REAL(wp), PUBLIC ::   rn_icechk_cel    !: rate of ice spuriously gained/lost (at any gridcell) 
     201   REAL(wp), PUBLIC ::   rn_icechk_glo    !: rate of ice spuriously gained/lost (globally) 
    198202   LOGICAL , PUBLIC ::   ln_icediahsb     !: flag for ice diag (T) or not (F) 
    199203   LOGICAL , PUBLIC ::   ln_icectl        !: flag for sea-ice points output (T) or not (F) 
    … …  
    330334 
    331335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond fraction 
     336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_ip      !: mean melt pond depth                     [m] 
    332337   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per unit area     [m] 
    333338 
    … …  
    351356   !! * Ice diagnostics 
    352357   !!---------------------------------------------------------------------- 
    353    ! thd refers to changes induced by thermodynamics 
    354    ! trp   ''         ''     ''       advection (transport of ice) 
    355    ! 
    356358   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi   !: transport of ice volume 
    357359   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs   !: transport of snw volume 
    … …  
    365367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw     !: snw volume variation   [m/s]  
    366368 
     369   !!---------------------------------------------------------------------- 
     370   !! * Ice conservation 
     371   !!---------------------------------------------------------------------- 
     372   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_v        !: conservation of ice volume 
     373   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_s        !: conservation of ice salt 
     374   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_t        !: conservation of ice heat 
     375   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fv       !: conservation of ice volume 
     376   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fs       !: conservation of ice salt 
     377   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_ft       !: conservation of ice heat 
    367378   ! 
    368379   !!---------------------------------------------------------------------- 
    … …  
    389400      INTEGER :: ice_alloc 
    390401      ! 
    391       INTEGER :: ierr(15), ii 
     402      INTEGER :: ierr(16), ii 
    392403      !!----------------------------------------------------------------- 
    393404      ierr(:) = 0 
    … …  
    440451 
    441452      ii = ii + 1 
    442       ALLOCATE( at_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) ) 
     453      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) ) 
    443454 
    444455      ! * Old values of global variables 
    … …  
    461472         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) ) 
    462473 
     474      ! * Ice conservation 
     475      ii = ii + 1 
     476      ALLOCATE( diag_v (jpi,jpj) , diag_s (jpi,jpj) , diag_t (jpi,jpj),   &  
     477         &      diag_fv(jpi,jpj) , diag_fs(jpi,jpj) , diag_ft(jpi,jpj), STAT=ierr(ii) ) 
     478       
    463479      ! * SIMIP diagnostics 
    464480      ii = ii + 1 

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icecor.F90

    r11413 r11586  
    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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icectl.F90

    r11413 r11586  
    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 year  spuriously gained/lost 
     49   REAL(wp), PARAMETER ::   zchk_s   = 1.e-4   ! g/m2/s  <=> 1mm of ice per year  spuriously gained/lost (considering s=10g/kg) 
     50   REAL(wp), PARAMETER ::   zchk_t   = 3.      ! W/m2    <=> 1mm of ice per year  spuriously gained/lost (considering Lf=3e5J/kg) 
     51    
    4352   !! * Substitutions 
    4453#  include "vectopt_loop_substitute.h90" 
    … …  
    5968      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    6069      !!              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 
     70      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine violations are set to 
    6271      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    6372      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
    … …  
    6877      REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
    6978      !! 
    70       REAL(wp) ::   zv, zs, zt, zfs, zfv, zft 
    71       REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin 
     79      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
     80         &          zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 
    7281      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 
     82      REAL(wp) ::   zarea 
    7583      !!------------------------------------------------------------------- 
    7684      ! 
    7785      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 
     86 
     87         pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) 
     88         pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) 
     89         pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 
     90 
     91         ! mass flux 
     92         pdiag_fv = glob_sum( 'icectl',  & 
     93            &                         ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
     94            &                           wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t ) 
     95         ! salt flux 
     96         pdiag_fs = glob_sum( 'icectl',  & 
     97            &                         ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + & 
     98            &                           sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) 
     99         ! heat flux 
     100         pdiag_ft = glob_sum( 'icectl',  & 
     101            &                         (   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  & 
     102            &                           - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) 
     103 
     104      ELSEIF( icount == 1 ) THEN 
     105 
     106         ! -- mass diag -- ! 
     107         zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_rdtice       & 
     108            &         + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn +       & 
     109            &                                 wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & 
     110            &                                 wfx_ice_sub + wfx_spr ) * e1e2t )                                           & 
     111            &         - pdiag_fv 
    85112         ! 
    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  
     113         ! -- salt diag -- ! 
     114         zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice  & 
     115            &         + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni +           & 
     116            &                                 sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) & 
     117            &         - pdiag_fs 
    91118         ! 
    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 
     119         ! -- heat diag -- ! 
     120         zdiag_heat = ( glob_sum( 'icectl', ( SUM(SUM(e_i, dim=4), dim=3) + SUM(SUM(e_s, dim=4), dim=3) ) * e1e2t ) - pdiag_t & 
     121            &         ) * r1_rdtice                                                                                           & 
     122            &         + glob_sum( 'icectl', (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                      & 
     123            &                                - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )                    & 
     124            &         - pdiag_ft 
     125 
     126         ! -- min/max diag -- ! 
     127         zdiag_amax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     128         zdiag_vmin  = glob_min( 'icectl', v_i ) 
     129         zdiag_amin  = glob_min( 'icectl', a_i ) 
     130         zdiag_smin  = glob_min( 'icectl', sv_i ) 
     131         zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
     132         zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
     133 
     134         ! -- advection scheme is conservative? -- ! 
     135         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 
     136         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0 
     137 
     138         ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     139         zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
     140 
     141         IF( lwp ) THEN 
    157142            ! 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 
     143            IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 
     144               &                   WRITE(numout,*)   cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 
     145            IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
     146               &                   WRITE(numout,*)   cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice 
     147            IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) & 
     148               &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice 
     149            ! check negative values 
     150            IF( zdiag_vmin  < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i < 0         = ',zdiag_vmin 
     151            IF( zdiag_amin  < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i < 0         = ',zdiag_amin 
     152            IF( zdiag_smin  < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i < 0         = ',zdiag_smin 
     153            IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i < 0         = ',zdiag_eimin 
     154            IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s < 0         = ',zdiag_esmin 
    161155            ! 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 
     156            IF( zdiag_amax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
     157               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
     158            ! check if advection scheme is conservative 
     159            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     160               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
    175161         ENDIF 
    176162         ! 
    … …  
    179165   END SUBROUTINE ice_cons_hsm 
    180166 
    181  
    182167   SUBROUTINE ice_cons_final( cd_routine ) 
    183168      !!------------------------------------------------------------------- 
    … …  
    188173      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
    189174      !!              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 
     175      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine the violation are set to 
    191176      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years. 
    192177      !!              For salt and heat thresholds, ice is considered to have a salinity of 10  
    193178      !!              and a heat content of 3e5 J/kg (=latent heat of fusion)  
    194179      !!------------------------------------------------------------------- 
    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 
     180      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine 
     181      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat 
     182      REAL(wp) ::   zarea 
    199183      !!------------------------------------------------------------------- 
    200184 
    201185      ! 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 
     186      ! -- mass diag -- ! 
     187      zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) 
     188 
     189      ! -- salt diag -- ! 
     190      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) 
     191 
     192      ! -- heat diag -- ! 
    208193      ! 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 
     194      !!zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  & 
     195      !!   &                              ) * e1e2t ) 
     196 
     197      ! ice area (+epsi10 to set a threshold > 0 when there is no ice)  
     198      zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) 
     199 
     200      IF( lwp ) THEN 
     201         IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) & 
     202            &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice 
     203         IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) & 
     204            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice 
     205         !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice 
    222206      ENDIF 
    223207      ! 
    224208   END SUBROUTINE ice_cons_final 
    225209 
     210   SUBROUTINE ice_cons2D( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft ) 
     211      !!------------------------------------------------------------------- 
     212      !!                       ***  ROUTINE ice_cons2D *** 
     213      !! 
     214      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine 
     215      !!                     + test if ice concentration and volume are > 0 
     216      !! 
     217      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true 
     218      !!              It stops the code if there is a violation of conservation at any gridcell 
     219      !!------------------------------------------------------------------- 
     220      INTEGER         , INTENT(in) ::   icount        ! called at: =0 the begining of the routine, =1  the end 
     221      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine 
     222      REAL(wp)        , DIMENSION(jpi,jpj), INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft 
     223      !! 
     224      REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
     225         &                              zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax   
     226      INTEGER ::   jl, jk 
     227      LOGICAL ::   ll_stop_m = .FALSE. 
     228      LOGICAL ::   ll_stop_s = .FALSE. 
     229      LOGICAL ::   ll_stop_t = .FALSE. 
     230      CHARACTER(len=120) ::   clnam   ! filename for the output 
     231      !!------------------------------------------------------------------- 
     232      ! 
     233      IF( icount == 0 ) THEN 
     234 
     235         pdiag_v = SUM( v_i  * rhoi + v_s * rhos, dim=3 ) 
     236         pdiag_s = SUM( sv_i * rhoi             , dim=3 ) 
     237         pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 
     238 
     239         ! mass flux 
     240         pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd  +  & 
     241            &       wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr 
     242         ! salt flux 
     243         pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam  
     244         ! heat flux 
     245         pdiag_ft =   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  &  
     246            &       - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr 
     247 
     248      ELSEIF( icount == 1 ) THEN 
     249 
     250         ! -- mass diag -- ! 
     251         zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos, dim=3 ) - pdiag_v ) * r1_rdtice                             & 
     252            &         + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
     253            &             wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr )           & 
     254            &         - pdiag_fv 
     255         IF( MAXVAL( ABS(zdiag_mass) ) > zchk_m * rn_icechk_cel )   ll_stop_m = .TRUE. 
     256         ! 
     257         ! -- salt diag -- ! 
     258         zdiag_salt =   ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice                                                  & 
     259            &         + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) & 
     260            &         - pdiag_fs 
     261         IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel )   ll_stop_s = .TRUE. 
     262         ! 
     263         ! -- heat diag -- ! 
     264         zdiag_heat =   ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_rdtice & 
     265            &         + (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                                &  
     266            &            - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr )                                        & 
     267            &         - pdiag_ft 
     268         IF( MAXVAL( ABS(zdiag_heat) ) > zchk_t * rn_icechk_cel )   ll_stop_t = .TRUE. 
     269         ! 
     270         ! -- other diags -- ! 
     271         ! a_i < 0 
     272         zdiag_amin(:,:) = 0._wp 
     273         DO jl = 1, jpl 
     274            WHERE( a_i(:,:,jl) < 0._wp )   zdiag_amin(:,:) = 1._wp 
     275         ENDDO 
     276         ! v_i < 0 
     277         zdiag_vmin(:,:) = 0._wp 
     278         DO jl = 1, jpl 
     279            WHERE( v_i(:,:,jl) < 0._wp )   zdiag_vmin(:,:) = 1._wp 
     280         ENDDO 
     281         ! s_i < 0 
     282         zdiag_smin(:,:) = 0._wp 
     283         DO jl = 1, jpl 
     284            WHERE( s_i(:,:,jl) < 0._wp )   zdiag_smin(:,:) = 1._wp 
     285         ENDDO 
     286         ! e_i < 0 
     287         zdiag_emin(:,:) = 0._wp 
     288         DO jl = 1, jpl 
     289            DO jk = 1, nlay_i 
     290               WHERE( e_i(:,:,jk,jl) < 0._wp )   zdiag_emin(:,:) = 1._wp 
     291            ENDDO 
     292         ENDDO 
     293         ! a_i > amax 
     294         !WHERE( SUM( a_i, dim=3 ) > ( MAX( rn_amax_n, rn_amax_s ) + epsi10 )   ;   zdiag_amax(:,:) = 1._wp 
     295         !ELSEWHERE                                                             ;   zdiag_amax(:,:) = 0._wp 
     296         !END WHERE 
     297 
     298         IF( ll_stop_m .OR. ll_stop_s .OR. ll_stop_t ) THEN 
     299            clnam = 'diag_ice_conservation_'//cd_routine 
     300            CALL ice_cons_wri( clnam, zdiag_mass, zdiag_salt, zdiag_heat, zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin ) 
     301         ENDIF 
     302 
     303         IF( ll_stop_m )   CALL ctl_stop( 'STOP', cd_routine//': ice mass conservation issue' ) 
     304         IF( ll_stop_s )   CALL ctl_stop( 'STOP', cd_routine//': ice salt conservation issue' ) 
     305         IF( ll_stop_t )   CALL ctl_stop( 'STOP', cd_routine//': ice heat conservation issue' ) 
     306          
     307      ENDIF 
     308 
     309   END SUBROUTINE ice_cons2D 
     310 
     311   SUBROUTINE ice_cons_wri( cdfile_name, pdiag_mass, pdiag_salt, pdiag_heat, pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin ) 
     312      !!--------------------------------------------------------------------- 
     313      !!                 ***  ROUTINE ice_cons_wri  *** 
     314      !!         
     315      !! ** Purpose :   create a NetCDF file named cdfile_name which contains  
     316      !!                the instantaneous fields when conservation issue occurs 
     317      !! 
     318      !! ** Method  :   NetCDF files using ioipsl 
     319      !!---------------------------------------------------------------------- 
     320      CHARACTER(len=*), INTENT( in ) ::   cdfile_name      ! name of the file created 
     321      REAL(wp), DIMENSION(:,:), INTENT( in ) ::   pdiag_mass, pdiag_salt, pdiag_heat, & 
     322         &                                        pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax   
     323      !! 
     324      INTEGER ::   inum 
     325      !!---------------------------------------------------------------------- 
     326      !  
     327      IF(lwp) WRITE(numout,*) 
     328      IF(lwp) WRITE(numout,*) 'ice_cons_wri : single instantaneous ice state' 
     329      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  named :', cdfile_name, '...nc' 
     330      IF(lwp) WRITE(numout,*)                 
     331 
     332      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
     333       
     334      CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 )    ! ice mass spurious lost/gain 
     335      CALL iom_rstput( 0, 0, inum, 'cons_salt', pdiag_salt(:,:) , ktype = jp_r8 )    ! ice salt spurious lost/gain 
     336      CALL iom_rstput( 0, 0, inum, 'cons_heat', pdiag_heat(:,:) , ktype = jp_r8 )    ! ice heat spurious lost/gain 
     337      ! other diags 
     338      CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 )    !  
     339      CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 )    !  
     340      CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 )    !  
     341      CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 )    !  
     342       
     343      CALL iom_close( inum ) 
     344 
     345   END SUBROUTINE ice_cons_wri 
    226346    
    227347   SUBROUTINE ice_ctl( kt ) 
    … …  
    246366      ialert_id = 2 ! reference number of this alert 
    247367      cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert 
    248  
    249368      DO jl = 1, jpl 
    250369         DO jj = 1, jpj 
    251370            DO ji = 1, jpi 
    252371               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) 
     372                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    258373                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    259374               ENDIF 
    … …  
    269384         DO ji = 1, jpi 
    270385            IF(   h_i(ji,jj,jl)  >  50._wp   ) THEN 
     386               WRITE(numout,*) ' ALERTE 3 :   Very thick ice' 
    271387               !CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    272388               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    … …  
    280396      DO jj = 1, jpj 
    281397         DO ji = 1, jpi 
    282             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
     398            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2.  .AND.  & 
    283399               &  at_i(ji,jj) > 0._wp   ) THEN 
     400               WRITE(numout,*) ' ALERTE 4 :   Very fast ice' 
    284401               !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,*)  
     402               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     403            ENDIF 
     404         END DO 
     405      END DO 
     406 
     407      ! Alert on salt flux 
     408      ialert_id = 5 ! reference number of this alert 
     409      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
     410      DO jj = 1, jpj 
     411         DO ji = 1, jpi 
     412            IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
     413               WRITE(numout,*) ' ALERTE 5 :   High salt flux' 
     414               !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    293415               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    294416            ENDIF 
    … …  
    302424         DO ji = 1, jpi 
    303425            IF(   tmask(ji,jj,1) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
     426               WRITE(numout,*) ' ALERTE 6 :   Ice on continents' 
    304427               !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                ! 
    314428               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    315429            ENDIF 
    … …  
    325439            DO ji = 1, jpi 
    326440               IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     441                  WRITE(numout,*) ' ALERTE 7 :   Very fresh ice' 
    327442!                 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,*)  
    331443                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    332444               ENDIF 
    … …  
    335447      END DO 
    336448! 
     449      ! Alert if qns very big 
     450      ialert_id = 8 ! reference number of this alert 
     451      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
     452      DO jj = 1, jpj 
     453         DO ji = 1, jpi 
     454            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     455               ! 
     456               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
     457               !CALL ice_prt( kt, ji, jj, 2, '   ') 
     458               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     459               ! 
     460            ENDIF 
     461         END DO 
     462      END DO 
     463      !+++++ 
    337464 
    338465!     ! Alert if too old ice 
    … …  
    345472                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    346473                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
     474                  WRITE(numout,*) ' ALERTE 9 :   Wrong ice age' 
    347475                  !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    348476                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    … …  
    351479         END DO 
    352480      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   
     481   
    393482      ! Alert if very warm ice 
    394483      ialert_id = 10 ! reference number of this alert 
    … …  
    400489               DO ji = 1, jpi 
    401490                  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 
     491                  IF( t_i(ji,jj,jk,jl) > ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
     492                     &                            .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
     493                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
     494                    inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    411495                  ENDIF 
    412496               END DO 
    … …  
    435519   END SUBROUTINE ice_ctl 
    436520  
    437     
    438521   SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 ) 
    439522      !!------------------------------------------------------------------- 
    … …  
    443526      !!                in ocean.ouput  
    444527      !!                3 possibilities exist  
    445       !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1) 
     528      !!                n = 1/-1 -> simple ice state 
    446529      !!                n = 2    -> exhaustive state 
    447530      !!                n = 3    -> ice/ocean salt fluxes 
    … …  
    482565               WRITE(numout,*) ' - Cell values ' 
    483566               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    484                WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj) 
    485567               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
     568               WRITE(numout,*) ' ato_i         : ', ato_i(ji,jj)        
    486569               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    487570               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
    … …  
    503586               END DO 
    504587            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              
    514588 
    515589            !-------------------- 
    … …  
    525599               WRITE(numout,*) ' - Cell values ' 
    526600               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    527                WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj) 
    528601               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    529602               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    … …  
    624697      !! 
    625698      !!------------------------------------------------------------------- 
    626       CHARACTER(len=*), INTENT(in)  :: cd_routine    ! name of the routine 
    627       INTEGER                       :: jk, jl        ! dummy loop indices 
     699      CHARACTER(len=*), INTENT(in) ::  cd_routine    ! name of the routine 
     700      INTEGER                      ::  jk, jl        ! dummy loop indices 
    628701       
    629702      CALL prt_ctl_info(' ========== ') 
    … …  
    684757       
    685758   END SUBROUTINE ice_prt3D 
    686  
     759       
    687760#else 
    688761   !!---------------------------------------------------------------------- 

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedia.F90

    r11413 r11586  
    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_cel, rn_icechk_glo, 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 (gridcell)    rn_icechk_cel = ', rn_icechk_cel 
     195         WRITE(numout,*) '         threshold for conservation (global)      rn_icechk_glo = ', rn_icechk_glo 
     196         WRITE(numout,*) '      Output          heat/mass/salt budget       ln_icediahsb  = ', ln_icediahsb 
     197         WRITE(numout,*) '      control prints for a given grid point       ln_icectl     = ', ln_icectl 
     198         WRITE(numout,*) '         chosen grid point position          (iiceprt,jiceprt)  = (', iiceprt,',', jiceprt,')' 
    197199      ENDIF 
    198200      !       

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rdgrft.F90

    r11348 r11586  
    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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg.F90

    r11413 r11586  
    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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icedyn_rhg_evp.F90

    r11413 r11586  
    121121      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    122122      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
    123       REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars 
     123      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    124124      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    125125      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
    … …  
    130130      REAL(wp) ::   zshear, zdum1, zdum2 
    131131      ! 
    132       REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors 
    133132      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    134133      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    … …  
    244243      CALL ice_strength 
    245244 
    246       ! scale factors 
    247       DO jj = 2, jpjm1 
    248          DO ji = fs_2, fs_jpim1 
    249             z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) ) 
    250             z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) ) 
    251          END DO 
    252       END DO 
    253  
    254245      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    255246      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
    … …  
    280271 
    281272            ! Ocean currents at U-V points 
    282             v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
    283                &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    284              
    285             u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
    286                &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     273            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 
     274            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 
    287275 
    288276            ! Coriolis at T points (m*f) 
    … …  
    459447                  &                  ) * r1_e1e2v(ji,jj) 
    460448               ! 
    461                !                !--- u_ice at V point 
    462                u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    463                   &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     449               !                !--- ice currents at U-V point 
     450               v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 
     451               u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 
    464452               ! 
    465                !                !--- v_ice at U point 
    466                v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    467                   &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    468453            END DO 
    469454         END DO 
    … …  
    494479                  ! 
    495480                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    496                   zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    497                   ! 
    498                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    499                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     481                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     482                  ! 
     483                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     484                  !                                         1 = sliding friction : TauB < RHS 
     485                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    500486                  ! 
    501487                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    502488                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    503                         &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     489                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    504490                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    505491                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    … …  
    508494                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    509495                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
    510                         &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     496                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    511497                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    512498                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    … …  
    544530                  ! 
    545531                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    546                   zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    547                   ! 
    548                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    549                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     532                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     533                  ! 
     534                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     535                  !                                         1 = sliding friction : TauB < RHS 
     536                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    550537                  ! 
    551538                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    552539                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    553                         &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     540                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    554541                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    555                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    556                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    557                         &            )   * zmsk00x(ji,jj) 
     542                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     543                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     544                        &           )   * zmsk00x(ji,jj) 
    558545                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    559546                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
    560                         &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     547                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    561548                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    562549                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    … …  
    596583                  ! 
    597584                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    598                   zTauE = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    599                   ! 
    600                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    601                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE + ztaux_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     585                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     586                  ! 
     587                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     588                  !                                         1 = sliding friction : TauB < RHS 
     589                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    602590                  ! 
    603591                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    604592                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    605                         &                                  + zTauE + zTauO * u_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     593                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    606594                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    607                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    608                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    609                         &            )   * zmsk00x(ji,jj) 
     595                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     596                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     597                        &           )   * zmsk00x(ji,jj) 
    610598                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    611599                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
    612                         &                                     + zTauE + zTauO * u_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     600                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    613601                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    614602                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    … …  
    646634                  ! 
    647635                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    648                   zTauE = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    649                   ! 
    650                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    651                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE + ztauy_base(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     636                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     637                  ! 
     638                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     639                  !                                         1 = sliding friction : TauB < RHS 
     640                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    652641                  ! 
    653642                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    654643                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    655                         &                                  + zTauE + zTauO * v_ice(ji,jj) )                                        & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     644                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    656645                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    657646                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    … …  
    660649                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    661650                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
    662                         &                                     + zTauE + zTauO * v_ice(ji,jj) )                                     & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     651                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    663652                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    664653                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceistate.F90

    r11413 r11586  
    101101      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays 
    102102      !! 
    103       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d 
     103      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 
    104104      !-------------------------------------------------------------------- 
    105105 
    … …  
    209209            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) 
    210210            ! 
    211             ! ponds 
     211            ! pond concentration 
    212212            IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) & 
    213                &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     213               &     si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc. 
     214               &                              * si(jp_ati)%fnow(:,:,1)  
    214215            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) 
     216            ! 
     217            ! pond depth 
    215218            IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 
    216219               &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    … …  
    238241               zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:) 
    239242               ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:) 
    240                zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) 
     243               zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.  
    241244               zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
    242245            ELSEWHERE 
    … …  
    248251               zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:) 
    249252               ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:) 
    250                zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) 
     253               zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
    251254               zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 
    252255            END WHERE 
    253256            ! 
    254257         ENDIF 
     258 
     259         ! make sure ponds = 0 if no ponds scheme 
     260         IF ( .NOT.ln_pnd ) THEN 
     261            zapnd_ini(:,:) = 0._wp 
     262            zhpnd_ini(:,:) = 0._wp 
     263         ENDIF 
     264          
    255265         !-------------! 
    256266         ! fill fields ! 
    … …  
    275285         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti)  , zt_su_ini ) 
    276286         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti)  , zsm_i_ini ) 
     287         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini ) 
     288         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini ) 
    277289 
    278290         ! allocate temporary arrays 
    279291         ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 
    280             &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl) ) 
     292            &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 
    281293          
    282294         ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 
    283          CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                   zhi_2d, zhs_2d, zai_2d , & 
    284             &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti),   zti_2d, zts_2d, ztsu_2d, zsi_2d ) 
     295         CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   & 
     296            &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   & 
     297            &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), & 
     298            &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d ) 
    285299 
    286300         ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 
    … …  
    289303            zts_3d(:,:,jl) = rt0 * tmask(:,:,1) 
    290304         END DO 
    291          CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d  , h_i  ) 
    292          CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d  , h_s  ) 
    293          CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d  , a_i  ) 
    294          CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d  , zti_3d ) 
    295          CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d  , zts_3d ) 
    296          CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su ) 
    297          CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d  , s_i  ) 
     305         CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d   , h_i    ) 
     306         CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d   , h_s    ) 
     307         CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d   , a_i    ) 
     308         CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d   , zti_3d ) 
     309         CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d   , zts_3d ) 
     310         CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d  , t_su   ) 
     311         CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d   , s_i    ) 
     312         CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   ) 
     313         CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   ) 
    298314 
    299315         ! deallocate temporary arrays 
    300316         DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 
    301             &        zti_2d, zts_2d, ztsu_2d, zsi_2d ) 
    302  
    303          ! Melt ponds: distribute uniformely over the categories 
    304          IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN 
    305             DO jl = 1, jpl 
    306                a_ip_frac(:,:,jl) = zapnd_ini(:,:) 
    307                h_ip     (:,:,jl) = zhpnd_ini(:,:) 
    308                a_ip     (:,:,jl) = a_ip_frac(:,:,jl) * a_i (:,:,jl)  
    309                v_ip     (:,:,jl) = h_ip     (:,:,jl) * a_ip(:,:,jl) 
    310             END DO 
    311          ENDIF 
    312            
     317            &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 
     318 
    313319         ! calculate extensive and intensive variables 
    314320         CALL ice_var_salprof ! for sz_i 
    … …  
    350356         END DO 
    351357 
     358         ! Melt ponds 
     359         WHERE( a_i > epsi10 ) 
     360            a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
     361         ELSEWHERE 
     362            a_ip_frac(:,:,:) = 0._wp 
     363         END WHERE 
     364         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
     365           
    352366         ! specific temperatures for coupled runs 
    353367         tn_ice(:,:,:) = t_su(:,:,:) 
    … …  
    512526      ENDIF 
    513527      ! 
     528      IF( .NOT.ln_pnd ) THEN 
     529         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0. 
     530         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0. 
     531         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' ) 
     532      ENDIF 
     533      ! 
    514534   END SUBROUTINE ice_istate_init 
    515535 

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/iceitd.F90

    r11348 r11586  
    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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd.F90

    r11348 r11586  
    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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_do.F90

    r11348 r11586  
    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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icethd_pnd.F90

    r11348 r11586  
    205205      INTEGER  ::   ios, ioptio   ! Local integer 
    206206      !! 
    207       NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
     207      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
    208208      !!------------------------------------------------------------------- 
    209209      ! 
    … …  
    221221         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    222222         WRITE(numout,*) '   Namelist namicethd_pnd:' 
    223          WRITE(numout,*) '      Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 
    224          WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST 
    225          WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
    226          WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd 
    227          WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb 
     223         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd 
     224         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 
     225         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST 
     226         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
     227         WRITE(numout,*) '            Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd 
     228         WRITE(numout,*) '         Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb 
    228229      ENDIF 
    229230      ! 
    230231      !                             !== set the choice of ice pond scheme ==! 
    231232      ioptio = 0 
    232                                                             nice_pnd = np_pndNO 
    233       IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
    234       IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF 
    235       IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' ) 
     233      IF( .NOT.ln_pnd ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndNO     ;   ENDIF 
     234      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
     235      IF( ln_pnd_H12  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF 
     236      IF( ioptio /= 1 )   & 
     237         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' ) 
    236238      ! 
    237239      SELECT CASE( nice_pnd ) 

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90

    r11413 r11586  
    4747   !!   ice_var_zapneg    : remove negative ice fields 
    4848   !!   ice_var_roundoff  : remove negative values arising from roundoff erros 
    49    !!   ice_var_itd       : convert N-cat to M-cat 
    5049   !!   ice_var_bv        : brine volume 
    5150   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature 
    5251   !!   ice_var_sshdyn    : compute equivalent ssh in lead 
     52   !!   ice_var_itd       : convert N-cat to M-cat 
    5353   !!---------------------------------------------------------------------- 
    5454   USE dom_oce        ! ocean space and time domain 
    … …  
    115115      ! 
    116116      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction   
    117       ! 
    118       ALLOCATE( z1_at_i(jpi,jpj) ) 
    119       WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:) 
    120       ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp 
    121       END WHERE 
    122       ! 
    123       tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    124       WHERE( at_i(:,:)<=epsi20 ); tm_su(:,:) = rt0; END WHERE       
    125       ! 
     117 
    126118      ! The following fields are calculated for diagnostics and outputs only 
    127119      ! ==> Do not use them for other purposes 
    128120      IF( kn > 1 ) THEN 
    129121         ! 
    130          ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 
     122         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) ) 
     123         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:) 
     124         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp 
     125         END WHERE 
    131126         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:) 
    132127         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp 
    … …  
    141136         !          
    142137         !                          ! mean temperature (K), salinity and age 
     138         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    143139         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    144140         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:) 
    … …  
    158154         !                           ! put rt0 where there is no ice 
    159155         WHERE( at_i(:,:)<=epsi20 ) 
     156            tm_su(:,:) = rt0 
    160157            tm_si(:,:) = rt0 
    161158            tm_i (:,:) = rt0 
    162159            tm_s (:,:) = rt0 
    163160         END WHERE 
    164  
    165          DEALLOCATE( z1_vt_i , z1_vt_s ) 
     161         ! 
     162         !                           ! mean melt pond depth 
     163         WHERE( at_ip(:,:) > epsi20 )   ;   hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 
     164         ELSEWHERE                      ;   hm_ip(:,:) = 0._wp 
     165         END WHERE          
     166         ! 
     167         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 
     168         ! 
    166169      ENDIF 
    167       ! 
    168       DEALLOCATE( z1_at_i ) 
    169170      ! 
    170171   END SUBROUTINE ice_var_agg 
    … …  
    664665      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    665666      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    666       IF ( ln_pnd_H12 ) THEN 
     667      IF( ln_pnd_H12 ) THEN 
    667668         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    668669         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
    … …  
    785786   !! ** Purpose :  converting N-cat ice to jpl ice categories 
    786787   !!------------------------------------------------------------------- 
    787    SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    788       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     788   SUBROUTINE ice_var_itd_1c1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     789      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    789790      !!------------------------------------------------------------------- 
    790791      !! ** Purpose :  converting 1-cat ice to 1 ice category 
    … …  
    792793      REAL(wp), DIMENSION(:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    793794      REAL(wp), DIMENSION(:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    794       REAL(wp), DIMENSION(:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    795       REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     795      REAL(wp), DIMENSION(:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     796      REAL(wp), DIMENSION(:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    796797      !!------------------------------------------------------------------- 
    797798      ! == thickness and concentration == ! 
    … …  
    800801      pa_i(:) = pati(:) 
    801802      ! 
    802       ! == temperature and salinity == ! 
    803       IF( PRESENT( pt_i  ) )   pt_i (:) = ptmi (:) 
    804       IF( PRESENT( pt_s  ) )   pt_s (:) = ptms (:) 
    805       IF( PRESENT( pt_su ) )   pt_su(:) = ptmsu(:) 
    806       IF( PRESENT( ps_i  ) )   ps_i (:) = psmi (:) 
     803      ! == temperature and salinity and ponds == ! 
     804      pt_i (:) = ptmi (:) 
     805      pt_s (:) = ptms (:) 
     806      pt_su(:) = ptmsu(:) 
     807      ps_i (:) = psmi (:) 
     808      pa_ip(:) = patip(:) 
     809      ph_ip(:) = phtip(:) 
    807810       
    808811   END SUBROUTINE ice_var_itd_1c1c 
    809812 
    810    SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    811       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     813   SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     814      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    812815      !!------------------------------------------------------------------- 
    813816      !! ** Purpose :  converting N-cat ice to 1 ice category 
    … …  
    815818      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    816819      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    817       REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    818       REAL(wp), DIMENSION(:)  , INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     820      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     821      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    819822      ! 
    820823      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z1_ai, z1_vi, z1_vs 
    … …  
    826829      ! 
    827830      ! == thickness and concentration == ! 
    828       ALLOCATE( z1_ai(idim) ) 
     831      ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) ) 
    829832      ! 
    830833      pa_i(:) = SUM( pati(:,:), dim=2 ) 
    … …  
    838841      ! 
    839842      ! == temperature and salinity == ! 
    840       IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 
    841          ! 
    842          ALLOCATE( z1_vi(idim), z1_vs(idim) ) 
    843          ! 
    844          WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 
    845          ELSEWHERE                                 ;   z1_vi(:) = 0._wp 
    846          END WHERE 
    847          WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 
    848          ELSEWHERE                                 ;   z1_vs(:) = 0._wp 
    849          END WHERE 
    850          ! 
    851          IF( PRESENT( pt_i  ) )   pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    852          IF( PRESENT( pt_s  ) )   pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
    853          IF( PRESENT( pt_su ) )   pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
    854          IF( PRESENT( ps_i  ) )   ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    855          ! 
    856          DEALLOCATE( z1_vi, z1_vs ) 
    857          ! 
    858       ENDIF 
    859       ! 
    860       DEALLOCATE( z1_ai ) 
     843      WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp )   ;   z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 
     844      ELSEWHERE                                 ;   z1_vi(:) = 0._wp 
     845      END WHERE 
     846      WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp )   ;   z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 
     847      ELSEWHERE                                 ;   z1_vs(:) = 0._wp 
     848      END WHERE 
     849      pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     850      pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
     851      pt_su(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
     852      ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     853 
     854      ! == ponds == ! 
     855      pa_ip(:) = SUM( patip(:,:), dim=2 ) 
     856      WHERE( pa_ip(:) /= 0._wp )   ;   ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 
     857      ELSEWHERE                    ;   ph_ip(:) = 0._wp 
     858      END WHERE 
     859      ! 
     860      DEALLOCATE( z1_ai, z1_vi, z1_vs ) 
    861861      ! 
    862862   END SUBROUTINE ice_var_itd_Nc1c 
    863863    
    864    SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    865       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     864   SUBROUTINE ice_var_itd_1cMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     865      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    866866      !!------------------------------------------------------------------- 
    867867      !! 
    … …  
    894894      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    895895      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    896       REAL(wp), DIMENSION(:)  , INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    897       REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     896      REAL(wp), DIMENSION(:)  , INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     897      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    898898      ! 
    899899      INTEGER , DIMENSION(4) ::   itest 
     900      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra 
    900901      INTEGER  ::   ji, jk, jl 
    901902      INTEGER  ::   idim, i_fill, jl0   
    … …  
    10101011      ! 
    10111012      ! == temperature and salinity == ! 
    1012       IF( PRESENT( pt_i  ) ) THEN 
    1013          DO jl = 1, jpl 
    1014             pt_i(:,jl) = ptmi (:) 
    1015          END DO 
    1016       ENDIF 
    1017       IF( PRESENT( pt_s  ) ) THEN 
    1018          DO jl = 1, jpl 
    1019             pt_s (:,jl) = ptms (:) 
    1020          END DO 
    1021       ENDIF 
    1022       IF( PRESENT( pt_su ) ) THEN 
    1023          DO jl = 1, jpl 
    1024             pt_su(:,jl) = ptmsu(:) 
    1025          END DO 
    1026       ENDIF 
    1027       IF( PRESENT( ps_i  ) ) THEN 
    1028          DO jl = 1, jpl 
    1029             ps_i (:,jl) = psmi (:) 
    1030          END DO 
    1031       ENDIF 
     1013      DO jl = 1, jpl 
     1014         pt_i (:,jl) = ptmi (:) 
     1015         pt_s (:,jl) = ptms (:) 
     1016         pt_su(:,jl) = ptmsu(:) 
     1017         ps_i (:,jl) = psmi (:) 
     1018         ps_i (:,jl) = psmi (:)          
     1019      END DO 
     1020      ! 
     1021      ! == ponds == ! 
     1022      ALLOCATE( zfra(idim) ) 
     1023      ! keep the same pond fraction atip/ati for each category 
     1024      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:) 
     1025      ELSEWHERE                   ;   zfra(:) = 0._wp 
     1026      END WHERE 
     1027      DO jl = 1, jpl 
     1028         pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     1029      END DO 
     1030      ! keep the same v_ip/v_i ratio for each category 
     1031      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) ) 
     1032      ELSEWHERE                                 ;   zfra(:) = 0._wp 
     1033      END WHERE 
     1034      DO jl = 1, jpl 
     1035         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1036         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
     1037         END WHERE 
     1038      END DO 
     1039      DEALLOCATE( zfra ) 
    10321040      ! 
    10331041   END SUBROUTINE ice_var_itd_1cMc 
    10341042 
    1035    SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,       ph_i, ph_s, pa_i, & 
    1036       &                         ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 
     1043   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                       ph_i, ph_s, pa_i, & 
     1044      &                         ptmi, ptms, ptmsu, psmi, patip, phtip,   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 
    10371045      !!------------------------------------------------------------------- 
    10381046      !! 
    … …  
    10651073      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
    10661074      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables 
    1067       REAL(wp), DIMENSION(:,:), INTENT(in)   , OPTIONAL ::   ptmi, ptms, ptmsu, psmi    ! input  ice/snow temp & sal 
    1068       REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL ::   pt_i, pt_s, pt_su, ps_i    ! output ice/snow temp & sal 
     1075      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip    ! input  ice/snow temp & sal & ponds 
     1076      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    10691077      ! 
    10701078      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2 
    10711079      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin 
    1072       REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp 
     1080      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra 
    10731081      ! 
    10741082      REAL(wp), PARAMETER ::   ztrans = 0.25_wp 
    … …  
    10881096         pa_i(:,:) = pati(:,:) 
    10891097         ! 
    1090          ! == temperature and salinity == ! 
    1091          IF( PRESENT( pt_i  ) )   pt_i (:,:) = ptmi (:,:) 
    1092          IF( PRESENT( pt_s  ) )   pt_s (:,:) = ptms (:,:) 
    1093          IF( PRESENT( pt_su ) )   pt_su(:,:) = ptmsu(:,:) 
    1094          IF( PRESENT( ps_i  ) )   ps_i (:,:) = psmi (:,:) 
     1098         ! == temperature and salinity and ponds == ! 
     1099         pt_i (:,:) = ptmi (:,:) 
     1100         pt_s (:,:) = ptms (:,:) 
     1101         pt_su(:,:) = ptmsu(:,:) 
     1102         ps_i (:,:) = psmi (:,:) 
     1103         pa_ip(:,:) = patip(:,:) 
     1104         ph_ip(:,:) = phtip(:,:) 
    10951105         !                              ! ---------------------- ! 
    10961106      ELSEIF( icat == 1 ) THEN          ! input cat = 1          ! 
    10971107         !                              ! ---------------------- ! 
    1098          CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:) ) 
    1099 !!         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1),            ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
    1100 !!            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) ) 
     1108         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 
     1109            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), & 
     1110            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 
     1111            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:)  ) 
    11011112         !                              ! ---------------------- ! 
    11021113      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         ! 
    11031114         !                              ! ---------------------- ! 
    1104          CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1) ) 
    1105 !!         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:),            ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
    1106 !!            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )          
     1115         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 
     1116            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), & 
     1117            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 
     1118            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1)  ) 
    11071119         !                              ! ----------------------- ! 
    11081120      ELSE                              ! input cat /= output cat ! 
    … …  
    11941206         ! == temperature and salinity == ! 
    11951207         ! 
    1196          IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 
    1197             ! 
    1198             ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 
    1199             ! 
    1200             WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 
    1201             ELSEWHERE                                               ;   z1_ai(:) = 0._wp 
     1208         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 
     1209         ! 
     1210         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 
     1211         ELSEWHERE                                               ;   z1_ai(:) = 0._wp 
     1212         END WHERE 
     1213         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 
     1214         ELSEWHERE                                               ;   z1_vi(:) = 0._wp 
     1215         END WHERE 
     1216         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 
     1217         ELSEWHERE                                               ;   z1_vs(:) = 0._wp 
     1218         END WHERE 
     1219         ! 
     1220         ! fill all the categories with the same value 
     1221         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1222         DO jl = 1, jpl 
     1223            pt_i (:,jl) = ztmp(:) 
     1224         END DO 
     1225         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
     1226         DO jl = 1, jpl 
     1227            pt_s (:,jl) = ztmp(:) 
     1228         END DO 
     1229         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
     1230         DO jl = 1, jpl 
     1231            pt_su(:,jl) = ztmp(:) 
     1232         END DO 
     1233         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
     1234         DO jl = 1, jpl 
     1235            ps_i (:,jl) = ztmp(:) 
     1236         END DO 
     1237         ! 
     1238         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 
     1239         ! 
     1240         ! == ponds == ! 
     1241         ALLOCATE( zfra(idim) ) 
     1242         ! keep the same pond fraction atip/ati for each category 
     1243         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 ) 
     1244         ELSEWHERE                                   ;   zfra(:) = 0._wp 
     1245         END WHERE 
     1246         DO jl = 1, jpl 
     1247            pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 
     1248         END DO 
     1249         ! keep the same v_ip/v_i ratio for each category 
     1250         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 
     1251            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 
     1252         ELSEWHERE 
     1253            zfra(:) = 0._wp 
     1254         END WHERE 
     1255         DO jl = 1, jpl 
     1256            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 
     1257            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp 
    12021258            END WHERE 
    1203             WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 
    1204             ELSEWHERE                                               ;   z1_vi(:) = 0._wp 
    1205             END WHERE 
    1206             WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 
    1207             ELSEWHERE                                               ;   z1_vs(:) = 0._wp 
    1208             END WHERE 
    1209             ! 
    1210             ! fill all the categories with the same value 
    1211             IF( PRESENT( pt_i  ) ) THEN 
    1212                ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    1213                DO jl = 1, jpl 
    1214                   pt_i (:,jl) = ztmp(:) 
    1215                END DO 
    1216             ENDIF 
    1217             IF( PRESENT( pt_s  ) ) THEN 
    1218                ztmp(:) =  SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 
    1219                DO jl = 1, jpl 
    1220                   pt_s (:,jl) = ztmp(:) 
    1221                END DO 
    1222             ENDIF 
    1223             IF( PRESENT( pt_su ) ) THEN 
    1224                ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:) 
    1225                DO jl = 1, jpl 
    1226                   pt_su(:,jl) = ztmp(:) 
    1227                END DO 
    1228             ENDIF 
    1229             IF( PRESENT( ps_i  ) ) THEN 
    1230                ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 
    1231                DO jl = 1, jpl 
    1232                   ps_i (:,jl) = ztmp(:) 
    1233                END DO 
    1234             ENDIF 
    1235             ! 
    1236             DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 
    1237             ! 
    1238          ENDIF 
     1259         END DO 
     1260         DEALLOCATE( zfra ) 
    12391261         ! 
    12401262      ENDIF 

  • NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icewri.F90

    r11413 r11586  
    114114      ! melt ponds 
    115115      IF( iom_use('iceapnd' ) )   CALL iom_put( 'iceapnd', at_ip  * zmsk00      )                                           ! melt pond total fraction 
     116      IF( iom_use('icehpnd' ) )   CALL iom_put( 'icehpnd', hm_ip  * zmsk00      )                                           ! melt pond depth 
    116117      IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area 
    117118      ! salt 
Note: See TracChangeset for help on using the changeset viewer.