New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13957 – NEMO

Changeset 13957


Ignore:
Timestamp:
2020-12-01T23:04:18+01:00 (3 years ago)
Author:
clem
Message:

correct bugs in level ponds. I do not think they were working anymore. Plus transform virtual ponds into real ponds (i.e. proper ice-pond mass exchanges). The system is now: ice+snow+ponds. Plus make sure conservation is still ensured. Plus add a ctl_stop in topo melt ponds since work is still ongoing

Location:
NEMO/branches/2020/SI3_martin_ponds
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/SI3_martin_ponds/cfgs/SHARED/field_def_nemo-ice.xml

    r13908 r13957  
    5151          <field id="icehlid"      long_name="melt pond lid depth"                                     standard_name="sea_ice_meltpondlid_depth"                 unit="m" />  
    5252          <field id="icevlid"      long_name="melt pond lid volume"                                    standard_name="sea_ice_meltpondlid_volume"                unit="m" />  
    53           <field id="dvpn_mlt"     long_name="pond volume tendency due to surface melt"                standard_name="sea_ice_pondvolume_tendency_melt"          unit="cm/d" /> 
    54           <field id="dvpn_lid"     long_name="pond volume tendency due to exchanges with lid"          standard_name="sea_ice_pondvolume_tendency_lids"          unit="cm/d" /> 
    55           <field id="dvpn_rnf"     long_name="pond volume tendency due to runoff"                      standard_name="sea_ice_pondvolume_tendency_runoff"        unit="cm/d" /> 
    56           <field id="dvpn_drn"     long_name="pond volume tendency due to drainage"                    standard_name="sea_ice_pondvolume_tendency_drainage"      unit="cm/d" />   
     53          <field id="dvpn_mlt"     long_name="pond volume tendency due to surface melt"                standard_name="sea_ice_pondvolume_tendency_melt"          unit="kg/m2/s" /> 
     54          <field id="dvpn_lid"     long_name="pond volume tendency due to exchanges with lid"          standard_name="sea_ice_pondvolume_tendency_lids"          unit="kg/m2/s" /> 
     55          <field id="dvpn_rnf"     long_name="pond volume tendency due to runoff"                      standard_name="sea_ice_pondvolume_tendency_runoff"        unit="kg/m2/s" /> 
     56          <field id="dvpn_drn"     long_name="pond volume tendency due to drainage"                    standard_name="sea_ice_pondvolume_tendency_drainage"      unit="kg/m2/s" />   
    5757      
    5858     <!-- heat --> 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/ice.F90

    r13908 r13957  
    311311   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice         !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1] 
    312312 
    313    ! meltwater arrays to save for melt ponds (mv - could be grouped in a single meltwater volume array) 
    314    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m] 
    315    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m] 
    316  
    317313   !!---------------------------------------------------------------------- 
    318314   !! * Ice global state variables 
    … …  
    368364   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   vt_il         !: total melt pond lid volume per gridcell area [m] 
    369365 
     366   ! meltwater arrays to save for melt ponds (mv - could be grouped in a single meltwater volume array) 
     367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m] 
     368   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m] 
     369 
    370370   !!---------------------------------------------------------------------- 
    371371   !! * Global variables at before time step 
    372372   !!---------------------------------------------------------------------- 
    373373   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 
     374   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_ip_b, v_il_b             !: ponds and lids volumes 
    374375   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b              !: 
    375376   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                      !: snow heat content 
    … …  
    398399   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s]  
    399400   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_aice         !: ice conc.  variation   [s-1]  
     401   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vpnd         !: pond volume variation  [m/s]  
    400402   ! 
    401403   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_mass     !: advection of mass (kg/m2/s) 
    … …  
    464466      ii = ii + 1 
    465467      ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) ,  & 
    466          &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) ,                    & 
    467468         &      h_i        (jpi,jpj,jpl) , a_i    (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,  & 
    468469         &      v_s        (jpi,jpj,jpl) , h_s    (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,  & 
    … …  
    485486      ii = ii + 1 
    486487      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  & 
    487          &      v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) ) 
     488         &      v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) ,                     & 
     489         &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , STAT = ierr(ii) ) 
    488490 
    489491      ii = ii + 1 
    … …  
    493495      ii = ii + 1 
    494496      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),         & 
     497         &      v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) ,                                                         & 
    495498         &      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) , & 
    496499         &      STAT=ierr(ii) ) 
    … …  
    507510      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),                      &  
    508511         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),                      & 
    509          &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj), & 
     512         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), & 
    510513         &      diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 
    511514 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icectl.F90

    r13601 r13957  
    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_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 
    9596         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 ) 
    … …  
    112113 
    113114         ! -- mass diag -- ! 
    114          zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_Dt_ice       & 
     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_Dt_ice                                                                          & 
    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' ) & 
    167                &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice 
     175               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 
    168176            IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 
    169                &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rdt_ice 
     177               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rDt_ice 
    170178         ENDIF 
    171179         ! 
    … …  
    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_Dt_ice                             & 
     271         zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice    & 
    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 ) 
    … …  
    776790      ! -- mass diag -- ! 
    777791      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub & 
    778          &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice 
     792         &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rDt_ice 
    779793      zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 
    780794 
    781795      ! -- salt diag -- ! 
    782       zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3 
     796      zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rDt_ice * 1.e-3 
    783797      zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 
    784798 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_pra.F90

    r13908 r13957  
    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 ) & 
    … …  
    320321         ! 
    321322         ! --- diagnostics --- ! 
    322          diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     323         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     324            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 
    323325            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
    324326         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_umx.F90

    r13911 r13957  
    182182 
    183183         ! diagnostics 
    184          zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos 
     184         zdiag_adv_mass(:,:) =   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     185            &                  + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 
    185186         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi 
    186187         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & 
    … …  
    379380         ! 
    380381         ! --- diagnostics --- ! 
    381          diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 
     382         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 
     383            &                                        + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 
    382384            &                                        - zdiag_adv_mass(:,:) ) * z1_dt 
    383385         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi & 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceistate.F90

    r13472 r13957  
    426426      ! 4) Snow-ice mass (case ice is fully embedded) 
    427427      !---------------------------------------------- 
    428       snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass 
     428      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3  )   ! snow+ice mass 
    429429      snwice_mass_b(:,:) = snwice_mass(:,:) 
    430430      ! 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceitd.F90

    r13908 r13957  
    2929   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    3030   USE prtctl         ! Print control 
     31   USE timing         ! Timing 
    3132 
    3233   IMPLICIT NONE 
    … …  
    8788      REAL(wp), DIMENSION(jpij,0:jpl) ::   zhbnew          ! new boundaries of ice categories 
    8889      !!------------------------------------------------------------------ 
     90      IF( ln_timing )   CALL timing_start('iceitd_rem') 
    8991 
    9092      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution'  
    … …  
    328330      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    329331      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     332      IF( ln_timing    )   CALL timing_stop ('iceitd_rem') 
    330333      ! 
    331334   END SUBROUTINE ice_itd_rem 
    … …  
    491494                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 
    492495                  !                                               
    493                   ztrans          = v_ip_2d(ji,jl1) * zworka(ji)     ! Pond volume (also proportional to da/a) 
     496                  ztrans          = v_ip_2d(ji,jl1) * zworkv(ji)     ! Pond volume 
    494497                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    495498                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
    496499                  ! 
    497500                  IF ( ln_pnd_lids ) THEN                            ! Pond lid volume 
    498                      ztrans          = v_il_2d(ji,jl1) * zworka(ji) 
     501                     ztrans          = v_il_2d(ji,jl1) * zworkv(ji) 
    499502                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 
    500503                     v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans 
    … …  
    606609      REAL(wp), DIMENSION(jpij,jpl-1) ::   zdaice, zdvice   ! ice area and volume transferred 
    607610      !!------------------------------------------------------------------ 
     611      IF( ln_timing )   CALL timing_start('iceitd_reb') 
    608612      ! 
    609613      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
    … …  
    635639               jdonor(ji,jl)  = jl  
    636640               ! how much of a_i you send in cat sup is somewhat arbitrary 
    637                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    638                !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)   
    639                !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 
    640                !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 
    641                !!          zdaice(ji,jl)  = a_i_1d(ji) 
    642                !!          zdvice(ji,jl)  = v_i_1d(ji) 
    643                !!clem: these are from UCL and work ok 
    644                zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp 
    645                zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     641               ! these are from CICE => transfer everything 
     642               !!zdaice(ji,jl)  = a_i_1d(ji) 
     643               !!zdvice(ji,jl)  = v_i_1d(ji) 
     644               ! these are from LLN => transfer only half of the category 
     645               zdaice(ji,jl)  =                       0.5_wp  * a_i_1d(ji) 
     646               zdvice(ji,jl)  = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl) 
    646647            END DO 
    647648            ! 
    … …  
    686687      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 
    687688      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
     689      IF( ln_timing    )   CALL timing_stop ('iceitd_reb') 
    688690      ! 
    689691   END SUBROUTINE ice_itd_reb 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icesbc.F90

    r13472 r13957  
    6262      !!------------------------------------------------------------------- 
    6363      ! 
    64       IF( ln_timing )   CALL timing_start('ice_sbc') 
     64      IF( ln_timing )   CALL timing_start('icesbc') 
    6565      ! 
    6666      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    8989      ENDIF 
    9090      ! 
    91       IF( ln_timing )   CALL timing_stop('ice_sbc') 
     91      IF( ln_timing )   CALL timing_stop('icesbc') 
    9292      ! 
    9393   END SUBROUTINE ice_sbc_tau 
    … …  
    122122      !!-------------------------------------------------------------------- 
    123123      ! 
    124       IF( ln_timing )   CALL timing_start('ice_sbc_flx') 
     124      IF( ln_timing )   CALL timing_start('icesbc') 
    125125 
    126126      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    176176      ENDIF 
    177177      ! 
    178       IF( ln_timing )   CALL timing_stop('ice_sbc_flx') 
     178      IF( ln_timing )   CALL timing_stop('icesbc') 
    179179      ! 
    180180   END SUBROUTINE ice_sbc_flx 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icestp.F90

    r13908 r13957  
    121121      !!---------------------------------------------------------------------- 
    122122      ! 
    123       IF( ln_timing )   CALL timing_start('ice_stp') 
     123      IF( ln_timing )   CALL timing_start('icestp') 
    124124      ! 
    125125      !                                      !-----------------------! 
    … …  
    215215!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    216216      ! 
    217       IF( ln_timing )   CALL timing_stop('ice_stp') 
     217      IF( ln_timing )   CALL timing_stop('icestp') 
    218218      ! 
    219219   END SUBROUTINE ice_stp 
    … …  
    370370      v_i_b (:,:,:)   = v_i (:,:,:)     ! ice volume 
    371371      v_s_b (:,:,:)   = v_s (:,:,:)     ! snow volume 
     372      v_ip_b(:,:,:)   = v_ip(:,:,:)     ! pond volume 
     373      v_il_b(:,:,:)   = v_il(:,:,:)     ! pond lid volume 
    372374      sv_i_b(:,:,:)   = sv_i(:,:,:)     ! salt content 
    373375      e_s_b (:,:,:,:) = e_s (:,:,:,:)   ! snow thermal energy 
    … …  
    429431         diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp 
    430432         diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp 
     433         diag_aice(ji,jj) = 0._wp ;   diag_vpnd(ji,jj) = 0._wp 
    431434 
    432435         tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here) 
    … …  
    485488         diag_vsnw(:,:) = diag_vsnw(:,:) & 
    486489            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos 
     490         diag_vpnd(:,:) = diag_vpnd(:,:) & 
     491            &             + SUM(     v_ip + v_il          - v_ip_b - v_il_b                , dim=3 ) * r1_Dt_ice * rhow 
    487492         ! 
    488493         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd.F90

    r13908 r13957  
    259259      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    260260      !                    
    261       IF ( ln_pnd  )          CALL ice_thd_pnd                      ! --- Melt ponds  
     261      IF ( ln_pnd .AND. ln_icedH ) & 
     262         &                    CALL ice_thd_pnd                      ! --- Melt ponds  
    262263      ! 
    263264      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
    … …  
    267268                              CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
    268269      ! 
    269       oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation      
     270      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! ice natural aging incrementation      
    270271      ! 
    271272      ! convergence tests 
    … …  
    378379            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
    379380         END DO 
    380          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
    381          CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    382          CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) ) 
    383381         ! 
    384382         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            ) 
    … …  
    410408         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          ) 
    411409         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          ) 
    412          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          ) 
    413410         ! 
    414411         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          ) 
    … …  
    465462         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti) 
    466463         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 
    467          v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
    468          v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 
    469464         oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti) 
    470465          
    … …  
    484479            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
    485480         END DO 
    486          CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
    487          CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    488          CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,kl) ) 
    489481         ! 
    490482         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
    … …  
    502494         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        ) 
    503495         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        ) 
    504          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    505496         ! 
    506497         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        ) 
    … …  
    541532         CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 
    542533         CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 
    543          CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) ) 
    544          CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) ) 
    545534         CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) ) 
    546535         ! check convergence of heat diffusion scheme 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_pnd.F90

    r13931 r13957  
    3636   INTEGER ::              nice_pnd    ! choice of the type of pond scheme 
    3737   !                                   ! associated indices: 
    38    INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme 
    39    INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant ice pond scheme 
    40    INTEGER, PARAMETER ::   np_pndLEV = 2   ! Level ice pond scheme 
     38   INTEGER, PARAMETER ::   np_pndNO   = 0   ! No pond scheme 
     39   INTEGER, PARAMETER ::   np_pndCST  = 1   ! Constant ice pond scheme 
     40   INTEGER, PARAMETER ::   np_pndLEV  = 2   ! Level ice pond scheme 
    4141   INTEGER, PARAMETER ::   np_pndTOPO = 3   ! Level ice pond scheme 
    4242 
    43    REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds 
    44       zhi_min = 0.1_wp        , & ! minimum ice thickness with ponds (m)  
    45       zTd     = 0.15_wp       , & ! temperature difference for freeze-up (C) 
    46       zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
    47  
    48    !-------------------------------------------------------------------------- 
    49    !  
    50    ! Pond volume per area budget diags 
    51    !   
    52    ! The idea of diags is the volume of ponds per grid cell area is 
     43   !--------------------------------------------------------------------------  
     44   ! Diagnostics for pond volume per area 
    5345   ! 
    5446   ! dV/dt = mlt + drn + lid + rnf 
    … …  
    5951   ! 
    6052   ! In topo mode, the pond water lost because it is in the snow is not included in the budget 
     53   ! In level mode, all terms are incorporated 
    6154   ! 
    62    ! In level mode, all terms are incorporated 
    63  
    64    REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & ! pond volume budget diagnostics 
    65       diag_dvpn_mlt,  &     !! meltwater pond volume input      [m/day] 
    66       diag_dvpn_drn,  &     !! pond volume lost by drainage     [m/day] 
    67       diag_dvpn_lid,  &     !! exchange with lid / refreezing   [m/day] 
    68       diag_dvpn_rnf         !! meltwater pond lost to runoff    [m/day] 
    69        
    70    REAL(wp), ALLOCATABLE, DIMENSION(:) ::   & ! pond volume budget diagnostics (1d) 
    71       diag_dvpn_mlt_1d,  &  !! meltwater pond volume input      [m/day] 
    72       diag_dvpn_drn_1d,  &  !! pond volume lost by drainage     [m/day] 
    73       diag_dvpn_lid_1d,  &  !! exchange with lid / refreezing   [m/day] 
    74       diag_dvpn_rnf_1d      !! meltwater pond lost to runoff    [m/day] 
     55   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_mlt       ! meltwater pond volume input      [kg/m2/s] 
     56   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_drn       ! pond volume lost by drainage     [-] 
     57   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_lid       ! exchange with lid / refreezing   [-] 
     58   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   diag_dvpn_rnf       ! meltwater pond lost to runoff    [-]       
     59   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_mlt_1d    ! meltwater pond volume input      [-] 
     60   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_drn_1d    ! pond volume lost by drainage     [-] 
     61   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_lid_1d    ! exchange with lid / refreezing   [-] 
     62   REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   diag_dvpn_rnf_1d    ! meltwater pond lost to runoff    [-] 
    7563 
    7664   !! * Substitutions 
    … …  
    9078      !! ** Purpose :   change melt pond fraction and thickness 
    9179      !! 
    92       !! Note: Melt ponds affect only radiative transfer for now 
    93       !! 
    94       !!       No heat, no salt. 
    95       !!       The melt water they carry is collected but  
    96       !!       not removed from fw budget or released to the ocean 
    97       !! 
    98       !!       A wfx_pnd has been coded for diagnostic purposes 
    99       !!       It is not fully consistent yet. 
    100       !! 
    101       !!       The current diagnostic lacks a contribution from drainage 
    102       !! 
     80      !! ** Note    : Melt ponds affect only radiative transfer for now 
     81      !!              No heat, no salt. 
     82      !!              The current diagnostics lacks a contribution from drainage 
    10383      !!------------------------------------------------------------------- 
    104       !! 
     84      INTEGER ::   ji, jj, jl        ! loop indices 
     85      !!------------------------------------------------------------------- 
    10586       
    10687      ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 
    10788      ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 
    108  
    109       diag_dvpn_mlt(:,:) = 0._wp    ;   diag_dvpn_drn(:,:) = 0._wp 
    110       diag_dvpn_lid(:,:) = 0._wp    ;   diag_dvpn_rnf(:,:) = 0._wp 
     89      ! 
     90      diag_dvpn_mlt (:,:) = 0._wp   ;   diag_dvpn_drn (:,:) = 0._wp 
     91      diag_dvpn_lid (:,:) = 0._wp   ;   diag_dvpn_rnf (:,:) = 0._wp 
    11192      diag_dvpn_mlt_1d(:) = 0._wp   ;   diag_dvpn_drn_1d(:) = 0._wp 
    11293      diag_dvpn_lid_1d(:) = 0._wp   ;   diag_dvpn_rnf_1d(:) = 0._wp 
    11394 
    114       SELECT CASE ( nice_pnd ) 
     95      !------------------------------------- 
     96      !  Remove ponds where ice has vanished 
     97      !------------------------------------- 
     98      at_i(:,:) = SUM( a_i, dim=3 ) 
    11599      ! 
    116       CASE (np_pndCST)   ;   CALL pnd_CST                              !==  Constant melt ponds  ==! 
    117          ! 
    118       CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==! 
    119          ! 
    120       CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==! 
    121          ! 
    122       END SELECT 
     100      DO jl = 1, jpl 
     101         DO_2D( 1, 1, 1, 1 ) 
     102            IF( v_i(ji,jj,jl) < epsi10 .OR. at_i(ji,jj) < epsi10 ) THEN 
     103               wfx_pnd  (ji,jj)    = wfx_pnd(ji,jj) + ( v_ip(ji,jj,jl) + v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 
     104               a_ip     (ji,jj,jl) = 0._wp 
     105               v_ip     (ji,jj,jl) = 0._wp 
     106               v_il     (ji,jj,jl) = 0._wp 
     107               h_ip     (ji,jj,jl) = 0._wp 
     108               h_il     (ji,jj,jl) = 0._wp 
     109               a_ip_frac(ji,jj,jl) = 0._wp 
     110            ENDIF 
     111         END_2D 
     112      END DO 
     113       
     114      !------------------------------ 
     115      !  Identify grid cells with ice 
     116      !------------------------------ 
     117      npti = 0   ;   nptidx(:) = 0 
     118      DO_2D( 1, 1, 1, 1 ) 
     119         IF( at_i(ji,jj) >= epsi10 ) THEN 
     120            npti = npti + 1 
     121            nptidx( npti ) = (jj - 1) * jpi + ji 
     122         ENDIF 
     123      END_2D 
     124 
     125      !------------------------------------ 
     126      !  Select melt pond scheme to be used 
     127      !------------------------------------ 
     128      IF( npti > 0 ) THEN 
     129         SELECT CASE ( nice_pnd ) 
     130            ! 
     131         CASE (np_pndCST)   ;   CALL pnd_CST                              !==  Constant melt ponds  ==! 
     132            ! 
     133         CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==! 
     134            ! 
     135         CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==! 
     136            ! 
     137         END SELECT 
     138      ENDIF 
     139 
     140      !------------------------------------ 
     141      !  Diagnostics 
     142      !------------------------------------ 
     143      CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting 
     144      CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid 
     145      CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage 
     146      CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow 
    123147      ! 
    124  
    125       IF( iom_use('dvpn_mlt'  ) )   CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting 
    126       IF( iom_use('dvpn_lid'  ) )   CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid 
    127       IF( iom_use('dvpn_drn'  ) )   CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage 
    128       IF( iom_use('dvpn_rnf'  ) )   CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow 
    129  
    130       DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf ) 
     148      DEALLOCATE( diag_dvpn_mlt   , diag_dvpn_lid   , diag_dvpn_drn   , diag_dvpn_rnf    ) 
    131149      DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 
    132150       
    133151   END SUBROUTINE ice_thd_pnd  
    134  
    135152 
    136153 
    … …  
    151168      !! ** References : Bush, G.W., and Trump, D.J. (2017) 
    152169      !!------------------------------------------------------------------- 
    153       INTEGER  ::   ji        ! loop indices 
     170      INTEGER  ::   ji, jl    ! loop indices 
     171      REAL(wp) ::   zdv_pnd   ! Amount of water going into the ponds & lids 
    154172      !!------------------------------------------------------------------- 
    155       DO ji = 1, npti 
    156          ! 
    157          IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 
    158             h_ip_1d(ji)      = rn_hpnd     
    159             a_ip_1d(ji)      = rn_apnd * a_i_1d(ji) 
    160             h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever 
    161          ELSE 
    162             h_ip_1d(ji)      = 0._wp     
    163             a_ip_1d(ji)      = 0._wp 
    164             h_il_1d(ji)      = 0._wp 
    165          ENDIF 
    166          ! 
     173      DO jl = 1, jpl 
     174          
     175         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d    (1:npti), a_i    (:,:,jl) ) 
     176         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d   (1:npti), t_su   (:,:,jl) ) 
     177         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d   (1:npti), a_ip   (:,:,jl) ) 
     178         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d   (1:npti), h_ip   (:,:,jl) ) 
     179         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d   (1:npti), h_il   (:,:,jl) ) 
     180         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:)    ) 
     181 
     182         DO ji = 1, npti 
     183            ! 
     184            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 
     185            ! 
     186            IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 
     187               h_ip_1d(ji)      = rn_hpnd     
     188               a_ip_1d(ji)      = rn_apnd * a_i_1d(ji) 
     189               h_il_1d(ji)      = 0._wp    ! no pond lids whatsoever 
     190            ELSE 
     191               h_ip_1d(ji)      = 0._wp     
     192               a_ip_1d(ji)      = 0._wp 
     193               h_il_1d(ji)      = 0._wp 
     194            ENDIF 
     195            ! 
     196            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     197            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     198            ! 
     199            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 
     200            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 
     201            ! 
     202         END DO 
     203 
     204         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d   (1:npti), a_ip   (:,:,jl) ) 
     205         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d   (1:npti), h_ip   (:,:,jl) ) 
     206         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d   (1:npti), h_il   (:,:,jl) ) 
     207         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d   (1:npti), v_ip   (:,:,jl) ) 
     208         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d   (1:npti), v_il   (:,:,jl) ) 
     209         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:)    ) 
     210 
    167211      END DO 
    168212      ! 
    … …  
    203247      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 --- 
    204248      !! 
    205       !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 --- 
    206       !!                                     perm = permability of sea-ice 
     249      !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi * flush         --- from Flocco et al 2007 --- 
     250      !!                                     perm = permability of sea-ice                              + correction from Hunke et al 2012 (flush) 
    207251      !!                                     visc = water viscosity 
    208252      !!                                     Hp   = height of top of the pond above sea-level 
    209253      !!                                     Hi   = ice thickness thru which there is flushing 
     254      !!                                     flush= correction otherwise flushing is excessive 
    210255      !! 
    211256      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness 
    … …  
    218263      !! ** Note       :   Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds.  
    219264      !! 
    220       !! ** References :   Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 
    221       !!                    
    222       !!------------------------------------------------------------------- 
    223        
     265      !! ** References :   Flocco and Feltham (JGR, 2007) 
     266      !!                   Flocco et al       (JGR, 2010) 
     267      !!                   Holland et al      (J. Clim, 2012) 
     268      !!                   Hunke et al        (OM 2012) 
     269      !!-------------------------------------------------------------------   
    224270      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
    225271      !! 
    … …  
    228274      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity 
    229275      !! 
    230       REAL(wp) ::   zfr_mlt, zdv_mlt                  ! fraction and volume of available meltwater retained for melt ponding 
     276      REAL(wp) ::   zfr_mlt, zdv_mlt, zdv_avail       ! fraction and volume of available meltwater retained for melt ponding 
    231277      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes 
     278      REAL(wp) ::   zdv_pnd                           ! Amount of water going into the ponds & lids 
    232279      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh 
    233280      REAL(wp) ::   zv_ip_max                         ! max pond volume allowed 
    234281      REAL(wp) ::   zdT                               ! zTp-t_su 
    235       REAL(wp) ::   zsbr                              ! Brine salinity 
     282      REAL(wp) ::   zsbr, ztmelts                     ! Brine salinity 
    236283      REAL(wp) ::   zperm                             ! permeability of sea ice 
    237284      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
    238285      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse 
    239       REAL(wp) ::   zvold                             ! 
    240       !! 
    241       INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
    242        
     286      !! 
     287      INTEGER  ::   ji, jk, jl                        ! loop indices 
    243288      !!------------------------------------------------------------------- 
    244        
    245289      z1_rhow   = 1._wp / rhow  
    246290      z1_aspect = 1._wp / zaspect 
    247291      z1_Tp     = 1._wp / zTp  
    248292       
    249       !----------------------------------------------------------------------------------------------- 
    250       !  Identify grid cells with ice 
    251       !----------------------------------------------------------------------------------------------- 
    252       at_i(:,:) = SUM( a_i, dim=3 ) 
     293      CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d          (1:npti), at_i          ) 
     294      CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d       (1:npti), wfx_pnd       ) 
     295       
     296      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 
     297      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 
     298      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 
     299      CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 
     300 
     301      DO jl = 1, jpl 
     302 
     303         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i (:,:,jl) ) 
     304         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d  (1:npti), h_i (:,:,jl) ) 
     305         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su(:,:,jl) ) 
     306         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,jl) ) 
     307         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip(:,:,jl) ) 
     308         CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il(:,:,jl) ) 
     309 
     310         CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum(1:npti), dh_i_sum_2d(:,:,jl) ) 
     311         CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt(1:npti), dh_s_mlt_2d(:,:,jl) ) 
     312 
     313         DO jk = 1, nlay_i 
     314            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 
     315            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,jl) ) 
     316         END DO 
     317          
     318         !----------------------- 
     319         ! Melt pond calculations 
     320         !----------------------- 
     321         DO ji = 1, npti 
     322            ! 
     323            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 
     324            !                                                            !----------------------------------------------------! 
     325            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 ! 
     326               !                                                         !----------------------------------------------------! 
     327               !--- Remove ponds on thin ice or tiny ice fractions 
     328               a_ip_1d(ji) = 0._wp 
     329               h_ip_1d(ji) = 0._wp 
     330               h_il_1d(ji) = 0._wp 
     331               !                                                         !--------------------------------! 
     332            ELSE                                                         ! Case ice thickness >= rn_himin ! 
     333               !                                                         !--------------------------------! 
     334               v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
     335               v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     336               ! 
     337               !------------------! 
     338               ! case ice melting ! 
     339               !------------------! 
     340               ! 
     341               !--- available meltwater for melt ponding (zdv_avail) ---! 
     342               zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 
     343               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 
     344               zdv_mlt   = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors?  
     345               ! 
     346               !--- overflow ---! 
     347               ! 
     348               ! area driven overflow 
     349               !    If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
     350               !       a_ip_max = zfr_mlt * a_i 
     351               !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     352               zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
     353               zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     354 
     355               ! depth driven overflow 
     356               !    If pond depth exceeds half the ice thickness then reduce the pond volume 
     357               !       h_ip_max = 0.5 * h_i 
     358               !       => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     359               zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
     360               zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     361 
     362               !--- Pond growing ---! 
     363               v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     364               ! 
     365               !--- Lid melting ---! 
     366               IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
     367               ! 
     368               !-------------------! 
     369               ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
     370               !-------------------! 
     371               ! 
     372               zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
     373               ! 
     374               !--- Pond contraction (due to refreezing) ---! 
     375               IF( ln_pnd_lids ) THEN 
     376                  ! 
     377                  !--- Lid growing and subsequent pond shrinking ---!  
     378                  zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
     379                     &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
     380 
     381                  ! Lid growing 
     382                  v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) 
     383 
     384                  ! Pond shrinking 
     385                  v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 
     386 
     387               ELSE 
     388                  zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp )  ! Holland 2012 (eq. 6)  
     389                  ! Pond shrinking 
     390                  v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 
     391               ENDIF 
     392               ! 
     393               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     394               ! v_ip     = h_ip * a_ip 
     395               ! 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) 
     396               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 
     397               h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     398               ! 
     399 
     400               !------------------------------------------------!             
     401               ! Pond drainage through brine network (flushing) ! 
     402               !------------------------------------------------! 
     403               ! height of top of the pond above sea-level 
     404               zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 
     405 
     406               ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
     407               DO jk = 1, nlay_i 
     408                  ! MV Assur is inconsistent with SI3 
     409                  !!zsbr = - 1.2_wp                                  & 
     410                  !!   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
     411                  !!   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
     412                  !!   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
     413                  !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     414                  ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 
     415                  ztmelts  = -rTmlt * sz_i_1d(ji,jk) 
     416                  ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 ) 
     417               END DO 
     418               zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
     419 
     420               ! Do the drainage using Darcy's law 
     421               zdv_flush   = -zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! zflush comes from Hunke et al. (2012) 
     422               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0  
     423               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
     424 
     425               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     426               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 
     427               h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     428 
     429               !--- Corrections and lid thickness ---! 
     430               IF( ln_pnd_lids ) THEN 
     431                  !--- retrieve lid thickness from volume ---! 
     432                  IF( a_ip_1d(ji) > 0.01_wp ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
     433                  ELSE                               ;   h_il_1d(ji) = 0._wp 
     434                  ENDIF 
     435                  !--- remove ponds if lids are much larger than ponds ---! 
     436                  IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     437                     a_ip_1d(ji) = 0._wp 
     438                     h_ip_1d(ji) = 0._wp 
     439                     h_il_1d(ji) = 0._wp 
     440                  ENDIF 
     441               ENDIF 
     442 
     443               ! diagnostics: dvpnd = mlt+rnf+lid+drn 
     444               diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow *   zdv_avail             * r1_Dt_ice   ! > 0, surface melt input 
     445               diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_Dt_ice   ! < 0, runoff 
     446               diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + rhow *   zdv_frz               * r1_Dt_ice   ! < 0, shrinking 
     447               diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow *   zdv_flush             * r1_Dt_ice   ! < 0, drainage 
     448               ! 
     449            ENDIF 
     450            ! 
     451            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     452            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     453            ! 
     454            zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 
     455            wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 
     456            ! 
     457         END DO 
     458 
     459         !-------------------------------------------------------------------- 
     460         ! Retrieve 2D arrays 
     461         !-------------------------------------------------------------------- 
     462         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,jl) ) 
     463         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d(1:npti), h_ip(:,:,jl) ) 
     464         CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d(1:npti), h_il(:,:,jl) ) 
     465         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,jl) ) 
     466         CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,jl) ) 
     467         ! 
     468      END DO 
    253469      ! 
    254       npti = 0   ;   nptidx(:) = 0 
    255       DO_2D( 1, 1, 1, 1 ) 
    256          IF ( at_i(ji,jj) > epsi10 ) THEN 
    257             npti = npti + 1 
    258             nptidx( npti ) = (jj - 1) * jpi + ji 
    259          ENDIF 
    260       END_2D 
    261        
    262       !----------------------------------------------------------------------------------------------- 
    263       ! Prepare 1D arrays 
    264       !----------------------------------------------------------------------------------------------- 
    265  
    266       IF( npti > 0 ) THEN 
    267        
    268          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum      ) 
    269          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti)   , wfx_sum          ) 
    270          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti)   , wfx_pnd          ) 
    271  
    272          CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d (1:npti)      , at_i             ) 
    273           
    274          CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 
    275          CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 
    276          CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 
    277          CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 
    278        
    279          DO jl = 1, jpl 
    280           
    281             CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d      (1:npti), a_i      (:,:,jl) ) 
    282             CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d      (1:npti), h_i      (:,:,jl) ) 
    283             CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d     (1:npti), t_su     (:,:,jl) ) 
    284             CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
    285             CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
    286             CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
    287  
    288             CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum    (1:npti), dh_i_sum_2d     (:,:,jl) ) 
    289             CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt    (1:npti), dh_s_mlt_2d     (:,:,jl) ) 
    290              
    291             DO jk = 1, nlay_i 
    292                CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  ) 
    293             END DO 
    294              
    295             !----------------------------------------------------------------------------------------------- 
    296             ! Go for ponds 
    297             !----------------------------------------------------------------------------------------------- 
    298  
    299             DO ji = 1, npti 
    300                !                                                            !----------------------------------------------------! 
    301                IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
    302                   !                                                         !----------------------------------------------------! 
    303                   !--- Remove ponds on thin ice or tiny ice fractions 
    304                   a_ip_1d(ji)      = 0._wp 
    305                   h_ip_1d(ji)      = 0._wp 
    306                   h_il_1d(ji)      = 0._wp 
    307                   !                                                         !--------------------------------! 
    308                ELSE                                                         ! Case ice thickness >= rn_himin ! 
    309                   !                                                         !--------------------------------! 
    310                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
    311                   v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
    312                   ! 
    313                   !------------------! 
    314                   ! case ice melting ! 
    315                   !------------------! 
    316                   ! 
    317                   !--- available meltwater for melt ponding ---! 
    318                   zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    319                   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 
    320                   zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
    321                    
    322                   diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_Dt_ice                      ! surface melt input diag 
    323                   diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_Dt_ice   ! runoff diag 
    324                   ! 
    325                   !--- overflow ---! 
    326                   ! 
    327                   ! 1) area driven overflow 
    328                   ! 
    329                   ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 
    330                   !    a_ip_max = zfr_mlt * a_i 
    331                   !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    332                   zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    333                   zvold     = zdv_mlt 
    334                   zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )                      
    335                    
    336                   !  
    337                   ! 2) depth driven overflow 
    338                   ! 
    339                   ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    340                   !    h_ip_max = 0.5 * h_i 
    341                   !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    342                   zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear 
    343                   zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    344                   diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_Dt_ice       ! runoff diag - overflow contribution 
    345              
    346                   !--- Pond growing ---! 
    347                   v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
    348                   ! 
    349                   !--- Lid melting ---! 
    350                   IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
    351                   ! 
    352                   !--- mass flux ---! 
    353                   ! MV I would recommend to remove that 
    354                   ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 
    355                    
    356                   IF( zdv_mlt > 0._wp ) THEN 
    357                      zfac = zdv_mlt * rhow * r1_Dt_ice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
    358                      wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    359                      ! 
    360                      zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
    361                      wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    362                      wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    363                   ENDIF 
    364  
    365                   !-------------------! 
    366                   ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
    367                   !-------------------! 
    368                   ! 
    369                   zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
    370                   ! 
    371                   !--- Pond contraction (due to refreezing) ---! 
    372                   zvold       = v_ip_1d(ji) ! for diag 
    373  
    374                   IF( ln_pnd_lids ) THEN 
    375                      ! 
    376                      !--- Lid growing and subsequent pond shrinking ---!  
    377                      zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
    378                         &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
    379                 
    380                      ! Lid growing 
    381                      v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
    382                 
    383                      ! Pond shrinking 
    384                      v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
    385                       
    386                   ELSE 
    387                    
    388                      ! Pond shrinking 
    389                      v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
    390                       
    391                   ENDIF 
    392  
    393                   diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_Dt_ice   ! shrinking counted as lid diagnostic 
    394  
    395                   ! 
    396                   !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    397                   ! v_ip     = h_ip * a_ip 
    398                   ! 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) 
    399                   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 
    400                   h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    401  
    402                   !------------------------------------------------!             
    403                   ! Pond drainage through brine network (flushing) ! 
    404                   !------------------------------------------------! 
    405                   ! height of top of the pond above sea-level 
    406                   zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 
    407              
    408                   ! Calculate the permeability of the ice  
    409                   ! mv- note here that a linear expression for permeability is fine,  
    410                   ! simpler and more consistent with the rest of SI3 code 
    411                   DO jk = 1, nlay_i 
    412                      ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
    413                   END DO 
    414                   zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
    415              
    416                   ! Do the drainage using Darcy's law 
    417                   zdv_flush   = - zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 
    418                   zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
    419                   ! zdv_flush   = 0._wp ! MV remove pond drainage for now 
    420                   v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
    421                    
    422                   diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_Dt_ice   ! shrinking counted as lid diagnostic 
    423              
    424                   ! MV --- why pond drainage does not give back water into freshwater flux ?             
    425                   !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    426                   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 
    427                   h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    428  
    429                   !--- Corrections and lid thickness ---! 
    430                   IF( ln_pnd_lids ) THEN 
    431                      !--- retrieve lid thickness from volume ---! 
    432                      IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
    433                      ELSE                              ;   h_il_1d(ji) = 0._wp 
    434                      ENDIF 
    435                      !--- remove ponds if lids are much larger than ponds ---! 
    436                      IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
    437                         a_ip_1d(ji)      = 0._wp 
    438                         h_ip_1d(ji)      = 0._wp 
    439                         h_il_1d(ji)      = 0._wp 
    440                      ENDIF 
    441                   ENDIF 
    442                   ! 
    443                ENDIF 
    444                 
    445             END DO ! ji 
    446  
    447             !----------------------------------------------------------------------------------------------- 
    448             ! Retrieve 2D arrays 
    449             !----------------------------------------------------------------------------------------------- 
    450              
    451             v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
    452             v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 
    453              
    454             CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
    455             CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
    456             CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
    457             CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d     (1:npti), v_ip     (:,:,jl) ) 
    458             CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d     (1:npti), v_il     (:,:,jl) ) 
    459     
    460          END DO ! ji 
    461                    
    462          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 
    463          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        ) 
    464          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    465           
    466          CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        ) 
    467          CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        ) 
    468          CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        ) 
    469          CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        )                   
    470                             
     470      CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd ) 
    471471      ! 
    472       ENDIF 
    473        
     472      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        ) 
     473      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        ) 
     474      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        ) 
     475      CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        ) 
     476      ! 
    474477   END SUBROUTINE pnd_LEV 
    475478 
    … …  
    507510      !! 
    508511      !!------------------------------------------------------------------- 
     512      REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds 
     513         zTd     = 0.15_wp       , & ! temperature difference for freeze-up (C) 
     514         zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
     515 
    509516 
    510517      ! local variables 
    … …  
    536543      ! a_ip      -> apond 
    537544      ! a_ip_frac -> apnd 
     545 
     546      CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) 
    538547       
    539548      !--------------------------------------------------------------- 
    … …  
    628637         DO_2D( 1, 1, 1, 1 ) 
    629638 
    630             IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
     639            IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
    631640                   
    632641               !-------------------------- 
    … …  
    774783 
    775784      END DO   ! jl 
     785 
    776786 
    777787   END SUBROUTINE pnd_TOPO 
    … …  
    848858          viscosity = 1.79e-3_wp     ! kinematic water viscosity in kg/m/s 
    849859 
    850        INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     860      REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds 
     861         zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
     862 
     863      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
    851864 
    852865       a_ip(:,:,:) = 0._wp 
    … …  
    855868       DO_2D( 1, 1, 1, 1 ) 
    856869  
    857              IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
     870             IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
    858871  
    859872        !------------------------------------------------------------------- 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/iceupdate.F90

    r13643 r13957  
    9494      REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace 
    9595      !!--------------------------------------------------------------------- 
    96       IF( ln_timing )   CALL timing_start('ice_update') 
     96      IF( ln_timing )   CALL timing_start('iceupdate') 
    9797 
    9898      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    154154         ! ice-ocean  mass flux 
    155155         wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    156             &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 
     156            &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
    157157          
    158158         ! snw-ocean mass flux 
    … …  
    160160          
    161161         ! total mass flux at the ocean/ice interface 
    162          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 
    163          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 
     162         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 
     163         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 
    164164 
    165165         ! Salt flux at the ocean surface       
    … …  
    172172         snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step 
    173173         !                                               ! new mass per unit area 
    174          snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  
     174         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)) )  
    175175         !                                               ! time evolution of snow+ice mass 
    176176         snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice 
    … …  
    286286      IF( ln_icectl         )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 
    287287      IF( sn_cfctl%l_prtctl )   CALL ice_prt3D     ('iceupdate')                                       ! prints 
    288       IF( ln_timing         )   CALL timing_stop   ('ice_update')                                      ! timing 
     288      IF( ln_timing         )   CALL timing_stop   ('iceupdate')                                       ! timing 
    289289      ! 
    290290   END SUBROUTINE ice_update_flx 
    … …  
    324324      REAL(wp) ::   zflagi                          !   -      - 
    325325      !!--------------------------------------------------------------------- 
    326       IF( ln_timing )   CALL timing_start('ice_update_tau') 
     326      IF( ln_timing )   CALL timing_start('ice_update') 
    327327 
    328328      IF( kt == nit000 .AND. lwp ) THEN 
    … …  
    376376      CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition 
    377377      ! 
    378       IF( ln_timing )   CALL timing_stop('ice_update_tau') 
     378      IF( ln_timing )   CALL timing_stop('ice_update') 
    379379      !   
    380380   END SUBROUTINE ice_update_tau 

  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icevar.F90

    r13911 r13957  
    236236      z1_zhmax =  1._wp / hi_max(jpl)                
    237237      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 
    238          h_i  (:,:,jpl) = zhmax 
     238         h_i   (:,:,jpl) = zhmax 
    239239         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax  
    240240         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) 
    … …  
    252252      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow 
    253253      ELSEWHERE                            ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * &   ! lid is in between. Expose part of the pond 
    254          &                                                       ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min ) 
     254         &                                                       ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 
    255255      END WHERE 
    256256      ! 
    … …  
    534534         DO_2D( 1, 1, 1, 1 ) 
    535535            ! update exchanges with ocean 
    536             sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice 
    537             wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_Dt_ice 
    538             wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_Dt_ice 
     536            sfx_res(ji,jj)  = sfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice 
     537            wfx_res(ji,jj)  = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_Dt_ice 
     538            wfx_res(ji,jj)  = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_Dt_ice 
     539            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_Dt_ice 
    539540            ! 
    540541            a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 
    … …  
    551552            v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
    552553            v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 
     554            h_ip (ji,jj,jl) = h_ip (ji,jj,jl) * zswitch(ji,jj) 
     555            h_il (ji,jj,jl) = h_il (ji,jj,jl) * zswitch(ji,jj) 
    553556            ! 
    554557         END_2D 
    555558         ! 
    556          !----------------------------------------------------------------- 
    557          ! zap small ponds 
    558          !----------------------------------------------------------------- 
    559          IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    560             DO_2D( 1, 1, 1, 1 ) 
    561                IF ( v_ip(ji,jj,jl) <= epsi10 ) THEN 
    562                   a_ip(ji,jj,jl)      = 0._wp 
    563                   a_ip_frac(ji,jj,jl) = 0._wp 
    564                   v_ip(ji,jj,jl)      = 0._wp 
    565                   h_ip(ji,jj,jl)      = 0._wp 
    566                   v_il(ji,jj,jl)      = 0._wp 
    567                   h_il(ji,jj,jl)      = 0._wp 
    568                ENDIF 
    569             END_2D 
    570          ENDIF 
    571559      END DO  
    572560 
    … …  
    650638               psv_i  (ji,jj,jl) = 0._wp 
    651639            ENDIF 
     640            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 
     641               wfx_pnd(ji,jj)    = wfx_pnd(ji,jj) + pv_il(ji,jj,jl) * rhow * z1_dt 
     642               pv_il  (ji,jj,jl) = 0._wp 
     643            ENDIF 
     644            IF( pv_ip(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 
     645               wfx_pnd(ji,jj)    = wfx_pnd(ji,jj) + pv_ip(ji,jj,jl) * rhow * z1_dt 
     646               pv_ip  (ji,jj,jl) = 0._wp 
     647            ENDIF 
    652648         END_2D 
    653649         ! 
    … …  
    658654      WHERE( pa_i  (:,:,:) < 0._wp )   pa_i  (:,:,:) = 0._wp 
    659655      WHERE( pa_ip (:,:,:) < 0._wp )   pa_ip (:,:,:) = 0._wp 
    660       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    661       WHERE( pv_il (:,:,:) < 0._wp )   pv_il (:,:,:) = 0._wp !    but it does not change conservation, so keep it this way is ok 
    662656      ! 
    663657   END SUBROUTINE ice_var_zapneg 
Note: See TracChangeset for help on using the changeset viewer.