Changeset 12636


Ignore:
Timestamp:
2020-04-01T09:44:12+02:00 (7 months ago)
Author:
dancopsey
Message:

Add:

  • Melt pond lids
  • Melt pond maximum area and thickness
  • Melt pond vertical flushing
  • Area contributing to melt ponds depends on total ice fraction
Location:
NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/field_def_nemo-ice.xml

    r11575 r12636  
    253253          <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" />  
    254254          <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        />  
     255          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        />  
     256          <field id="icelhpnd_cat" long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" />  
    255257          <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        /> 
    256258          <field id="iceage_cat"   long_name="Ice age per category"                              unit="days"    detect_missing_value="true" /> 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/cfgs/SHARED/namelist_ice_ref

    r11649 r12636  
    177177   ln_pnd           = .false.         !  activate melt ponds or not 
    178178     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
     179       rn_pnd_min   = 0.15            !  Minimum ice fraction that contributes to melt ponds 
     180       rn_pnd_max   = 0.7             !  Maximum ice fraction that contributes to melt ponds 
     181       ln_pnd_totfrac = .false.       !  Use total ice fraction instead of category ice fraction in melt pond contribution fracton 
     182       ln_pnd_overflow = .false.      !  Allow ponds to overflow and have vertical flushing 
     183       ln_pnd_lids  = .false.         !  Melt ponds can have frozen lids 
     184       ln_use_pndmass = .true.        !  Use melt pond mass flux diagnostic, passing it to the ocean for emp 
    179185     ln_pnd_CST     = .false.         !  activate constant  melt ponds 
    180186       rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/doc/namelists/namthd_pnd

    r11536 r12636  
    44   ln_pnd           = .false.         !  activate melt ponds or not 
    55     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
     6       rn_pnd_min   = 0.15            !  Minimum ice fraction that contributes to melt ponds 
     7       rn_pnd_max   = 0.7             !  Maximum ice fraction that contributes to melt ponds 
     8       ln_pnd_totfrac = .false.       !  Use total ice fraction instead of category ice fraction in melt pond contribution fracton 
     9       ln_pnd_overflow = .false.      !  Allow ponds to overflow and have vertical flushing 
     10       ln_pnd_lids  = .false.         !  Melt ponds can have frozen lids 
     11       ln_use_pndmass = .true.        !  Use melt pond mass flux diagnostic, passing it to the ocean for emp 
    612     ln_pnd_CST     = .false.         !  activate constant  melt ponds 
    713       rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice.F90

    r11715 r12636  
    7070   !! a_ip        |      -      |    Ice pond concentration       |       | 
    7171   !! v_ip        |      -      |    Ice pond volume per unit area| m     | 
     72   !! lh_ip       !    lh_ip_1d !    Ice pond lid thickness       ! m     ! 
    7273   !!                                                                     | 
    7374   !!-------------|-------------|---------------------------------|-------| 
     
    195196   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1) 
    196197   LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo 
     198   REAL(wp), PUBLIC ::   rn_pnd_min       !: Minimum ice fraction that contributes to melt ponds 
     199   REAL(wp), PUBLIC ::   rn_pnd_max       !: Maximum ice fraction that contributes to melt ponds 
     200   LOGICAL,  PUBLIC ::   ln_pnd_totfrac   !: Use total ice fraction instead of category ice fraction 
     201                                          !: when determining ice fraction that contributes to melt ponds 
     202   LOGICAL,  PUBLIC ::   ln_pnd_overflow  !: Allow ponds to overflow and have vertical flushing 
     203   LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Melt ponds can have frozen lids 
     204   LOGICAL,  PUBLIC ::   ln_use_pndmass   !: Use melt pond mass flux diagnostic, passing it to the ocean for emp 
    197205 
    198206   !                                     !!** ice-diagnostics namelist (namdia) ** 
     
    331339   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area      [m] 
    332340   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i) 
     341   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_eff   !: melt pond effective fraction (not covered up by lid) (a_ip/a_i) 
    333342   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m] 
     343   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   lh_ip      !: melt pond lid thickness                  [m] 
    334344 
    335345   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration 
     
    448458 
    449459      ii = ii + 1 
    450       ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 
     460      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  & 
     461         &      lh_ip(jpi,jpj,jpl), a_ip_eff(jpi,jpj,jpl) , STAT = ierr(ii) ) 
    451462 
    452463      ii = ii + 1 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice1d.F90

    r11715 r12636  
    128128   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !: 
    129129   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !: 
     130   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_eff_1d   !: 
     131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   lh_ip_1d      !: Ice pond lid thickness   [m] 
    130132 
    131133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     
    157159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d 
    158160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d  
     161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lh_ip_2d  
    159162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d  
    160163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_2d 
     
    208211         &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &     
    209212         &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  & 
    210          &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) ,                   & 
    211          &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   & 
     213         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , lh_ip_1d(jpij) ,  & 
     214         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) , a_ip_eff_1d(jpij) ,                               & 
    212215         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
    213216      ! 
     
    226229      ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  & 
    227230         &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  & 
    228          &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                      & 
     231         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl),  & 
    229232         &      STAT=ierr(ii) ) 
    230233 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icealb.F90

    r11715 r12636  
    9696      LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo 
    9797      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area) 
     98                                                                  !  This is the effective fraction not covered up by a pond lid 
    9899      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth 
    99100      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn.F90

    r11715 r12636  
    101101      ELSEWHERE 
    102102         h_ip(:,:,:) = 0._wp 
     103         lh_ip(:,:,:) = 0._wp 
    103104      END WHERE 
    104105      ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv.F90

    r11715 r12636  
    8484         !                             !-----------------------! 
    8585         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 
    86             &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     86            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i ) 
    8787         !                             !-----------------------! 
    8888      CASE( np_advPRA )                ! PRATHER scheme        ! 
    8989         !                             !-----------------------! 
    9090         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, & 
    91             &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     91            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i ) 
    9292      END SELECT 
    9393 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_pra.F90

    r11715 r12636  
    4444   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction 
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxlhp, sylhp, sxxlhp, syylhp, sxylhp   ! melt pond lid thickness 
    4647 
    4748   !! * Substitutions 
     
    5556 
    5657   SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  & 
    57       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    5859      !!---------------------------------------------------------------------- 
    5960      !!                **  routine ice_dyn_adv_pra  ** 
     
    7879      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    7980      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     81      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    8082      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8183      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    8890      REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw 
    8991      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    90       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     92      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0lhp 
    9193      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    9294      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     
    129131            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    130132            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     133            z0lhp(:,:,jl)  = plh_ip(:,:,jl) * e1e2t(:,:)   ! Melt pond lid thickness 
    131134         ENDIF 
    132135      END DO 
     
    167170               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    168171               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     172               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness 
     173               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 
    169174            ENDIF 
    170175         END DO 
     
    202207               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    203208               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
     209               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness 
     210               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)  
    204211            ENDIF 
    205212         END DO 
     
    225232            pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    226233            pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     234            plh_ip(:,:,jl) = z0lhp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    227235         ENDIF 
    228236      END DO 
     
    231239      !     Remove negative values (conservation is ensured) 
    232240      !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    233       CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     241      CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    234242      ! 
    235243      ! --- Ensure snow load is not too big --- ! 
     
    653661         &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
    654662         &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     663         &      sxlhp(jpi,jpj,jpl) , sylhp (jpi,jpj,jpl), sxxlhp (jpi,jpj,jpl), syylhp (jpi,jpj,jpl), sxylhp (jpi,jpj,jpl),   & 
    655664         ! 
    656665         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & 
     
    765774               CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp ) 
    766775               CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp ) 
     776               !                                                     ! melt pond lid thickness 
     777               CALL iom_get( numrir, jpdom_autoglo, 'sxlhp' , sxlhp  ) 
     778               CALL iom_get( numrir, jpdom_autoglo, 'sylhp' , sylhp  ) 
     779               CALL iom_get( numrir, jpdom_autoglo, 'sxxlhp', sxxlhp ) 
     780               CALL iom_get( numrir, jpdom_autoglo, 'syylhp', syylhp ) 
     781               CALL iom_get( numrir, jpdom_autoglo, 'sxylhp', sxylhp ) 
    767782            ENDIF 
    768783            ! 
     
    782797               sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction 
    783798               sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume 
     799               sxlhp  = 0._wp  ;   sylhp  = 0._wp  ;   sxxlhp  = 0._wp  ;   syylhp  = 0._wp  ;   sxylhp  = 0._wp  ! melt pond lid thickness 
    784800            ENDIF 
    785801         ENDIF 
     
    862878            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 
    863879            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 
     880            !                                                        ! melt pond lid thickness 
     881            CALL iom_rstput( iter, nitrst, numriw, 'sxlhp' , sxlhp  ) 
     882            CALL iom_rstput( iter, nitrst, numriw, 'sylhp' , sylhp  ) 
     883            CALL iom_rstput( iter, nitrst, numriw, 'sxxlhp', sxxlhp ) 
     884            CALL iom_rstput( iter, nitrst, numriw, 'syylhp', syylhp ) 
     885            CALL iom_rstput( iter, nitrst, numriw, 'sxylhp', sxylhp ) 
    864886         ENDIF 
    865887         ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_umx.F90

    r11715 r12636  
    6060 
    6161   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    62       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     62      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    6363      !!---------------------------------------------------------------------- 
    6464      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
     
    8585      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration 
    8686      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     87      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    8788      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8889      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    334335            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
    335336               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
     337            ! lid thickness 
     338            zamsk = 0._wp 
     339            zhvar(:,:,:) = plh_ip(:,:,:) * z1_aip(:,:,:) 
     340            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     341               &                                      zhvar, plh_ip, zua_ups, zva_ups ) 
     342             
    336343         ENDIF 
    337344         ! 
     
    350357         ! Remove negative values (conservation is ensured) 
    351358         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    352          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     359         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    353360         ! 
    354361         ! Make sure ice thickness is not too big 
     
    15191526      !! 
    15201527      !! ** Purpose : Thickness correction in case advection scheme creates 
    1521       !!              abnormally tick ice or snow 
     1528      !!              abnormally thick ice or snow 
    15221529      !! 
    15231530      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_rdgrft.F90

    r11715 r12636  
    503503      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1 
    504504      REAL(wp)                  ::   airft1, oirft1, aprft1 
    505       REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges 
    506       REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice 
     505      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, lhprdg  ! area etc of new ridges 
     506      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft  ! area etc of rafted ice 
    507507      ! 
    508508      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges 
     
    578578                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
    579579                  vprdg (ji) = v_ip_2d(ji,jl1) * afrdg 
     580                  lhprdg(ji) = lh_ip_2d(ji,jl1) * afrdg 
    580581                  aprft1     = a_ip_2d(ji,jl1) * afrft 
    581582                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 
    582583                  vprft (ji) = v_ip_2d(ji,jl1) * afrft 
     584                  lhprft(ji) = lh_ip_2d(ji,jl1) * afrft 
    583585               ENDIF 
    584586 
     
    610612                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    611613                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
     614                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - lhprdg(ji) - lhprft(ji) 
    612615               ENDIF 
    613616            ENDIF 
     
    706709                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         &  
    707710                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   ) 
     711                     lh_ip_2d (ji,jl2) = lh_ip_2d(ji,jl2) + (   lhprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
     712                        &                                   + lhprft (ji) * rn_fpndrft * zswitch(ji)   ) 
    708713                  ENDIF 
    709714                   
     
    736741      !---------------- 
    737742      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    738       CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     743      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d ) 
    739744      ! 
    740745   END SUBROUTINE rdgrft_shift 
     
    848853         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    849854         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     855         CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 
    850856         DO jl = 1, jpl 
    851857            DO jk = 1, nlay_s 
     
    874880         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    875881         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     882         CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) ) 
    876883         DO jl = 1, jpl 
    877884            DO jk = 1, nlay_s 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceistate.F90

    r11715 r12636  
    155155      a_ip     (:,:,:) = 0._wp 
    156156      v_ip     (:,:,:) = 0._wp 
     157      lh_ip    (:,:,:) = 0._wp 
    157158      a_ip_frac(:,:,:) = 0._wp 
     159      a_ip_eff (:,:,:) = 0._wp 
    158160      h_ip     (:,:,:) = 0._wp 
    159161      ! 
     
    362364            a_ip_frac(:,:,:) = 0._wp 
    363365         END WHERE 
     366         a_ip_eff(:,:,:) = a_ip_frac(:,:,:) 
    364367         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    365368           
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceitd.F90

    r11715 r12636  
    411411      CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 
    412412      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
     413      CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 
    413414      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
    414415      DO jl = 1, jpl 
     
    483484                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 
    484485                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
     486                  !                                               
     487                  ztrans          = lh_ip_2d(ji,jl1) * zworka(ji)     ! Pond lid thickness 
     488                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - ztrans 
     489                  lh_ip_2d(ji,jl2) = lh_ip_2d(ji,jl2) + ztrans 
    485490               ENDIF 
    486491               ! 
     
    527532      ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20) 
    528533      !       because of truncation error ( i.e. 1. - 1. /= 0 ) 
    529       CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     534      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d ) 
    530535 
    531536      ! at_i must be <= rn_amax 
     
    555560      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip ) 
    556561      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip ) 
     562      CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip ) 
    557563      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
    558564      DO jl = 1, jpl 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icerst.F90

    r11715 r12636  
    132132      CALL iom_rstput( iter, nitrst, numriw, 'a_ip' , a_ip  ) 
    133133      CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip  ) 
     134      CALL iom_rstput( iter, nitrst, numriw, 'lh_ip', lh_ip ) 
    134135      ! Snow enthalpy 
    135136      DO jk = 1, nlay_s  
     
    171172      INTEGER           ::   jk 
    172173      LOGICAL           ::   llok 
    173       INTEGER           ::   id0, id1, id2, id3, id4   ! local integer 
     174      INTEGER           ::   id0, id1, id2, id3, id4, id5   ! local integer 
    174175      CHARACTER(len=25) ::   znam 
    175176      CHARACTER(len=2)  ::   zchar, zchar1 
     
    250251            v_ip(:,:,:) = 0._wp 
    251252         ENDIF 
     253         a_ip_eff(:,:,:) = a_ip(:,:,:) 
     254         ! melt pond lids 
     255         id5 = iom_varid( numrir, 'lh_ip' , ldstop = .FALSE. ) 
     256         IF( id5 > 0 ) THEN 
     257            CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip) 
     258         ELSE 
     259            lh_ip(:,:,:) = 0._wp 
     260         ENDIF 
    252261         ! fields needed for Met Office (Jules) coupling 
    253262         IF( ln_cpl ) THEN 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icesbc.F90

    r11715 r12636  
    132132 
    133133      ! --- cloud-sky and overcast-sky ice albedos --- ! 
    134       CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) 
     134      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os ) 
    135135 
    136136      ! albedo depends on cloud fraction because of non-linear spectral effects 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd.F90

    r11715 r12636  
    355355         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    356356         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
     357         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 
     358         CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d     (1:npti), lh_ip     (:,:,kl) ) 
    357359         ! 
    358360         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            ) 
     
    461463         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
    462464         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
     465         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) ) 
     466         CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d    (1:npti), lh_ip    (:,:,kl) ) 
    463467         ! 
    464468         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90

    r11715 r12636  
    3838   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012) 
    3939 
     40   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp 
     41 
    4042   !! * Substitutions 
    4143#  include "vectopt_loop_substitute.h90" 
     
    8991         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 
    9092            a_ip_frac_1d(ji) = rn_apnd 
     93            a_ip_eff_1d(ji)  = rn_apnd 
    9194            h_ip_1d(ji)      = rn_hpnd     
    9295            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji) 
    9396         ELSE 
    9497            a_ip_frac_1d(ji) = 0._wp 
     98            a_ip_eff_1d(ji)  = 0._wp 
    9599            h_ip_1d(ji)      = 0._wp     
    96100            a_ip_1d(ji)      = 0._wp 
     
    129133      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio 
    130134      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
    131       ! 
     135      REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6 ! Maximum meltpond depth change due to leaking or overflow (m s-1) 
     136      REAL(wp), PARAMETER ::   pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 
     137      REAL(wp), PARAMETER ::   pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation 
     138      ! 
     139      REAL(wp) ::   tot_mlt          ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 
    132140      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding 
    133       REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding 
     141      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change) 
    134142      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    135143      REAL(wp) ::   z1_rhow          ! inverse freshwater density 
    136144      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio 
     145      REAL(wp) ::   z1_rhoi          ! inverse ice density 
    137146      REAL(wp) ::   zfac, zdum 
    138       ! 
    139       INTEGER  ::   ji   ! loop indices 
     147      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing 
     148      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing 
     149      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable) 
     150      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m) 
     151      REAL(wp) ::   v_ip_old         ! Pond volume before leaking back to the ocean 
     152      REAL(wp) ::   dh_ip_over       ! Pond thickness change due to overflow or leaking 
     153      REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean 
     154      REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond 
     155      REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice 
     156      REAL(wp) ::   Sbr              ! Brine salinity 
     157      REAL(wp), DIMENSION(nlay_i) ::   phi     ! liquid fraction 
     158      REAL(wp) ::   perm             ! Permeability of the sea ice 
     159      REAL(wp) ::   za_ip            ! Temporary pond fraction 
     160      REAL(wp) ::   max_h_diff_ts    ! Maximum meltpond depth change due to leaking or overflow (m per ts) 
     161      REAL(wp) ::   lfrac_pnd        ! The fraction of the meltpond exposed (not inder a frozen lid) 
     162       
     163      ! 
     164      INTEGER  ::   ji, jk   ! loop indices 
    140165      !!------------------------------------------------------------------- 
    141166      z1_rhow        = 1._wp / rhow  
    142167      z1_zpnd_aspect = 1._wp / zpnd_aspect 
    143168      z1_Tp          = 1._wp / zTp  
     169      z1_rhoi        = 1._wp / rhoi 
     170      max_h_diff_ts  = max_h_diff_s * rdt_ice 
     171 
     172      ! Define time-independent field for use in refreezing 
     173      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow) 
    144174 
    145175      DO ji = 1, npti 
    146          !                                                        !--------------------------------! 
    147          IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  ! 
    148             !                                                     !--------------------------------! 
    149             !--- Remove ponds on thin ice 
     176 
     177         !                                                            !----------------------------------------------------! 
     178         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
     179            !                                                         !----------------------------------------------------! 
     180            !--- Remove ponds on thin ice or tiny ice fractions 
    150181            a_ip_1d(ji)      = 0._wp 
    151182            a_ip_frac_1d(ji) = 0._wp 
    152183            h_ip_1d(ji)      = 0._wp 
     184            lh_ip_1d(ji)     = 0._wp 
    153185            !                                                     !--------------------------------! 
    154186         ELSE                                                     ! Case ice thickness >= rn_himin ! 
    155187            !                                                     !--------------------------------! 
    156188            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step 
    157             ! 
    158             ! available meltwater for melt ponding [m, >0] and fraction 
    159             zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    160             zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    161             !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
     189 
     190            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 
     191            za_ip = a_ip_1d(ji) 
     192            IF ( za_ip < epsi06 ) za_ip = epsi06 
     193            ! 
     194            ! available meltwater for melt ponding [m, >0] and fraction  
     195            tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow  
     196            IF ( ln_pnd_totfrac ) THEN 
     197               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction 
     198            ELSE 
     199               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji)  ! Use category ice fraction  
     200            ENDIF  
     201            zdv_mlt = zfr_mlt * tot_mlt  
    162202            ! 
    163203            !--- Pond gowth ---! 
    164             ! v_ip should never be negative, otherwise code crashes 
    165             v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 
     204            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     205            ! 
     206            !--- Lid shrinking. ---! 
     207            IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 
    166208            ! 
    167209            ! melt pond mass flux (<0) 
    168             IF( zdv_mlt > 0._wp ) THEN 
    169                zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice 
     210            IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN 
     211               zfac = zdv_mlt * rhow * r1_rdtice 
    170212               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    171213               ! 
     
    177219            ! 
    178220            !--- Pond contraction (due to refreezing) ---! 
    179             v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 
     221            IF ( ln_pnd_lids ) THEN 
     222 
     223               ! Code to use if using melt pond lids 
     224               IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 
     225                  t_grad = (zTp+rt0) - t_su_1d(ji) 
     226 
     227                  ! The following equation is a rearranged form of: 
     228                  ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 
     229                  ! where: lid_thickness_start = lh_ip_1d(ji) 
     230                  !        lid_thickness_end = lh_ip_end 
     231                  !        omega_dt is a bunch of terms in the equation that do not change 
     232                  ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 
     233                  ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density).  
     234 
     235                  lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
     236                  zdh_frz = lh_ip_end - lh_ip_1d(ji) 
     237 
     238                  ! Pond shrinking 
     239                  v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 
     240 
     241                  ! Lid growing 
     242                  IF ( ln_pnd_lids )  lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 
     243               ELSE 
     244                  zdh_frz = 0._wp 
     245               END IF 
     246 
     247            ELSE 
     248               ! Code to use if not using melt pond lids 
     249               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 
     250            ENDIF 
     251 
     252            ! 
     253            ! Make sure pond volume or lid thickness has not gone negative 
     254            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp  
     255            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 
    180256            ! 
    181257            ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 
     
    184260            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    185261            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji) 
     262             
     263            !--- Pond overflow and vertical flushing ---! 
     264            IF ( ln_pnd_overflow ) THEN 
     265 
     266               ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
     267               IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN 
     268                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     269                   dh_ip_over = zpnd_aspect * zfr_mlt - h_ip_1d(ji)   ! This will be a negative number 
     270                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
     271                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
     272                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     273                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     274                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     275               ENDIF 
     276 
     277               ! If pond depth exceeds half the ice thickness then reduce the pond volume 
     278               IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 
     279                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     280                   dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji)                ! This will be a negative number 
     281                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
     282                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
     283                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     284                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     285                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     286               ENDIF 
     287 
     288               !-- Vertical flushing of pond water --! 
     289               ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth. 
     290               ! This assumes the pond is sitting on top of the ice. 
     291               h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 
     292 
     293               ! The depth through which water percolates is the distance under the melt pond 
     294               h_percolation = h_i_1d(ji) - h_ip_1d(ji) 
     295 
     296               ! Calculate the permeability of the ice (Assur 1958) 
     297               DO jk = 1, nlay_i 
     298                   Sbr = - 1.2_wp                         & 
     299                         - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    & 
     300                         - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 & 
     301                         - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3 
     302                   phi(jk) = sz_i_1d(ji,jk)/Sbr 
     303               END DO 
     304               perm = 3.0e-08_wp * (minval(phi))**3 
     305 
     306               ! Do the drainage using Darcy's law 
     307               IF ( perm > 0._wp ) THEN 
     308                   dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number 
     309                   dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative 
     310 
     311                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     312                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
     313                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
     314                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     315                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     316                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     317               ENDIF 
     318            ENDIF 
     319 
     320            ! If lid thickness is ten times greater than pond thickness then remove pond  
     321            IF ( ln_pnd_lids ) THEN            
     322               IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     323                   a_ip_1d(ji)      = 0._wp 
     324                   a_ip_frac_1d(ji) = 0._wp 
     325                   h_ip_1d(ji)      = 0._wp 
     326                   lh_ip_1d(ji)     = 0._wp 
     327                   v_ip_1d(ji)      = 0._wp 
     328               ENDIF 
     329            ENDIF 
     330 
     331            ! If any of the previous changes has removed all the ice thickness then remove ice area. 
     332            IF ( h_i_1d(ji) == 0._wp ) THEN 
     333                a_i_1d(ji) = 0._wp 
     334                h_s_1d(ji) = 0._wp 
     335            ENDIF 
     336 
    186337            ! 
    187338         ENDIF 
     339 
     340         ! Calculate how much meltpond is exposed to the atmosphere (not under a frozen lid)         
     341         IF ( lh_ip_1d(ji) < pnd_lid_min ) THEN       ! Pond lid is very thin. Expose all the pond. 
     342            lfrac_pnd = 1.0 
     343         ELSE 
     344            IF ( lh_ip_1d(ji) > pnd_lid_max ) THEN    ! Pond lid is very thick. Cover all the pond up with ice and nosw. 
     345              lfrac_pnd = 0.0 
     346            ELSE                                      ! Pond lid is of intermediate thickness. Expose part of the melt pond. 
     347              lfrac_pnd = ( lh_ip_1d(ji) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 
     348            ENDIF 
     349         ENDIF 
     350          
     351         ! Calculate the melt pond effective area 
     352         a_ip_eff_1d(ji) = a_ip_frac_1d(ji) * lfrac_pnd 
     353 
    188354      END DO 
    189355      ! 
     
    205371      INTEGER  ::   ios, ioptio   ! Local integer 
    206372      !! 
    207       NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
     373      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb,          & 
     374                            rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac,  & 
     375                            ln_use_pndmass 
    208376      !!------------------------------------------------------------------- 
    209377      ! 
     
    223391         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd 
    224392         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 
     393         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds       rn_pnd_min = ', rn_pnd_min 
     394         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds       rn_pnd_max = ', rn_pnd_max 
     395         WRITE(numout,*) '            Use total ice fraction instead of category ice fraction   ln_pnd_totfrac  = ',ln_pnd_totfrac 
     396         WRITE(numout,*) '            Allow ponds to overflow and have vertical flushing        ln_pnd_overflow = ',ln_pnd_overflow 
     397         WRITE(numout,*) '            Melt ponds can have frozen lids                           ln_pnd_lids     = ',ln_pnd_lids 
     398         WRITE(numout,*) '            Use melt pond mass flux diagnostic, passing to ocean      ln_use_pndmass  = ',ln_use_pndmass 
    225399         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST 
    226400         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceupdate.F90

    r11715 r12636  
    185185      ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 
    186186      !------------------------------------------------------------------ 
    187       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 
     187      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    188188      ! 
    189189      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90

    r11715 r12636  
    533533               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
    534534               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     535               lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj) 
    535536               ! 
    536537            END DO 
     
    555556 
    556557 
    557    SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     558   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    558559      !!------------------------------------------------------------------- 
    559560      !!                   ***  ROUTINE ice_var_zapneg *** 
     
    570571      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    571572      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     573      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    572574      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    573575      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    637639      WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+) 
    638640      !                                                        but it does not change conservation, so keep it this way is ok 
     641      WHERE( plh_ip (:,:,:) < 0._wp )   plh_ip (:,:,:) = 0._wp 
    639642      ! 
    640643   END SUBROUTINE ice_var_zapneg 
    641644 
    642645 
    643    SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     646   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i ) 
    644647      !!------------------------------------------------------------------- 
    645648      !!                   ***  ROUTINE ice_var_roundoff *** 
     
    654657      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    655658      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     659      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness 
    656660      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    657661      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    668672         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    669673         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
     674         WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 )    plh_ip(1:npti,:)   = 0._wp   ! lh_ip must be >= 0 
    670675      ENDIF 
    671676      ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icewri.F90

    r11715 r12636  
    149149 
    150150      ! --- category-dependent fields --- ! 
     151      IF( iom_use('icelhpnd_cat') )   CALL iom_put( 'icelhpnd_cat', lh_ip          * zmsk00l                                   ) ! melt pond lid thickness for categories 
    151152      IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0% 
    152153      IF( iom_use('iceconc_cat' ) )   CALL iom_put( 'iceconc_cat' , a_i            * zmsk00l                                   ) ! area for categories 
     
    165166      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories 
    166167      IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories 
     168 
     169      IF( iom_use('icevol_cat'  ) )   CALL iom_put( "icevol_cat"  , v_i            * zmsk00l                                   ) ! volume for categories 
     170      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories 
    167171 
    168172      !------------------ 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/BDY/bdyice.F90

    r11715 r12636  
    9898                 &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                & 
    9999                 &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      ) 
     100            CALL lbc_lnk_multi( 'bdyice', lh_ip,'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      ) 
    100101            ! exchange 4d arrays :   third dimension = 1   and then   third dimension = jpk 
    101102            CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
     
    163164            a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration 
    164165            h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth 
     166            lh_ip(ji,jj, jl) = ( lh_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond lid thickness 
    165167            ! 
    166168            sz_i(ji,jj,:,jl) = s_i(ji,jj,jl) 
     
    170172               a_ip(ji,jj,jl) = 0._wp 
    171173               h_ip(ji,jj,jl) = 0._wp 
     174               lh_ip(ji,jj,jl) = 0._wp 
    172175            ENDIF 
    173176            ! 
     
    231234               a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl) 
    232235               h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl) 
     236               lh_ip(ji,jj, jl) = lh_ip(ib,jb, jl) 
    233237               ! 
    234238               sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) 
     
    280284               a_ip(ji,jj,  jl) = 0._wp 
    281285               v_ip(ji,jj,  jl) = 0._wp 
     286               lh_ip(ji,jj, jl) = 0._wp 
    282287               t_su(ji,jj,  jl) = rt0 
    283288               t_s (ji,jj,:,jl) = rt0 
  • NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/OCE/SBC/sbccpl.F90

    r11715 r12636  
    20472047      INTEGER, INTENT(in) ::   kt 
    20482048      ! 
     2049      ! 
    20492050      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    20502051      INTEGER ::   isec, info   ! local integer 
    20512052      REAL(wp) ::   zumax, zvmax 
    20522053      REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 
    2053       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4    
     2054      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4 
     2055 
    20542056      !!---------------------------------------------------------------------- 
    20552057      ! 
     
    22562258            SELECT CASE( sn_snd_mpnd%clcat )   
    22572259            CASE( 'yes' )   
    2258                ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl) 
    2259                ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)   
     2260 
     2261               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl) 
     2262               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl) 
     2263 
    22602264            CASE( 'no' )   
    22612265               ztmp3(:,:,:) = 0.0   
Note: See TracChangeset for help on using the changeset viewer.