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 8597 – NEMO

Changeset 8597


Ignore:
Timestamp:
2017-10-05T16:44:46+02:00 (7 years ago)
Author:
clem
Message:

first step to make melt ponds compliant with the new code

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml

    r8562 r8597  
    132132         <field id="icevmp"       long_name="melt pond volume"                                             standard_name="sea_ice_meltpond_volume"                            unit="m"            />  
    133133 
    134          <field id="iceconc_cat"  long_name="Sea-ice area fractions in thickness categories"               unit=""   grid_ref="grid_T_3D_ncatice" /> 
     134         <field id="iceconc_cat"  long_name="Sea-ice concentration in thickness categories"                unit=""   grid_ref="grid_T_3D_ncatice" /> 
    135135         <field id="icethic_cat"  long_name="Sea-ice thickness in thickness categories"                    unit="m"  grid_ref="grid_T_3D_ncatice" /> 
    136136         <field id="snowthic_cat" long_name="Snow thickness in thickness categories"                       unit="m"  grid_ref="grid_T_3D_ncatice" /> 
    … …  
    140140         <field id="icetemp_cat"  long_name="Ice temperature for categories"                               unit="degC"   grid_ref="grid_T_3D_ncatice" /> 
    141141         <field id="snwtemp_cat"  long_name="Snow temperature for categories"                              unit="degC"   grid_ref="grid_T_3D_ncatice" /> 
    142          <field id="iceamp_cat"   long_name="Ice melt pond fraction for categories"                        unit="%"      grid_ref="grid_T_3D_ncatice" />  
     142         <field id="iceamp_cat"   long_name="Ice melt pond concentration for categories"                   unit="%"      grid_ref="grid_T_3D_ncatice" />  
    143143         <field id="icevmp_cat"   long_name="Ice melt pond volume for categories"                          unit="m"      grid_ref="grid_T_3D_ncatice" />  
     144         <field id="icehmp_cat"   long_name="Ice melt pond thickness for categories"                       unit="m"      grid_ref="grid_T_3D_ncatice" />  
     145         <field id="iceafp_cat"   long_name="Ice melt pond fraction for categories"                        unit="m"      grid_ref="grid_T_3D_ncatice" />  
    144146 
    145147         <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     /> 
    … …  
    186188         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="kg/m2/s"   /> 
    187189         <field id="vfxthin"      long_name="thermo ice prod. for thin ice(20cm) + open water"             unit="kg/m2/s"   /> 
     190         <field id="vfxpnd"       long_name="melt pond water flux to ocean"                                unit="kg/m2/s"   /> 
    188191 
    189192         <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"   /> 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r8592 r8597  
    125125   ln_icedA         = .true.          !  activate lateral melting param. (T) or not (F) 
    126126   ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F) 
    127    ln_icedS         = .true.          !  activate gravity drainage and flushing (T) or not (F) 
     127   ln_icedS         = .true.          !  activate brine drainage (T) or not (F) 
    128128/ 
    129129!------------------------------------------------------------------------------ 
    … …  
    176176&nammp      !   Melt ponds 
    177177!------------------------------------------------------------------------------ 
    178    ln_pnd           = .false.         !  active melt ponds 
    179    ln_pnd_rad       = .false.         !  active melt ponds radiative coupling 
    180    ln_pnd_fw        = .false.         !  active melt ponds freshwater coupling 
    181    nn_pnd_scheme    =   0             !  type of melt pond scheme  : =0 prescribed ( Tsu=0 ), =1 empirical, =2 topographic 
    182    rn_apnd          =   0.2           !  prescribed pond fraction, at Tsu=0  : (0<rn_apnd<1, nn_pnd_scheme = 0) 
    183    rn_hpnd          =   0.05          !  prescribed pond depth, at Tsu=0     : (0<rn_apnd<1, nn_pnd_scheme = 0) 
     178   ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
     179      ln_pnd_fwb    = .false.         !     melt ponds store freshwater or not 
     180   ln_pnd_CST       = .false.         !  activate constant melt ponds 
     181      rn_apnd       =   0.2           !     prescribed pond fraction, at Tsu=0 
     182      rn_hpnd       =   0.05          !     prescribed pond depth, at Tsu=0 
     183   ln_pnd_alb       = .false.         !  melt ponds affect albedo or not 
    184184/ 
    185185!------------------------------------------------------------------------------ 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8569 r8597  
    204204   REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU] 
    205205 
    206    ! MV MP 2016 
    207    !                                     !!** melt pond namelist (nammp) 
    208    LOGICAL , PUBLIC ::   ln_pnd           !: activate ponds or not 
    209    LOGICAL , PUBLIC ::   ln_pnd_rad       !: ponds radiatively active or not 
    210    LOGICAL , PUBLIC ::   ln_pnd_fw        !: ponds active wrt meltwater or not 
    211    INTEGER , PUBLIC ::   nn_pnd_scheme    !: type of melt pond scheme:   =0 prescribed, =1 empirical, =2 topographic 
    212    REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1), only if nn_pnd_scheme = 0 
    213    REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1), only if nn_pnd_scheme = 0 
    214    ! END MV MP 2016 
     206   !                                     !!** namelist namthd_pnd 
     207   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012 
     208   LOGICAL , PUBLIC ::   ln_pnd_fwb       !: melt ponds store freshwater 
     209   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
     210   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1) 
     211   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1) 
     212   LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo 
     213 
    215214   !                                     !!** ice-diagnostics namelist (namdia) ** 
    216215   LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90

    r8592 r8597  
    4949CONTAINS 
    5050 
    51    SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, pafrac_pnd, ph_pnd, palb_cs, palb_os ) 
     51   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os ) 
    5252      !!---------------------------------------------------------------------- 
    5353      !!               ***  ROUTINE ice_alb  *** 
    … …  
    9898      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice       !  sea-ice thickness 
    9999      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw       !  snow depth 
     100      LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo 
    100101      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area) 
    101102      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth 
    … …  
    125126               IF( ph_snw(ji,jj,jl) == 0._wp ) THEN 
    126127                  zafrac_snw = 0._wp 
    127                   zafrac_pnd = pafrac_pnd(ji,jj,jl) 
     128                  IF( ld_pnd_alb ) THEN 
     129                     zafrac_pnd = pafrac_pnd(ji,jj,jl) 
     130                  ELSE 
     131                     zafrac_pnd = 0._wp 
     132                  ENDIF 
    128133                  zafrac_ice = 1._wp - zafrac_pnd 
    129134               ELSE 
    … …  
    134139               ! 
    135140               !                       !--- Bare ice albedo (for hi > 150cm) 
    136                IF( zafrac_pnd > 0._wp ) THEN 
     141               IF( ld_pnd_alb ) THEN 
    137142                  zalb_ice = rn_alb_idry 
    138143               ELSE 
    … …  
    154159               ENDIF 
    155160               !                       !--- Ponded ice albedo 
    156                IF( zafrac_pnd > 0._wp ) THEN 
     161               IF( ld_pnd_alb ) THEN 
    157162                  zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd )  
    158163               ELSE 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90

    r8563 r8597  
    172172            END DO 
    173173         END DO 
    174       END DO 
    175              
    176       IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i 
    177          WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    178       ENDIF 
     174      END DO  
     175      !                                           !-- correct pond fraction to avoid a_ip > a_i 
     176      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    179177      ! 
    180178   END SUBROUTINE Hbig 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90

    r8569 r8597  
    120120            z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    121121         END DO 
    122          IF ( nn_pnd_scheme > 0 ) THEN 
     122         IF ( ln_pnd_H12 ) THEN 
    123123            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    124124            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
    … …  
    167167                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    168168               END DO 
    169                IF ( nn_pnd_scheme > 0 ) THEN 
     169               IF ( ln_pnd_H12 ) THEN 
    170170                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction -- 
    171171                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    … …  
    220220                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    221221               END DO 
    222                IF ( nn_pnd_scheme > 0 ) THEN 
     222               IF ( ln_pnd_H12 ) THEN 
    223223                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction --- 
    224224                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    … …  
    249249         END DO 
    250250         ! MV MP 2016 
    251          IF ( nn_pnd_scheme > 0 ) THEN 
     251         IF ( ln_pnd_H12 ) THEN 
    252252            pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:) 
    253253            pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:) 
    … …  
    755755                  sxyage(:,:,jl)= z2d(:,:) 
    756756               END DO 
    757                IF ( nn_pnd_scheme > 0 ) THEN 
     757               IF ( ln_pnd_H12 ) THEN 
    758758                  DO jl = 1, jpl  
    759759                     WRITE(zchar,'(I2.2)') jl 
    … …  
    833833               syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp 
    834834               sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp 
    835                IF ( nn_pnd_scheme > 0 ) THEN 
     835               IF ( ln_pnd_H12 ) THEN 
    836836                  sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp  
    837837                  syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp  
    … …  
    854854            syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp 
    855855            sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp 
    856             IF ( nn_pnd_scheme > 0 ) THEN 
     856            IF ( ln_pnd_H12 ) THEN 
    857857               sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp  
    858858               syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp  
    … …  
    989989            END DO 
    990990         END DO 
    991          IF ( nn_pnd_scheme > 0 ) THEN 
     991         IF ( ln_pnd_H12 ) THEN 
    992992            DO jl = 1, jpl  
    993993               WRITE(zchar,'(I2.2)') jl 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90

    r8580 r8597  
    124124            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
    125125            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content 
    126             IF ( nn_pnd_scheme > 0 ) THEN 
     126            IF ( ln_pnd_H12 ) THEN 
    127127               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
    128128               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90

    r8565 r8597  
    253253         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
    254254         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
    255          IF ( nn_pnd_scheme > 0 ) THEN 
    256             CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    257             CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    258          ENDIF 
     255         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
     256         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    259257         DO jl = 1, jpl 
    260258            DO jk = 1, nlay_s 
    … …  
    270268         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
    271269         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    272          IF ( nn_pnd_scheme > 0 ) THEN 
    273             CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    274          ENDIF 
     270         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    275271          
    276272         !-----------------------------------------------------------------------------! 
    … …  
    320316         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
    321317         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
    322          IF ( nn_pnd_scheme > 0 ) THEN 
    323             CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    324             CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    325          ENDIF 
     318         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
     319         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    326320         DO jl = 1, jpl 
    327321            DO jk = 1, nlay_s 
    … …  
    337331         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
    338332         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    339          IF ( nn_pnd_scheme > 0 ) THEN 
    340             CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    341          ENDIF 
     333         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    342334 
    343335      ENDIF ! npti > 0 
    … …  
    651643 
    652644               !MV MP 2016 
    653                IF ( nn_pnd_scheme > 0 ) THEN 
     645               IF ( ln_pnd_H12 ) THEN 
    654646                  aprdg1     = a_ip_2d(ji,jl1) * afrdg 
    655647                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
    … …  
    679671               ! Put the melt pond water into the ocean 
    680672               !------------------------------------------             
    681                IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN 
     673               IF ( ln_pnd_fwb ) THEN 
    682674                  wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
    683675                     &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice 
    … …  
    702694               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1 
    703695               ! MV MP 2016 
    704                IF ( nn_pnd_scheme > 0 ) THEN 
     696               IF ( ln_pnd_H12 ) THEN 
    705697                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    706698                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
    … …  
    766758                     &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 
    767759                  ! MV MP 2016 
    768                   IF ( nn_pnd_scheme > 0 ) THEN 
     760                  IF ( ln_pnd_H12 ) THEN 
    769761                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
    770762                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   ) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

    r8592 r8597  
    130130 
    131131      ! --- cloud-sky and overcast-sky ice albedos --- ! 
    132       CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, zalb_cs, zalb_os ) 
     132      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) 
    133133 
    134134      ! albedo depends on cloud fraction because of non-linear spectral effects 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90

    r8564 r8597  
    9494      REAL(wp) ::   ztmelts, zdh 
    9595      INTEGER  ::   i_hemis, i_fill, jl0 
    96       REAL(wp) ::   zarg, zV, zconv, zdv 
     96      REAL(wp) ::   zarg, zV, zconv, zdv, zfac 
    9797      INTEGER , DIMENSION(4)           ::   itest 
    9898      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d 
    … …  
    358358         tn_ice (:,:,:) = t_su (:,:,:) 
    359359 
    360          ! MV MP 2016 
    361360         ! Melt pond volume and fraction 
    362          IF ( ln_pnd ) THEN 
    363             DO jl = 1, jpl 
    364                a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:) 
    365                h_ip     (:,:,jl) = 0.05 * zswitch(:,:) 
    366                a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl)  
    367                v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl) 
    368             END DO 
    369          ELSE 
    370             a_ip(:,:,:)      = 0._wp 
    371             v_ip(:,:,:)      = 0._wp 
    372             a_ip_frac(:,:,:) = 0._wp 
    373             h_ip     (:,:,:) = 0._wp 
    374          ENDIF 
    375          ! END MV MP 2016 
     361         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp 
     362         ELSE                                     ;   zfac = 0._wp   ;   ENDIF  
     363         DO jl = 1, jpl 
     364            a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac 
     365            h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac 
     366         END DO 
     367         a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:)  
     368         v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:) 
    376369 
    377370      ELSE ! if ln_iceini=false 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

    r8565 r8597  
    293293            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
    294294               a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
    295                ! MV MP 2016 
    296                IF ( nn_pnd_scheme > 0 ) THEN 
    297                   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    298                ENDIF 
    299                ! END MV MP 2016 
     295               IF ( ln_pnd_H12 )    a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    300296               h_i_1d(ji) = rn_himin 
    301297            ENDIF 
    … …  
    457453               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    458454               !   
    459                ! MV MP 2016  
    460                IF ( nn_pnd_scheme > 0 ) THEN 
     455               IF ( ln_pnd_H12 ) THEN 
    461456                  !                                                  ! Pond fraction 
    462457                  ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
    … …  
    468463                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
    469464               ENDIF 
    470                ! END MV MP 2016 
    471465               ! 
    472466            ENDIF   ! jl1 >0 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90

    r8592 r8597  
    189189      !! ** purpose  :   read restart file 
    190190      !!---------------------------------------------------------------------- 
    191       INTEGER  :: jk, jl 
    192       REAL(wp) :: zfice, ziter 
     191      INTEGER  ::   jk, jl 
     192      INTEGER  ::   id1            ! local integer 
     193      REAL(wp) ::   zfice, ziter 
    193194      REAL(wp), DIMENSION(jpi,jpj) ::   z2d 
    194195      CHARACTER(len=25) ::   znam 
    … …  
    244245      END DO 
    245246 
    246       DO jl = 1, jpl  
    247          WRITE(zchar,'(I2.2)') jl 
    248          znam = 'a_ip'//'_htc'//zchar 
    249          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    250          a_ip(:,:,jl) = z2d(:,:) 
    251          znam = 'v_ip'//'_htc'//zchar 
    252          CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    253          v_ip(:,:,jl) = z2d(:,:) 
    254       END DO 
     247      id1 = iom_varid( numrir, 'a_ip_htc01' , ldstop = .FALSE. ) 
     248      IF( id1 > 0 ) THEN                       ! fields exist (melt ponds) 
     249         DO jl = 1, jpl  
     250            WRITE(zchar,'(I2.2)') jl 
     251            znam = 'a_ip'//'_htc'//zchar 
     252            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     253            a_ip(:,:,jl) = z2d(:,:) 
     254            znam = 'v_ip'//'_htc'//zchar 
     255            CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
     256            v_ip(:,:,jl) = z2d(:,:) 
     257         END DO 
     258      ELSE                                     ! start from rest 
     259         IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it' 
     260         a_ip(:,:,:) = 0._wp 
     261         v_ip(:,:,:) = 0._wp 
     262      ENDIF 
    255263 
    256264      DO jl = 1, jpl  

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8564 r8597  
    170170         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics       
    171171         ! 
    172          IF ( ln_pnd )                  CALL lim_mp( kt )             ! -- Melt ponds 
     172                                        CALL lim_mp( kt )             ! -- Melt ponds 
    173173         ! 
    174174         IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90

    r8592 r8597  
    166166               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj)  
    167167 
    168             IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj) 
     168            IF ( ln_pnd_fwb )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj) 
    169169 
    170170            ! add the snow melt water to snow mass flux to the ocean 
    … …  
    199199      ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 
    200200      !------------------------------------------------------------------ 
    201       CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
     201      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    202202      ! 
    203203      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    … …  
    246246      CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt  
    247247      CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt  
    248       IF ( ln_pnd )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux 
     248      IF ( ln_pnd_fwb )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux 
    249249 
    250250      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8592 r8597  
    9797      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 
    9898 
    99       ! MV MP 2016 
    100       IF ( ln_pnd ) THEN                     ! Melt pond 
    101          at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) 
    102          vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 
    103       ENDIF 
    104       ! END MP 2016 
     99      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 
     100      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 
    105101 
    106102      ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction   
    … …  
    479475               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice 
    480476               hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0 
     477               IF( ln_pnd_fwb ) THEN  
     478                  wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 
     479               ENDIF 
    481480               !----------------------------------------------------------------- 
    482481               ! Zap snow energy  
    … …  
    499498               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj) 
    500499 
    501             END DO 
    502          END DO 
    503  
    504          IF( ln_pnd ) THEN  
    505             DO jj = 1 , jpj 
    506                DO ji = 1 , jpi 
    507                   IF( ln_pnd_fw )   & 
    508                      &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 
    509                   a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
    510                   v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
    511                END DO 
    512             END DO 
    513          ENDIF 
    514           
     500               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
     501               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     502 
     503            END DO 
     504         END DO 
     505         
    515506      END DO  
    516507 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90

    r8564 r8597  
    128128      CALL iom_put( "snowvol"     , vt_s  * zswi      )        ! snow volume 
    129129       
    130       IF ( ln_pnd ) THEN 
    131          CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction 
    132          CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area 
    133       ENDIF 
     130      CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction 
     131      CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area 
    134132 
    135133      !---------------------------------- 
    … …  
    145143      IF ( iom_use('brinevol_cat') )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )          ! brine volume 
    146144 
    147       IF ( ln_pnd ) THEN 
    148          IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories 
    149          IF ( iom_use('icevmp_cat') )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories 
    150          IF ( iom_use('icehmp_cat') )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories 
    151          IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories 
    152       ENDIF 
     145      IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories 
     146      IF ( iom_use('icevmp_cat') )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories 
     147      IF ( iom_use('icehmp_cat') )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories 
     148      IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories 
    153149 
    154150      !-------------------------------- 
    … …  
    316312      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   & 
    317313      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    318  
    319       ! MV MP 2016 
    320       IF ( ln_pnd ) THEN 
    321          CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   & 
    322       &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    323          CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   & 
    324       &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    325       ENDIF 
    326       ! END MV MP 2016 
    327  
     314      CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   & 
     315      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     316      CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   & 
     317      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    328318      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   & 
    329319      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    … …  
    370360      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
    371361 
    372       ! MV MP 2016 
    373       IF ( ln_pnd ) THEN 
    374          CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) ) 
    375          CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) ) 
    376       ENDIF 
    377       ! END MV MP 2016 
     362      CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) ) 
     363      CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) ) 
    378364 
    379365      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) ) 
    … …  
    384370      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) ) 
    385371      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) ) 
    386       IF ( ln_pnd ) & 
    387          CALL histwrite( kid, "vfxpnd", kt, wfx_pnd     , jpi*jpj, (/1/) ) 
     372      CALL histwrite( kid, "vfxpnd", kt, wfx_pnd        , jpi*jpj, (/1/) ) 
    388373 
    389374      CALL histwrite( kid, "sithicat", kt, h_i         , jpi*jpj*jpl, (/1/) )     

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90

    r8564 r8597  
    77   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1) 
    88   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6 
    9    !!                 ! NB: Only lim_mp_cstt and lim_mp_cesm work in this 
    10    !!                       version 
    119   !!---------------------------------------------------------------------- 
    1210#if defined key_lim3 
    … …  
    1614   !!   lim_mp_init      : some initialization and namelist read 
    1715   !!   lim_mp           : main calling routine 
    18    !!   lim_mp_topo      : main melt pond routine for the "topographic" formulation (FloccoFeltham) 
    19    !!   lim_mp_area      : ??? compute melt pond fraction per category 
    20    !!   lim_mp_perm      : computes permeability (should be a FUNCTION!) 
    21    !!   calc_hpond       : computes melt pond depth 
    22    !!   permeability_phy : computes permeability 
    23  
    2416   !!---------------------------------------------------------------------- 
    2517   USE phycst           ! physical constants 
    2618   USE dom_oce          ! ocean space and time domain 
    27 !  USE sbc_ice          ! Surface boundary condition: ice   fields 
    2819   USE ice              ! LIM-3 variables 
    29 !  USE icecons          ! conservation tests 
    30 !  USE limctl           ! control prints 
    31 !  USE limvar 
    3220   ! 
    3321   USE lbclnk           ! lateral boundary conditions - MPP exchanges 
    … …  
    3725   USE timing           ! Timing 
    3826 
    39 !OLI_CODE    USE ice_oce, ONLY: rdt_ice, tatm_ice 
    40 !OLI_CODE    USE phycst 
    41 !OLI_CODE    USE dom_ice 
    42 !OLI_CODE    USE dom_oce 
    43 !OLI_CODE    USE sbc_oce 
    44 !OLI_CODE    USE sbc_ice 
    45 !OLI_CODE    USE par_ice 
    46 !OLI_CODE    USE par_oce 
    47 !OLI_CODE    USE ice 
    48 !OLI_CODE    USE thd_ice 
    49 !OLI_CODE    USE in_out_manager 
    50 !OLI_CODE    USE lbclnk 
    51 !OLI_CODE    USE lib_mpp 
    52 !OLI_CODE  
    53 !OLI_CODE    IMPLICIT NONE 
    54 !OLI_CODE    PRIVATE 
    55 !OLI_CODE  
    56 !OLI_CODE    PUBLIC   lim_mp_init 
    57 !OLI_CODE    PUBLIC   lim_mp 
    58  
    5927   IMPLICIT NONE 
    6028   PRIVATE 
    … …  
    6230   PUBLIC   lim_mp_init    ! routine called by icestp.F90 
    6331   PUBLIC   lim_mp         ! routine called by icestp.F90 
     32 
     33   INTEGER ::              nice_pnd   ! choice of the type of pond scheme 
     34   !                                               ! associated indices: 
     35   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme 
     36   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme 
     37   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012) 
    6438 
    6539   !! * Substitutions 
    … …  
    7145   !!---------------------------------------------------------------------- 
    7246CONTAINS 
    73  
    74    SUBROUTINE lim_mp_init  
    75       !!------------------------------------------------------------------- 
    76       !!                  ***  ROUTINE lim_mp_init   *** 
    77       !! 
    78       !! ** Purpose : Physical constants and parameters linked to melt ponds 
    79       !!      over sea ice 
    80       !! 
    81       !! ** Method  :  Read the nammp  namelist and check the melt pond   
    82       !!       parameter values called at the first timestep (nit000) 
    83       !! 
    84       !! ** input   :   Namelist nammp   
    85       !!------------------------------------------------------------------- 
    86       INTEGER  ::   ios                 ! Local integer output status for namelist read 
    87       NAMELIST/nammp/  ln_pnd, ln_pnd_rad, ln_pnd_fw, nn_pnd_scheme, rn_apnd, rn_hpnd 
    88       !!------------------------------------------------------------------- 
    89  
    90       REWIND( numnam_ice_ref )              ! Namelist nammp  in reference namelist : Melt Ponds   
    91       READ  ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901) 
    92 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp  in reference namelist', lwp ) 
    93  
    94       REWIND( numnam_ice_cfg )              ! Namelist nammp  in configuration namelist : Melt Ponds 
    95       READ  ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 ) 
    96 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp ) 
    97       IF(lwm) WRITE ( numoni, nammp ) 
    98        
    99       IF(lwp) THEN                        ! control print 
    100          WRITE(numout,*) 
    101          WRITE(numout,*) 'lim_mp_init : ice parameters for melt ponds' 
    102          WRITE(numout,*) '~~~~~~~~~~~~' 
    103          WRITE(numout,*) '   Active melt ponds                                           ln_pnd        = ', ln_pnd 
    104          WRITE(numout,*) '   Active melt ponds radiative coupling                        ln_pnd_rad    = ', ln_pnd_rad 
    105          WRITE(numout,*) '   Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw 
    106          WRITE(numout,*) '   Type of melt pond scheme =0 presc, =1 empirical = 2 topo    nn_pnd_scheme = ', nn_pnd_scheme 
    107          WRITE(numout,*) '   Prescribed pond fraction                                    rn_apnd       = ', rn_apnd 
    108          WRITE(numout,*) '   Prescribed pond depth                                       rn_hpnd       = ', rn_hpnd 
    109       ENDIF 
    110  
    111       IF ( .NOT. ln_pnd ) THEN 
    112          IF(lwp) WRITE(numout,*) 
    113          IF(lwp) WRITE(numout,*) ' Melt ponds are not activated ' 
    114          IF(lwp) WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. ' 
    115          IF(lwp) WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero ' 
    116          ln_pnd_rad    = .FALSE. 
    117          ln_pnd_fw     = .FALSE. 
    118          nn_pnd_scheme = 0 
    119          rn_apnd       = 0._wp 
    120          rn_hpnd       = 0._wp 
    121  
    122          IF(lwp) THEN                     ! control print 
    123             WRITE(numout,*) '    Active melt ponds radiative coupling                        ln_pnd_rad    = ', ln_pnd_rad 
    124             WRITE(numout,*) '    Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw 
    125             WRITE(numout,*) '    Type of melt pond scheme =0 presc, =1 empirical = 2 topo    nn_pnd_scheme = ', nn_pnd_scheme 
    126             WRITE(numout,*) '    Prescribed pond fraction                                    rn_apnd       = ', rn_apnd 
    127             WRITE(numout,*) '    Prescribed pond depth                                       rn_hpnd       = ', rn_hpnd 
    128          ENDIF 
    129       ENDIF 
    130  
    131       IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN 
    132          IF(lwp) WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false '  
    133          ln_pnd_fw     = .FALSE. 
    134          IF(lwp) THEN                     ! control print 
    135             WRITE(numout,*) '    Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw 
    136          ENDIF 
    137       ENDIF 
    138  
    139       IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN 
    140          IF(lwp) WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 ' 
    141          IF(lwp) WRITE(numout,*) ' Run aborted ' 
    142          CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' ) 
    143       ENDIF 
    144  
    145       ! 
    146    END SUBROUTINE lim_mp_init 
    147     
    14847 
    14948   SUBROUTINE lim_mp( kt ) 
    … …  
    15857      !!              -  
    15958      !!------------------------------------------------------------------------------------ 
    160  
    161       INTEGER, INTENT(in) :: kt    ! number of iteration 
    162       INTEGER  ::   ji, jj, jl     ! dummy loop indices 
    163  
     59      INTEGER, INTENT(in) ::   kt    ! number of iteration 
     60      INTEGER ::   ji, jj, jl     ! dummy loop indices 
    16461      !!------------------------------------------------------------------- 
    16562 
    16663      IF( nn_timing == 1 )  CALL timing_start('lim_mp') 
    16764 
    168       SELECT CASE ( nn_pnd_scheme ) 
    169  
    170          CASE (0) 
    171  
    172             CALL lim_mp_cstt       ! staircase melt ponds 
    173  
    174          CASE (1) 
    175  
    176             CALL lim_mp_cesm       ! empirical melt ponds 
    177  
    178          CASE (2) 
    179  
    180             CALL lim_mp_topo   &   ! topographic melt ponds  
    181                       &          (at_i, a_i,                                       & 
    182                       &          vt_i, v_i, v_s,            t_i, sz_i, a_ip_frac,   & 
    183                       &          h_ip,     t_su) 
    184  
     65      SELECT CASE ( nice_pnd ) 
     66 
     67      CASE (np_pndCST) 
     68         !                             !------------------------------! 
     69         CALL lim_mp_CST               ! Constant melt ponds          ! 
     70         !                             !------------------------------! 
     71      CASE (np_pndH12) 
     72         !                             !------------------------------! 
     73         CALL lim_mp_H12               ! Holland et al 2012 melt ponds! 
     74         !                             !------------------------------! 
    18575      END SELECT 
    18676 
    … …  
    18979   END SUBROUTINE lim_mp  
    19080 
    191    SUBROUTINE lim_mp_cstt  
    192        !!------------------------------------------------------------------- 
    193        !!                ***  ROUTINE lim_mp_cstt  *** 
     81   SUBROUTINE lim_mp_CST  
     82       !!------------------------------------------------------------------- 
     83       !!                ***  ROUTINE lim_mp_CST  *** 
    19484       !! 
    19585       !! ** Purpose    : Compute melt pond evolution 
    … …  
    20696       !!     
    20797       !!------------------------------------------------------------------- 
    208        INTEGER                             :: ji, jj, jl 
    209        REAL(wp)                            :: z1_jpl            ! 1/jpl 
    210        !!------------------------------------------------------------------- 
    211  
    212        z1_jpl     = 1. / FLOAT(jpl) 
    21398 
    21499       WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) )  
     100!!       WHERE ( ( a_i > 0._wp ) .AND. ( t_su >= rt0 ) )  
    215101          a_ip_frac = rn_apnd 
    216102          h_ip      = rn_hpnd     
    … …  
    226112       wfx_pnd(:,:) = 0._wp 
    227113 
    228    END SUBROUTINE lim_mp_cstt 
    229  
    230    SUBROUTINE lim_mp_cesm 
    231        !!------------------------------------------------------------------- 
    232        !!                ***  ROUTINE lim_mp_cesm  *** 
     114   END SUBROUTINE lim_mp_CST 
     115 
     116   SUBROUTINE lim_mp_H12 
     117       !!------------------------------------------------------------------- 
     118       !!                ***  ROUTINE lim_mp_H12  *** 
    233119       !! 
    234120       !! ** Purpose    : Compute melt pond evolution 
    … …  
    249135       !!------------------------------------------------------------------- 
    250136 
    251        INTEGER, DIMENSION(jpij)         :: indxi             ! compressed indices for cells with ice melting 
    252        INTEGER, DIMENSION(jpij)         :: indxj             ! 
    253  
    254        REAL(wp), DIMENSION(jpi,jpj)     :: zwfx_mlw          ! available meltwater for melt ponding 
    255        REAL(wp), DIMENSION(jpi,jpj,jpl) :: zrfrac            ! fraction of available meltwater retained for melt ponding 
    256  
    257        REAL(wp), PARAMETER :: zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding 
    258        REAL(wp), PARAMETER :: zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            '' 
    259        REAL(wp), PARAMETER :: zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds 
    260        REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 
    261  
    262        REAL(wp) :: zhi               ! dummy ice thickness 
    263        REAL(wp) :: zhs               ! dummy snow depth 
    264        REAL(wp) :: zTp               ! reference temperature 
    265        REAL(wp) :: zdTs              ! dummy temperature difference 
    266        REAL(wp) :: z1_rhofw          ! inverse freshwater density 
    267        REAL(wp) :: z1_zpnd_aspect    ! inverse pond aspect ratio 
    268        REAL(wp) :: zvpold            ! dummy pond volume 
    269  
    270        INTEGER  :: ji, jj, jl, ij    ! loop indices 
    271        INTEGER  :: icells            ! size of dummy array 
     137       INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting 
     138       INTEGER, DIMENSION(jpij)         ::   indxj             ! 
     139 
     140       REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding 
     141       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zrfrac            ! fraction of available meltwater retained for melt ponding 
     142 
     143       REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding 
     144       REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            '' 
     145       REAL(wp), PARAMETER ::   zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds 
     146       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio 
     147 
     148       REAL(wp) ::   zhi               ! dummy ice thickness 
     149       REAL(wp) ::   zhs               ! dummy snow depth 
     150       REAL(wp) ::   zTp               ! reference temperature 
     151       REAL(wp) ::   zdTs              ! dummy temperature difference 
     152       REAL(wp) ::   z1_rhofw          ! inverse freshwater density 
     153       REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio 
     154       REAL(wp) ::   zvpold            ! dummy pond volume 
     155 
     156       INTEGER  ::   ji, jj, jl, ij    ! loop indices 
     157       INTEGER  ::   icells            ! size of dummy array 
    272158       !!------------------------------------------------------------------- 
    273159        z1_rhofw       = 1. / rhofw  
    … …  
    301187           ! Identify grid cells where ponds should be updated (can probably be improved) 
    302188           !------------------------------------------------------------------------------ 
    303   
    304189           indxi(:) = 0 
    305190           indxj(:) = 0 
    306191           icells   = 0 
    307   
    308192           DO jj = 1, jpj 
    309193             DO ji = 1, jpi 
    … …  
    313197                   indxj(icells) = jj 
    314198                ENDIF 
    315              END DO                 ! ji 
    316           END DO                    ! jj 
     199             END DO 
     200          END DO 
    317201  
    318202          DO ij = 1, icells 
    … …  
    331215                h_ip(ji,jj,jl)      = 0._wp 
    332216  
    333                 IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean 
     217                IF ( ln_pnd_fwb ) & !--- Give freshwater to the ocean 
    334218                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)  
    335219  
    … …  
    352236                !--- Dump meltwater due to refreezing ( of course this is wrong 
    353237                !--- but this parameterization is too simple ) 
    354                 IF ( ln_pnd_fw ) & 
     238                IF ( ln_pnd_fwb ) & 
    355239                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice 
    356240  
    … …  
    361245  
    362246                a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 
    363   
    364              !----------------------------------------------------------- 
    365              ! Limit pond depth 
    366              !----------------------------------------------------------- 
    367              ! The original version has pond depth limitation, which I did not 
    368              ! keep here. Maybe we need it later on 
    369              ! 
    370247              
    371248             ENDIF 
    … …  
    377254        !--- Remove retained meltwater from surface fluxes  
    378255  
    379         IF ( ln_pnd_fw ) THEN 
    380   
    381             wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )  
    382  
    383             wfx_sum(:,:)     = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
     256        IF ( ln_pnd_fwb ) THEN 
     257  
     258           wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )  
     259           wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
    384260 
    385261        ENDIF 
    386262 
    387    END SUBROUTINE lim_mp_cesm 
    388  
    389    SUBROUTINE lim_mp_topo    (aice,      aicen,     & 
    390                               vice,      vicen,     & 
    391                               vsnon,                & 
    392                               ticen,     salin,     & 
    393                               a_ip_frac, h_ip,      & 
    394                                          Tsfc ) 
    395        !!------------------------------------------------------------------- 
    396        !!                ***  ROUTINE lim_mp_topo  *** 
    397        !! 
    398        !! ** Purpose :   Compute melt pond evolution based on the ice  
    399        !!                topography as inferred from the ice thickness  
    400        !!                distribution.   
    401        !! 
    402        !! ** Method  :   This code is initially based on Flocco and Feltham  
    403        !!                (2007) and Flocco et al. (2010). More to come... 
    404        !! 
    405        !! ** Tunable parameters : 
    406        !!  
    407        !! ** Note : 
    408        !! 
    409        !! ** References 
    410        !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond  
    411        !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi:  
    412        !!    10.1029/2006JC003836. 
    413        !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of 
    414        !!    a physically based melt pond scheme into the sea ice component of a 
    415        !!    climate model.  J. Geophys. Res. 115, C08012,  
    416        !!    doi: 10.1029/2009JC005568. 
    417        !!     
    418        !!------------------------------------------------------------------- 
    419   
    420        REAL (wp), DIMENSION (jpi,jpj), & 
    421           INTENT(IN) :: & 
    422           aice, &    ! total ice area fraction 
    423           vice       ! total ice volume (m) 
    424   
    425        REAL (wp), DIMENSION (jpi,jpj,jpl), & 
    426           INTENT(IN) :: & 
    427           aicen, &   ! ice area fraction, per category 
    428           vsnon, &   ! snow volume, per category (m) 
    429           vicen      ! ice volume, per category (m) 
    430   
    431        REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & 
    432           INTENT(IN) :: & 
    433           ticen, &   ! ice enthalpy, per category 
    434           salin 
    435   
    436        REAL (wp), DIMENSION (jpi,jpj,jpl), & 
    437           INTENT(INOUT) :: & 
    438           a_ip_frac , &   ! pond area fraction of ice, per ice category 
    439           h_ip       ! pond depth, per ice category (m) 
    440   
    441        REAL (wp), DIMENSION (jpi,jpj,jpl), & 
    442           INTENT(IN) :: & 
    443           Tsfc       ! snow/sea ice surface temperature 
    444   
    445        ! local variables 
    446        REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
    447           zTsfcn, & ! ice/snow surface temperature (C) 
    448           zvolpn, & ! pond volume per unit area, per category (m) 
    449           zvuin     ! water-equivalent volume of ice lid on melt pond ('upper ice', m) 
    450   
    451        REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
    452           zapondn,& ! pond area fraction, per category 
    453           zhpondn   ! pond depth, per category (m) 
    454   
    455        REAL (wp), DIMENSION (jpi,jpj) :: & 
    456           zvolp       ! total volume of pond, per unit area of pond (m) 
    457   
    458        REAL (wp) :: & 
    459           zhi,    & ! ice thickness (m) 
    460           zdHui,  & ! change in thickness of ice lid (m) 
    461           zomega, & ! conduction 
    462           zdTice, & ! temperature difference across ice lid (C) 
    463           zdvice, & ! change in ice volume (m) 
    464           zTavg,  & ! mean surface temperature across categories (C) 
    465           zTp,    & ! pond freezing temperature (C) 
    466           zdvn      ! change in melt pond volume for fresh water budget 
    467        INTEGER, DIMENSION (jpi*jpj) :: & 
    468           indxi, indxj    ! compressed indices for cells with ice melting 
    469   
    470        INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices 
    471   
    472        INTEGER, DIMENSION (jpl) :: & 
    473           kcells          ! cells where ice lid combines with vice 
    474   
    475        INTEGER, DIMENSION (jpi*jpj,jpl) :: & 
    476           indxii, indxjj  ! i,j indices for kcells loop 
    477   
    478        REAL (wp), parameter :: & 
    479           zhicemin  = 0.1_wp , & ! minimum ice thickness with ponds (m)  
    480           zTd       = 0.15_wp, & ! temperature difference for freeze-up (C) 
    481           zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) 
    482           zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) 
    483           z0       = 0._wp,    &  
    484           zTimelt   = 0._wp,    & 
    485           z01      = 0.01_wp,  & 
    486           z25      = 0.25_wp,  & 
    487           z5       = 0.5_wp 
    488  
    489        !--------------------------------------------------------------- 
    490        ! Initialization 
    491        !--------------------------------------------------------------- 
    492   
    493        zhpondn(:,:,:) = 0._wp 
    494        zapondn(:,:,:) = 0._wp 
    495        indxii(:,:) = 0 
    496        indxjj(:,:) = 0 
    497        kcells(:)   = 0 
    498  
    499        zvolp(:,:) = wfx_sum(:,:) + wfx_snw_sum(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds 
    500        zTsfcn(:,:,:) = zTsfcn(:,:,:) - rt0                   ! Convert in Celsius 
    501   
    502        ! The freezing temperature for meltponds is assumed slightly below 0C, 
    503        ! as if meltponds had a little salt in them.  The salt budget is not 
    504        ! altered for meltponds, but if it were then an actual pond freezing  
    505        ! temperature could be computed. 
    506   
    507        ! zTp = zTimelt - zTd  ---> for lids 
    508   
    509        !----------------------------------------------------------------- 
    510        ! Identify grid cells with ponds 
    511        !----------------------------------------------------------------- 
    512   
    513        icells = 0 
    514        DO j = 1, jpj 
    515        DO i = 1, jpi 
    516           zhi = z0 
    517           IF (aice(i,j) > epsi10 ) zhi = vice(i,j)/aice(i,j) 
    518           IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. & 
    519              zvolp(i,j) > zmin_volp*aice(i,j)) THEN 
    520              icells = icells + 1 
    521              indxi(icells) = i 
    522              indxj(icells) = j 
    523           ELSE  ! remove ponds on thin ice 
    524              !fpond(i,j) = fpond(i,j) - zvolp(i,j) 
    525              zvolpn(i,j,:) = z0 
    526              zvuin (i,j,:) = z0 
    527              zvolp (i,j) = z0 
    528           END IF 
    529        END DO                     ! i 
    530        END DO                     ! j 
    531   
    532        DO ij = 1, icells 
    533           i = indxi(ij) 
    534           j = indxj(ij) 
    535   
    536           !-------------------------------------------------------------- 
    537           ! calculate pond area and depth 
    538           !-------------------------------------------------------------- 
    539           CALL lim_mp_area(aice(i,j),vice(i,j), & 
    540                     aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), & 
    541                     ticen(i,j,:,:), salin(i,j,:,:), & 
    542                     zvolpn(i,j,:), zvolp(i,j), & 
    543                     zapondn(i,j,:),zhpondn(i,j,:), zdvn) 
    544          ! outputs are 
    545          ! - zdvn 
    546          ! - zvolpn 
    547          ! - zvolp 
    548          ! - zapondn 
    549          ! - zhpondn 
    550  
    551           wfx_pnd(i,j) = wfx_pnd(i,j) + zdvn ! update flux from ponds to ocean 
    552   
    553           ! mean surface temperature MV - why do we need that ? --> for the lid 
    554  
    555           ! zTavg = z0 
    556           ! DO n = 1, jpl 
    557           !   zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n) 
    558           ! END DO 
    559           ! zTavg = zTavg / aice(i,j) 
    560   
    561        END DO ! ij 
    562   
    563        !--------------------------------------------------------------- 
    564        ! Update pond volume and fraction 
    565        !--------------------------------------------------------------- 
    566   
    567        a_ip(:,:,:) = zapondn(:,:,:) 
    568        v_ip(:,:,:) = zapondn(:,:,:) * zhpondn(:,:,:) 
    569        a_ip_frac(:,:,:) = 0._wp 
    570        h_ip     (:,:,:) = 0._wp 
    571  
    572     END SUBROUTINE lim_mp_topo 
    573  
    574     SUBROUTINE lim_mp_area(aice,vice, & 
    575                          aicen, vicen, vsnon, ticen, & 
    576                          salin, zvolpn, zvolp,         & 
    577                          zapondn,zhpondn,dvolp) 
    578  
    579        !!------------------------------------------------------------------- 
    580        !!                ***  ROUTINE lim_mp_area *** 
    581        !! 
    582        !! ** Purpose : Given the total volume of meltwater, update 
    583        !!              pond fraction (a_ip) and depth (should be volume) 
    584        !! 
    585        !! ** 
    586        !!               
    587        !!------------------------------------------------------------------ 
    588  
    589        REAL (wp), INTENT(IN) :: & 
    590           aice,vice 
    591   
    592        REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
    593           aicen, vicen, vsnon 
    594   
    595        REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 
    596           ticen, salin 
    597   
    598        REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 
    599           zvolpn 
    600   
    601        REAL (wp), INTENT(INOUT) :: & 
    602           zvolp, dvolp 
    603   
    604        REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 
    605           zapondn, zhpondn 
    606   
    607        INTEGER :: & 
    608           n, ns,   & 
    609           m_index, & 
    610           permflag 
    611   
    612        REAL (wp), DIMENSION(jpl) :: & 
    613           hicen, & 
    614           hsnon, & 
    615           asnon, & 
    616           alfan, & 
    617           betan, & 
    618           cum_max_vol, & 
    619           reduced_aicen         
    620   
    621        REAL (wp), DIMENSION(0:jpl) :: & 
    622           cum_max_vol_tmp 
    623   
    624        REAL (wp) :: & 
    625           hpond, & 
    626           drain, & 
    627           floe_weight, & 
    628           pressure_head, & 
    629           hsl_rel, & 
    630           deltah, & 
    631           perm, & 
    632           msno 
    633   
    634        REAL (wp), parameter :: &  
    635           viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s 
    636           z0        = 0.0_wp    , & 
    637           c1        = 1.0_wp    , & 
    638           p4        = 0.4_wp    , & 
    639           p6        = 0.6_wp    , & 
    640           epsi10      = 1.0e-11_wp 
    641            
    642       !-----------| 
    643       !           | 
    644       !           |-----------| 
    645       !___________|___________|______________________________________sea-level 
    646       !           |           | 
    647       !           |           |---^--------| 
    648       !           |           |   |        | 
    649       !           |           |   |        |-----------|              |------- 
    650       !           |           |   |alfan(n)|           |              | 
    651       !           |           |   |        |           |--------------| 
    652       !           |           |   |        |           |              | 
    653       !---------------------------v------------------------------------------- 
    654       !           |           |   ^        |           |              | 
    655       !           |           |   |        |           |--------------| 
    656       !           |           |   |betan(n)|           |              | 
    657       !           |           |   |        |-----------|              |------- 
    658       !           |           |   |        | 
    659       !           |           |---v------- | 
    660       !           |           | 
    661       !           |-----------| 
    662       !           | 
    663       !-----------| 
    664       
    665        !------------------------------------------------------------------- 
    666        ! initialize 
    667        !------------------------------------------------------------------- 
    668   
    669        DO n = 1, jpl 
    670   
    671           zapondn(n) = z0 
    672           zhpondn(n) = z0 
    673   
    674           !---------------------------------------- 
    675           ! X) compute the effective snow fraction 
    676           !---------------------------------------- 
    677           IF (aicen(n) < epsi10)  THEN 
    678              hicen(n) =  z0  
    679              hsnon(n) = z0 
    680              reduced_aicen(n) = z0 
    681           ELSE 
    682              hicen(n) = vicen(n) / aicen(n) 
    683              hsnon(n) = vsnon(n) / aicen(n) 
    684              reduced_aicen(n) = c1 ! n=jpl 
    685              IF (n < jpl) reduced_aicen(n) = aicen(n) & 
    686                                   * (-0.024_wp*hicen(n) + 0.832_wp)  
    687              asnon(n) = reduced_aicen(n)  ! effective snow fraction (empirical) 
    688              ! MV should check whether this makes sense to have the same effective snow fraction in here 
    689           END IF 
    690   
    691  ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
    692  ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
    693  ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
    694  ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
    695  ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
    696  ! Here, alfa = 60% of the ice area (and since hice is constant in a category,  
    697  ! alfan = 60% of the ice volume) in each category lies above the reference line,  
    698  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
    699  
    700  ! MV:   
    701  ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 
    702  ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 
    703  
    704  ! Where does that choice come from 
    705   
    706           alfan(n) = 0.6 * hicen(n) 
    707           betan(n) = 0.4 * hicen(n) 
    708          
    709           cum_max_vol(n)     = z0 
    710           cum_max_vol_tmp(n) = z0 
    711       
    712        END DO ! jpl 
    713   
    714        cum_max_vol_tmp(0) = z0 
    715        drain = z0 
    716        dvolp = z0 
    717  
    718        !---------------------------------------------------------- 
    719        ! x) Drain overflow water, update pond fraction and volume 
    720        !---------------------------------------------------------- 
    721         
    722        !-------------------------------------------------------------------------- 
    723        ! the maximum amount of water that can be contained up to each ice category 
    724        !-------------------------------------------------------------------------- 
    725  
    726        ! MV 
    727        ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 
    728        ! Then the excess volume cum_max_vol(jl) drains out of the system 
    729        ! It should be added to wfx_pnd 
    730        ! END MV 
    731       
    732        DO n = 1, jpl-1 ! last category can not hold any volume 
    733   
    734           IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN 
    735   
    736              ! total volume in level including snow 
    737              cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 
    738                 (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))  
    739   
    740              ! subtract snow solid volumes from lower categories in current level 
    741              DO ns = 1, n  
    742                 cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 
    743                    - rhosn/rhofw * &   ! free air fraction that can be filled by water 
    744                      asnon(ns)  * &    ! effective areal fraction of snow in that category 
    745                      max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- & 
    746                                    alfan(n)), z0)   
    747              END DO 
    748  
    749           ELSE ! assume higher categories unoccupied 
    750              cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 
    751           END IF 
    752           !IF (cum_max_vol_tmp(n) < z0) THEN 
    753           !   call abort_ice('negative melt pond volume') 
    754           !END IF 
    755        END DO 
    756        cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
    757        cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
    758       
    759        !---------------------------------------------------------------- 
    760        ! is there more meltwater than can be held in the floe? 
    761        !---------------------------------------------------------------- 
    762        IF (zvolp >= cum_max_vol(jpl)) THEN 
    763           drain = zvolp - cum_max_vol(jpl) + epsi10 
    764           zvolp = zvolp - drain ! update meltwater volume available 
    765           dvolp = drain         ! this is the drained water 
    766           IF (zvolp < epsi10) THEN 
    767              dvolp = dvolp + zvolp 
    768              zvolp = z0 
    769           END IF 
    770        END IF 
    771       
    772        ! height and area corresponding to the remaining volume 
    773   
    774 !      call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 
    775 !           zvolp, cum_max_vol, hpond, m_index) 
    776       
    777        DO n=1, m_index 
    778           zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update 
    779                                                    !  volume instead, no ? 
    780           zapondn(n) = reduced_aicen(n)  
    781           ! in practise, pond fraction depends on the empirical snow fraction 
    782           ! so in turn on ice thickness 
    783        END DO 
    784       
    785        !------------------------------------------------------------------------ 
    786        ! Drainage through brine network (permeability) 
    787        !------------------------------------------------------------------------ 
    788        !!! drainage due to ice permeability - Darcy's law 
    789       
    790        ! sea water level  
    791        msno = z0 
    792        DO n=1,jpl 
    793          msno = msno + vsnon(n) * rhosn 
    794        END DO 
    795        floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice 
    796        hsl_rel = floe_weight / rau0 & 
    797                - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 
    798       
    799        deltah = hpond - hsl_rel 
    800        pressure_head = grav * rau0 * max(deltah, z0) 
    801   
    802        ! drain IF ice is permeable     
    803        permflag = 0 
    804        IF (pressure_head > z0) THEN 
    805        DO n = 1, jpl-1 
    806           IF (hicen(n) /= z0) THEN 
    807              perm = 0. ! MV ugly dummy patch 
    808              CALL lim_mp_perm(ticen(:,n), salin(:,n), vicen(n), perm) 
    809              IF (perm > z0) permflag = 1 
    810  
    811              drain = perm*zapondn(n)*pressure_head*rdt_ice / & 
    812                                       (viscosity*hicen(n)) 
    813              dvolp = dvolp + min(drain, zvolp) 
    814              zvolp = max(zvolp - drain, z0) 
    815              IF (zvolp < epsi10) THEN 
    816                 dvolp = dvolp + zvolp 
    817                 zvolp = z0 
    818              END IF 
    819           END IF 
    820        END DO 
    821    
    822        ! adjust melt pond DIMENSIONs 
    823        IF (permflag > 0) THEN 
    824           ! recompute pond depth     
    825 !         CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 
    826 !                         zvolp, cum_max_vol, hpond, m_index) 
    827           DO n=1, m_index 
    828              zhpondn(n) = hpond - alfan(n) + alfan(1) 
    829              zapondn(n) = reduced_aicen(n)  
    830           END DO 
    831        END IF 
    832        END IF ! pressure_head 
    833   
    834        !------------------------------- 
    835        ! X) remove water from the snow 
    836        !------------------------------- 
    837        !------------------------------------------------------------------------ 
    838        ! total melt pond volume in category DOes not include snow volume 
    839        ! snow in melt ponds is not melted 
    840        !------------------------------------------------------------------------ 
    841   
    842        ! Calculate pond volume for lower categories 
    843        DO n=1,m_index-1 
    844           zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow 
    845                    - (rhosn/rhofw) * asnon(n) * min(hsnon(n), zhpondn(n)) 
    846        END DO 
    847   
    848        ! Calculate pond volume for highest category = remaining pond volume 
    849  
    850        ! The following is completely unclear to Martin at least 
    851        ! Could we redefine properly and recode in a more readable way ? 
    852  
    853        ! m_index = last category with melt pond 
    854  
    855        IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water 
    856  
    857        IF (m_index > 1) THEN 
    858          IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 
    859            zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) !  
    860          ELSE 
    861            zvolpn(m_index) = z0 
    862            zhpondn(m_index) = z0 
    863            zapondn(m_index) = z0 
    864            ! If remaining pond volume is negative reduce pond volume of  
    865            ! lower category 
    866            IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &  
    867              zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))& 
    868                                 + zvolp 
    869          END IF 
    870        END IF 
    871   
    872        DO n=1,m_index 
    873           IF (zapondn(n) > epsi10) THEN 
    874               zhpondn(n) = zvolpn(n) / zapondn(n) 
    875           ELSE 
    876              dvolp = dvolp + zvolpn(n) 
    877              zhpondn(n) = z0 
    878              zvolpn(n) = z0 
    879              zapondn(n) = z0 
    880           end IF 
    881        END DO 
    882        DO n = m_index+1, jpl 
    883           zhpondn(n) = z0 
    884           zapondn(n) = z0 
    885           zvolpn (n) = z0 
    886        END DO 
    887   
    888     END SUBROUTINE lim_mp_area 
    889  
    890 !OLI_CODE    
    891 !OLI_CODE  
    892 !OLI_CODE    SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, & 
    893 !OLI_CODE                          zvolp, cum_max_vol,                & 
    894 !OLI_CODE                          hpond, m_index) 
    895 !OLI_CODE       !!------------------------------------------------------------------- 
    896 !OLI_CODE       !!                ***  ROUTINE calc_hpond  *** 
    897 !OLI_CODE       !! 
    898 !OLI_CODE       !! ** Purpose :   Compute melt pond depth 
    899 !OLI_CODE       !!------------------------------------------------------------------- 
    900 !OLI_CODE      
    901 !OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
    902 !OLI_CODE          aicen, & 
    903 !OLI_CODE          asnon, & 
    904 !OLI_CODE          hsnon, & 
    905 !OLI_CODE          rhos,  & 
    906 !OLI_CODE          alfan, & 
    907 !OLI_CODE          cum_max_vol 
    908 !OLI_CODE      
    909 !OLI_CODE       REAL (wp), INTENT(IN) :: & 
    910 !OLI_CODE          zvolp 
    911 !OLI_CODE      
    912 !OLI_CODE       REAL (wp), INTENT(OUT) :: & 
    913 !OLI_CODE          hpond 
    914 !OLI_CODE      
    915 !OLI_CODE       INTEGER, INTENT(OUT) :: & 
    916 !OLI_CODE          m_index 
    917 !OLI_CODE      
    918 !OLI_CODE       INTEGER :: n, ns 
    919 !OLI_CODE      
    920 !OLI_CODE       REAL (wp), DIMENSION(0:jpl+1) :: & 
    921 !OLI_CODE          hitl, & 
    922 !OLI_CODE          aicetl 
    923 !OLI_CODE      
    924 !OLI_CODE       REAL (wp) :: & 
    925 !OLI_CODE          rem_vol, & 
    926 !OLI_CODE          area, & 
    927 !OLI_CODE          vol, & 
    928 !OLI_CODE          tmp, & 
    929 !OLI_CODE          z0   = 0.0_wp,    &    
    930 !OLI_CODE          epsi10 = 1.0e-11_wp 
    931 !OLI_CODE      
    932 !OLI_CODE       !---------------------------------------------------------------- 
    933 !OLI_CODE       ! hpond is zero if zvolp is zero - have we fully drained?  
    934 !OLI_CODE       !---------------------------------------------------------------- 
    935 !OLI_CODE      
    936 !OLI_CODE       IF (zvolp < epsi10) THEN 
    937 !OLI_CODE        hpond = z0 
    938 !OLI_CODE        m_index = 0 
    939 !OLI_CODE       ELSE 
    940 !OLI_CODE         
    941 !OLI_CODE        !---------------------------------------------------------------- 
    942 !OLI_CODE        ! Calculate the category where water fills up to  
    943 !OLI_CODE        !---------------------------------------------------------------- 
    944 !OLI_CODE         
    945 !OLI_CODE        !----------| 
    946 !OLI_CODE        !          | 
    947 !OLI_CODE        !          | 
    948 !OLI_CODE        !          |----------|                                     -- -- 
    949 !OLI_CODE        !__________|__________|_________________________________________ ^ 
    950 !OLI_CODE        !          |          |             rem_vol     ^                | Semi-filled 
    951 !OLI_CODE        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
    952 !OLI_CODE        !          |          |          |              | 
    953 !OLI_CODE        !          |          |          |              |hpond 
    954 !OLI_CODE        !          |          |          |----------|   |     |------- 
    955 !OLI_CODE        !          |          |          |          |   |     | 
    956 !OLI_CODE        !          |          |          |          |---v-----| 
    957 !OLI_CODE        !          |          | m_index  |          |         | 
    958 !OLI_CODE        !------------------------------------------------------------- 
    959 !OLI_CODE         
    960 !OLI_CODE        m_index = 0  ! 1:m_index categories have water in them 
    961 !OLI_CODE        DO n = 1, jpl 
    962 !OLI_CODE           IF (zvolp <= cum_max_vol(n)) THEN 
    963 !OLI_CODE              m_index = n 
    964 !OLI_CODE              IF (n == 1) THEN 
    965 !OLI_CODE                 rem_vol = zvolp 
    966 !OLI_CODE              ELSE 
    967 !OLI_CODE                 rem_vol = zvolp - cum_max_vol(n-1) 
    968 !OLI_CODE              END IF 
    969 !OLI_CODE              exit ! to break out of the loop 
    970 !OLI_CODE           END IF 
    971 !OLI_CODE        END DO 
    972 !OLI_CODE        m_index = min(jpl-1, m_index) 
    973 !OLI_CODE         
    974 !OLI_CODE        !---------------------------------------------------------------- 
    975 !OLI_CODE        ! semi-filled layer may have m_index different snow in it 
    976 !OLI_CODE        !---------------------------------------------------------------- 
    977 !OLI_CODE         
    978 !OLI_CODE        !-----------------------------------------------------------  ^ 
    979 !OLI_CODE        !                                                             |  alfan(m_index+1) 
    980 !OLI_CODE        !                                                             | 
    981 !OLI_CODE        !hitl(3)-->                             |----------|          | 
    982 !OLI_CODE        !hitl(2)-->                |------------| * * * * *|          | 
    983 !OLI_CODE        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
    984 !OLI_CODE        !hitl(0)-->-------------------------------------------------  |  ^ 
    985 !OLI_CODE        !                various snow from lower categories          |  |alfa(m_index) 
    986 !OLI_CODE         
    987 !OLI_CODE        ! hitl - heights of the snow layers from thinner and current categories 
    988 !OLI_CODE        ! aicetl - area of each snow depth in this layer 
    989 !OLI_CODE         
    990 !OLI_CODE        hitl(:) = z0 
    991 !OLI_CODE        aicetl(:) = z0 
    992 !OLI_CODE        DO n = 1, m_index 
    993 !OLI_CODE           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
    994 !OLI_CODE                                  alfan(m_index+1) - alfan(m_index)), z0) 
    995 !OLI_CODE           aicetl(n) = asnon(n) 
    996 !OLI_CODE            
    997 !OLI_CODE           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
    998 !OLI_CODE        END DO 
    999 !OLI_CODE        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
    1000 !OLI_CODE        aicetl(m_index+1) = z0 
    1001 !OLI_CODE         
    1002 !OLI_CODE        !---------------------------------------------------------------- 
    1003 !OLI_CODE        ! reorder array according to hitl  
    1004 !OLI_CODE        ! snow heights not necessarily in height order 
    1005 !OLI_CODE        !---------------------------------------------------------------- 
    1006 !OLI_CODE         
    1007 !OLI_CODE        DO ns = 1, m_index+1 
    1008 !OLI_CODE           DO n = 0, m_index - ns + 1 
    1009 !OLI_CODE              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
    1010 !OLI_CODE                 tmp = hitl(n) 
    1011 !OLI_CODE                 hitl(n) = hitl(n+1) 
    1012 !OLI_CODE                 hitl(n+1) = tmp 
    1013 !OLI_CODE                 tmp = aicetl(n) 
    1014 !OLI_CODE                 aicetl(n) = aicetl(n+1) 
    1015 !OLI_CODE                 aicetl(n+1) = tmp 
    1016 !OLI_CODE              END IF 
    1017 !OLI_CODE           END DO 
    1018 !OLI_CODE        END DO 
    1019 !OLI_CODE         
    1020 !OLI_CODE        !---------------------------------------------------------------- 
    1021 !OLI_CODE        ! divide semi-filled layer into set of sublayers each vertically homogenous 
    1022 !OLI_CODE        !---------------------------------------------------------------- 
    1023 !OLI_CODE         
    1024 !OLI_CODE        !hitl(3)---------------------------------------------------------------- 
    1025 !OLI_CODE        !                                                       | * * * * * * * *   
    1026 !OLI_CODE        !                                                       |* * * * * * * * *  
    1027 !OLI_CODE        !hitl(2)---------------------------------------------------------------- 
    1028 !OLI_CODE        !                                    | * * * * * * * *  | * * * * * * * *   
    1029 !OLI_CODE        !                                    |* * * * * * * * * |* * * * * * * * *  
    1030 !OLI_CODE        !hitl(1)---------------------------------------------------------------- 
    1031 !OLI_CODE        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *   
    1032 !OLI_CODE        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *  
    1033 !OLI_CODE        !hitl(0)---------------------------------------------------------------- 
    1034 !OLI_CODE        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)             
    1035 !OLI_CODE         
    1036 !OLI_CODE        ! move up over layers incrementing volume 
    1037 !OLI_CODE        DO n = 1, m_index+1 
    1038 !OLI_CODE            
    1039 !OLI_CODE           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
    1040 !OLI_CODE                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
    1041 !OLI_CODE            
    1042 !OLI_CODE           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
    1043 !OLI_CODE            
    1044 !OLI_CODE           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
    1045 !OLI_CODE              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - & 
    1046 !OLI_CODE                           alfan(1) 
    1047 !OLI_CODE              exit 
    1048 !OLI_CODE           ELSE  ! still in sub-layer below the sub-layer with the depth 
    1049 !OLI_CODE              rem_vol = rem_vol - vol 
    1050 !OLI_CODE           END IF 
    1051 !OLI_CODE            
    1052 !OLI_CODE        END DO 
    1053 !OLI_CODE         
    1054 !OLI_CODE       END IF 
    1055 !OLI_CODE      
    1056 !OLI_CODE    END SUBROUTINE calc_hpond 
    1057 !OLI_CODE    
    1058 !OLI_CODE  
    1059     SUBROUTINE lim_mp_perm(ticen, salin, vicen, perm) 
    1060        !!------------------------------------------------------------------- 
    1061        !!                ***  ROUTINE lim_mp_perm *** 
    1062        !! 
    1063        !! ** Purpose :   Determine the liquid fraction of brine in the ice 
    1064        !!                and its permeability 
    1065        !!------------------------------------------------------------------- 
    1066        REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
    1067           ticen,  & ! energy of melting for each ice layer (J/m2) 
    1068           salin 
    1069   
    1070        REAL (wp), INTENT(IN) :: & 
    1071           vicen     ! ice volume 
    1072       
    1073        REAL (wp), INTENT(OUT) :: & 
    1074           perm      ! permeability 
    1075   
    1076        REAL (wp) ::   & 
    1077           Sbr        ! brine salinity 
    1078   
    1079        REAL (wp), DIMENSION(nlay_i) ::   & 
    1080           Tin, &    ! ice temperature 
    1081           phi       ! liquid fraction 
    1082   
    1083        INTEGER :: k 
    1084       
    1085        REAL (wp) :: & 
    1086           c2    = 2.0_wp 
    1087    
    1088        !----------------------------------------------------------------- 
    1089        ! Compute ice temperatures from enthalpies using quadratic formula 
    1090        !----------------------------------------------------------------- 
    1091   
    1092        DO k = 1,nlay_i 
    1093           Tin(k) = ticen(k)  
    1094        END DO 
    1095   
    1096        !----------------------------------------------------------------- 
    1097        ! brine salinity and liquid fraction 
    1098        !----------------------------------------------------------------- 
    1099   
    1100        IF (maxval(Tin-rtt) <= -c2) THEN 
    1101   
    1102           DO k = 1,nlay_i 
    1103              Sbr = - 1.2_wp                 & 
    1104                    -21.8_wp     * (Tin(k)-rtt)    & 
    1105                    - 0.919_wp   * (Tin(k)-rtt)**2 & 
    1106                    - 0.01878_wp * (Tin(k)-rtt)**3 
    1107              phi(k) = salin(k)/Sbr ! liquid fraction 
    1108           END DO ! k 
    1109          
    1110        ELSE 
    1111   
    1112           DO k = 1,nlay_i 
    1113              Sbr = -17.6_wp    * (Tin(k)-rtt)    & 
    1114                    - 0.389_wp  * (Tin(k)-rtt)**2 & 
    1115                    - 0.00362_wp* (Tin(k)-rtt)**3 
    1116              phi(k) = salin(k)/Sbr ! liquid fraction 
    1117           END DO 
    1118   
    1119        END IF 
    1120       
    1121        !----------------------------------------------------------------- 
    1122        ! permeability 
    1123        !----------------------------------------------------------------- 
    1124   
    1125        perm = 3.0e-08_wp * (minval(phi))**3 ! REFERENCE PLEASE (this fucking 
    1126                                             ! bastard of Golden) 
    1127       
    1128     END SUBROUTINE lim_mp_perm 
    1129 !OLI_CODE    
    1130 !OLI_CODE #else 
    1131 !OLI_CODE    !!---------------------------------------------------------------------- 
    1132 !OLI_CODE    !!   Default option         Dummy Module          No LIM-3 sea-ice model 
    1133 !OLI_CODE    !!---------------------------------------------------------------------- 
    1134 !OLI_CODE CONTAINS 
    1135 !OLI_CODE    SUBROUTINE lim_mp_init          ! Empty routine 
    1136 !OLI_CODE    END SUBROUTINE lim_mp_init    
    1137 !OLI_CODE    SUBROUTINE lim_mp               ! Empty routine 
    1138 !OLI_CODE    END SUBROUTINE lim_mp 
    1139 !OLI_CODE    SUBROUTINE compute_mp_topo      ! Empty routine 
    1140 !OLI_CODE    END SUBROUTINE compute_mp_topo         
    1141 !OLI_CODE    SUBROUTINE pond_area            ! Empty routine 
    1142 !OLI_CODE    END SUBROUTINE pond_area    
    1143 !OLI_CODE    SUBROUTINE calc_hpond           ! Empty routine 
    1144 !OLI_CODE    END SUBROUTINE calc_hpond    
    1145 !OLI_CODE    SUBROUTINE permeability_phy     ! Empty routine 
    1146 !OLI_CODE    END SUBROUTINE permeability_phy    
    1147 !OLI_CODE #endif 
    1148 !OLI_CODE    !!====================================================================== 
    1149 !OLI_CODE END MODULE limmp_topo 
    1150  
    1151  
    1152 !OLI_CODE MODULE limmp_topo 
    1153 !OLI_CODE    !!====================================================================== 
    1154 !OLI_CODE    !!                       ***  MODULE limmp_topo *** 
    1155 !OLI_CODE    !! LIM-3 sea-ice :  computation of melt ponds' properties 
    1156 !OLI_CODE    !!====================================================================== 
    1157 !OLI_CODE    !! History :  Original code by Daniela Flocco and Adrian Turner 
    1158 !OLI_CODE    !!            ! 2012-09 (O. Lecomte) Adaptation for routine inclusion in  
    1159 !OLI_CODE    !!                      NEMO-LIM3.1 
    1160 !OLI_CODE    !!            ! 2016-11 (O. Lecomte, C. Rousset, M. Vancoppenolle)  
    1161 !OLI_CODE    !!                      Adaptation for merge with NEMO-LIM3.6 
    1162 !OLI_CODE    !!---------------------------------------------------------------------- 
    1163 !OLI_CODE #if defined key_lim3 
    1164 !OLI_CODE    !!---------------------------------------------------------------------- 
    1165 !OLI_CODE    !!   'key_lim3'                                      LIM-3 sea-ice model 
    1166 !OLI_CODE    !!---------------------------------------------------------------------- 
    1167 !OLI_CODE    !!   lim_mp_init     : melt pond properties initialization 
    1168 !OLI_CODE    !!   lim_mp          : melt pond routine caller 
    1169 !OLI_CODE    !!   compute_mp_topo : Actual melt pond routine 
    1170 !OLI_CODE    !!---------------------------------------------------------------------- 
    1171 !OLI_CODE    USE ice_oce, ONLY: rdt_ice, tatm_ice 
    1172 !OLI_CODE    USE phycst 
    1173 !OLI_CODE    USE dom_ice 
    1174 !OLI_CODE    USE dom_oce 
    1175 !OLI_CODE    USE sbc_oce 
    1176 !OLI_CODE    USE sbc_ice 
    1177 !OLI_CODE    USE par_ice 
    1178 !OLI_CODE    USE par_oce 
    1179 !OLI_CODE    USE ice 
    1180 !OLI_CODE    USE thd_ice 
    1181 !OLI_CODE    USE in_out_manager 
    1182 !OLI_CODE    USE lbclnk 
    1183 !OLI_CODE    USE lib_mpp 
    1184 !OLI_CODE  
    1185 !OLI_CODE    IMPLICIT NONE 
    1186 !OLI_CODE    PRIVATE 
    1187 !OLI_CODE  
    1188 !OLI_CODE    PUBLIC   lim_mp_init 
    1189 !OLI_CODE    PUBLIC   lim_mp 
    1190 !OLI_CODE  
    1191 !OLI_CODE CONTAINS 
    1192 !OLI_CODE  
    1193 !OLI_CODE    SUBROUTINE lim_mp_init 
    1194 !OLI_CODE       !!------------------------------------------------------------------- 
    1195 !OLI_CODE       !!                ***  ROUTINE lim_mp_init  *** 
    1196 !OLI_CODE       !! 
    1197 !OLI_CODE       !! ** Purpose :   Initialize melt ponds 
    1198 !OLI_CODE       !!------------------------------------------------------------------- 
    1199 !OLI_CODE       a_ip_frac(:,:,:)    = 0._wp 
    1200 !OLI_CODE       a_ip(:,:,:)         = 0._wp 
    1201 !OLI_CODE       h_ip(:,:,:)         = 0._wp 
    1202 !OLI_CODE       v_ip(:,:,:)         = 0._wp 
    1203 !OLI_CODE       h_il(:,:,:)         = 0._wp 
    1204 !OLI_CODE       v_il(:,:,:)         = 0._wp 
    1205 !OLI_CODE          
    1206 !OLI_CODE    END SUBROUTINE lim_mp_init 
    1207 !OLI_CODE  
    1208 !OLI_CODE  
    1209 !OLI_CODE    SUBROUTINE lim_mp 
    1210 !OLI_CODE       !!------------------------------------------------------------------- 
    1211 !OLI_CODE       !!                ***  ROUTINE lim_mp  *** 
    1212 !OLI_CODE       !! 
    1213 !OLI_CODE       !! ** Purpose :   Compute surface heat flux and call main melt pond 
    1214 !OLI_CODE       !!                routine 
    1215 !OLI_CODE       !!------------------------------------------------------------------- 
    1216 !OLI_CODE  
    1217 !OLI_CODE       INTEGER  ::   ji, jj, jl   ! dummy loop indices 
    1218 !OLI_CODE  
    1219 !OLI_CODE       fsurf(:,:) = 0.e0 
    1220 !OLI_CODE       DO jl = 1, jpl 
    1221 !OLI_CODE          DO jj = 1, jpj 
    1222 !OLI_CODE             DO ji = 1, jpi 
    1223 !OLI_CODE                fsurf(ji,jj) = fsurf(ji,jj) + a_i(ji,jj,jl) * & 
    1224 !OLI_CODE                      (qns_ice(ji,jj,jl) + (1.0 - izero(ji,jj,jl)) & 
    1225 !OLI_CODE                        * qsr_ice(ji,jj,jl)) 
    1226 !OLI_CODE             END DO 
    1227 !OLI_CODE          END DO 
    1228 !OLI_CODE       END DO 
    1229 !OLI_CODE  
    1230 !OLI_CODE       CALL compute_mp_topo(at_i, a_i,                               & 
    1231 !OLI_CODE                          vt_i, v_i, v_s, rhosn_glo, t_i, sz_i, a_ip_frac,  & 
    1232 !OLI_CODE                       h_ip, h_il, t_su, tatm_ice, diag_sur_me*rdt_ice, & 
    1233 !OLI_CODE                          fsurf, fwoc) 
    1234 !OLI_CODE  
    1235 !OLI_CODE       at_ip(:,:) = 0.0 
    1236 !OLI_CODE       vt_ip(:,:) = 0.0 
    1237 !OLI_CODE       vt_il(:,:) = 0.0 
    1238 !OLI_CODE       DO jl = 1, jpl 
    1239 !OLI_CODE         DO jj = 1, jpj 
    1240 !OLI_CODE            DO ji = 1, jpi 
    1241 !OLI_CODE               a_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * & 
    1242 !OLI_CODE                                                   a_i(ji,jj,jl)) 
    1243 !OLI_CODE               v_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * & 
    1244 !OLI_CODE                                   a_i(ji,jj,jl) * h_ip(ji,jj,jl)) 
    1245 !OLI_CODE               v_il(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * & 
    1246 !OLI_CODE                                   a_i(ji,jj,jl) * h_il(ji,jj,jl)) 
    1247 !OLI_CODE               at_ip(ji,jj)   = at_ip(ji,jj) + a_ip(ji,jj,jl) 
    1248 !OLI_CODE               vt_ip(ji,jj)   = vt_ip(ji,jj) + v_ip(ji,jj,jl) 
    1249 !OLI_CODE               vt_il(ji,jj)   = vt_il(ji,jj) + v_il(ji,jj,jl) 
    1250 !OLI_CODE            END DO 
    1251 !OLI_CODE         END DO 
    1252 !OLI_CODE       END DO 
    1253 !OLI_CODE  
    1254 !OLI_CODE    END SUBROUTINE lim_mp 
    1255 !OLI_CODE  
    1256 !OLI_CODE  
    1257 !OLI_CODE    SUBROUTINE compute_mp_topo(aice,      aicen,     & 
    1258 !OLI_CODE                               vice,      vicen,     & 
    1259 !OLI_CODE                               vsnon,     rhos,      & 
    1260 !OLI_CODE                               ticen,     salin,     & 
    1261 !OLI_CODE                               a_ip_frac, h_ip,      & 
    1262 !OLI_CODE                               h_il,      Tsfc,      & 
    1263 !OLI_CODE                               potT,      meltt,     & 
    1264 !OLI_CODE                               fsurf,     fwoc) 
    1265 !OLI_CODE       !!------------------------------------------------------------------- 
    1266 !OLI_CODE       !!                ***  ROUTINE compute_mp_topo  *** 
    1267 !OLI_CODE       !! 
    1268 !OLI_CODE       !! ** Purpose :   Compute melt pond evolution based on the ice  
    1269 !OLI_CODE       !!                topography as inferred from the ice thickness  
    1270 !OLI_CODE       !!                distribution.   
    1271 !OLI_CODE       !! 
    1272 !OLI_CODE       !! ** Method  :   This code is initially based on Flocco and Feltham  
    1273 !OLI_CODE       !!                (2007) and Flocco et al. (2010). More to come... 
    1274 !OLI_CODE       !! 
    1275 !OLI_CODE       !! ** Tunable parameters : 
    1276 !OLI_CODE       !!  
    1277 !OLI_CODE       !! ** Note : 
    1278 !OLI_CODE       !! 
    1279 !OLI_CODE       !! ** References 
    1280 !OLI_CODE       !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond  
    1281 !OLI_CODE       !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi:  
    1282 !OLI_CODE       !!    10.1029/2006JC003836. 
    1283 !OLI_CODE       !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of 
    1284 !OLI_CODE       !!    a physically based melt pond scheme into the sea ice component of a 
    1285 !OLI_CODE       !!    climate model.  J. Geophys. Res. 115, C08012,  
    1286 !OLI_CODE       !!    doi: 10.1029/2009JC005568. 
    1287 !OLI_CODE       !!     
    1288 !OLI_CODE       !!------------------------------------------------------------------- 
    1289 !OLI_CODE  
    1290 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), & 
    1291 !OLI_CODE          INTENT(IN) :: & 
    1292 !OLI_CODE          aice, &    ! total ice area fraction 
    1293 !OLI_CODE          vice       ! total ice volume (m) 
    1294 !OLI_CODE  
    1295 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
    1296 !OLI_CODE          INTENT(IN) :: & 
    1297 !OLI_CODE          aicen, &   ! ice area fraction, per category 
    1298 !OLI_CODE          vsnon, &   ! snow volume, per category (m) 
    1299 !OLI_CODE          rhos,  &   ! equivalent snow density, per category (kg/m^3) 
    1300 !OLI_CODE          vicen      ! ice volume, per category (m) 
    1301 !OLI_CODE  
    1302 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & 
    1303 !OLI_CODE          INTENT(IN) :: & 
    1304 !OLI_CODE          ticen, &   ! ice enthalpy, per category 
    1305 !OLI_CODE          salin 
    1306 !OLI_CODE  
    1307 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
    1308 !OLI_CODE          INTENT(INOUT) :: & 
    1309 !OLI_CODE          a_ip_frac , &   ! pond area fraction of ice, per ice category 
    1310 !OLI_CODE          h_ip , &   ! pond depth, per ice category (m) 
    1311 !OLI_CODE          h_il       ! Refrozen ice lid thickness, per ice category (m) 
    1312 !OLI_CODE  
    1313 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), & 
    1314 !OLI_CODE          INTENT(IN) :: & 
    1315 !OLI_CODE          potT,  &   ! air potential temperature 
    1316 !OLI_CODE          meltt, &   ! total surface meltwater flux 
    1317 !OLI_CODE          fsurf      ! thermodynamic heat flux at ice/snow surface (W/m^2) 
    1318 !OLI_CODE  
    1319 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), & 
    1320 !OLI_CODE          INTENT(INOUT) :: & 
    1321 !OLI_CODE          fwoc      ! fresh water flux to the ocean (from draining and other pond volume adjustments) 
    1322 !OLI_CODE                    ! (m) 
    1323 !OLI_CODE  
    1324 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), & 
    1325 !OLI_CODE          INTENT(IN) :: & 
    1326 !OLI_CODE          Tsfc       ! snow/sea ice surface temperature 
    1327 !OLI_CODE  
    1328 !OLI_CODE       ! local variables 
    1329 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
    1330 !OLI_CODE          zTsfcn, & ! ice/snow surface temperature (C) 
    1331 !OLI_CODE          zvolpn, & ! pond volume per unit area, per category (m) 
    1332 !OLI_CODE          zvuin     ! water-equivalent volume of ice lid on melt pond ('upper ice', m) 
    1333 !OLI_CODE  
    1334 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 
    1335 !OLI_CODE          zapondn,& ! pond area fraction, per category 
    1336 !OLI_CODE          zhpondn   ! pond depth, per category (m) 
    1337 !OLI_CODE  
    1338 !OLI_CODE       REAL (wp), DIMENSION (jpi,jpj) :: & 
    1339 !OLI_CODE          zvolp       ! total volume of pond, per unit area of pond (m) 
    1340 !OLI_CODE  
    1341 !OLI_CODE       REAL (wp) :: & 
    1342 !OLI_CODE          zhi,    & ! ice thickness (m) 
    1343 !OLI_CODE          zdHui,  & ! change in thickness of ice lid (m) 
    1344 !OLI_CODE          zomega, & ! conduction 
    1345 !OLI_CODE          zdTice, & ! temperature difference across ice lid (C) 
    1346 !OLI_CODE          zdvice, & ! change in ice volume (m) 
    1347 !OLI_CODE          zTavg,  & ! mean surface temperature across categories (C) 
    1348 !OLI_CODE          zTp,    & ! pond freezing temperature (C) 
    1349 !OLI_CODE          zdvn      ! change in melt pond volume for fresh water budget 
    1350 !OLI_CODE       INTEGER, DIMENSION (jpi*jpj) :: & 
    1351 !OLI_CODE          indxi, indxj    ! compressed indices for cells with ice melting 
    1352 !OLI_CODE  
    1353 !OLI_CODE       INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices 
    1354 !OLI_CODE  
    1355 !OLI_CODE       INTEGER, DIMENSION (jpl) :: & 
    1356 !OLI_CODE          kcells          ! cells where ice lid combines with vice 
    1357 !OLI_CODE  
    1358 !OLI_CODE       INTEGER, DIMENSION (jpi*jpj,jpl) :: & 
    1359 !OLI_CODE          indxii, indxjj  ! i,j indices for kcells loop 
    1360 !OLI_CODE  
    1361 !OLI_CODE       REAL (wp), parameter :: & 
    1362 !OLI_CODE          zhicemin  = 0.1_wp , & ! minimum ice thickness with ponds (m)  
    1363 !OLI_CODE          zTd       = 0.15_wp, & ! temperature difference for freeze-up (C) 
    1364 !OLI_CODE          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) 
    1365 !OLI_CODE          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) 
    1366 !OLI_CODE          z0       = 0._wp,    &  
    1367 !OLI_CODE          zTimelt   = 0._wp,    & 
    1368 !OLI_CODE          z01      = 0.01_wp,  & 
    1369 !OLI_CODE          z25      = 0.25_wp,  & 
    1370 !OLI_CODE          z5       = 0.5_wp,   & 
    1371 !OLI_CODE          epsi10     = 1.0e-11_wp 
    1372 !OLI_CODE       !--------------------------------------------------------------- 
    1373 !OLI_CODE       ! initialize 
    1374 !OLI_CODE       !--------------------------------------------------------------- 
    1375 !OLI_CODE  
    1376 !OLI_CODE       DO j = 1, jpj 
    1377 !OLI_CODE          DO i = 1, jpi 
    1378 !OLI_CODE             zvolp(i,j) = z0 
    1379 !OLI_CODE          END DO 
    1380 !OLI_CODE       END DO 
    1381 !OLI_CODE       DO n = 1, jpl 
    1382 !OLI_CODE          DO j = 1, jpj 
    1383 !OLI_CODE             DO i = 1, jpi 
    1384 !OLI_CODE                ! load tracers 
    1385 !OLI_CODE                zvolp(i,j) = zvolp(i,j) + h_ip(i,j,n) & 
    1386 !OLI_CODE                                      * a_ip_frac(i,j,n) * aicen(i,j,n) 
    1387 !OLI_CODE                zTsfcn(i,j,n) = Tsfc(i,j,n) - rtt ! convert in Celsius - Oli 
    1388 !OLI_CODE                zvuin (i,j,n) = h_il(i,j,n) & 
    1389 !OLI_CODE                             * a_ip_frac(i,j,n) * aicen(i,j,n) 
    1390 !OLI_CODE  
    1391 !OLI_CODE                zhpondn(i,j,n) = z0     ! pond depth, per category 
    1392 !OLI_CODE                zapondn(i,j,n) = z0     ! pond area,  per category 
    1393 !OLI_CODE             END DO 
    1394 !OLI_CODE          END DO 
    1395 !OLI_CODE          indxii(:,n) = 0 
    1396 !OLI_CODE          indxjj(:,n) = 0 
    1397 !OLI_CODE          kcells  (n) = 0 
    1398 !OLI_CODE       END DO 
    1399 !OLI_CODE  
    1400 !OLI_CODE       ! The freezing temperature for meltponds is assumed slightly below 0C, 
    1401 !OLI_CODE       ! as if meltponds had a little salt in them.  The salt budget is not 
    1402 !OLI_CODE       ! altered for meltponds, but if it were then an actual pond freezing  
    1403 !OLI_CODE       ! temperature could be computed. 
    1404 !OLI_CODE  
    1405 !OLI_CODE       zTp = zTimelt - zTd 
    1406 !OLI_CODE  
    1407 !OLI_CODE       !----------------------------------------------------------------- 
    1408 !OLI_CODE       ! Identify grid cells with ponds 
    1409 !OLI_CODE       !----------------------------------------------------------------- 
    1410 !OLI_CODE  
    1411 !OLI_CODE       icells = 0 
    1412 !OLI_CODE       DO j = 1, jpj 
    1413 !OLI_CODE       DO i = 1, jpi 
    1414 !OLI_CODE          zhi = z0 
    1415 !OLI_CODE          IF (aice(i,j) > epsi10) zhi = vice(i,j)/aice(i,j) 
    1416 !OLI_CODE          IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. & 
    1417 !OLI_CODE             zvolp(i,j) > zmin_volp*aice(i,j)) THEN 
    1418 !OLI_CODE             icells = icells + 1 
    1419 !OLI_CODE             indxi(icells) = i 
    1420 !OLI_CODE             indxj(icells) = j 
    1421 !OLI_CODE          ELSE  ! remove ponds on thin ice 
    1422 !OLI_CODE             !fpond(i,j) = fpond(i,j) - zvolp(i,j) 
    1423 !OLI_CODE             zvolpn(i,j,:) = z0 
    1424 !OLI_CODE             zvuin (i,j,:) = z0 
    1425 !OLI_CODE             zvolp (i,j) = z0 
    1426 !OLI_CODE          END IF 
    1427 !OLI_CODE       END DO                     ! i 
    1428 !OLI_CODE       END DO                     ! j 
    1429 !OLI_CODE  
    1430 !OLI_CODE       DO ij = 1, icells 
    1431 !OLI_CODE          i = indxi(ij) 
    1432 !OLI_CODE          j = indxj(ij) 
    1433 !OLI_CODE  
    1434 !OLI_CODE          !-------------------------------------------------------------- 
    1435 !OLI_CODE          ! calculate pond area and depth 
    1436 !OLI_CODE          !-------------------------------------------------------------- 
    1437 !OLI_CODE          CALL pond_area(aice(i,j),vice(i,j),rhos(i,j,:), & 
    1438 !OLI_CODE                    aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), & 
    1439 !OLI_CODE                    ticen(i,j,:,:), salin(i,j,:,:), & 
    1440 !OLI_CODE                    zvolpn(i,j,:), zvolp(i,j), & 
    1441 !OLI_CODE                    zapondn(i,j,:),zhpondn(i,j,:), zdvn) 
    1442 !OLI_CODE  
    1443 !OLI_CODE          fwoc(i,j) = fwoc(i,j) + zdvn ! -> Goes to fresh water budget 
    1444 !OLI_CODE  
    1445 !OLI_CODE          ! mean surface temperature 
    1446 !OLI_CODE          zTavg = z0 
    1447 !OLI_CODE          DO n = 1, jpl 
    1448 !OLI_CODE             zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n) 
    1449 !OLI_CODE          END DO 
    1450 !OLI_CODE          zTavg = zTavg / aice(i,j) 
    1451 !OLI_CODE  
    1452 !OLI_CODE          DO n = 1, jpl-1 
    1453 !OLI_CODE                         
    1454 !OLI_CODE             IF (zvuin(i,j,n) > epsi10) THEN 
    1455 !OLI_CODE  
    1456 !OLI_CODE          !---------------------------------------------------------------- 
    1457 !OLI_CODE          ! melting: floating upper ice layer melts in whole or part 
    1458 !OLI_CODE          !---------------------------------------------------------------- 
    1459 !OLI_CODE !               IF (zTsfcn(i,j,n) > zTp) THEN 
    1460 !OLI_CODE                IF (zTavg > zTp) THEN 
    1461 !OLI_CODE  
    1462 !OLI_CODE                   zdvice = min(meltt(i,j)*zapondn(i,j,n), zvuin(i,j,n)) 
    1463 !OLI_CODE                   IF (zdvice > epsi10) THEN 
    1464 !OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) - zdvice 
    1465 !OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) + zdvice 
    1466 !OLI_CODE                      zvolp (i,j)   = zvolp (i,j)   + zdvice 
    1467 !OLI_CODE                      !fwoc(i,j)   = fwoc(i,j)   + zdvice 
    1468 !OLI_CODE                         
    1469 !OLI_CODE                      IF (zvuin(i,j,n) < epsi10 .and. zvolpn(i,j,n) > puny) THEN 
    1470 !OLI_CODE                         ! ice lid melted and category is pond covered 
    1471 !OLI_CODE                         zvolpn(i,j,n) = zvolpn(i,j,n) + zvuin(i,j,n) 
    1472 !OLI_CODE                         !fwoc(i,j)   = fwoc(i,j)   + zvuin(i,j,n) 
    1473 !OLI_CODE                         zvuin(i,j,n)  = z0 
    1474 !OLI_CODE                      END IF 
    1475 !OLI_CODE                      zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n) 
    1476 !OLI_CODE                   END IF 
    1477 !OLI_CODE  
    1478 !OLI_CODE          !---------------------------------------------------------------- 
    1479 !OLI_CODE          ! freezing: existing upper ice layer grows 
    1480 !OLI_CODE          !---------------------------------------------------------------- 
    1481 !OLI_CODE                ELSE IF (zvolpn(i,j,n) > epsi10) THEN ! zTavg <= zTp 
    1482 !OLI_CODE  
    1483 !OLI_CODE                 ! dIFferential growth of base of surface floating ice layer 
    1484 !OLI_CODE                   zdTice = max(-zTavg, z0) ! > 0 
    1485 !OLI_CODE                   zomega = rcdic*zdTice * zr1_rlfus  
    1486 !OLI_CODE                   zdHui = sqrt(zomega*rdt_ice + z25*(zvuin(i,j,n)/  & 
    1487 !OLI_CODE                         aicen(i,j,n))**2)- z5 * zvuin(i,j,n)/aicen(i,j,n) 
    1488 !OLI_CODE  
    1489 !OLI_CODE                   zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))    
    1490 !OLI_CODE                   IF (zdvice > epsi10) THEN 
    1491 !OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) + zdvice 
    1492 !OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice 
    1493 !OLI_CODE                      zvolp (i,j)   = zvolp (i,j)   - zdvice 
    1494 !OLI_CODE                      !fwoc(i,j)   = fwoc(i,j)   - zdvice 
    1495 !OLI_CODE                      zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n) 
    1496 !OLI_CODE                   END IF 
    1497 !OLI_CODE  
    1498 !OLI_CODE                END IF ! zTavg 
    1499 !OLI_CODE  
    1500 !OLI_CODE          !---------------------------------------------------------------- 
    1501 !OLI_CODE          ! freezing: upper ice layer begins to form 
    1502 !OLI_CODE          ! note: albedo does not change 
    1503 !OLI_CODE          !---------------------------------------------------------------- 
    1504 !OLI_CODE             ELSE ! zvuin < epsi10 
    1505 !OLI_CODE                      
    1506 !OLI_CODE                ! thickness of newly formed ice 
    1507 !OLI_CODE                ! the surface temperature of a meltpond is the same as that 
    1508 !OLI_CODE                ! of the ice underneath (0C), and the thermodynamic surface  
    1509 !OLI_CODE                ! flux is the same 
    1510 !OLI_CODE                zdHui = max(-fsurf(i,j)*rdt_ice*zr1_rlfus, z0) 
    1511 !OLI_CODE                zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))   
    1512 !OLI_CODE                IF (zdvice > epsi10) THEN 
    1513 !OLI_CODE                   zvuin (i,j,n) = zdvice 
    1514 !OLI_CODE                   zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice 
    1515 !OLI_CODE                   zvolp (i,j)   = zvolp (i,j)   - zdvice 
    1516 !OLI_CODE                   !fwoc(i,j)   = fwoc(i,j)   - zdvice 
    1517 !OLI_CODE                   zhpondn(i,j,n)= zvolpn(i,j,n) / zapondn(i,j,n) 
    1518 !OLI_CODE                END IF 
    1519 !OLI_CODE                      
    1520 !OLI_CODE             END IF  ! zvuin 
    1521 !OLI_CODE  
    1522 !OLI_CODE          END DO ! jpl 
    1523 !OLI_CODE  
    1524 !OLI_CODE       END DO ! ij 
    1525 !OLI_CODE  
    1526 !OLI_CODE       !--------------------------------------------------------------- 
    1527 !OLI_CODE       ! remove ice lid if there is no liquid pond 
    1528 !OLI_CODE       ! zvuin may be nonzero on category jpl due to dynamics 
    1529 !OLI_CODE       !--------------------------------------------------------------- 
    1530 !OLI_CODE  
    1531 !OLI_CODE       DO j = 1, jpj 
    1532 !OLI_CODE       DO i = 1, jpi 
    1533 !OLI_CODE          DO n = 1, jpl 
    1534 !OLI_CODE             IF (aicen(i,j,n) > epsi10 .and. zvolpn(i,j,n) < puny & 
    1535 !OLI_CODE                                     .and. zvuin (i,j,n) > epsi10) THEN 
    1536 !OLI_CODE                kcells(n) = kcells(n) + 1 
    1537 !OLI_CODE                indxij    = kcells(n) 
    1538 !OLI_CODE                indxii(indxij,n) = i 
    1539 !OLI_CODE                indxjj(indxij,n) = j 
    1540 !OLI_CODE             END IF 
    1541 !OLI_CODE          END DO 
    1542 !OLI_CODE       END DO                     ! i 
    1543 !OLI_CODE       END DO                     ! j 
    1544 !OLI_CODE  
    1545 !OLI_CODE       DO n = 1, jpl 
    1546 !OLI_CODE  
    1547 !OLI_CODE          IF (kcells(n) > 0) THEN 
    1548 !OLI_CODE          DO ij = 1, kcells(n) 
    1549 !OLI_CODE             i = indxii(ij,n) 
    1550 !OLI_CODE             j = indxjj(ij,n) 
    1551 !OLI_CODE             fwoc(i,j) = fwoc(i,j) + rhoic/rauw * zvuin(i,j,n) ! Completely refrozen lid goes into ocean (to be changed) 
    1552 !OLI_CODE             zvuin(i,j,n) = z0  
    1553 !OLI_CODE          END DO    ! ij 
    1554 !OLI_CODE          END IF 
    1555 !OLI_CODE  
    1556 !OLI_CODE          ! reload tracers 
    1557 !OLI_CODE          DO j = 1, jpj 
    1558 !OLI_CODE             DO i = 1, jpi 
    1559 !OLI_CODE                IF (zapondn(i,j,n) > epsi10) THEN 
    1560 !OLI_CODE                   h_il(i,j,n) = zvuin(i,j,n) / zapondn(i,j,n) 
    1561 !OLI_CODE                ELSE 
    1562 !OLI_CODE                   zvuin(i,j,n) = z0 
    1563 !OLI_CODE                   h_il(i,j,n) = z0 
    1564 !OLI_CODE                END IF 
    1565 !OLI_CODE                IF (aicen(i,j,n) > epsi10) THEN 
    1566 !OLI_CODE                   a_ip_frac(i,j,n) = zapondn(i,j,n) / aicen(i,j,n) * & 
    1567 !OLI_CODE                         (1.0_wp - MAX(z0, SIGN(1.0_wp, -zvolpn(i,j,n)))) 
    1568 !OLI_CODE                   h_ip(i,j,n) = zhpondn(i,j,n) 
    1569 !OLI_CODE                ELSE 
    1570 !OLI_CODE                   a_ip_frac(i,j,n) = z0 
    1571 !OLI_CODE                   h_ip(i,j,n) = z0 
    1572 !OLI_CODE                   h_il(i,j,n) = z0 
    1573 !OLI_CODE                END IF 
    1574 !OLI_CODE             END DO ! i 
    1575 !OLI_CODE          END DO    ! j 
    1576 !OLI_CODE  
    1577 !OLI_CODE       END DO       ! n 
    1578 !OLI_CODE  
    1579 !OLI_CODE    END SUBROUTINE compute_mp_topo 
    1580 !OLI_CODE  
    1581 !OLI_CODE  
    1582 !OLI_CODE    SUBROUTINE pond_area(aice,vice,rhos,             & 
    1583 !OLI_CODE                         aicen, vicen, vsnon, ticen, & 
    1584 !OLI_CODE                         salin, zvolpn, zvolp,         & 
    1585 !OLI_CODE                         zapondn,zhpondn,dvolp) 
    1586 !OLI_CODE       !!------------------------------------------------------------------- 
    1587 !OLI_CODE       !!                ***  ROUTINE pond_area  *** 
    1588 !OLI_CODE       !! 
    1589 !OLI_CODE       !! ** Purpose :   Compute melt pond area, depth and melting rates 
    1590 !OLI_CODE       !!------------------------------------------------------------------ 
    1591 !OLI_CODE       REAL (wp), INTENT(IN) :: & 
    1592 !OLI_CODE          aice,vice 
    1593 !OLI_CODE  
    1594 !OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
    1595 !OLI_CODE          aicen, vicen, vsnon, rhos 
    1596 !OLI_CODE  
    1597 !OLI_CODE       REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 
    1598 !OLI_CODE          ticen, salin 
    1599 !OLI_CODE  
    1600 !OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 
    1601 !OLI_CODE          zvolpn 
    1602 !OLI_CODE  
    1603 !OLI_CODE       REAL (wp), INTENT(INOUT) :: & 
    1604 !OLI_CODE          zvolp, dvolp 
    1605 !OLI_CODE  
    1606 !OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 
    1607 !OLI_CODE          zapondn, zhpondn 
    1608 !OLI_CODE  
    1609 !OLI_CODE       INTEGER :: & 
    1610 !OLI_CODE          n, ns,   & 
    1611 !OLI_CODE          m_index, & 
    1612 !OLI_CODE          permflag 
    1613 !OLI_CODE  
    1614 !OLI_CODE       REAL (wp), DIMENSION(jpl) :: & 
    1615 !OLI_CODE          hicen, & 
    1616 !OLI_CODE          hsnon, & 
    1617 !OLI_CODE          asnon, & 
    1618 !OLI_CODE          alfan, & 
    1619 !OLI_CODE          betan, & 
    1620 !OLI_CODE          cum_max_vol, & 
    1621 !OLI_CODE          reduced_aicen         
    1622 !OLI_CODE  
    1623 !OLI_CODE       REAL (wp), DIMENSION(0:jpl) :: & 
    1624 !OLI_CODE          cum_max_vol_tmp 
    1625 !OLI_CODE  
    1626 !OLI_CODE       REAL (wp) :: & 
    1627 !OLI_CODE          hpond, & 
    1628 !OLI_CODE          drain, & 
    1629 !OLI_CODE          floe_weight, & 
    1630 !OLI_CODE          pressure_head, & 
    1631 !OLI_CODE          hsl_rel, & 
    1632 !OLI_CODE          deltah, & 
    1633 !OLI_CODE          perm, & 
    1634 !OLI_CODE          apond, & 
    1635 !OLI_CODE          msno 
    1636 !OLI_CODE  
    1637 !OLI_CODE       REAL (wp), parameter :: &  
    1638 !OLI_CODE          viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s 
    1639 !OLI_CODE          z0        = 0.0_wp    , & 
    1640 !OLI_CODE          c1        = 1.0_wp    , & 
    1641 !OLI_CODE          p4        = 0.4_wp    , & 
    1642 !OLI_CODE          p6        = 0.6_wp    , & 
    1643 !OLI_CODE          epsi10      = 1.0e-11_wp 
    1644 !OLI_CODE           
    1645 !OLI_CODE      !-----------| 
    1646 !OLI_CODE      !           | 
    1647 !OLI_CODE      !           |-----------| 
    1648 !OLI_CODE      !___________|___________|______________________________________sea-level 
    1649 !OLI_CODE      !           |           | 
    1650 !OLI_CODE      !           |           |---^--------| 
    1651 !OLI_CODE      !           |           |   |        | 
    1652 !OLI_CODE      !           |           |   |        |-----------|              |------- 
    1653 !OLI_CODE      !           |           |   |alfan(n)|           |              | 
    1654 !OLI_CODE      !           |           |   |        |           |--------------| 
    1655 !OLI_CODE      !           |           |   |        |           |              | 
    1656 !OLI_CODE      !---------------------------v------------------------------------------- 
    1657 !OLI_CODE      !           |           |   ^        |           |              | 
    1658 !OLI_CODE      !           |           |   |        |           |--------------| 
    1659 !OLI_CODE      !           |           |   |betan(n)|           |              | 
    1660 !OLI_CODE      !           |           |   |        |-----------|              |------- 
    1661 !OLI_CODE      !           |           |   |        | 
    1662 !OLI_CODE      !           |           |---v------- | 
    1663 !OLI_CODE      !           |           | 
    1664 !OLI_CODE      !           |-----------| 
    1665 !OLI_CODE      !           | 
    1666 !OLI_CODE      !-----------| 
    1667 !OLI_CODE      
    1668 !OLI_CODE       !------------------------------------------------------------------- 
    1669 !OLI_CODE       ! initialize 
    1670 !OLI_CODE       !------------------------------------------------------------------- 
    1671 !OLI_CODE  
    1672 !OLI_CODE       DO n = 1, jpl 
    1673 !OLI_CODE  
    1674 !OLI_CODE          zapondn(n) = z0 
    1675 !OLI_CODE          zhpondn(n) = z0 
    1676 !OLI_CODE  
    1677 !OLI_CODE          IF (aicen(n) < epsi10)  THEN 
    1678 !OLI_CODE             hicen(n) =  z0  
    1679 !OLI_CODE             hsnon(n) = z0 
    1680 !OLI_CODE             reduced_aicen(n) = z0 
    1681 !OLI_CODE          ELSE 
    1682 !OLI_CODE             hicen(n) = vicen(n) / aicen(n) 
    1683 !OLI_CODE             hsnon(n) = vsnon(n) / aicen(n) 
    1684 !OLI_CODE             reduced_aicen(n) = c1 ! n=jpl 
    1685 !OLI_CODE             IF (n < jpl) reduced_aicen(n) = aicen(n) & 
    1686 !OLI_CODE                                  * (-0.024_wp*hicen(n) + 0.832_wp)  
    1687 !OLI_CODE             asnon(n) = reduced_aicen(n)  
    1688 !OLI_CODE          END IF 
    1689 !OLI_CODE  
    1690 !OLI_CODE ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
    1691 !OLI_CODE ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
    1692 !OLI_CODE ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
    1693 !OLI_CODE ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
    1694 !OLI_CODE ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
    1695 !OLI_CODE ! Here, alfa = 60% of the ice area (and since hice is constant in a category,  
    1696 !OLI_CODE ! alfan = 60% of the ice volume) in each category lies above the reference line,  
    1697 !OLI_CODE ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
    1698 !OLI_CODE  
    1699 !OLI_CODE          alfan(n) = p6 * hicen(n) 
    1700 !OLI_CODE          betan(n) = p4 * hicen(n) 
    1701 !OLI_CODE         
    1702 !OLI_CODE          cum_max_vol(n)     = z0 
    1703 !OLI_CODE          cum_max_vol_tmp(n) = z0 
    1704 !OLI_CODE      
    1705 !OLI_CODE       END DO ! jpl 
    1706 !OLI_CODE  
    1707 !OLI_CODE       cum_max_vol_tmp(0) = z0 
    1708 !OLI_CODE       drain = z0 
    1709 !OLI_CODE       dvolp = z0 
    1710 !OLI_CODE      
    1711 !OLI_CODE       !-------------------------------------------------------------------------- 
    1712 !OLI_CODE       ! the maximum amount of water that can be contained up to each ice category 
    1713 !OLI_CODE       !-------------------------------------------------------------------------- 
    1714 !OLI_CODE      
    1715 !OLI_CODE       DO n = 1, jpl-1 ! last category can not hold any volume 
    1716 !OLI_CODE  
    1717 !OLI_CODE          IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN 
    1718 !OLI_CODE  
    1719 !OLI_CODE             ! total volume in level including snow 
    1720 !OLI_CODE             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 
    1721 !OLI_CODE                (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n))  
    1722 !OLI_CODE  
    1723 !OLI_CODE  
    1724 !OLI_CODE             ! subtract snow solid volumes from lower categories in current level 
    1725 !OLI_CODE             DO ns = 1, n  
    1726 !OLI_CODE                cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 
    1727 !OLI_CODE                   - rhos(ns)/rauw  * &    ! fraction of snow that is occupied by solid ??rauw 
    1728 !OLI_CODE                     asnon(ns)  * &    ! area of snow from that category 
    1729 !OLI_CODE                     max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- & 
    1730 !OLI_CODE                                   alfan(n)), z0)   
    1731 !OLI_CODE                                       ! thickness of snow from ns layer in n layer 
    1732 !OLI_CODE             END DO 
    1733 !OLI_CODE  
    1734 !OLI_CODE          ELSE ! assume higher categories unoccupied 
    1735 !OLI_CODE             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 
    1736 !OLI_CODE          END IF 
    1737 !OLI_CODE          !IF (cum_max_vol_tmp(n) < z0) THEN 
    1738 !OLI_CODE          !   call abort_ice('negative melt pond volume') 
    1739 !OLI_CODE          !END IF 
    1740 !OLI_CODE       END DO 
    1741 !OLI_CODE       cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
    1742 !OLI_CODE       cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
    1743 !OLI_CODE      
    1744 !OLI_CODE       !---------------------------------------------------------------- 
    1745 !OLI_CODE       ! is there more meltwater than can be held in the floe? 
    1746 !OLI_CODE       !---------------------------------------------------------------- 
    1747 !OLI_CODE       IF (zvolp >= cum_max_vol(jpl)) THEN 
    1748 !OLI_CODE          drain = zvolp - cum_max_vol(jpl) + epsi10 
    1749 !OLI_CODE          zvolp = zvolp - drain 
    1750 !OLI_CODE          dvolp = drain 
    1751 !OLI_CODE          IF (zvolp < epsi10) THEN 
    1752 !OLI_CODE             dvolp = dvolp + zvolp 
    1753 !OLI_CODE             zvolp = z0 
    1754 !OLI_CODE          END IF 
    1755 !OLI_CODE       END IF 
    1756 !OLI_CODE      
    1757 !OLI_CODE       ! height and area corresponding to the remaining volume 
    1758 !OLI_CODE  
    1759 !OLI_CODE       call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 
    1760 !OLI_CODE            zvolp, cum_max_vol, hpond, m_index) 
    1761 !OLI_CODE      
    1762 !OLI_CODE       DO n=1, m_index 
    1763 !OLI_CODE          zhpondn(n) = hpond - alfan(n) + alfan(1) 
    1764 !OLI_CODE          zapondn(n) = reduced_aicen(n)  
    1765 !OLI_CODE       END DO 
    1766 !OLI_CODE       apond = sum(zapondn(1:m_index)) 
    1767 !OLI_CODE      
    1768 !OLI_CODE       !------------------------------------------------------------------------ 
    1769 !OLI_CODE       ! drainage due to ice permeability - Darcy's law 
    1770 !OLI_CODE       !------------------------------------------------------------------------ 
    1771 !OLI_CODE      
    1772 !OLI_CODE       ! sea water level  
    1773 !OLI_CODE       msno = z0 
    1774 !OLI_CODE       DO n=1,jpl 
    1775 !OLI_CODE         msno = msno + vsnon(n) * rhos(n) 
    1776 !OLI_CODE       END DO 
    1777 !OLI_CODE       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice 
    1778 !OLI_CODE       hsl_rel = floe_weight / rau0 & 
    1779 !OLI_CODE               - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 
    1780 !OLI_CODE      
    1781 !OLI_CODE       deltah = hpond - hsl_rel 
    1782 !OLI_CODE       pressure_head = grav * rau0 * max(deltah, z0) 
    1783 !OLI_CODE  
    1784 !OLI_CODE       ! drain IF ice is permeable     
    1785 !OLI_CODE       permflag = 0 
    1786 !OLI_CODE       IF (pressure_head > z0) THEN 
    1787 !OLI_CODE       DO n = 1, jpl-1 
    1788 !OLI_CODE          IF (hicen(n) /= z0) THEN 
    1789 !OLI_CODE             CALL permeability_phi(ticen(:,n), salin(:,n), vicen(n), perm) 
    1790 !OLI_CODE             IF (perm > z0) permflag = 1 
    1791 !OLI_CODE             drain = perm*zapondn(n)*pressure_head*rdt_ice / & 
    1792 !OLI_CODE                                      (viscosity*hicen(n)) 
    1793 !OLI_CODE             dvolp = dvolp + min(drain, zvolp) 
    1794 !OLI_CODE             zvolp = max(zvolp - drain, z0) 
    1795 !OLI_CODE             IF (zvolp < epsi10) THEN 
    1796 !OLI_CODE                dvolp = dvolp + zvolp 
    1797 !OLI_CODE                zvolp = z0 
    1798 !OLI_CODE             END IF 
    1799 !OLI_CODE          END IF 
    1800 !OLI_CODE       END DO 
    1801 !OLI_CODE   
    1802 !OLI_CODE       ! adjust melt pond DIMENSIONs 
    1803 !OLI_CODE       IF (permflag > 0) THEN 
    1804 !OLI_CODE          ! recompute pond depth     
    1805 !OLI_CODE          CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 
    1806 !OLI_CODE                          zvolp, cum_max_vol, hpond, m_index) 
    1807 !OLI_CODE          DO n=1, m_index 
    1808 !OLI_CODE             zhpondn(n) = hpond - alfan(n) + alfan(1) 
    1809 !OLI_CODE             zapondn(n) = reduced_aicen(n)  
    1810 !OLI_CODE          END DO 
    1811 !OLI_CODE          apond = sum(zapondn(1:m_index)) 
    1812 !OLI_CODE       END IF 
    1813 !OLI_CODE       END IF ! pressure_head 
    1814 !OLI_CODE  
    1815 !OLI_CODE       !------------------------------------------------------------------------ 
    1816 !OLI_CODE       ! total melt pond volume in category DOes not include snow volume 
    1817 !OLI_CODE       ! snow in melt ponds is not melted 
    1818 !OLI_CODE       !------------------------------------------------------------------------ 
    1819 !OLI_CODE  
    1820 !OLI_CODE       ! Calculate pond volume for lower categories 
    1821 !OLI_CODE       DO n=1,m_index-1 
    1822 !OLI_CODE          zvolpn(n) = zapondn(n) * zhpondn(n) & 
    1823 !OLI_CODE                   - (rhos(n)/rauw) * asnon(n) * min(hsnon(n), zhpondn(n))!  
    1824 !OLI_CODE       END DO 
    1825 !OLI_CODE  
    1826 !OLI_CODE       ! Calculate pond volume for highest category = remaining pond volume 
    1827 !OLI_CODE       IF (m_index == 1) zvolpn(m_index) = zvolp 
    1828 !OLI_CODE       IF (m_index > 1) THEN 
    1829 !OLI_CODE         IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 
    1830 !OLI_CODE           zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) 
    1831 !OLI_CODE         ELSE 
    1832 !OLI_CODE           zvolpn(m_index) = z0 
    1833 !OLI_CODE           zhpondn(m_index) = z0 
    1834 !OLI_CODE           zapondn(m_index) = z0 
    1835 !OLI_CODE           ! If remaining pond volume is negative reduce pond volume of  
    1836 !OLI_CODE           ! lower category 
    1837 !OLI_CODE           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) &  
    1838 !OLI_CODE             zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))& 
    1839 !OLI_CODE                                + zvolp 
    1840 !OLI_CODE         END IF 
    1841 !OLI_CODE       END IF 
    1842 !OLI_CODE  
    1843 !OLI_CODE       DO n=1,m_index 
    1844 !OLI_CODE          IF (zapondn(n) > epsi10) THEN 
    1845 !OLI_CODE              zhpondn(n) = zvolpn(n) / zapondn(n) 
    1846 !OLI_CODE          ELSE 
    1847 !OLI_CODE             dvolp = dvolp + zvolpn(n) 
    1848 !OLI_CODE             zhpondn(n) = z0 
    1849 !OLI_CODE             zvolpn(n) = z0 
    1850 !OLI_CODE             zapondn(n) = z0 
    1851 !OLI_CODE          end IF 
    1852 !OLI_CODE       END DO 
    1853 !OLI_CODE       DO n = m_index+1, jpl 
    1854 !OLI_CODE          zhpondn(n) = z0 
    1855 !OLI_CODE          zapondn(n) = z0 
    1856 !OLI_CODE          zvolpn (n) = z0 
    1857 !OLI_CODE       END DO 
    1858 !OLI_CODE  
    1859 !OLI_CODE    END SUBROUTINE pond_area 
    1860 !OLI_CODE    
    1861 !OLI_CODE  
    1862 !OLI_CODE    SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, & 
    1863 !OLI_CODE                          zvolp, cum_max_vol,                & 
    1864 !OLI_CODE                          hpond, m_index) 
    1865 !OLI_CODE       !!------------------------------------------------------------------- 
    1866 !OLI_CODE       !!                ***  ROUTINE calc_hpond  *** 
    1867 !OLI_CODE       !! 
    1868 !OLI_CODE       !! ** Purpose :   Compute melt pond depth 
    1869 !OLI_CODE       !!------------------------------------------------------------------- 
    1870 !OLI_CODE      
    1871 !OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
    1872 !OLI_CODE          aicen, & 
    1873 !OLI_CODE          asnon, & 
    1874 !OLI_CODE          hsnon, & 
    1875 !OLI_CODE          rhos,  & 
    1876 !OLI_CODE          alfan, & 
    1877 !OLI_CODE          cum_max_vol 
    1878 !OLI_CODE      
    1879 !OLI_CODE       REAL (wp), INTENT(IN) :: & 
    1880 !OLI_CODE          zvolp 
    1881 !OLI_CODE      
    1882 !OLI_CODE       REAL (wp), INTENT(OUT) :: & 
    1883 !OLI_CODE          hpond 
    1884 !OLI_CODE      
    1885 !OLI_CODE       INTEGER, INTENT(OUT) :: & 
    1886 !OLI_CODE          m_index 
    1887 !OLI_CODE      
    1888 !OLI_CODE       INTEGER :: n, ns 
    1889 !OLI_CODE      
    1890 !OLI_CODE       REAL (wp), DIMENSION(0:jpl+1) :: & 
    1891 !OLI_CODE          hitl, & 
    1892 !OLI_CODE          aicetl 
    1893 !OLI_CODE      
    1894 !OLI_CODE       REAL (wp) :: & 
    1895 !OLI_CODE          rem_vol, & 
    1896 !OLI_CODE          area, & 
    1897 !OLI_CODE          vol, & 
    1898 !OLI_CODE          tmp, & 
    1899 !OLI_CODE          z0   = 0.0_wp,    &    
    1900 !OLI_CODE          epsi10 = 1.0e-11_wp 
    1901 !OLI_CODE      
    1902 !OLI_CODE       !---------------------------------------------------------------- 
    1903 !OLI_CODE       ! hpond is zero if zvolp is zero - have we fully drained?  
    1904 !OLI_CODE       !---------------------------------------------------------------- 
    1905 !OLI_CODE      
    1906 !OLI_CODE       IF (zvolp < epsi10) THEN 
    1907 !OLI_CODE        hpond = z0 
    1908 !OLI_CODE        m_index = 0 
    1909 !OLI_CODE       ELSE 
    1910 !OLI_CODE         
    1911 !OLI_CODE        !---------------------------------------------------------------- 
    1912 !OLI_CODE        ! Calculate the category where water fills up to  
    1913 !OLI_CODE        !---------------------------------------------------------------- 
    1914 !OLI_CODE         
    1915 !OLI_CODE        !----------| 
    1916 !OLI_CODE        !          | 
    1917 !OLI_CODE        !          | 
    1918 !OLI_CODE        !          |----------|                                     -- -- 
    1919 !OLI_CODE        !__________|__________|_________________________________________ ^ 
    1920 !OLI_CODE        !          |          |             rem_vol     ^                | Semi-filled 
    1921 !OLI_CODE        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
    1922 !OLI_CODE        !          |          |          |              | 
    1923 !OLI_CODE        !          |          |          |              |hpond 
    1924 !OLI_CODE        !          |          |          |----------|   |     |------- 
    1925 !OLI_CODE        !          |          |          |          |   |     | 
    1926 !OLI_CODE        !          |          |          |          |---v-----| 
    1927 !OLI_CODE        !          |          | m_index  |          |         | 
    1928 !OLI_CODE        !------------------------------------------------------------- 
    1929 !OLI_CODE         
    1930 !OLI_CODE        m_index = 0  ! 1:m_index categories have water in them 
    1931 !OLI_CODE        DO n = 1, jpl 
    1932 !OLI_CODE           IF (zvolp <= cum_max_vol(n)) THEN 
    1933 !OLI_CODE              m_index = n 
    1934 !OLI_CODE              IF (n == 1) THEN 
    1935 !OLI_CODE                 rem_vol = zvolp 
    1936 !OLI_CODE              ELSE 
    1937 !OLI_CODE                 rem_vol = zvolp - cum_max_vol(n-1) 
    1938 !OLI_CODE              END IF 
    1939 !OLI_CODE              exit ! to break out of the loop 
    1940 !OLI_CODE           END IF 
    1941 !OLI_CODE        END DO 
    1942 !OLI_CODE        m_index = min(jpl-1, m_index) 
    1943 !OLI_CODE         
    1944 !OLI_CODE        !---------------------------------------------------------------- 
    1945 !OLI_CODE        ! semi-filled layer may have m_index different snow in it 
    1946 !OLI_CODE        !---------------------------------------------------------------- 
    1947 !OLI_CODE         
    1948 !OLI_CODE        !-----------------------------------------------------------  ^ 
    1949 !OLI_CODE        !                                                             |  alfan(m_index+1) 
    1950 !OLI_CODE        !                                                             | 
    1951 !OLI_CODE        !hitl(3)-->                             |----------|          | 
    1952 !OLI_CODE        !hitl(2)-->                |------------| * * * * *|          | 
    1953 !OLI_CODE        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
    1954 !OLI_CODE        !hitl(0)-->-------------------------------------------------  |  ^ 
    1955 !OLI_CODE        !                various snow from lower categories          |  |alfa(m_index) 
    1956 !OLI_CODE         
    1957 !OLI_CODE        ! hitl - heights of the snow layers from thinner and current categories 
    1958 !OLI_CODE        ! aicetl - area of each snow depth in this layer 
    1959 !OLI_CODE         
    1960 !OLI_CODE        hitl(:) = z0 
    1961 !OLI_CODE        aicetl(:) = z0 
    1962 !OLI_CODE        DO n = 1, m_index 
    1963 !OLI_CODE           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
    1964 !OLI_CODE                                  alfan(m_index+1) - alfan(m_index)), z0) 
    1965 !OLI_CODE           aicetl(n) = asnon(n) 
    1966 !OLI_CODE            
    1967 !OLI_CODE           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
    1968 !OLI_CODE        END DO 
    1969 !OLI_CODE        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
    1970 !OLI_CODE        aicetl(m_index+1) = z0 
    1971 !OLI_CODE         
    1972 !OLI_CODE        !---------------------------------------------------------------- 
    1973 !OLI_CODE        ! reorder array according to hitl  
    1974 !OLI_CODE        ! snow heights not necessarily in height order 
    1975 !OLI_CODE        !---------------------------------------------------------------- 
    1976 !OLI_CODE         
    1977 !OLI_CODE        DO ns = 1, m_index+1 
    1978 !OLI_CODE           DO n = 0, m_index - ns + 1 
    1979 !OLI_CODE              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
    1980 !OLI_CODE                 tmp = hitl(n) 
    1981 !OLI_CODE                 hitl(n) = hitl(n+1) 
    1982 !OLI_CODE                 hitl(n+1) = tmp 
    1983 !OLI_CODE                 tmp = aicetl(n) 
    1984 !OLI_CODE                 aicetl(n) = aicetl(n+1) 
    1985 !OLI_CODE                 aicetl(n+1) = tmp 
    1986 !OLI_CODE              END IF 
    1987 !OLI_CODE           END DO 
    1988 !OLI_CODE        END DO 
    1989 !OLI_CODE         
    1990 !OLI_CODE        !---------------------------------------------------------------- 
    1991 !OLI_CODE        ! divide semi-filled layer into set of sublayers each vertically homogenous 
    1992 !OLI_CODE        !---------------------------------------------------------------- 
    1993 !OLI_CODE         
    1994 !OLI_CODE        !hitl(3)---------------------------------------------------------------- 
    1995 !OLI_CODE        !                                                       | * * * * * * * *   
    1996 !OLI_CODE        !                                                       |* * * * * * * * *  
    1997 !OLI_CODE        !hitl(2)---------------------------------------------------------------- 
    1998 !OLI_CODE        !                                    | * * * * * * * *  | * * * * * * * *   
    1999 !OLI_CODE        !                                    |* * * * * * * * * |* * * * * * * * *  
    2000 !OLI_CODE        !hitl(1)---------------------------------------------------------------- 
    2001 !OLI_CODE        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *   
    2002 !OLI_CODE        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *  
    2003 !OLI_CODE        !hitl(0)---------------------------------------------------------------- 
    2004 !OLI_CODE        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)             
    2005 !OLI_CODE         
    2006 !OLI_CODE        ! move up over layers incrementing volume 
    2007 !OLI_CODE        DO n = 1, m_index+1 
    2008 !OLI_CODE            
    2009 !OLI_CODE           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
    2010 !OLI_CODE                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
    2011 !OLI_CODE            
    2012 !OLI_CODE           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
    2013 !OLI_CODE            
    2014 !OLI_CODE           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
    2015 !OLI_CODE              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - & 
    2016 !OLI_CODE                           alfan(1) 
    2017 !OLI_CODE              exit 
    2018 !OLI_CODE           ELSE  ! still in sub-layer below the sub-layer with the depth 
    2019 !OLI_CODE              rem_vol = rem_vol - vol 
    2020 !OLI_CODE           END IF 
    2021 !OLI_CODE            
    2022 !OLI_CODE        END DO 
    2023 !OLI_CODE         
    2024 !OLI_CODE       END IF 
    2025 !OLI_CODE      
    2026 !OLI_CODE    END SUBROUTINE calc_hpond 
    2027 !OLI_CODE    
    2028 !OLI_CODE  
    2029 !OLI_CODE    SUBROUTINE permeability_phi(ticen, salin, vicen, perm) 
    2030 !OLI_CODE       !!------------------------------------------------------------------- 
    2031 !OLI_CODE       !!                ***  ROUTINE permeability_phi  *** 
    2032 !OLI_CODE       !! 
    2033 !OLI_CODE       !! ** Purpose :   Determine the liquid fraction of brine in the ice 
    2034 !OLI_CODE       !!                and its permeability 
    2035 !OLI_CODE       !!------------------------------------------------------------------- 
    2036 !OLI_CODE       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
    2037 !OLI_CODE          ticen,  & ! energy of melting for each ice layer (J/m2) 
    2038 !OLI_CODE          salin 
    2039 !OLI_CODE  
    2040 !OLI_CODE       REAL (wp), INTENT(IN) :: & 
    2041 !OLI_CODE          vicen     ! ice volume 
    2042 !OLI_CODE      
    2043 !OLI_CODE       REAL (wp), INTENT(OUT) :: & 
    2044 !OLI_CODE          perm      ! permeability 
    2045 !OLI_CODE  
    2046 !OLI_CODE       REAL (wp) ::   & 
    2047 !OLI_CODE          Sbr        ! brine salinity 
    2048 !OLI_CODE  
    2049 !OLI_CODE       REAL (wp), DIMENSION(nlay_i) ::   & 
    2050 !OLI_CODE          Tin, &    ! ice temperature 
    2051 !OLI_CODE          phi       ! liquid fraction 
    2052 !OLI_CODE  
    2053 !OLI_CODE       INTEGER :: k 
    2054 !OLI_CODE      
    2055 !OLI_CODE       REAL (wp) :: & 
    2056 !OLI_CODE          c2    = 2.0_wp 
    2057 !OLI_CODE   
    2058 !OLI_CODE       !----------------------------------------------------------------- 
    2059 !OLI_CODE       ! Compute ice temperatures from enthalpies using quadratic formula 
    2060 !OLI_CODE       !----------------------------------------------------------------- 
    2061 !OLI_CODE  
    2062 !OLI_CODE       DO k = 1,nlay_i 
    2063 !OLI_CODE          Tin(k) = ticen(k)  
    2064 !OLI_CODE       END DO 
    2065 !OLI_CODE  
    2066 !OLI_CODE       !----------------------------------------------------------------- 
    2067 !OLI_CODE       ! brine salinity and liquid fraction 
    2068 !OLI_CODE       !----------------------------------------------------------------- 
    2069 !OLI_CODE  
    2070 !OLI_CODE       IF (maxval(Tin-rtt) <= -c2) THEN 
    2071 !OLI_CODE  
    2072 !OLI_CODE          DO k = 1,nlay_i 
    2073 !OLI_CODE             Sbr = - 1.2_wp                 & 
    2074 !OLI_CODE                   -21.8_wp     * (Tin(k)-rtt)    & 
    2075 !OLI_CODE                   - 0.919_wp   * (Tin(k)-rtt)**2 & 
    2076 !OLI_CODE                   - 0.01878_wp * (Tin(k)-rtt)**3 
    2077 !OLI_CODE             phi(k) = salin(k)/Sbr ! liquid fraction 
    2078 !OLI_CODE          END DO ! k 
    2079 !OLI_CODE         
    2080 !OLI_CODE       ELSE 
    2081 !OLI_CODE  
    2082 !OLI_CODE          DO k = 1,nlay_i 
    2083 !OLI_CODE             Sbr = -17.6_wp    * (Tin(k)-rtt)    & 
    2084 !OLI_CODE                   - 0.389_wp  * (Tin(k)-rtt)**2 & 
    2085 !OLI_CODE                   - 0.00362_wp* (Tin(k)-rtt)**3 
    2086 !OLI_CODE             phi(k) = salin(k)/Sbr ! liquid fraction 
    2087 !OLI_CODE          END DO 
    2088 !OLI_CODE  
    2089 !OLI_CODE       END IF 
    2090 !OLI_CODE      
    2091 !OLI_CODE       !----------------------------------------------------------------- 
    2092 !OLI_CODE       ! permeability 
    2093 !OLI_CODE       !----------------------------------------------------------------- 
    2094 !OLI_CODE  
    2095 !OLI_CODE       perm = 3.0e-08_wp * (minval(phi))**3 
    2096 !OLI_CODE      
    2097 !OLI_CODE    END SUBROUTINE permeability_phi 
    2098 !OLI_CODE    
    2099 !OLI_CODE #else 
    2100 !OLI_CODE    !!---------------------------------------------------------------------- 
    2101 !OLI_CODE    !!   Default option         Dummy Module          No LIM-3 sea-ice model 
    2102 !OLI_CODE    !!---------------------------------------------------------------------- 
    2103 !OLI_CODE CONTAINS 
    2104 !OLI_CODE    SUBROUTINE lim_mp_init          ! Empty routine 
    2105 !OLI_CODE    END SUBROUTINE lim_mp_init    
    2106 !OLI_CODE    SUBROUTINE lim_mp               ! Empty routine 
    2107 !OLI_CODE    END SUBROUTINE lim_mp 
    2108 !OLI_CODE    SUBROUTINE compute_mp_topo      ! Empty routine 
    2109 !OLI_CODE    END SUBROUTINE compute_mp_topo         
    2110 !OLI_CODE    SUBROUTINE pond_area            ! Empty routine 
    2111 !OLI_CODE    END SUBROUTINE pond_area    
    2112 !OLI_CODE    SUBROUTINE calc_hpond           ! Empty routine 
    2113 !OLI_CODE    END SUBROUTINE calc_hpond    
    2114 !OLI_CODE    SUBROUTINE permeability_phy     ! Empty routine 
    2115 !OLI_CODE    END SUBROUTINE permeability_phy    
    2116 !OLI_CODE #endif 
    2117 !OLI_CODE    !!====================================================================== 
    2118 !OLI_CODE END MODULE limmp_topo 
    2119 !OLI_CODE  
     263   END SUBROUTINE lim_mp_H12 
     264 
     265   SUBROUTINE lim_mp_init  
     266      !!------------------------------------------------------------------- 
     267      !!                  ***  ROUTINE lim_mp_init   *** 
     268      !! 
     269      !! ** Purpose : Physical constants and parameters linked to melt ponds 
     270      !!      over sea ice 
     271      !! 
     272      !! ** Method  :  Read the nammp  namelist and check the melt pond   
     273      !!       parameter values called at the first timestep (nit000) 
     274      !! 
     275      !! ** input   :   Namelist nammp   
     276      !!------------------------------------------------------------------- 
     277      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read 
     278      NAMELIST/nammp/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
     279      !!------------------------------------------------------------------- 
     280 
     281      REWIND( numnam_ice_ref )              ! Namelist nammp  in reference namelist : Melt Ponds   
     282      READ  ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901) 
     283901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp  in reference namelist', lwp ) 
     284 
     285      REWIND( numnam_ice_cfg )              ! Namelist nammp  in configuration namelist : Melt Ponds 
     286      READ  ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 ) 
     287902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp ) 
     288      IF(lwm) WRITE ( numoni, nammp ) 
     289       
     290      IF(lwp) THEN                        ! control print 
     291         WRITE(numout,*) 
     292         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds' 
     293         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
     294         WRITE(numout,*) '   Namelist namicethd_pnd:' 
     295         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12 
     296         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb 
     297         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST 
     298         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
     299         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd 
     300         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb 
     301      ENDIF 
     302      ! 
     303      !                             !== set the choice of ice pond scheme ==! 
     304      ioptio = 0 
     305                                                            nice_pnd = np_pndNO 
     306      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
     307      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF 
     308      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' ) 
     309 
     310      SELECT CASE( nice_pnd ) 
     311      CASE( np_pndNO )          
     312         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF 
     313         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF 
     314      CASE( np_pndCST) 
     315         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF 
     316      END SELECT 
     317      ! 
     318   END SUBROUTINE lim_mp_init 
     319    
    2120320#else 
    2121321   !!---------------------------------------------------------------------- 
    2122    !!   Default option          Empty module           NO LIM sea-ice model 
    2123    !!---------------------------------------------------------------------- 
    2124 CONTAINS 
    2125    SUBROUTINE lim_mp_init     ! Empty routine 
    2126    END SUBROUTINE lim_mp_init 
    2127    SUBROUTINE lim_mp          ! Empty routine 
    2128    END SUBROUTINE lim_mp      
    2129    SUBROUTINE lim_mp_topo     ! Empty routine 
    2130    END SUBROUTINE lim_mp_topo 
    2131    SUBROUTINE lim_mp_cesm     ! Empty routine 
    2132    END SUBROUTINE lim_mp_cesm 
    2133    SUBROUTINE lim_mp_area     ! Empty routine 
    2134    END SUBROUTINE lim_mp_area 
    2135    SUBROUTINE lim_mp_perm     ! Empty routine 
    2136    END SUBROUTINE lim_mp_perm 
     322   !!   Default option          Empty module          NO ESIM sea-ice model 
     323   !!---------------------------------------------------------------------- 
    2137324#endif  
    2138325 
Note: See TracChangeset for help on using the changeset viewer.