Changeset 12449


Ignore:
Timestamp:
2020-02-24T18:04:40+01:00 (8 months ago)
Author:
dancopsey
Message:
  • Snow into meltponds contributes to snow to ice diagnostics affecting thickness and salinity of the ice.
  • Avoid divide by zero errors by:
    • Removing meltponds from small ice areas.
    • Removing divide by pond fractions apart from lid shrinking there a minimum pond fraction is imnposed
  • Ponds leaking into the ocean have a salinity flux to the ocean (also a new diagnostic).
  • Add vertical fluxing of pond water.
Location:
NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints
Files:
10 edited

Legend:

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

    r12439 r12449  
    115115          <field id="sfxres"       long_name="ice-ocean salt flux from undiagnosed processes"               unit="kg/m2/s" /> 
    116116          <field id="sfxsub"       long_name="ice-ocean salt flux from ice sublimation"                     unit="kg/m2/s" /> 
     117          <field id="sfxpnd"       long_name="ice-ocean salt flux from meltponds leaking into the ocean"    unit="kg/m2/s" /> 
    117118 
    118119     <!-- mass fluxes --> 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90

    r12439 r12449  
    266266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: salt flux due to correction on ice thick. (residual) [pss.kg.m-2.s-1 => g.m-2.s-1] 
    267267   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
     268   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_pnd     !: salt flux due to melt ponds leaking into the ocean   [pss.kg.m-2.s-1 => g.m-2.s-1] 
    268269 
    269270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth           [W.m-2] 
     
    418419         &      sfx_res    (jpi,jpj) , sfx_bri   (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  & 
    419420         &      sfx_bog    (jpi,jpj) , sfx_bom   (jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,  & 
     421         &      sfx_pnd    (jpi,jpj) ,                                                                 & 
    420422         &      hfx_res    (jpi,jpj) , hfx_snw   (jpi,jpj) , hfx_sub(jpi,jpj) ,                        &  
    421423         &      qt_atm_oi  (jpi,jpj) , qt_oce_ai (jpi,jpj) , fhld   (jpi,jpj) ,                        & 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90

    r12439 r12449  
    9393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_lam_1d 
    9494   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_dyn_1d 
     95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_pnd_1d  
    9596 
    9697   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d 
     
    207208         &      sfx_bri_1d    (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  & 
    208209         &      sfx_sni_1d    (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  & 
    209          &      sfx_lam_1d    (jpij) , sfx_dyn_1d(jpij)  , STAT=ierr(ii) ) 
     210         &      sfx_lam_1d    (jpij) , sfx_dyn_1d(jpij)  , sfx_pnd_1d (jpij) , STAT=ierr(ii) ) 
    210211      ! 
    211212      ii = ii + 1 
     
    216217         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , lh_ip_1d(jpij) ,  & 
    217218         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) , dh_i_pnd(jpij) , a_pnd_avail_1d(jpij) ,           & 
     219         &      dh_s_pnd(jpij) ,                                                                        & 
    218220         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
    219221      ! 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icectl.F90

    r11081 r12449  
    8787         pdiag_fs = glob_sum( 'icectl',                                                                     & 
    8888            &                  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    89             &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
    90             &                  ) *  e1e2t(:,:) ) * zconv  
     89            &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) +  & 
     90            &                    sfx_pnd(:,:) ) *  e1e2t(:,:) ) * zconv  
    9191         ! 
    9292         !                          ! heat flux 
     
    116116         zfs = glob_sum( 'icectl',                                                                       & 
    117117            &              ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
    118             &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &  
    119             &              ) * e1e2t(:,:) ) * zconv - pdiag_fs 
     118            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) +  &  
     119            &                sfx_pnd(:,:) ) * e1e2t(:,:) ) * zconv - pdiag_fs 
    120120 
    121121         ! heat flux 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90

    r12369 r12449  
    486486      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    487487      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp 
     488      sfx_pnd(:,:) = 0._wp 
    488489      ! 
    489490      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90

    r12439 r12449  
    455455         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub          ) 
    456456         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam          ) 
     457         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd          ) 
    457458         ! 
    458459         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd       ) 
     
    548549         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub        ) 
    549550         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam        ) 
     551         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd        ) 
    550552         ! 
    551553         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd     ) 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90

    r12439 r12449  
    9999      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness 
    100100      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah 
     101      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah_pnd  ! Change in thickness due to melting into ponds 
    101102      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting  
    102103 
     
    281282               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting 
    282283               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
     284               zdeltah_pnd(ji,jk) = zdeltah  (ji,jk) * a_pnd_avail_1d(ji) 
    283285                
    284286               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) 
    285                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
     287               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * (zdeltah(ji,jk) - zdeltah_pnd(ji,jk)) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
    286288                
    287289               ! updates available heat + thickness 
    288290               dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk) 
    289                dh_s_pnd(ji)    = dh_s_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji)                   ! Cumulate surface melt on areas contributing to melt ponds 
     291               dh_s_pnd(ji)    = dh_s_pnd(ji) + zdeltah_pnd(ji,jk)                ! Cumulate surface melt on areas contributing to melt ponds 
    290292               zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    291293               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
    292                h_i_1d  (ji)    = MAX( 0._wp , h_i_1d(ji) - zdeltah(ji,jk) * a_pnd_avail_1d(ji) * rhos * r1_rhoi )   ! Snow melt trapped in ponds contrbutes to ice thickness 
    293294               zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 
     295 
     296               ! Snow melt trapped in ponds contrbutes to ice thickness and snow to ice conversion 
     297               h_i_1d  (ji)    = MAX( 0._wp , h_i_1d(ji) - zdeltah_pnd(ji,jk) * rhos * r1_rhoi )  
     298               wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice 
     299               wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice 
    294300               ! 
    295301            ENDIF 
     
    400406               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
    401407                
    402                zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0] 
    403                 
    404                zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
     408               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean or meltponds [kg/m2, >0] 
     409                
     410               zdeltah(ji,jk) = - zfmdt * r1_rhoi * (1._wp - a_pnd_avail_1d(ji))                    ! Melt of layer jk [m, <0] into the ocean 
     411               zdeltah_pnd(ji,jk) = - zfmdt * r1_rhoi * a_pnd_avail_1d(ji)                          ! Melt of layer jk [m, <0] into meltponds 
    405412                
    406413               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] 
     
    408415               zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 
    409416                
    410                dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk) * (1._wp - a_pnd_avail_1d(ji))         ! Cumulate surface melt on areas not contributing to melt ponds 
    411                dh_i_pnd(ji)   = dh_i_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji)                   ! Cumulate surface melt on areas contributing to melt ponds 
     417               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)          ! Cumulate surface melt on areas not contributing to melt ponds 
     418               dh_i_pnd(ji)   = dh_i_pnd(ji) + zdeltah_pnd(ji,jk)      ! Cumulate surface melt on areas contributing to melt ponds 
    412419                
    413420               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90

    r12439 r12449  
    3939 
    4040   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    !  pond lid thickness above which the ponds disappear from the albedo calculation 
    42    REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005    !  pond lid thickness below which the full pond area is used in the albedo calculation 
    43  
    44    REAL(wp), PARAMETER :: snow_lid_min = 0.15  ! Snow thickness above which form a lid of size pnd_lid_min on the melt ponds 
    45    REAL(wp), PARAMETER :: snow_lid_max = 0.25  ! Snow thickness above which form a lid of size pnd_lid_max on the melt ponds 
     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 
    4645 
    4746   !! * Substitutions 
     
    138137      ! 
    139138      REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding 
    140       REAL(wp) ::   zdh_mlt          ! available meltwater for melt ponding (equivalent thickness change) 
     139      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change) 
    141140      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    142141      REAL(wp) ::   z1_rhow          ! inverse freshwater density 
     
    151150      REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean 
    152151      REAL(wp) ::   weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min 
    153       ! 
    154       INTEGER  ::   ji   ! loop indices 
     152      REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond 
     153      REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice 
     154      REAL(wp) ::   Sbr              ! Brine salinity 
     155      REAL(wp), DIMENSION(nlay_i) ::   phi     ! liquid fraction 
     156      REAL(wp) ::   perm             ! Permeability of the sea ice 
     157      REAL(wp) ::   drain            ! Amount of melt pond draining the sea ice per m2 
     158      REAL(wp) ::   za_ip            ! Temporary pond fraction 
     159  
     160      ! 
     161      INTEGER  ::   ji, jk   ! loop indices 
    155162      !!------------------------------------------------------------------- 
    156163      z1_rhow        = 1._wp / rhow  
     
    168175         END IF 
    169176 
    170          !                                                        !--------------------------------! 
    171          IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin ! 
    172             !                                                     !--------------------------------! 
    173             !--- Remove ponds on thin ice 
     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 
    174181            a_ip_1d(ji)      = 0._wp 
    175182            a_ip_frac_1d(ji) = 0._wp 
     
    177184            lh_ip_1d(ji)     = 0._wp 
    178185 
    179             zdh_mlt = 0._wp 
     186            zdv_mlt = 0._wp 
    180187            zdh_frz = 0._wp 
    181188            !                                                     !--------------------------------! 
     
    183190            !                                                     !--------------------------------! 
    184191            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step 
     192 
     193            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 
     194            za_ip = a_ip_1d(ji) 
     195            IF ( za_ip < epsi06 ) za_ip = epsi06 
    185196            ! 
    186197            ! available meltwater for melt ponding 
    187             ! This is the change in ice thickness due to melt scaled up by the realive areas of the meltpond 
    188             ! and the area of sea ice feeding the melt ponds. 
    189             zdh_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) / a_ip_1d(ji) 
     198            ! This is the change in ice volume due to melt 
     199            zdv_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) 
    190200            ! 
    191201            !--- Pond gowth ---! 
    192             v_ip_1d(ji) = v_ip_1d(ji) + zdh_mlt * a_ip_1d(ji) 
     202            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
    193203            ! 
    194204            !--- Lid shrinking. ---! 
    195             lh_ip_1d(ji) = lh_ip_1d(ji) - zdh_mlt 
     205            lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 
    196206            ! 
    197207            ! 
     
    231241            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji) 
    232242 
     243            ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume 
     244            IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN 
     245                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     246                h_ip_1d(ji) = zpnd_aspect * a_pnd_avail_1d(ji) 
     247                a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji) 
     248                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     249                a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
     250 
     251                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean. 
     252                ! Scale this up to make water depth meaned over sea ice. 
     253                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number 
     254 
     255                ! Output this water loss as a mass flux diagnostic 
     256                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice 
     257 
     258                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference 
     259                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi ) 
     260 
     261                ! Pass this as a salinity flux to the ocean 
     262                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice 
     263            ENDIF 
     264 
     265            ! If pond depth exceeds half the ice thickness then reduce the pond volume 
     266            IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 
     267                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     268                h_ip_1d(ji) = 0.5_wp * h_i_1d(ji) 
     269                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     270                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     271                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     272 
     273                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean. 
     274                ! Scale this up to make water depth meaned over sea ice. 
     275                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number 
     276 
     277                ! Output this water loss as a mass flux diagnostic 
     278                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice 
     279 
     280                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference 
     281                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi ) 
     282 
     283                ! Pass this as a salinity flux to the ocean 
     284                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice 
     285            ENDIF 
     286 
     287            !-- Vertical flushing of pond water --! 
     288            ! 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. 
     289            ! This assumes the pond is sitting on top of the ice. 
     290            h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 
     291 
     292            ! The depth through which water percolates is the distance under the melt pond 
     293            h_percolation = h_i_1d(ji) - h_ip_1d(ji) 
     294 
     295            ! Calculate the permeability of the ice (Assur 1958) 
     296            DO jk = 1, nlay_i 
     297                Sbr = - 1.2_wp                         & 
     298                      - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    & 
     299                      - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 & 
     300                      - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3 
     301                phi(jk) = sz_i_1d(ji,jk)/Sbr 
     302            END DO 
     303            perm = 3.0e-08_wp * (minval(phi))**3 
     304 
     305            ! Do the drainage using Darcy's law 
     306            IF ( perm > 0._wp ) THEN 
     307                drain = perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) 
     308 
     309                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
     310                h_ip_1d(ji) = h_ip_1d(ji) - drain 
     311                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
     312                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
     313                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
     314 
     315                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean. 
     316                ! Scale this up to make water depth meaned over sea ice. 
     317                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number 
     318 
     319                ! Output this water loss as a mass flux diagnostic 
     320                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice 
     321 
     322                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference 
     323                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi ) 
     324 
     325                ! Pass this as a salinity flux to the ocean 
     326                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice 
     327            ENDIF 
     328 
    233329            ! If lid thickness is ten times greater than pond thickness then remove pond             
    234330            IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     
    240336            END IF 
    241337 
    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                 h_ip_1d(ji) = zpnd_aspect * a_pnd_avail_1d(ji) 
    246                 a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji) 
    247                 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
    248                 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    249  
    250                 ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean. 
    251                 ! Scale this up to make water depth meaned over sea ice. 
    252                 dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number 
    253  
    254                 ! Output this water loss as a mass flux diagnostic 
    255                 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice 
    256  
    257                 ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference 
    258                 h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi ) 
    259             ENDIF 
    260  
    261             ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    262             IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 
    263                 v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use 
    264                 h_ip_1d(ji) = 0.5_wp * h_i_1d(ji) 
    265                 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 
    266                 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 
    267                 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 
    268  
    269                 ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean. 
    270                 ! Scale this up to make water depth meaned over sea ice. 
    271                 dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number 
    272  
    273                 ! Output this water loss as a mass flux diagnostic 
    274                 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice 
    275  
    276                 ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference 
    277                 h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi ) 
    278             ENDIF 
    279  
    280             ! If any of the previous two IF tests has removed all the ice thickness then remove ice area. 
     338            ! If any of the previous changes has removed all the ice thickness then remove ice area. 
    281339            IF ( h_i_1d(ji) == 0._wp ) THEN 
    282340                a_i_1d(ji) = 0._wp 
     
    284342            ENDIF 
    285343 
    286             ! If snow thickness exceeds snow_lid_min then form a very thin lid (so snow can go over the top) 
    287             IF ( h_s_1d(ji) > snow_lid_min .AND. h_s_1d(ji) < snow_lid_max ) THEN 
    288                weighted_lid_snow = (snow_lid_max - h_s_1d(ji))/(snow_lid_max-snow_lid_min) * pnd_lid_min +             & 
    289                                    (h_s_1d(ji) - snow_lid_min)/(snow_lid_max-snow_lid_min) * pnd_lid_max 
    290                lh_ip_1d(ji) = MAX(lh_ip_1d(ji), weighted_lid_snow) 
    291             ENDIF 
    292             IF ( h_s_1d(ji) >= snow_lid_max ) THEN 
    293                lh_ip_1d(ji) = MAX(lh_ip_1d(ji), pnd_lid_max) 
    294             ENDIF 
    295  
    296344            ! 
    297345         ENDIF 
     
    299347         IF (to_print(ji) == 10) THEN 
    300348            write(numout,*)'icethd_pnd: h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji) 
    301             write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdh_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdh_mlt 
     349            write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdv_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdv_mlt 
    302350            write(numout,*)'icethd_pnd: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi 
    303             write(numout,*)'icethd_pnd: lh_ip_1d(ji), zdh_mlt, zdh_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', zdh_mlt, ' ', zdh_frz, ' ', t_su_1d(ji) 
     351            write(numout,*)'icethd_pnd: lh_ip_1d(ji), sfx_pnd_1d(ji), zdh_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', sfx_pnd_1d(ji), ' ', zdh_frz, ' ', t_su_1d(ji) 
    304352            write(numout,*)'icethd_pnd: a_pnd_avail_1d(ji), at_i_1d(ji), wfx_pnd_1d(ji), h_i_1d(ji) = ', a_pnd_avail_1d(ji), ' ', at_i_1d(ji), ' ', wfx_pnd_1d(ji), ' ', h_i_1d(ji) 
     353            write(numout,*)'icethd_pnd: drain, perm, h_percolation, h_gravity_head = ',drain, ' ', perm, ' ', h_percolation, ' ', h_gravity_head 
     354            write(numout,*)'icethd_pnd: t_i_1d(ji,1), sz_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', sz_i_1d(ji,1) 
    305355         END IF 
    306356 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_sal.F90

    r10888 r12449  
    5959      INTEGER  ::   ji, jk                       ! dummy loop indices  
    6060      REAL(wp) ::   iflush, igravdr              ! local scalars 
    61       REAL(wp) ::   zs_sni, zs_i_gd, zs_i_fl, zs_i_si, zs_i_bg   ! local scalars 
     61      REAL(wp) ::   zs_sni, zs_i_gd, zs_i_fl, zs_i_si, zs_i_bg, zs_i_sp   ! local scalars 
    6262      REAL(wp) ::   z1_time_gd, z1_time_fl 
    6363      !!--------------------------------------------------------------------- 
     
    7979               zs_sni  = sss_1d(ji) * ( rhoi - rhos ) * r1_rhoi                     ! Salinity of snow ice 
    8080               zs_i_si = ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice     
     81               zs_i_sp = ( 0._wp       - s_i_1d(ji) ) * dh_s_pnd(ji)   / h_i_1d(ji) ! snow-pond 
    8182               zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth 
    8283               ! Update salinity (nb: salt flux already included in icethd_dh) 
    83                s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si 
     84               s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si + zs_i_sp 
    8485            ENDIF 
    8586            ! 
     
    100101               sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_rdtice 
    101102            ENDIF 
     103            ! 
    102104         END DO 
    103105         ! 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90

    r12394 r12449  
    165165            !------------------------------------------ 
    166166            sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   & 
    167                &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) 
     167               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) + sfx_pnd(ji,jj) 
    168168             
    169169            ! Mass of snow and ice per unit area    
     
    197197      ! 
    198198      ! --- salt fluxes [kg/m2/s] --- ! 
    199       !                           ! sfxice =  sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam 
     199      !                           ! sfxice =  sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam + sfxpnd 
    200200      IF( iom_use('sfxice'  ) )   CALL iom_put( "sfxice", sfx     * 1.e-03 )   ! salt flux from total ice growth/melt 
    201201      IF( iom_use('sfxbog'  ) )   CALL iom_put( "sfxbog", sfx_bog * 1.e-03 )   ! salt flux from bottom growth 
     
    209209      IF( iom_use('sfxres'  ) )   CALL iom_put( "sfxres", sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes 
    210210      IF( iom_use('sfxsub'  ) )   CALL iom_put( "sfxsub", sfx_sub * 1.e-03 )   ! salt flux from sublimation 
     211      IF( iom_use('sfxpnd'  ) )   CALL iom_put( "sfxpnd", sfx_pnd * 1.e-03 )   ! salt flux from melt ponds leaking into the ocean 
    211212 
    212213      ! --- mass fluxes [kg/m2/s] --- ! 
Note: See TracChangeset for help on using the changeset viewer.