Changeset 12500


Ignore:
Timestamp:
2020-03-02T17:57:51+01:00 (8 months ago)
Author:
dancopsey
Message:

Add namelist variables to control new functionality.

Location:
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/cfgs/SHARED/namelist_ice_ref

    r11918 r12500  
    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_fix_cpl/doc/namelists/namthd_pnd

    r11536 r12500  
    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_fix_cpl/src/ICE/ice.F90

    r12402 r12500  
    196196   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1) 
    197197   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 
    198205 
    199206   !                                     !!** ice-diagnostics namelist (namdia) ** 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice1d.F90

    r12484 r12500  
    111111   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot      !: Snow accretion/ablation        [m] 
    112112   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sum      !: Ice surface ablation [m] 
    113    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_pnd      !: Ice surface adlation into ponds [m] 
    114    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_pnd      !: Snow surface adlation into ponds [m] 
    115113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_itm      !: Ice internal ablation [m] 
    116114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bom      !: Ice bottom ablation  [m] 
     
    131129   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !: 
    132130   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   lh_ip_1d      !: Ice pond lid thickness   [m] 
    133    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_pnd_avail_1d !: Fraction of sea ice available for melt ponding 
    134131 
    135132   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     
    214211         &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  & 
    215212         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , lh_ip_1d(jpij) ,  & 
    216          &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) , dh_i_pnd(jpij) , a_pnd_avail_1d(jpij) ,           & 
    217          &      dh_s_pnd(jpij) ,                                                                        & 
     213         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   & 
    218214         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
    219215      ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icerst.F90

    r12402 r12500  
    172172      INTEGER           ::   jk 
    173173      LOGICAL           ::   llok 
    174       INTEGER           ::   id0, id1, id2, id3, id4   ! local integer 
     174      INTEGER           ::   id0, id1, id2, id3, id4, id5   ! local integer 
    175175      CHARACTER(len=25) ::   znam 
    176176      CHARACTER(len=2)  ::   zchar, zchar1 
     
    246246            CALL iom_get( numrir, jpdom_autoglo, 'a_ip' , a_ip ) 
    247247            CALL iom_get( numrir, jpdom_autoglo, 'v_ip' , v_ip ) 
    248             CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip) 
    249248         ELSE                                     ! start from rest 
    250249            IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it to zero' 
    251250            a_ip(:,:,:) = 0._wp 
    252251            v_ip(:,:,:) = 0._wp 
     252         ENDIF 
     253         ! melt pond lids 
     254         id5 = iom_varid( numrir, 'lh_ip' , ldstop = .FALSE. ) 
     255         IF( id5 > 0 ) THEN 
     256            CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip) 
     257         ELSE 
    253258            lh_ip(:,:,:) = 0._wp 
    254259         ENDIF 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90

    r12484 r12500  
    220220            dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 
    221221            dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
    222             dh_i_pnd  (1:npti) = 0._wp 
    223             dh_s_pnd  (1:npti) = 0._wp 
    224222            !                                       
    225223                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90

    r12484 r12500  
    2424   USE lib_mpp        ! MPP library 
    2525   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    26    USE icethd_pnd, ONLY : a_pnd_avail 
    2726    
    2827   IMPLICIT NONE 
     
    9998      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness 
    10099      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah 
    101       REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah_pnd  ! Change in thickness due to melting into ponds 
    102100      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting  
    103101 
     
    114112         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile 
    115113      END SELECT 
    116  
    117       ! Initialise fraction of sea ice available for melt ponding 
    118       DO ji = 1, npti 
    119          a_pnd_avail_1d(ji) = a_pnd_avail * at_i_1d(ji) 
    120       END DO 
    121114 
    122115      ! initialize layer thicknesses and enthalpies 
     
    249242               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting 
    250243               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
    251                zdeltah_pnd(ji,jk) = zdeltah  (ji,jk) * a_pnd_avail_1d(ji) 
    252244                
    253245               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0) 
     
    256248               ! updates available heat + thickness 
    257249               dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk) 
    258                dh_s_pnd(ji)    = dh_s_pnd(ji) + zdeltah_pnd(ji,jk)                ! Cumulate surface melt on areas contributing to melt ponds 
    259250               zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    260251               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
     
    352343                
    353344               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
    354                zdeltah_pnd(ji,jk) = - zfmdt * r1_rhoi * a_pnd_avail_1d(ji)                          ! Melt of layer jk [m, <0] into meltponds 
    355345                
    356346               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     
    359349                
    360350               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt 
    361                dh_i_pnd(ji)   = dh_i_pnd(ji) + zdeltah_pnd(ji,jk)      ! Cumulate surface melt on areas contributing to melt ponds 
    362351                
    363352               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90

    r12497 r12500  
    3838   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012) 
    3939 
    40    REAL(wp), PUBLIC, PARAMETER :: a_pnd_avail = 0.7_wp   ! Fraction of sea ice available for melt ponding 
    4140   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp !  pond lid thickness above which the ponds disappear from the albedo calculation 
    4241   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005_wp !  pond lid thickness below which the full pond area is used in the albedo calculation 
     
    137136      REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6  ! Maximum meltpond depth change due to leaking or overflow (m s-1) 
    138137      ! 
     138      REAL(wp) ::   tot_mlt          ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 
    139139      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding 
    140140      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change) 
     
    151151      REAL(wp) ::   dh_ip_over       ! Pond thickness change due to overflow or leaking 
    152152      REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean 
    153       REAL(wp) ::   weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min 
    154153      REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond 
    155154      REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice 
     
    182181            h_ip_1d(ji)      = 0._wp 
    183182            lh_ip_1d(ji)     = 0._wp 
    184  
    185             zdv_mlt = 0._wp 
    186             zdh_frz = 0._wp 
    187183            !                                                     !--------------------------------! 
    188184         ELSE                                                     ! Case ice thickness >= rn_himin ! 
     
    194190            IF ( za_ip < epsi06 ) za_ip = epsi06 
    195191            ! 
    196             ! available meltwater for melt ponding 
    197             ! This is the change in ice volume due to melt 
    198             zdv_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) 
     192            ! available meltwater for melt ponding [m, >0] and fraction  
     193            tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow  
     194            IF ( ln_pnd_totfrac ) THEN 
     195               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction 
     196            ELSE 
     197               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji)  ! Use category ice fraction  
     198            ENDIF  
     199            zdv_mlt = zfr_mlt * tot_mlt  
    199200            ! 
    200201            !--- Pond gowth ---! 
     
    202203            ! 
    203204            !--- Lid shrinking. ---! 
    204             lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 
    205             ! 
     205            IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 
     206            ! 
     207            ! melt pond mass flux (<0) 
     208            IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN 
     209               zfac = zdv_mlt * rhow * r1_rdtice 
     210               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
     211               ! 
     212               ! adjust ice/snow melting flux to balance melt pond flux (>0) 
     213               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 
     214               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
     215               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
     216            ENDIF 
    206217            ! 
    207218            !--- Pond contraction (due to refreezing) ---! 
    208             IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 
    209                t_grad = (zTp+rt0) - t_su_1d(ji) 
    210                 
    211                ! The following equation is a rearranged form of: 
    212                ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 
    213                ! where: lid_thickness_start = lh_ip_1d(ji) 
    214                !        lid_thickness_end = lh_ip_end 
    215                !        omega_dt is a bunch of terms in the equation that do not change 
    216                ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same 
    217                ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification. 
    218                               
    219                lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
    220                zdh_frz = lh_ip_end - lh_ip_1d(ji) 
    221  
    222                ! Pond shrinking 
    223                v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 
    224  
    225                ! Lid growing 
    226                lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 
     219            IF ( ln_pnd_lids ) THEN 
     220 
     221               ! Code to use if using melt pond lids 
     222               IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 
     223                  t_grad = (zTp+rt0) - t_su_1d(ji) 
     224 
     225                  ! The following equation is a rearranged form of: 
     226                  ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 
     227                  ! where: lid_thickness_start = lh_ip_1d(ji) 
     228                  !        lid_thickness_end = lh_ip_end 
     229                  !        omega_dt is a bunch of terms in the equation that do not change 
     230                  ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 
     231                  ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density).  
     232 
     233                  lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
     234                  zdh_frz = lh_ip_end - lh_ip_1d(ji) 
     235 
     236                  ! Pond shrinking 
     237                  v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 
     238 
     239                  ! Lid growing 
     240                  F ( ln_pnd_lids )  lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 
     241               ELSE 
     242                  zdh_frz = 0._wp 
     243               END IF 
     244 
    227245            ELSE 
    228                zdh_frz = 0._wp 
    229             END IF 
     246               ! Code to use if not using melt pond lids 
     247               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 
     248            ENDIF 
    230249 
    231250            ! 
     
    239258            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    240259            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji) 
    241  
    242             ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume 
    243             IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN 
    244                 v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
    245                 dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji)   ! This will be a negative number 
    246                 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
    247                 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
    248                 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
    249                 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
    250                 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
    251  
     260             
     261            !--- Pond overflow and vertical flushing ---! 
     262            IF ( ln_pnd_overflow ) THEN 
     263 
     264               ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume 
     265               IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN 
     266                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     267                   dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji)   ! This will be a negative number 
     268                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
     269                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
     270                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     271                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     272                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     273               ENDIF 
     274 
     275               ! If pond depth exceeds half the ice thickness then reduce the pond volume 
     276               IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 
     277                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     278                   dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji)                ! This will be a negative number 
     279                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
     280                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
     281                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     282                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     283                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     284               ENDIF 
     285 
     286               !-- Vertical flushing of pond water --! 
     287               ! 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. 
     288               ! This assumes the pond is sitting on top of the ice. 
     289               h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 
     290 
     291               ! The depth through which water percolates is the distance under the melt pond 
     292               h_percolation = h_i_1d(ji) - h_ip_1d(ji) 
     293 
     294               ! Calculate the permeability of the ice (Assur 1958) 
     295               DO jk = 1, nlay_i 
     296                   Sbr = - 1.2_wp                         & 
     297                         - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    & 
     298                         - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 & 
     299                         - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3 
     300                   phi(jk) = sz_i_1d(ji,jk)/Sbr 
     301               END DO 
     302               perm = 3.0e-08_wp * (minval(phi))**3 
     303 
     304               ! Do the drainage using Darcy's law 
     305               IF ( perm > 0._wp ) THEN 
     306                   dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number 
     307                   dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative 
     308 
     309                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     310                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
     311                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
     312                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     313                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     314                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     315               ENDIF 
    252316            ENDIF 
    253317 
    254             ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    255             IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 
    256                 v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
    257                 dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji)                ! This will be a negative number 
    258                 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
    259                 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
    260                 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
    261                 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
    262                 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
    263  
     318            ! If lid thickness is ten times greater than pond thickness then remove pond  
     319            IF ( ln_pnd_lids ) THEN            
     320               IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     321                   a_ip_1d(ji)      = 0._wp 
     322                   a_ip_frac_1d(ji) = 0._wp 
     323                   h_ip_1d(ji)      = 0._wp 
     324                   lh_ip_1d(ji)     = 0._wp 
     325                   v_ip_1d(ji)      = 0._wp 
     326               ENDIF 
    264327            ENDIF 
    265  
    266             !-- Vertical flushing of pond water --! 
    267             ! 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. 
    268             ! This assumes the pond is sitting on top of the ice. 
    269             h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 
    270  
    271             ! The depth through which water percolates is the distance under the melt pond 
    272             h_percolation = h_i_1d(ji) - h_ip_1d(ji) 
    273  
    274             ! Calculate the permeability of the ice (Assur 1958) 
    275             DO jk = 1, nlay_i 
    276                 Sbr = - 1.2_wp                         & 
    277                       - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    & 
    278                       - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 & 
    279                       - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3 
    280                 phi(jk) = sz_i_1d(ji,jk)/Sbr 
    281             END DO 
    282             perm = 3.0e-08_wp * (minval(phi))**3 
    283  
    284             ! Do the drainage using Darcy's law 
    285             IF ( perm > 0._wp ) THEN 
    286                 dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number 
    287                 dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative 
    288  
    289                 v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
    290                 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit 
    291                 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 
    292                 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
    293                 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
    294                 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
    295  
    296             ENDIF 
    297  
    298             ! If lid thickness is ten times greater than pond thickness then remove pond             
    299             IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
    300                 a_ip_1d(ji)      = 0._wp 
    301                 a_ip_frac_1d(ji) = 0._wp 
    302                 h_ip_1d(ji)      = 0._wp 
    303                 lh_ip_1d(ji)     = 0._wp 
    304                 v_ip_1d(ji)      = 0._wp 
    305             END IF 
    306328 
    307329            ! If any of the previous changes has removed all the ice thickness then remove ice area. 
     
    333355      INTEGER  ::   ios, ioptio   ! Local integer 
    334356      !! 
    335       NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 
     357      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb,          & 
     358                            rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac,  & 
     359                            ln_use_pndmass 
    336360      !!------------------------------------------------------------------- 
    337361      ! 
     
    351375         WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd 
    352376         WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 
     377         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds       rn_pnd_min = ', rn_pnd_min 
     378         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds       rn_pnd_max = ', rn_pnd_max 
     379         WRITE(numout,*) '            Use total ice fraction instead of category ice fraction   ln_pnd_totfrac  = ',ln_pnd_totfrac 
     380         WRITE(numout,*) '            Allow ponds to overflow and have vertical flushing        ln_pnd_overflow = ',ln_pnd_overflow 
     381         WRITE(numout,*) '            Melt ponds can have frozen lids                           ln_pnd_lids     = ',ln_pnd_lids 
     382         WRITE(numout,*) '            Use melt pond mass flux diagnostic, passing to ocean      ln_use_pndmass  = ',ln_use_pndmass 
    353383         WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST 
    354384         WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd 
Note: See TracChangeset for help on using the changeset viewer.