Changeset 12484


Ignore:
Timestamp:
2020-02-28T11:57:10+01:00 (8 months ago)
Author:
dancopsey
Message:
  • Make meltponds have a maximum size.
  • Make meltponds leak into the ocean.
  • Make meltpond lids melt at the same rate as the pond underneath.
  • Make meltpond available area proportional to the total sea ice fraction (not category fraction).
Location:
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice1d.F90

    r12402 r12484  
    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] 
    113115   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_itm      !: Ice internal ablation [m] 
    114116   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bom      !: Ice bottom ablation  [m] 
     
    129131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !: 
    130132   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 
    131134 
    132135   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     
    211214         &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  & 
    212215         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , lh_ip_1d(jpij) ,  & 
    213          &      h_ip_1d (jpij) , a_ip_frac_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) ,                                                                        & 
    214218         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
    215219      ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90

    r12402 r12484  
    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 
    222224            !                                       
    223225                              CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- ! 
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90

    r11715 r12484  
    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 
    2627    
    2728   IMPLICIT NONE 
     
    9899      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness 
    99100      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah 
     101      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah_pnd  ! Change in thickness due to melting into ponds 
    100102      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting  
    101103 
     
    112114         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile 
    113115      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 
    114121 
    115122      ! initialize layer thicknesses and enthalpies 
     
    242249               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting 
    243250               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
     251               zdeltah_pnd(ji,jk) = zdeltah  (ji,jk) * a_pnd_avail_1d(ji) 
    244252                
    245253               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) 
     
    248256               ! updates available heat + thickness 
    249257               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 
    250259               zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    251260               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
     
    343352                
    344353               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 
    345355                
    346356               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] 
     
    349359                
    350360               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 
    351362                
    352363               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

    r12402 r12484  
    3737   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme 
    3838   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012) 
     39 
     40   REAL(wp), PUBLIC, PARAMETER :: a_pnd_avail = 0.7_wp   ! Fraction of sea ice available for melt ponding 
     41   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp !  pond lid thickness above which the ponds disappear from the albedo calculation 
     42   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 
     43 
     44   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp 
    3945 
    4046   !! * Substitutions 
     
    129135      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio 
    130136      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
    131       REAL(wp), PARAMETER ::   freezing_t  = 273.0_wp ! temperature below which refreezing occurs 
     137      REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6  ! Maximum meltpond depth change due to leaking or overflow (m s-1) 
    132138      ! 
    133139      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding 
    134       REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding 
    135       REAL(wp) ::   actual_mlt       ! actual melt water used to make melt ponds (per m2). 
    136                                      ! Need to multiply this by ice area to work on volumes. 
     140      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change) 
    137141      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    138142      REAL(wp) ::   z1_rhow          ! inverse freshwater density 
    139143      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio 
     144      REAL(wp) ::   z1_rhoi          ! inverse ice density 
    140145      REAL(wp) ::   zfac, zdum 
    141146      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing 
    142147      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing 
    143148      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable) 
    144       REAL(wp) ::   actual_frz       ! Amount of melt pond that freezes 
    145       ! 
    146       INTEGER  ::   ji   ! loop indices 
     149      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m) 
     150      REAL(wp) ::   v_ip_old         ! Pond volume before leaking back to the ocean 
     151      REAL(wp) ::   dh_ip_over       ! Pond thickness change due to overflow or leaking 
     152      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 
     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       
     162      ! 
     163      INTEGER  ::   ji, jk   ! loop indices 
    147164      !!------------------------------------------------------------------- 
    148165      z1_rhow        = 1._wp / rhow  
    149166      z1_zpnd_aspect = 1._wp / zpnd_aspect 
    150167      z1_Tp          = 1._wp / zTp  
     168      z1_rhoi        = 1._wp / rhoi 
     169      max_h_diff_ts  = max_h_diff_s * rdt_ice 
    151170 
    152171      ! Define time-independent field for use in refreezing 
     
    154173 
    155174      DO ji = 1, npti 
    156          !                                                        !--------------------------------! 
    157          IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  ! 
    158             !                                                     !--------------------------------! 
    159             !--- Remove ponds on thin ice 
     175 
     176         !                                                            !----------------------------------------------------! 
     177         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
     178            !                                                         !----------------------------------------------------! 
     179            !--- Remove ponds on thin ice or tiny ice fractions 
    160180            a_ip_1d(ji)      = 0._wp 
    161181            a_ip_frac_1d(ji) = 0._wp 
     
    163183            lh_ip_1d(ji)     = 0._wp 
    164184 
    165             actual_mlt = 0._wp 
    166             actual_frz = 0._wp 
     185            zdv_mlt = 0._wp 
     186            zdh_frz = 0._wp 
    167187            !                                                     !--------------------------------! 
    168188         ELSE                                                     ! Case ice thickness >= rn_himin ! 
    169189            !                                                     !--------------------------------! 
    170190            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step 
    171             ! 
    172             ! available meltwater for melt ponding [m, >0] and fraction 
    173             zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 
    174             zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    175             !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
    176             actual_mlt = zfr_mlt * zdv_mlt 
     191 
     192            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 
     193            za_ip = a_ip_1d(ji) 
     194            IF ( za_ip < epsi06 ) za_ip = epsi06 
     195            ! 
     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) 
    177199            ! 
    178200            !--- Pond gowth ---! 
    179             v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji) 
    180             ! 
    181             !--- Lid shrinking ---! 
    182             lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt 
     201            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     202            ! 
     203            !--- Lid shrinking. ---! 
     204            lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 
    183205            ! 
    184206            ! 
    185207            !--- Pond contraction (due to refreezing) ---! 
    186             IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN 
    187                t_grad = freezing_t - t_su_1d(ji) 
     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) 
    188210                
    189211               ! The following equation is a rearranged form of: 
     
    196218                              
    197219               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
    198                actual_frz = lh_ip_end - lh_ip_1d(ji) 
     220               zdh_frz = lh_ip_end - lh_ip_1d(ji) 
    199221 
    200222               ! Pond shrinking 
    201                v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji) 
     223               v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 
    202224 
    203225               ! Lid growing 
    204                lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz 
     226               lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 
    205227            ELSE 
    206                actual_frz = 0._wp 
     228               zdh_frz = 0._wp 
    207229            END IF 
    208230 
    209             ! melt pond mass flux (<0) 
    210             IF( zdv_mlt > 0._wp ) THEN 
    211                zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice 
    212                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    213                ! 
    214                ! adjust ice/snow melting flux to balance melt pond flux (>0) 
    215                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 
    216                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    217                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    218             ENDIF 
    219231            ! 
    220232            ! Make sure pond volume or lid thickness has not gone negative 
     
    227239            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    228240            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) = h_ip_1d(ji) + dh_ip_over 
     248                a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji) 
     249                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     250                a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
     251 
     252            ENDIF 
     253 
     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) = 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 
     264            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) = 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 
    229297 
    230298            ! If lid thickness is ten times greater than pond thickness then remove pond             
     
    236304                v_ip_1d(ji)      = 0._wp 
    237305            END IF 
     306 
     307            ! If any of the previous changes has removed all the ice thickness then remove ice area. 
     308            IF ( h_i_1d(ji) == 0._wp ) THEN 
     309                a_i_1d(ji) = 0._wp 
     310                h_s_1d(ji) = 0._wp 
     311            ENDIF 
     312 
    238313            ! 
    239314         ENDIF 
     315 
     316      END DO 
    240317      END DO 
    241318      ! 
Note: See TracChangeset for help on using the changeset viewer.