Changeset 12423


Ignore:
Timestamp:
2020-02-20T15:29:25+01:00 (9 months ago)
Author:
dancopsey
Message:

Correct pond melting.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90

    r12394 r12423  
    129129      REAL(wp), PARAMETER ::   zpnd_aspect = 0.174_wp   ! pond aspect ratio 
    130130      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
    131       REAL(wp), PARAMETER ::   freezing_t  = 273.0_wp ! temperature below which refreezing occurs 
    132131      ! 
    133132      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. 
     133      REAL(wp) ::   zdh_mlt          ! available meltwater for melt ponding (equivalent thickness change) 
    137134      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    138135      REAL(wp) ::   z1_rhow          ! inverse freshwater density 
     
    142139      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing 
    143140      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable) 
    144       REAL(wp) ::   actual_frz       ! Amount of melt pond that freezes 
     141      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m) 
    145142      ! 
    146143      INTEGER  ::   ji   ! loop indices 
     
    168165            lh_ip_1d(ji)     = 0._wp 
    169166 
    170             actual_mlt = 0._wp 
    171             actual_frz = 0._wp 
     167            zdh_mlt = 0._wp 
     168            zdh_frz = 0._wp 
    172169            !                                                     !--------------------------------! 
    173170         ELSE                                                     ! Case ice thickness >= rn_himin ! 
     
    175172            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step 
    176173            ! 
    177             ! available meltwater for melt ponding [m, >0] and fraction 
    178             zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 
    179             zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    180             !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
    181             actual_mlt = zfr_mlt * zdv_mlt 
     174            ! available meltwater for melt ponding 
     175            ! This is the change in ice thickness due to melt scaled up by the realive areas of the meltpond 
     176            ! and the area of sea ice feeding the melt ponds. 
     177            zdh_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * (zrmax * a_i_1d(ji)) / a_ip_1d(ji) 
    182178            ! 
    183179            !--- Pond gowth ---! 
    184             v_ip_1d(ji) = v_ip_1d(ji) + actual_mlt * a_i_1d(ji) 
    185             ! 
    186             !--- Lid shrinking ---! 
    187             lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt 
     180            v_ip_1d(ji) = v_ip_1d(ji) + zdh_mlt * a_ip_1d(ji) 
     181            ! 
     182            !--- Lid shrinking. ---! 
     183            lh_ip_1d(ji) = lh_ip_1d(ji) - zdh_mlt 
    188184            ! 
    189185            ! 
    190186            !--- Pond contraction (due to refreezing) ---! 
    191             IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN 
    192                t_grad = freezing_t - t_su_1d(ji) 
     187            IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 
     188               t_grad = (zTp+rt0) - t_su_1d(ji) 
    193189                
    194190               ! The following equation is a rearranged form of: 
     
    201197                              
    202198               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 
    203                actual_frz = lh_ip_end - lh_ip_1d(ji) 
     199               zdh_frz = lh_ip_end - lh_ip_1d(ji) 
    204200 
    205201               ! Pond shrinking 
    206                v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji) 
     202               v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 
    207203 
    208204               ! Lid growing 
    209                lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz 
     205               lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 
    210206            ELSE 
    211                actual_frz = 0._wp 
     207               zdh_frz = 0._wp 
    212208            END IF 
    213209 
    214210            ! melt pond mass flux (<0) 
    215             IF( zdv_mlt > 0._wp ) THEN 
    216                zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice 
     211            IF( zdh_mlt > 0._wp ) THEN 
     212               zfac = zdh_mlt * zrmax * a_i_1d(ji) * rhow * r1_rdtice 
    217213               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    218214               ! 
     
    246242         IF (to_print(ji) == 10) THEN 
    247243            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) 
    248             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 
     244            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 
    249245            write(numout,*)'icethd_pnd: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi 
    250             write(numout,*)'icethd_pnd: lh_ip_1d(ji), actual_mlt, actual_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', actual_mlt, ' ', actual_frz, ' ', t_su_1d(ji) 
     246            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) 
    251247         END IF 
    252248 
Note: See TracChangeset for help on using the changeset viewer.