Changeset 14026


Ignore:
Timestamp:
2020-12-03T09:48:10+01:00 (10 months ago)
Author:
clem
Message:

4.0-HEAD: fix bugs and defects related to tickets #2573 #2575 #2576 #2578. Sette passed and those fixes are now in the trunk, so unless there is a tricky trick somewhere, everything should be fine.

Location:
NEMO/releases/r4.0/r4.0-HEAD/src/ICE
Files:
15 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/ice.F90

    r13640 r14026  
    366366   !!---------------------------------------------------------------------- 
    367367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 
     368   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_ip_b, v_il_b             !: ponds and lids volumes 
    368369   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b              !: 
    369370   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                      !: snow heat content 
     
    392393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s]  
    393394   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_aice         !: ice conc.  variation   [s-1]  
     395   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vpnd         !: pond volume variation  [m/s]  
    394396   ! 
    395397   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_mass     !: advection of mass (kg/m2/s) 
     
    468470         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) ,  & 
    469471         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) ,  & 
    470          &      om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj)            , STAT=ierr(ii) ) 
     472         &      om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj), STAT=ierr(ii) ) 
    471473 
    472474      ii = ii + 1 
     
    486488      ii = ii + 1 
    487489      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl),         & 
     490         &      v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) ,                                                         & 
    488491         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 
    489492         &      STAT=ierr(ii) ) 
     
    500503      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),                      &  
    501504         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),                      & 
    502          &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj), & 
     505         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), & 
    503506         &      diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 
    504507 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icectl.F90

    r13589 r14026  
    8585      !! 
    8686      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, & 
    87          &          zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 
     87         &          zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, & 
     88         &          zdiag_eimin, zdiag_esmin, zdiag_simin 
    8889      REAL(wp) ::   zvtrp, zetrp 
    8990      REAL(wp) ::   zarea 
     
    9293      IF( icount == 0 ) THEN 
    9394 
    94          pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) 
    95          pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) 
     95         pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 
     96         pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi , dim=3 ) * e1e2t ) 
    9697         pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) 
    9798 
     
    112113 
    113114         ! -- mass diag -- ! 
    114          zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_rdtice       & 
     115            zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t )   & 
     116            &            - pdiag_v ) * r1_rdtice                                                                          & 
    115117            &         + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn +       & 
    116118            &                                 wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & 
     
    132134 
    133135         ! -- min/max diag -- ! 
    134          zdiag_amax  = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
    135          zdiag_vmin  = glob_min( 'icectl', v_i ) 
    136          zdiag_amin  = glob_min( 'icectl', a_i ) 
    137          zdiag_smin  = glob_min( 'icectl', sv_i ) 
     136         zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 
     137         zdiag_vimin = glob_min( 'icectl', v_i  ) 
     138         zdiag_vsmin = glob_min( 'icectl', v_s  ) 
     139         zdiag_vpmin = glob_min( 'icectl', v_ip ) 
     140         zdiag_vlmin = glob_min( 'icectl', v_il ) 
     141         zdiag_aimin = glob_min( 'icectl', a_i  ) 
     142         zdiag_simin = glob_min( 'icectl', sv_i ) 
    138143         zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 
    139144         zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) 
     
    155160               &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice 
    156161            ! check negative values 
    157             IF( zdiag_vmin  < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i < 0         = ',zdiag_vmin 
    158             IF( zdiag_amin  < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i < 0         = ',zdiag_amin 
    159             IF( zdiag_smin  < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i < 0         = ',zdiag_smin 
    160             IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i < 0         = ',zdiag_eimin 
    161             IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s < 0         = ',zdiag_esmin 
     162            IF( zdiag_vimin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i  < 0        = ',zdiag_vimin 
     163            IF( zdiag_vsmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_s  < 0        = ',zdiag_vsmin 
     164            IF( zdiag_vpmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_ip < 0        = ',zdiag_vpmin 
     165            IF( zdiag_vlmin < 0. ) WRITE(numout,*)   cd_routine,' : violation v_il < 0        = ',zdiag_vlmin 
     166            IF( zdiag_aimin < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i  < 0        = ',zdiag_aimin 
     167            IF( zdiag_simin < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i  < 0        = ',zdiag_simin 
     168            IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i  < 0        = ',zdiag_eimin 
     169            IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s  < 0        = ',zdiag_esmin 
    162170            ! check maximum ice concentration 
    163             IF( zdiag_amax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
    164                &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax 
     171            IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 
     172               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_aimax 
    165173            ! check if advection scheme is conservative 
    166174            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
     
    193201      ! water flux 
    194202      ! -- mass diag -- ! 
    195       zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
    196          &                              + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) 
     203      zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr   + wfx_sub + wfx_pnd  & 
     204         &                              + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) 
    197205 
    198206      ! -- salt diag -- ! 
     
    200208 
    201209      ! -- heat diag -- ! 
    202       zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
     210      zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 
    203211      ! equivalent to this: 
    204212      !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & 
     
    245253      IF( icount == 0 ) THEN 
    246254 
    247          pdiag_v = SUM( v_i  * rhoi + v_s * rhos, dim=3 ) 
    248          pdiag_s = SUM( sv_i * rhoi             , dim=3 ) 
     255         pdiag_v = SUM( v_i  * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) 
     256         pdiag_s = SUM( sv_i * rhoi , dim=3 ) 
    249257         pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 
    250258 
     
    261269 
    262270         ! -- mass diag -- ! 
    263          zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos, dim=3 ) - pdiag_v ) * r1_rdtice                             & 
     271         zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_rdtice    & 
    264272            &         + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 
    265273            &             wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr )           & 
     
    352360      CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 )    !  
    353361      CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 )    !  
     362      ! mean state 
     363      CALL iom_rstput( 0, 0, inum, 'icecon'    , SUM(a_i ,dim=3) , ktype = jp_r8 )    ! 
     364      CALL iom_rstput( 0, 0, inum, 'icevol'    , SUM(v_i ,dim=3) , ktype = jp_r8 )    ! 
     365      CALL iom_rstput( 0, 0, inum, 'snwvol'    , SUM(v_s ,dim=3) , ktype = jp_r8 )    ! 
     366      CALL iom_rstput( 0, 0, inum, 'pndvol'    , SUM(v_ip,dim=3) , ktype = jp_r8 )    ! 
     367      CALL iom_rstput( 0, 0, inum, 'lidvol'    , SUM(v_il,dim=3) , ktype = jp_r8 )    ! 
    354368       
    355369      CALL iom_close( inum ) 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn.F90

    r13284 r14026  
    8080      ! 
    8181      ! controls 
    82       IF( ln_timing )   CALL timing_start('icedyn') 
     82      IF( ln_timing )   CALL timing_start('ice_dyn') 
    8383      ! 
    8484      IF( kt == nit000 .AND. lwp ) THEN 
     
    175175      ! 
    176176      ! controls 
    177       IF( ln_timing )   CALL timing_stop ('icedyn') 
     177      IF( ln_timing )   CALL timing_stop ('ice_dyn') 
    178178      ! 
    179179   END SUBROUTINE ice_dyn 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_pra.F90

    r13634 r14026  
    156156 
    157157         ! diagnostics 
    158          zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     158         zdiag_adv_mass(:,:) =   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     159            &                  + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 
    159160         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
    160161         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     
    322323         ! 
    323324         ! --- diagnostics --- ! 
    324          diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     325         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     326            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 
    325327            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
    326328         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_umx.F90

    r13617 r14026  
    183183 
    184184         ! diagnostics 
    185          zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     185         zdiag_adv_mass(:,:) =   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     186            &                  + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 
    186187         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
    187188         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
     
    384385         ! 
    385386         ! --- diagnostics --- ! 
    386          diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     387         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     388            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 
    387389            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
    388390         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceistate.F90

    r13284 r14026  
    403403      ! 4) Snow-ice mass (case ice is fully embedded) 
    404404      !---------------------------------------------- 
    405       snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass 
     405      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3  )   ! snow+ice mass 
    406406      snwice_mass_b(:,:) = snwice_mass(:,:) 
    407407      ! 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceitd.F90

    r13617 r14026  
    2929   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    3030   USE prtctl         ! Print control 
     31   USE timing         ! Timing 
    3132 
    3233   IMPLICIT NONE 
     
    8586      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories 
    8687      !!------------------------------------------------------------------ 
     88      IF( ln_timing )   CALL timing_start('iceitd_rem') 
    8789 
    8890      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution'  
     
    318320      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    319321      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     322      IF( ln_timing    )   CALL timing_stop ('iceitd_rem') 
    320323      ! 
    321324   END SUBROUTINE ice_itd_rem 
     
    481484                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 
    482485                  !                                               
    483                   ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a) 
     486!!$                  ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a) 
     487                  ztrans          = v_ip_2d(ji,jl1) * zworkv(ji)     ! Pond volume 
    484488                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    485489                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
    486490                  ! 
    487491                  IF ( ln_pnd_lids ) THEN                            ! Pond lid volume 
    488                      ztrans          = v_il_2d(ji,jl1) * zworka(ji) 
     492!!$                     ztrans          = v_il_2d(ji,jl1) * zworka(ji) 
     493                     ztrans          = v_il_2d(ji,jl1) * zworkv(ji) 
    489494                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 
    490495                     v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans 
     
    592597      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred 
    593598      !!------------------------------------------------------------------ 
     599      IF( ln_timing )   CALL timing_start('iceitd_reb') 
    594600      ! 
    595601      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
     
    623629               jdonor(ji,jl)  = jl  
    624630               ! how much of a_i you send in cat sup is somewhat arbitrary 
    625                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    626                !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    627                !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    628                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    629                !!          zdaice(ji,jl)  = a_i_1d(ji) 
    630                !!          zdvice(ji,jl)  = v_i_1d(ji) 
    631                !!clem: these are from UCL and work ok 
    632                zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    633                zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     631               ! these are from CICE => transfer everything 
     632               !!zdaice(ji,jl)  = a_i_1d(ji) 
     633               !!zdvice(ji,jl)  = v_i_1d(ji) 
     634               ! these are from LLN => transfer only half of the category 
     635               zdaice(ji,jl)  =                       0.5_wp  * a_i_1d(ji) 
     636               zdvice(ji,jl)  = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl) 
    634637            END DO 
    635638            ! 
     
    677680      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    678681      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     682      IF( ln_timing    )   CALL timing_stop ('iceitd_reb') 
    679683      ! 
    680684   END SUBROUTINE ice_itd_reb 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icesbc.F90

    r13284 r14026  
    6161      !!------------------------------------------------------------------- 
    6262      ! 
    63       IF( ln_timing )   CALL timing_start('ice_sbc') 
     63      IF( ln_timing )   CALL timing_start('icesbc') 
    6464      ! 
    6565      IF( kt == nit000 .AND. lwp ) THEN 
     
    8686      ENDIF 
    8787      ! 
    88       IF( ln_timing )   CALL timing_stop('ice_sbc') 
     88      IF( ln_timing )   CALL timing_stop('icesbc') 
    8989      ! 
    9090   END SUBROUTINE ice_sbc_tau 
     
    119119      !!-------------------------------------------------------------------- 
    120120      ! 
    121       IF( ln_timing )   CALL timing_start('ice_sbc_flx') 
     121      IF( ln_timing )   CALL timing_start('icesbc') 
    122122 
    123123      IF( kt == nit000 .AND. lwp ) THEN 
     
    172172      ENDIF 
    173173      ! 
    174       IF( ln_timing )   CALL timing_stop('ice_sbc_flx') 
     174      IF( ln_timing )   CALL timing_stop('icesbc') 
    175175      ! 
    176176   END SUBROUTINE ice_sbc_flx 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icestp.F90

    r13720 r14026  
    120120      !!---------------------------------------------------------------------- 
    121121      ! 
    122       IF( ln_timing )   CALL timing_start('ice_stp') 
     122      IF( ln_timing )   CALL timing_start('icestp') 
    123123      ! 
    124124      !                                      !-----------------------! 
     
    214214!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    215215      ! 
    216       IF( ln_timing )   CALL timing_stop('ice_stp') 
     216      IF( ln_timing )   CALL timing_stop('icestp') 
    217217      ! 
    218218   END SUBROUTINE ice_stp 
     
    365365      v_i_b (:,:,:)   = v_i (:,:,:)     ! ice volume 
    366366      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume 
     367      v_ip_b(:,:,:)   = v_ip(:,:,:)     ! pond volume 
     368      v_il_b(:,:,:)   = v_il(:,:,:)     ! pond lid volume 
    367369      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content 
    368370      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy 
     
    425427            diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp 
    426428            diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp 
    427             diag_aice(ji,jj) = 0._wp 
     429            diag_aice(ji,jj) = 0._wp ;   diag_vpnd(ji,jj) = 0._wp 
    428430 
    429431            tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
     
    482484         diag_vsnw(:,:) = diag_vsnw(:,:) & 
    483485            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice * rhos 
     486         diag_vpnd(:,:) = diag_vpnd(:,:) & 
     487            &             + SUM(     v_ip + v_il          - v_ip_b - v_il_b                , dim=3 ) * r1_rdtice * rhow 
    484488         ! 
    485489         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd.F90

    r13642 r14026  
    171171            qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
    172172 
     173            ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 
     174            ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 
     175            IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 
     176               zqfr               = 0._wp 
     177               zqfr_pos           = 0._wp 
     178               qsb_ice_bot(ji,jj) = 0._wp 
     179            ENDIF 
     180            ! 
    173181            ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
    174182            !     qlead is the energy received from the atm. in the leads. 
     
    247255            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
    248256                              CALL ice_thd_dh                           ! Ice-Snow thickness    
    249                               CALL ice_thd_pnd                          ! Melt ponds formation 
    250257                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    251258            ENDIF 
     
    258265            ! 
    259266            IF( ln_icedA )    CALL ice_thd_da                       ! --- Lateral melting --- ! 
     267            ! 
     268            IF( ln_pnd .AND. ln_icedH ) & 
     269               &              CALL ice_thd_pnd                      ! --- Melt ponds formation --- ! 
    260270            ! 
    261271                              CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_dh.F90

    r13642 r14026  
    5555      !!                - Snow ice formation 
    5656      !! 
     57      !! ** Note     :  h=max(0,h+dh) are often used to ensure positivity of h. 
     58      !!                very small negative values can occur otherwise (e.g. -1.e-20) 
     59      !! 
    5760      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 
    5861      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646    
     
    7982      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean 
    8083 
    81       REAL(wp), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3) 
    8284      REAL(wp), DIMENSION(jpij) ::   zq_top      ! heat for surface ablation                   (J.m-2) 
    8385      REAL(wp), DIMENSION(jpij) ::   zq_bot      ! heat for bottom ablation                    (J.m-2) 
     
    8587      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    8688      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2) 
    87  
    88       REAL(wp), DIMENSION(jpij) ::   zdh_s_mel   ! snow melt  
    89       REAL(wp), DIMENSION(jpij) ::   zdh_s_pre   ! snow precipitation  
    90       REAL(wp), DIMENSION(jpij) ::   zdh_s_sub   ! snow sublimation 
    91  
    92       REAL(wp), DIMENSION(jpij,nlay_s) ::   zh_s      ! snw layer thickness 
    93       REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness 
    94       REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah 
    95       INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting  
    96  
     89      REAL(wp), DIMENSION(jpij) ::   zdeltah      
    9790      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing 
     91 
     92      INTEGER , DIMENSION(jpij,nlay_i)     ::   icount    ! number of layers vanishing by melting  
     93      REAL(wp), DIMENSION(jpij,0:nlay_i+1) ::   zh_i      ! ice layer thickness (m) 
     94      REAL(wp), DIMENSION(jpij,0:nlay_s  ) ::   zh_s      ! snw layer thickness (m) 
     95      REAL(wp), DIMENSION(jpij,0:nlay_s  ) ::   ze_s      ! snw layer enthalpy (J.m-3) 
    9896 
    9997      REAL(wp) ::   zswitch_sal 
     
    108106      END SELECT 
    109107 
    110       ! initialize layer thicknesses and enthalpies 
     108      ! initialize ice layer thicknesses and enthalpies 
     109      eh_i_old(1:npti,0:nlay_i+1) = 0._wp 
    111110      h_i_old (1:npti,0:nlay_i+1) = 0._wp 
    112       eh_i_old(1:npti,0:nlay_i+1) = 0._wp 
     111      zh_i    (1:npti,0:nlay_i+1) = 0._wp 
    113112      DO jk = 1, nlay_i 
    114113         DO ji = 1, npti 
     114            eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk) 
    115115            h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i 
    116             eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) 
     116            zh_i    (ji,jk) = h_i_1d(ji) * r1_nlay_i 
     117         END DO 
     118      END DO 
     119      ! 
     120      ! initialize snw layer thicknesses and enthalpies 
     121      zh_s(1:npti,0) = 0._wp 
     122      ze_s(1:npti,0) = 0._wp 
     123      DO jk = 1, nlay_s 
     124         DO ji = 1, npti 
     125            zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 
     126            ze_s(ji,jk) = e_s_1d(ji,jk) 
    117127         END DO 
    118128      END DO 
     
    141151         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji)  
    142152         zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    143       END DO 
    144  
    145       ! Ice and snow layer thicknesses 
    146       !------------------------------- 
    147       DO jk = 1, nlay_i 
    148          DO ji = 1, npti 
    149             zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i 
    150          END DO 
    151       END DO 
    152  
    153       DO jk = 1, nlay_s 
    154          DO ji = 1, npti 
    155             zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 
    156          END DO 
    157153      END DO 
    158154 
     
    167163         DO ji = 1, npti 
    168164            IF( t_s_1d(ji,jk) > rt0 ) THEN 
    169                hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! heat flux to the ocean [W.m-2], < 0 
    170                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux 
     165               hfx_res_1d    (ji) = hfx_res_1d    (ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! heat flux to the ocean [W.m-2], < 0 
     166               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos        * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux 
    171167               ! updates 
    172                dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk) 
    173                h_s_1d  (ji)    = h_s_1d(ji) - zh_s(ji,jk) 
     168               dh_s_mlt(ji)    =             dh_s_mlt(ji) - zh_s(ji,jk) 
     169               h_s_1d  (ji)    = MAX( 0._wp, h_s_1d  (ji) - zh_s(ji,jk) ) 
    174170               zh_s    (ji,jk) = 0._wp 
    175                e_s_1d  (ji,jk) = 0._wp 
    176                t_s_1d  (ji,jk) = rt0 
     171               ze_s    (ji,jk) = 0._wp 
    177172            END IF 
    178173         END DO 
     
    181176      ! Snow precipitation 
    182177      !------------------- 
    183       CALL ice_var_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing 
    184  
    185       zdeltah(1:npti,:) = 0._wp 
     178      CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing 
     179 
    186180      DO ji = 1, npti 
    187181         IF( sprecip_1d(ji) > 0._wp ) THEN 
     182            zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)   ! thickness of precip 
     183            ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) )                              ! enthalpy of the precip (>0, J.m-3) 
    188184            ! 
    189             ! --- precipitation --- 
    190             zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)   ! thickness change 
    191             zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3) 
     185            hfx_spr_1d(ji) = hfx_spr_1d(ji) + ze_s(ji,0) * zh_s(ji,0) * a_i_1d(ji) * r1_rdtice   ! heat flux from snow precip (>0, W.m-2) 
     186            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos       * zh_s(ji,0) * a_i_1d(ji) * r1_rdtice   ! mass flux, <0 
    192187            ! 
    193             hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat flux from snow precip (>0, W.m-2) 
    194             wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0 
    195              
    196             ! --- melt of falling snow --- 
    197             rswitch              = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
    198             zdeltah       (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change 
    199             zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                 ! bound melting  
    200             hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat used to melt snow (W.m-2, >0) 
    201             wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip), >0 
    202              
    203             ! updates available heat + precipitations after melting 
    204             dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1) 
    205             zq_top   (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )       
    206             zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
    207              
    208188            ! update thickness 
    209             h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 
    210             ! 
    211          ELSE 
    212             ! 
    213             zdh_s_pre(ji) = 0._wp 
    214             zqprec   (ji) = 0._wp 
    215             ! 
     189            h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0) 
    216190         ENDIF 
    217       END DO 
    218  
    219       ! recalculate snow layers 
    220       DO jk = 1, nlay_s 
    221          DO ji = 1, npti 
    222             zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 
    223          END DO 
    224191      END DO 
    225192 
    226193      ! Snow melting 
    227194      ! ------------ 
    228       ! If heat still available (zq_top > 0), then melt more snow 
    229       zdeltah(1:npti,:) = 0._wp 
    230       zdh_s_mel(1:npti) = 0._wp 
    231       DO jk = 1, nlay_s 
     195      ! If heat still available (zq_top > 0) 
     196      ! then all snw precip has been melted and we need to melt more snow 
     197      DO jk = 0, nlay_s 
    232198         DO ji = 1, npti 
    233199            IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN 
    234200               ! 
    235                rswitch          = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 
    236                zdeltah  (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change 
    237                zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting 
    238                zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
    239                 
    240                hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0) 
    241                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
     201               rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) ) 
     202               zdum    = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 )   ! thickness change 
     203               zdum    = MAX( zdum , - zh_s(ji,jk) )                           ! bound melting 
     204                
     205               hfx_snw_1d    (ji) = hfx_snw_1d    (ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_rdtice   ! heat used to melt snow(W.m-2, >0) 
     206               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos        * zdum * a_i_1d(ji) * r1_rdtice   ! snow melting only = water into the ocean 
    242207                
    243208               ! updates available heat + thickness 
    244                dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk) 
    245                zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    246                h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
    247                zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 
     209               dh_s_mlt(ji)    =              dh_s_mlt(ji)    + zdum 
     210               zq_top  (ji)    = MAX( 0._wp , zq_top  (ji)    + zdum * ze_s(ji,jk) ) 
     211               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d  (ji)    + zdum ) 
     212               zh_s    (ji,jk) = MAX( 0._wp , zh_s    (ji,jk) + zdum ) 
     213!!$               IF( zh_s(ji,jk) == 0._wp )   ze_s(ji,jk) = 0._wp 
    248214               ! 
    249215            ENDIF 
     
    255221      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    256222      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 
    257       zdeltah(1:npti,:) = 0._wp 
     223      zdeltah   (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 
     224      zevap_rema(1:npti) = 0._wp 
    258225      DO ji = 1, npti 
    259226         IF( evap_ice_1d(ji) > 0._wp ) THEN 
     227            zdeltah   (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rdt_ice, - h_s_1d(ji) )   ! amount of snw that sublimates, < 0             
     228            zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rdt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on) 
     229         ENDIF 
     230      END DO 
     231       
     232      DO jk = 0, nlay_s 
     233         DO ji = 1, npti 
     234            zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0 
    260235            ! 
    261             zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice ) 
    262             zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos   ! remaining evap in kg.m-2 (used for ice melting later on) 
    263             zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    264              
    265             hfx_sub_1d    (ji) = hfx_sub_1d(ji) + &   ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 
    266                &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  & 
    267                &                 * a_i_1d(ji) * r1_rdtice 
    268             wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice   ! Mass flux by sublimation 
    269              
    270             ! new snow thickness 
    271             h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 
    272             ! update precipitations after sublimation and correct sublimation 
    273             zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
    274             zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
    275             ! 
    276          ELSE 
    277             ! 
    278             zdh_s_sub (ji) = 0._wp 
    279             zevap_rema(ji) = 0._wp 
    280             ! 
    281          ENDIF 
    282       END DO 
    283        
    284       ! --- Update snow diags --- ! 
    285       DO ji = 1, npti 
    286          dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 
    287       END DO 
    288  
    289       ! Update temperature, energy 
    290       !--------------------------- 
    291       ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 
    292       DO jk = 1, nlay_s 
    293          DO ji = 1,npti 
    294             rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 
    295             e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            & 
    296               &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  & 
    297               &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 
    298          END DO 
    299       END DO 
    300        
     236            hfx_sub_1d    (ji) = hfx_sub_1d    (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_rdtice  ! Heat flux of snw that sublimates [W.m-2], < 0 
     237            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos        * zdum * a_i_1d(ji) * r1_rdtice  ! Mass flux by sublimation 
     238 
     239            ! update thickness 
     240            h_s_1d(ji)    = MAX( 0._wp , h_s_1d(ji)    + zdum ) 
     241            zh_s  (ji,jk) = MAX( 0._wp , zh_s  (ji,jk) + zdum ) 
     242!!$            IF( zh_s(ji,jk) == 0._wp )   ze_s(ji,jk) = 0._wp 
     243 
     244            ! update sublimation left 
     245            zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp ) 
     246         END DO 
     247      END DO 
     248 
     249      !       
    301250      !                       ! ============ ! 
    302251      !                       !     Ice      ! 
     
    305254      ! Surface ice melting  
    306255      !-------------------- 
    307       zdeltah(1:npti,:) = 0._wp ! important 
    308256      DO jk = 1, nlay_i 
    309257         DO ji = 1, npti 
     
    313261 
    314262               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]        
    315                zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0) 
    316                                                                       ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    317                zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
    318                                                                       ! this should normally not happen, but sometimes, heat diffusion leads to this 
    319                zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0 
    320                           
    321                dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting 
    322                 
    323                zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
    324  
    325                hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
    326                !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    327                sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux 
    328                !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    329                wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    330  
     263               zdE            =   0._wp                               ! Specific enthalpy difference (J/kg, <0) 
     264               !                                                          set up at 0 since no energy is needed to melt water...(it is already melted) 
     265               zdum           = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     266               !                                                          this should normally not happen, but sometimes, heat diffusion leads to this 
     267               zfmdt          = - zdum * rhoi                         ! Recompute mass flux [kg/m2, >0] 
     268               ! 
     269               dh_i_itm(ji)   = dh_i_itm(ji) + zdum                   ! Cumulate internal melting 
     270               ! 
     271               hfx_res_1d(ji) = hfx_res_1d(ji) + zEi  * zfmdt             * a_i_1d(ji) * r1_rdtice    ! Heat flux to the ocean [W.m-2], <0 
     272               !                                                                                          ice enthalpy zEi is "sent" to the ocean 
     273               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_rdtice    ! Mass flux 
     274               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice    ! Salt flux 
     275               !                                                                                          using s_i_1d and not sz_i_1d(jk) is ok 
    331276            ELSE                                        !-- Surface melting 
    332277                
     
    337282               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0] 
    338283                
    339                zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
    340                 
    341                zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
    342                 
    343                zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 
    344                 
    345                dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt 
    346                 
    347                zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
     284               zdum          = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
     285                
     286               zdum           = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     287 
     288               zq_top(ji)     = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat 
     289                
     290               dh_i_sum(ji)   = dh_i_sum(ji) + zdum                   ! Cumulate surface melt 
     291                
     292               zfmdt          = - rhoi * zdum                         ! Recompute mass flux [kg/m2, >0] 
    348293                
    349294               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    350295                
    351                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux >0 
    352                !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok) 
    353                hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0 
    354                hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], > 0   
    355                !  
    356                wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    357                 
     296               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt             * a_i_1d(ji) * r1_rdtice    ! Heat flux [W.m-2], < 0 
     297               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zdE  * zfmdt             * a_i_1d(ji) * r1_rdtice    ! Heat flux used in this process [W.m-2], > 0   
     298               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_rdtice    ! Mass flux 
     299               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice    ! Salt flux >0 
     300               !                                                                                          using s_i_1d and not sz_i_1d(jk) is ok)  
    358301            END IF 
    359              
     302            ! update thickness 
     303            zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 
     304            h_i_1d(ji)  = MAX( 0._wp, h_i_1d(ji)  + zdum ) 
     305            ! 
     306            ! update heat content (J.m-2) and layer thickness 
     307            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 
     308            h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 
     309            ! 
     310            ! 
    360311            ! Ice sublimation 
    361312            ! --------------- 
    362             zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 
    363             zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 
    364             dh_i_sub(ji)    = dh_i_sub(ji)    + zdum 
    365              
    366             sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice  ! Salt flux >0 
    367             !                                                                                          clem: flux is sent to the ocean for simplicity 
    368             !                                                                                                but salt should remain in the ice except 
    369             !                                                                                                if all ice is melted. => must be corrected 
    370             hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice      ! Heat flux [W.m-2], < 0 
    371  
    372             wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_rdtice           ! Mass flux > 0 
    373  
    374             ! update remaining mass flux 
    375             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi 
    376              
     313            zdum               = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi ) 
     314            ! 
     315            hfx_sub_1d(ji)     = hfx_sub_1d(ji)     + e_i_1d(ji,jk) * zdum              * a_i_1d(ji) * r1_rdtice ! Heat flux [W.m-2], < 0 
     316            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi          * zdum              * a_i_1d(ji) * r1_rdtice ! Mass flux > 0 
     317            sfx_sub_1d(ji)     = sfx_sub_1d(ji)     - rhoi          * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice ! Salt flux >0 
     318            !                                                                                                      clem: flux is sent to the ocean for simplicity 
     319            !                                                                                                            but salt should remain in the ice except 
     320            !                                                                                                            if all ice is melted. => must be corrected 
     321            ! update remaining mass flux and thickness 
     322            zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi             
     323            zh_i(ji,jk)    = MAX( 0._wp, zh_i(ji,jk) + zdum ) 
     324            h_i_1d(ji)     = MAX( 0._wp, h_i_1d(ji)  + zdum ) 
     325            dh_i_sub(ji)   = dh_i_sub(ji) + zdum 
     326 
     327            ! update heat content (J.m-2) and layer thickness 
     328            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 
     329            h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 
     330 
    377331            ! record which layers have disappeared (for bottom melting)  
    378332            !    => icount=0 : no layer has vanished 
    379333            !    => icount=5 : 5 layers have vanished 
    380             rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
     334            rswitch       = MAX( 0._wp , SIGN( 1._wp , - zh_i(ji,jk) ) )  
    381335            icount(ji,jk) = NINT( rswitch ) 
    382             zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    383336                         
    384             ! update heat content (J.m-2) and layer thickness 
    385             eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 
    386             h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    387          END DO 
    388       END DO 
    389        
    390       ! update ice thickness 
    391       DO ji = 1, npti 
    392          h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) ) 
    393       END DO 
    394  
     337         END DO 
     338      END DO 
     339       
    395340      ! remaining "potential" evap is sent to ocean 
    396341      DO ji = 1, npti 
     
    430375                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 ) 
    431376 
    432                s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity 
    433  
    434                ztmelts       = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C) 
    435  
    436                zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    437                 
    438                zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
    439                   &            - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 
    440  
    441                zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
    442  
    443                zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
    444  
    445                dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 
     377               s_i_new(ji)    = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity 
     378 
     379               ztmelts        = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C) 
     380 
     381               zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
     382                
     383               zEi            = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
     384                  &             - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts 
     385 
     386               zEw            = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
     387 
     388               zdE            = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
     389 
     390               dh_i_bog(ji)   = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 
    446391                
    447392            END DO 
    448393            ! Contribution to Energy and Salt Fluxes                                     
    449             zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
     394            zfmdt = - rhoi * dh_i_bog(ji)                                                              ! Mass flux x time step (kg/m2, < 0) 
    450395             
    451             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], >0 
    452             hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], <0 
    453              
    454             sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice     ! Salt flux, <0 
    455  
    456             wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice                   ! Mass flux, <0 
     396            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt                      * a_i_1d(ji) * r1_rdtice   ! Heat flux to the ocean [W.m-2], >0 
     397            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zdE  * zfmdt                      * a_i_1d(ji) * r1_rdtice   ! Heat flux used in this process [W.m-2], <0           
     398            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji)               * a_i_1d(ji) * r1_rdtice   ! Mass flux, <0 
     399            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * s_i_new(ji) * a_i_1d(ji) * r1_rdtice   ! Salt flux, <0 
     400 
     401            ! update thickness 
     402            zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji) 
     403            h_i_1d(ji)        = h_i_1d(ji)        + dh_i_bog(ji) 
    457404 
    458405            ! update heat content (J.m-2) and layer thickness 
     
    466413      ! Ice Basal melt 
    467414      !--------------- 
    468       zdeltah(1:npti,:) = 0._wp ! important 
    469415      DO jk = nlay_i, 1, -1 
    470416         DO ji = 1, npti 
     
    475421               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    476422 
    477                   zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0) 
    478                   zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    479                                                                     !    set up at 0 since no energy is needed to melt water...(it is already melted) 
    480                   zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
    481                                                                     ! this should normally not happen, but sometimes, heat diffusion leads to this 
    482  
    483                   dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk) 
    484  
    485                   zfmdt             = - zdeltah(ji,jk) * rhoi      ! Mass flux x time step > 0 
    486  
    487                   hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
    488                   !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    489                   sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux 
    490                   !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    491                   wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    492  
    493                   ! update heat content (J.m-2) and layer thickness 
    494                   eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 
    495                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    496  
     423                  zEi            = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0) 
     424                  zdE            = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
     425                  !                                                  set up at 0 since no energy is needed to melt water...(it is already melted) 
     426                  zdum           = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
     427                  !                                                  this should normally not happen, but sometimes, heat diffusion leads to this 
     428                  dh_i_itm (ji)  = dh_i_itm(ji) + zdum 
     429                  ! 
     430                  zfmdt          = - zdum * rhoi                 ! Mass flux x time step > 0 
     431                  ! 
     432                  hfx_res_1d(ji) = hfx_res_1d(ji) + zEi  * zfmdt             * a_i_1d(ji) * r1_rdtice   ! Heat flux to the ocean [W.m-2], <0 
     433                  !                                                                                         ice enthalpy zEi is "sent" to the ocean 
     434                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_rdtice   ! Mass flux 
     435                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice   ! Salt flux 
     436                  !                                                                                         using s_i_1d and not sz_i_1d(jk) is ok 
    497437               ELSE                                        !-- Basal melting 
    498438 
    499                   zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0) 
    500                   zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0) 
    501                   zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0) 
    502  
    503                   zfmdt           = - zq_bot(ji) / zdE                                        ! Mass flux x time step (kg/m2, >0) 
    504  
    505                   zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change 
    506  
    507                   zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change 
     439                  zEi            = - e_i_1d(ji,jk) * r1_rhoi                       ! Specific enthalpy of melting ice (J/kg, <0) 
     440                  zEw            = rcp * ztmelts                                   ! Specific enthalpy of meltwater (J/kg, <0) 
     441                  zdE            = zEi - zEw                                       ! Specific enthalpy difference   (J/kg, <0) 
     442 
     443                  zfmdt          = - zq_bot(ji) / zdE                              ! Mass flux x time step (kg/m2, >0) 
     444 
     445                  zdum           = - zfmdt * r1_rhoi                               ! Gross thickness change 
     446 
     447                  zdum           = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) )       ! bound thickness change 
    508448                   
    509                   zq_bot(ji)      = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors 
    510  
    511                   dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt 
    512  
    513                   zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0 
    514  
    515                   zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean 
    516  
    517                   hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0   
    518                   hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat used in this process [W.m-2], >0   
    519  
    520                   sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
    521                   !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok 
    522                    
    523                   wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    524  
    525                   ! update heat content (J.m-2) and layer thickness 
    526                   eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 
    527                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
     449                  zq_bot(ji)     = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors 
     450 
     451                  dh_i_bom(ji)   = dh_i_bom(ji) + zdum                             ! Update basal melt 
     452 
     453                  zfmdt          = - zdum * rhoi                                   ! Mass flux x time step > 0 
     454 
     455                  zQm            = zfmdt * zEw                                     ! Heat exchanged with ocean 
     456 
     457                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw  * zfmdt             * a_i_1d(ji) * r1_rdtice   ! Heat flux to the ocean [W.m-2], <0   
     458                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zdE  * zfmdt             * a_i_1d(ji) * r1_rdtice   ! Heat used in this process [W.m-2], >0 
     459                  wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum              * a_i_1d(ji) * r1_rdtice   ! Mass flux 
     460                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_rdtice   ! Salt flux 
     461                  !                                                                                         using s_i_1d and not sz_i_1d(jk) is ok 
    528462               ENDIF 
    529             
     463               ! update thickness 
     464               zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 
     465               h_i_1d(ji)  = MAX( 0._wp, h_i_1d(ji)  + zdum ) 
     466               ! 
     467               ! update heat content (J.m-2) and layer thickness 
     468               eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 
     469               h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 
    530470            ENDIF 
    531471         END DO 
    532472      END DO 
    533473 
    534       ! Update temperature, energy 
    535       ! -------------------------- 
    536       DO ji = 1, npti 
    537          h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) ) 
    538       END DO   
    539  
    540       ! If heat still available then melt more snow 
    541       !------------------------------------------- 
    542       zdeltah(1:npti,:) = 0._wp ! important 
    543       DO ji = 1, npti 
    544          zq_rema (ji)   = zq_top(ji) + zq_bot(ji)  
    545          rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow 
    546          rswitch        = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 
    547          zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 
    548          zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 
    549          dh_s_tot(ji)   = dh_s_tot(ji) + zdeltah(ji,1) 
    550          h_s_1d  (ji)   = h_s_1d  (ji) + zdeltah(ji,1) 
    551          
    552          zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2) 
    553          hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice   ! Heat used to melt snow, W.m-2 (>0) 
    554          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice       ! Mass flux 
    555          dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    556          !     
    557          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    558          !!!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
    559  
    560          IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
    561       END DO 
    562  
    563       ! 
     474      ! Remove snow if ice has melted entirely 
     475      ! -------------------------------------- 
     476      DO jk = 0, nlay_s 
     477         DO ji = 1,npti 
     478            IF( h_i_1d(ji) == 0._wp ) THEN 
     479               ! mass & energy loss to the ocean 
     480               hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice  ! heat flux to the ocean [W.m-2], < 0 
     481               wfx_res_1d(ji) = wfx_res_1d(ji) + rhos        * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice  ! mass flux 
     482 
     483               ! update thickness and energy 
     484               h_s_1d(ji)    = 0._wp 
     485               ze_s  (ji,jk) = 0._wp 
     486               zh_s  (ji,jk) = 0._wp 
     487            ENDIF 
     488         END DO 
     489      END DO 
     490       
     491      ! Snow load on ice 
     492      ! ----------------- 
     493      ! When snow load exceeds Archimede's limit and sst is positive, 
     494      ! snow-ice formation (next bloc) can lead to negative ice enthalpy. 
     495      ! Therefore we consider here that this excess of snow falls into the ocean 
     496      zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rau0) * r1_rhos 
     497      DO jk = 0, nlay_s 
     498         DO ji = 1, npti 
     499            IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN 
     500               ! snow layer thickness that falls into the ocean 
     501               zdum = MIN( zdeltah(ji) , zh_s(ji,jk) ) 
     502               ! mass & energy loss to the ocean 
     503               hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_rdtice  ! heat flux to the ocean [W.m-2], < 0 
     504               wfx_res_1d(ji) = wfx_res_1d(ji) + rhos        * zdum * a_i_1d(ji) * r1_rdtice  ! mass flux 
     505               ! update thickness and energy 
     506               h_s_1d(ji)    = MAX( 0._wp, h_s_1d(ji)  - zdum ) 
     507               zh_s  (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) 
     508               ! update snow thickness that still has to fall 
     509               zdeltah(ji)   = MAX( 0._wp, zdeltah(ji) - zdum ) 
     510            ENDIF 
     511         END DO 
     512      END DO 
     513       
    564514      ! Snow-Ice formation 
    565515      ! ------------------ 
    566       ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,  
    567       ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 
     516      ! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level,  
     517      ! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice) 
    568518      z1_rho = 1._wp / ( rhos+rau0-rhoi ) 
     519      zdeltah(1:npti) = 0._wp 
    569520      DO ji = 1, npti 
    570521         ! 
    571          dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho ) 
     522         dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho ) 
    572523 
    573524         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji) 
     
    579530         zQm            = zfmdt * zEw  
    580531          
    581          hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux 
    582  
    583          sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux 
     532         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw        * zfmdt * a_i_1d(ji) * r1_rdtice ! Heat flux 
     533         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_rdtice ! Salt flux 
    584534 
    585535         ! Case constant salinity in time: virtual salt flux to keep salinity constant 
    586536         IF( nn_icesal /= 2 )  THEN 
    587             sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean 
    588                &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice     ! and get  rn_icesal from the ocean  
     537            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d(ji) * zfmdt                 * a_i_1d(ji) * r1_rdtice  & ! put back sss_m     into the ocean 
     538               &                            - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_rdtice     ! and get  rn_icesal from the ocean  
    589539         ENDIF 
    590540 
    591541         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 
    592          wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice 
    593          wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_rdtice 
     542         wfx_sni_1d    (ji) = wfx_sni_1d    (ji) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_rdtice 
     543         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_rdtice 
     544 
     545         ! update thickness 
     546         zh_i(ji,0)  = zh_i(ji,0) + dh_snowice(ji) 
     547         zdeltah(ji) =              dh_snowice(ji) 
    594548 
    595549         ! update heat content (J.m-2) and layer thickness 
    596          eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw 
    597550         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 
    598           
    599       END DO 
    600  
    601       ! 
    602       ! Update temperature, energy 
    603       ! -------------------------- 
    604       DO ji = 1, npti 
    605          rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )  
    606          t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 
    607       END DO 
    608  
     551         eh_i_old(ji,0) = eh_i_old(ji,0) + zfmdt * zEw           ! 1st part (sea water enthalpy) 
     552 
     553      END DO 
     554      ! 
     555      DO jk = nlay_s, 0, -1   ! flooding of snow starts from the base 
     556         DO ji = 1, npti 
     557            zdum           = MIN( zdeltah(ji), zh_s(ji,jk) )     ! amount of snw that floods, > 0 
     558            zh_s(ji,jk)    = MAX( 0._wp, zh_s(ji,jk) - zdum )    ! remove some snow thickness 
     559            eh_i_old(ji,0) = eh_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy) 
     560            ! update dh_snowice 
     561            zdeltah(ji)    = MAX( 0._wp, zdeltah(ji) - zdum ) 
     562         END DO 
     563      END DO 
     564      ! 
     565      ! 
     566!!$      ! --- Update snow diags --- ! 
     567!!$      !!clem: this is wrong. dh_s_tot is not used anyway 
     568!!$      DO ji = 1, npti 
     569!!$         dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah(ji) + zdh_s_sub(ji) - dh_snowice(ji) 
     570!!$      END DO 
     571      ! 
     572      ! 
     573      ! Remapping of snw enthalpy on a regular grid 
     574      !-------------------------------------------- 
     575      CALL snw_ent( zh_s, ze_s, e_s_1d ) 
     576       
     577      ! recalculate t_s_1d from e_s_1d 
    609578      DO jk = 1, nlay_s 
    610579         DO ji = 1,npti 
    611             ! where there is no ice or no snow 
    612             rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) ) 
    613             ! mass & energy loss to the ocean 
    614             hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * & 
    615                &                              ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! heat flux to the ocean [W.m-2], < 0 
    616             wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * & 
    617                &                              ( rhos          * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_rdtice )  ! mass flux 
    618             ! update energy (mass is updated in the next loop) 
    619             e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 
    620             ! recalculate t_s_1d from e_s_1d 
    621             t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 
    622          END DO 
    623       END DO 
     580            IF( h_s_1d(ji) > 0._wp ) THEN 
     581               t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 
     582            ELSE 
     583               t_s_1d(ji,jk) = rt0 
     584            ENDIF 
     585         END DO 
     586      END DO 
     587 
     588      ! Note: remapping of ice enthalpy is done in icethd.F90 
    624589 
    625590      ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 --- 
    626591      WHERE( h_i_1d(1:npti) == 0._wp )    
    627          a_i_1d(1:npti) = 0._wp 
    628          h_s_1d(1:npti) = 0._wp 
     592         a_i_1d (1:npti) = 0._wp 
     593         h_s_1d (1:npti) = 0._wp 
     594         t_su_1d(1:npti) = rt0 
    629595      END WHERE 
    630       ! 
     596       
    631597   END SUBROUTINE ice_thd_dh 
    632598 
     599   SUBROUTINE snw_ent( ph_old, pe_old, pe_new ) 
     600      !!------------------------------------------------------------------- 
     601      !!               ***   ROUTINE snw_ent  *** 
     602      !! 
     603      !! ** Purpose : 
     604      !!           This routine computes new vertical grids in the snow,  
     605      !!           and consistently redistributes temperatures.  
     606      !!           Redistribution is made so as to ensure to energy conservation 
     607      !! 
     608      !! 
     609      !! ** Method  : linear conservative remapping 
     610      !!            
     611      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 
     612      !!            2) linear remapping on the new layers 
     613      !! 
     614      !! ------------ cum0(0)                        ------------- cum1(0) 
     615      !!                                    NEW      ------------- 
     616      !! ------------ cum0(1)               ==>      ------------- 
     617      !!     ...                                     ------------- 
     618      !! ------------                                ------------- 
     619      !! ------------ cum0(nlay_s+1)                 ------------- cum1(nlay_s) 
     620      !! 
     621      !! 
     622      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 
     623      !!------------------------------------------------------------------- 
     624      REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in   ) ::   ph_old             ! old thicknesses (m) 
     625      REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in   ) ::   pe_old             ! old enthlapies (J.m-3) 
     626      REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) ::   pe_new             ! new enthlapies (J.m-3, remapped) 
     627      ! 
     628      INTEGER  :: ji         !  dummy loop indices 
     629      INTEGER  :: jk0, jk1   !  old/new layer indices 
     630      ! 
     631      REAL(wp), DIMENSION(jpij,0:nlay_s+1) ::   zeh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces 
     632      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces 
     633      REAL(wp), DIMENSION(jpij)            ::   zhnew               ! new layers thicknesses 
     634      !!------------------------------------------------------------------- 
     635 
     636      !-------------------------------------------------------------------------- 
     637      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces 
     638      !-------------------------------------------------------------------------- 
     639      zeh_cum0(1:npti,0) = 0._wp  
     640      zh_cum0 (1:npti,0) = 0._wp 
     641      DO jk0 = 1, nlay_s+1 
     642         DO ji = 1, npti 
     643            zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + pe_old(ji,jk0-1) * ph_old(ji,jk0-1) 
     644            zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + ph_old(ji,jk0-1) 
     645         END DO 
     646      END DO 
     647 
     648      !------------------------------------ 
     649      !  2) Interpolation on the new layers 
     650      !------------------------------------ 
     651      ! new layer thickesses 
     652      DO ji = 1, npti 
     653         zhnew(ji) = SUM( ph_old(ji,0:nlay_s) ) * r1_nlay_s   
     654      END DO 
     655 
     656      ! new layers interfaces 
     657      zh_cum1(1:npti,0) = 0._wp 
     658      DO jk1 = 1, nlay_s 
     659         DO ji = 1, npti 
     660            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 
     661         END DO 
     662      END DO 
     663 
     664      zeh_cum1(1:npti,0:nlay_s) = 0._wp  
     665      ! new cumulative q*h => linear interpolation 
     666      DO jk0 = 1, nlay_s+1 
     667         DO jk1 = 1, nlay_s-1 
     668            DO ji = 1, npti 
     669               IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 
     670                  zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  & 
     671                     &                 zeh_cum0(ji,jk0  ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) )  & 
     672                     &             / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 
     673               ENDIF 
     674            END DO 
     675         END DO 
     676      END DO 
     677      ! to ensure that total heat content is strictly conserved, set: 
     678      zeh_cum1(1:npti,nlay_s) = zeh_cum0(1:npti,nlay_s+1)  
     679 
     680      ! new enthalpies 
     681      DO jk1 = 1, nlay_s 
     682         DO ji = 1, npti 
     683            rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) )  
     684            pe_new(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 
     685         END DO 
     686      END DO 
     687       
     688   END SUBROUTINE snw_ent 
     689 
     690    
    633691#else 
    634692   !!---------------------------------------------------------------------- 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_pnd.F90

    r13284 r14026  
    8383      !!------------------------------------------------------------------- 
    8484      INTEGER  ::   ji        ! loop indices 
     85      REAL(wp) ::   zdv_pnd   ! Amount of water going into the ponds & lids 
    8586      !!------------------------------------------------------------------- 
    8687      DO ji = 1, npti 
    8788         ! 
    88          IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 
     89         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 
     90         ! 
     91         IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 
    8992            h_ip_1d(ji)      = rn_hpnd     
    9093            a_ip_1d(ji)      = rn_apnd * a_i_1d(ji) 
     
    9699         ENDIF 
    97100         ! 
     101         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 
     102         wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_rdtice 
     103         ! 
    98104      END DO 
    99105      ! 
     
    134140      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 --- 
    135141      !! 
    136       !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 --- 
    137       !!                                     perm = permability of sea-ice 
     142      !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi * flush         --- from Flocco et al 2007 --- 
     143      !!                                     perm = permability of sea-ice                              + correction from Hunke et al 2012 (flush) 
    138144      !!                                     visc = water viscosity 
    139145      !!                                     Hp   = height of top of the pond above sea-level 
    140146      !!                                     Hi   = ice thickness thru which there is flushing 
     147      !!                                     flush= correction otherwise flushing is excessive 
    141148      !! 
    142149      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness 
     
    147154      !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 
    148155      !!  
    149       !! ** Note       :   mostly stolen from CICE 
     156      !! ** Note       :   mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 
    150157      !! 
    151158      !! ** References :   Flocco and Feltham (JGR, 2007) 
    152159      !!                   Flocco et al       (JGR, 2010) 
    153160      !!                   Holland et al      (J. Clim, 2012) 
     161      !!                   Hunke et al        (OM 2012) 
    154162      !!------------------------------------------------------------------- 
    155163      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
     
    158166      REAL(wp), PARAMETER ::   zTp     = -2._wp       ! reference temperature 
    159167      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity 
    160       !! 
    161       REAL(wp) ::   zfr_mlt, zdv_mlt                  ! fraction and volume of available meltwater retained for melt ponding 
     168      REAL(wp), PARAMETER ::   zflush  =  1.e-02_wp   ! tuning param to reduce flushing (from Hunke et al 2012) 
     169      !! 
     170      REAL(wp) ::   zfr_mlt, zdv_mlt, zdv_avail       ! fraction and volume of available meltwater retained for melt ponding 
    162171      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes 
     172      REAL(wp) ::   zdv_pnd                           ! Amount of water going into the ponds & lids 
    163173      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh 
    164174      REAL(wp) ::   zv_ip_max                         ! max pond volume allowed 
    165175      REAL(wp) ::   zdT                               ! zTp-t_su 
    166       REAL(wp) ::   zsbr                              ! Brine salinity 
     176      REAL(wp) ::   zsbr, ztmelts                     ! Brine salinity 
    167177      REAL(wp) ::   zperm                             ! permeability of sea ice 
    168178      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
     
    176186 
    177187      DO ji = 1, npti 
     188         ! 
     189         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 
    178190         !                                                            !----------------------------------------------------! 
    179          IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
     191         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < 0.01_wp ) THEN   ! Case ice thickness < rn_himin or tiny ice fraction ! 
    180192            !                                                         !----------------------------------------------------! 
    181193            !--- Remove ponds on thin ice or tiny ice fractions 
    182             a_ip_1d(ji)      = 0._wp 
    183             h_ip_1d(ji)      = 0._wp 
    184             h_il_1d(ji)      = 0._wp 
     194            a_ip_1d(ji) = 0._wp 
     195            h_ip_1d(ji) = 0._wp 
     196            h_il_1d(ji) = 0._wp 
    185197            !                                                         !--------------------------------! 
    186198         ELSE                                                         ! Case ice thickness >= rn_himin ! 
     
    193205            !------------------! 
    194206            ! 
    195             !--- available meltwater for melt ponding ---! 
    196             zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    197             zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed 
    198             zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
     207            !--- available meltwater for melt ponding (zdv_avail) ---! 
     208            zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 
     209            zfr_mlt   = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed 
     210            zdv_mlt   = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors?  
    199211            ! 
    200212            !--- overflow ---! 
    201             ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
    202             !    a_ip_max = zfr_mlt * a_i 
    203             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     213            ! 
     214            ! area driven overflow 
     215            !    If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
     216            !       a_ip_max = zfr_mlt * a_i 
     217            !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    204218            zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    205             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    206  
    207             ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    208             !    h_ip_max = 0.5 * h_i 
    209             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     219            zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     220 
     221            ! depth driven overflow 
     222            !    If pond depth exceeds half the ice thickness then reduce the pond volume 
     223            !       h_ip_max = 0.5 * h_i 
     224            !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    210225            zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
    211             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     226            zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    212227             
    213228            !--- Pond growing ---! 
     
    217232            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
    218233            ! 
    219             !--- mass flux ---! 
    220             IF( zdv_mlt > 0._wp ) THEN 
    221                zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
    222                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    223                ! 
    224                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
    225                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    226                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    227             ENDIF 
    228  
    229234            !-------------------! 
    230235            ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
     
    237242               ! 
    238243               !--- Lid growing and subsequent pond shrinking ---!  
    239                zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
     244               zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
    240245                  &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
    241246                
    242247               ! Lid growing 
    243                v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
     248               v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) 
    244249                
    245250               ! Pond shrinking 
    246                v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
     251               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 
    247252 
    248253            ELSE 
     254               zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp )  ! Holland 2012 (eq. 6)  
    249255               ! Pond shrinking 
    250                v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
     256               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 
    251257            ENDIF 
    252258            ! 
     
    254260            ! v_ip     = h_ip * a_ip 
    255261            ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 
    256             a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    257             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    258  
    259             !---------------!             
    260             ! Pond flushing ! 
    261             !---------------! 
     262            a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
     263            h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     264             
     265            !------------------------------------------------!             
     266            ! Pond drainage through brine network (flushing) ! 
     267            !------------------------------------------------! 
    262268            ! height of top of the pond above sea-level 
    263269            zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 
     
    265271            ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
    266272            DO jk = 1, nlay_i 
    267                zsbr = - 1.2_wp                                  & 
    268                   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
    269                   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
    270                   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
    271                ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     273               ! MV Assur is inconsistent with SI3 
     274               !!zsbr = - 1.2_wp                                  & 
     275               !!   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
     276               !!   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
     277               !!   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
     278               !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     279               ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 
     280               ztmelts  = -rTmlt * sz_i_1d(ji,jk) 
     281               ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 ) 
    272282            END DO 
    273283            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
    274284             
    275285            ! Do the drainage using Darcy's law 
    276             zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 
    277             zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
     286            zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * zflush ! zflush comes from Hunke et al. (2012) 
     287            zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0  
    278288            v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
    279289             
    280290            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    281             a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    282             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     291            a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
     292            h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    283293 
    284294            !--- Corrections and lid thickness ---! 
    285295            IF( ln_pnd_lids ) THEN 
    286296               !--- retrieve lid thickness from volume ---! 
    287                IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
    288                ELSE                              ;   h_il_1d(ji) = 0._wp 
     297               IF( a_ip_1d(ji) > 0.01_wp ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
     298               ELSE                               ;   h_il_1d(ji) = 0._wp 
    289299               ENDIF 
    290300               !--- remove ponds if lids are much larger than ponds ---! 
    291301               IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
    292                   a_ip_1d(ji)      = 0._wp 
    293                   h_ip_1d(ji)      = 0._wp 
    294                   h_il_1d(ji)      = 0._wp 
     302                  a_ip_1d(ji) = 0._wp 
     303                  h_ip_1d(ji) = 0._wp 
     304                  h_il_1d(ji) = 0._wp 
    295305               ENDIF 
    296306            ENDIF 
     307 
     308!!$            ! diagnostics: dvpnd = mlt+rnf+frz+drn 
     309!!$            diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow *   zdv_avail             * r1_rdtice   ! > 0, surface melt input 
     310!!$            diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_rdtice   ! < 0, runoff 
     311!!$            diag_dvpn_frz_1d(ji) = diag_dvpn_frz_1d(ji) + rhow *   zdv_frz               * r1_rdtice   ! < 0, shrinking 
     312!!$            diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow *   zdv_flush             * r1_rdtice   ! < 0, drainage 
    297313            ! 
    298314         ENDIF 
    299           
     315         ! 
     316         zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 
     317         wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_rdtice 
     318         ! 
    300319      END DO 
    301320      ! 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icethd_zdf_bl99.F90

    r13284 r14026  
    109109      REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function 
    110110      ! 
    111       REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    112       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    113       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
    114       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence 
    115       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    116       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    117       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i_cp ! copy 
    118       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    119       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
    120       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice 
    121       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice 
    122       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow 
    123       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow 
    124       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    125       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
    126       REAL(wp), DIMENSION(jpij)            ::   zkappa_comb ! Combined snow and ice surface conductivity 
    127       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    128       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
    129       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term 
    130       REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms 
    131       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    132       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    133       REAL(wp), DIMENSION(jpij)            ::   za_s_fra    ! ice fraction covered by snow  
    134       REAL(wp), DIMENSION(jpij)            ::   isnow       ! snow presence (1) or not (0)  
    135       REAL(wp), DIMENSION(jpij)            ::   isnow_comb  ! snow presence for met-office  
     111      REAL(wp), DIMENSION(jpij       )   ::   ztsuold     ! Old surface temperature in the ice 
     112      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztiold      ! Old temperature in the ice 
     113      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsold      ! Old temperature in the snow 
     114      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Temporary temperature in the ice to check the convergence 
     115      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow to check the convergence 
     116      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity 
     117      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   ztcond_i_cp ! copy 
     118      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice 
     119      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice 
     120      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice 
     121      REAL(wp), DIMENSION(jpij,0:nlay_i) ::   zeta_i      ! Eta factor in the ice 
     122      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow 
     123      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow 
     124      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow 
     125      REAL(wp), DIMENSION(jpij,0:nlay_s) ::   zeta_s      ! Eta factor in the snow 
     126      REAL(wp), DIMENSION(jpij)          ::   zkappa_comb ! Combined snow and ice surface conductivity 
     127      REAL(wp), DIMENSION(jpij)          ::   zq_ini      ! diag errors on heat 
     128      REAL(wp), DIMENSION(jpij)          ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
     129      REAL(wp), DIMENSION(jpij)          ::   za_s_fra    ! ice fraction covered by snow  
     130      REAL(wp), DIMENSION(jpij)          ::   isnow       ! snow presence (1) or not (0)  
     131      REAL(wp), DIMENSION(jpij)          ::   isnow_comb  ! snow presence for met-office  
     132      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindterm    ! 'Ind'ependent term 
     133      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zindtbis    ! Temporary 'ind'ependent term 
     134      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1)   ::   zdiagbis    ! Temporary 'dia'gonal term 
     135      REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) ::   ztrid       ! Tridiagonal system terms 
    136136      ! 
    137137      ! Mono-category 
     
    533533            ! Solve the tridiagonal system with Gauss elimination method. 
    534534            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    535             jm_maxt = 0 
    536             jm_mint = nlay_i+5 
    537             DO ji = 1, npti 
    538                jm_mint = MIN(jm_min(ji),jm_mint) 
    539                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    540             END DO 
    541  
    542             DO jk = jm_mint+1, jm_maxt 
    543                DO ji = 1, npti 
    544                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     535!!$            jm_maxt = 0 
     536!!$            jm_mint = nlay_i+5 
     537!!$            DO ji = 1, npti 
     538!!$               jm_mint = MIN(jm_min(ji),jm_mint) 
     539!!$               jm_maxt = MAX(jm_max(ji),jm_maxt) 
     540!!$            END DO 
     541!!$            !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 
     542!!$             
     543!!$            DO jk = jm_mint+1, jm_maxt 
     544!!$               DO ji = 1, npti 
     545!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     546!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
     547!!$                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
     548!!$               END DO 
     549!!$            END DO 
     550            ! clem: maybe one should find a way to reverse this loop for mpi performance 
     551            DO ji = 1, npti 
     552               jm_mint = jm_min(ji) 
     553               jm_maxt = jm_max(ji) 
     554               DO jm = jm_mint+1, jm_maxt 
    545555                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    546556                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
     
    564574            END DO 
    565575 
     576            ! snow temperatures       
    566577            DO ji = 1, npti 
    567578               ! Variables used after iterations 
    568579               ! Value must be frozen after convergence for MPP independance reason 
    569                IF ( .NOT. l_T_converged(ji) ) THEN 
    570                   ! snow temperatures       
    571                   IF( h_s_1d(ji) > 0._wp ) THEN 
    572                      t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
    573                   ENDIF 
    574                   ! surface temperature 
     580               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     581                  &   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     582            END DO 
     583            !!clem SNWLAY 
     584            DO jm = nlay_s, 2, -1 
     585               DO ji = 1, npti 
     586                  jk = jm - 1 
     587                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     588                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     589               END DO 
     590            END DO 
     591             
     592            ! surface temperature 
     593            DO ji = 1, npti 
     594               IF( .NOT. l_T_converged(ji) ) THEN 
    575595                  ztsub(ji) = t_su_1d(ji) 
    576596                  IF( t_su_1d(ji) < rt0 ) THEN 
    577                      t_su_1d(ji) = (  zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    578                         &           ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     597                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
     598                        &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    579599                  ENDIF 
    580600               ENDIF 
    581601            END DO 
    582             !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
    583602            ! 
    584603            !-------------------------------------------------------------- 
     
    727746            ! Solve the tridiagonal system with Gauss elimination method. 
    728747            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    729             jm_maxt = 0 
    730             jm_mint = nlay_i+5 
    731             DO ji = 1, npti 
    732                jm_mint = MIN(jm_min(ji),jm_mint) 
    733                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    734             END DO 
    735              
    736             DO jk = jm_mint+1, jm_maxt 
    737                DO ji = 1, npti 
    738                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     748!!$            jm_maxt = 0 
     749!!$            jm_mint = nlay_i+5 
     750!!$            DO ji = 1, npti 
     751!!$               jm_mint = MIN(jm_min(ji),jm_mint) 
     752!!$               jm_maxt = MAX(jm_max(ji),jm_maxt) 
     753!!$            END DO 
     754!!$             
     755!!$            DO jk = jm_mint+1, jm_maxt 
     756!!$               DO ji = 1, npti 
     757!!$                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
     758!!$                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
     759!!$                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1) 
     760!!$               END DO 
     761!!$            END DO 
     762            ! clem: maybe one should find a way to reverse this loop for mpi performance 
     763            DO ji = 1, npti 
     764               jm_mint = jm_min(ji) 
     765               jm_maxt = jm_max(ji) 
     766               DO jm = jm_mint+1, jm_maxt 
    739767                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    740                   zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)  / zdiagbis(ji,jm-1) 
     768                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
    741769               END DO 
    742770            END DO 
    743              
     771 
    744772            ! ice temperatures 
    745773            DO ji = 1, npti 
     
    761789            ! snow temperatures       
    762790            DO ji = 1, npti 
    763                ! Variable used after iterations 
     791               ! Variables used after iterations 
    764792               ! Value must be frozen after convergence for MPP independance reason 
    765                IF ( .NOT. l_T_converged(ji) ) THEN 
    766                   IF( h_s_1d(ji) > 0._wp ) THEN 
    767                      t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
    768                   ENDIF 
    769                ENDIF 
    770             END DO 
    771             !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
     793               IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     794                  &   t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 
     795            END DO 
     796            !!clem SNWLAY 
     797            DO jm = nlay_s, 2, -1 
     798               DO ji = 1, npti 
     799                  jk = jm - 1 
     800                  IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 
     801                     &   t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     802               END DO 
     803            END DO 
    772804            ! 
    773805            !-------------------------------------------------------------- 
     
    923955         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    924956         IF( h_s_1d(ji) >= zhs_ssl ) THEN 
    925             t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,1)   & 
    926                &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 
     957            t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s)   & 
     958               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1)      ) & 
    927959               &          / ( rn_cnd_s       * h_i_1d(ji) * r1_nlay_i & 
    928960               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceupdate.F90

    r13642 r14026  
    9595      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace 
    9696      !!--------------------------------------------------------------------- 
    97       IF( ln_timing )   CALL timing_start('ice_update') 
     97      IF( ln_timing )   CALL timing_start('iceupdate') 
    9898 
    9999      IF( kt == nit000 .AND. lwp ) THEN 
     
    156156            ! ice-ocean  mass flux 
    157157            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    158                &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 
     158               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
    159159 
    160160            ! snw-ocean mass flux 
     
    162162 
    163163            ! total mass flux at the ocean/ice interface 
    164             fmmflx(ji,jj) =                - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! ice-ocean mass flux saved at least for biogeochemical model 
    165             emp   (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! atm-ocean + ice-ocean mass flux 
     164            fmmflx(ji,jj) =                - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj)   ! ice-ocean mass flux saved at least for biogeochemical model 
     165            emp   (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj)   ! atm-ocean + ice-ocean mass flux 
    166166 
    167167            ! Salt flux at the ocean surface       
     
    174174            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step 
    175175            !                                               ! new mass per unit area 
    176             snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  
     176            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * ( vt_ip(ji,jj) + vt_il(ji,jj) ) )  
    177177            !                                               ! time evolution of snow+ice mass 
    178178            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice 
     
    289289      IF( ln_icectl    )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
    290290      IF( ln_ctl       )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
    291       IF( ln_timing    )   CALL timing_stop   ('ice_update')                                      ! timing 
     291      IF( ln_timing    )   CALL timing_stop   ('iceupdate')                                       ! timing 
    292292      ! 
    293293   END SUBROUTINE ice_update_flx 
     
    327327      REAL(wp) ::   zflagi                          !   -      - 
    328328      !!--------------------------------------------------------------------- 
    329       IF( ln_timing )   CALL timing_start('ice_update_tau') 
     329      IF( ln_timing )   CALL timing_start('iceupdate') 
    330330 
    331331      IF( kt == nit000 .AND. lwp ) THEN 
     
    383383      CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1., vtau, 'V', -1. )   ! lateral boundary condition 
    384384      ! 
    385       IF( ln_timing )   CALL timing_stop('ice_update_tau') 
     385      IF( ln_timing )   CALL timing_stop('iceupdate') 
    386386      !   
    387387   END SUBROUTINE ice_update_tau 
  • NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icevar.F90

    r13284 r14026  
    228228      z1_zhmax =  1._wp / hi_max(jpl)                
    229229      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 
    230          h_i  (:,:,jpl) = zhmax 
     230         h_i   (:,:,jpl) = zhmax 
    231231         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax  
    232232         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) 
     
    244244      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow 
    245245      ELSEWHERE                            ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * &   ! lid is in between. Expose part of the pond 
    246          &                                                       ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min ) 
     246         &                                                       ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 
    247247      END WHERE 
    248248      ! 
     
    545545            DO ji = 1 , jpi 
    546546               ! update exchanges with ocean 
    547                sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice 
    548                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice 
    549                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice 
     547               sfx_res(ji,jj)  = sfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice 
     548               wfx_res(ji,jj)  = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice 
     549               wfx_res(ji,jj)  = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice 
     550               wfx_pnd(ji,jj)  = wfx_pnd(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_rdtice 
    550551               ! 
    551552               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 
     
    562563               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
    563564               v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 
     565               h_ip (ji,jj,jl) = h_ip (ji,jj,jl) * zswitch(ji,jj) 
     566               h_il (ji,jj,jl) = h_il (ji,jj,jl) * zswitch(ji,jj) 
    564567               ! 
    565568            END DO 
     
    656659                  psv_i  (ji,jj,jl) = 0._wp 
    657660               ENDIF 
     661               IF( pv_ip(ji,jj,jl) < 0._wp .OR. pv_il(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 
     662                  wfx_pnd(ji,jj)    = wfx_pnd(ji,jj) + pv_il(ji,jj,jl) * rhow * z1_dt 
     663                  pv_il  (ji,jj,jl) = 0._wp 
     664               ENDIF 
     665               IF( pv_ip(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 
     666                  wfx_pnd(ji,jj)    = wfx_pnd(ji,jj) + pv_ip(ji,jj,jl) * rhow * z1_dt 
     667                  pv_ip  (ji,jj,jl) = 0._wp 
     668               ENDIF 
    658669            END DO 
    659670         END DO 
     
    665676      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp 
    666677      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp 
    667       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    668       WHERE( pv_il (:,:,:) < 0._wp )   pv_il (:,:,:) = 0._wp !    but it does not change conservation, so keep it this way is ok 
    669678      ! 
    670679   END SUBROUTINE ice_var_zapneg 
Note: See TracChangeset for help on using the changeset viewer.