New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12725 for NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_pnd.F90 – NEMO

Ignore:
Timestamp:
2020-04-09T15:48:28+02:00 (4 years ago)
Author:
clem
Message:

fix previous revision. Make sure a_ip and a_ip_eff are correct

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_pnd.F90

    r12720 r12725  
    8989         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 
    9090            a_ip_frac_1d(ji) = rn_apnd 
    91             a_ip_eff_1d(ji)  = rn_apnd 
    9291            h_ip_1d(ji)      = rn_hpnd     
    9392            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji) 
     
    9594         ELSE 
    9695            a_ip_frac_1d(ji) = 0._wp 
    97             a_ip_eff_1d(ji)  = 0._wp 
    9896            h_ip_1d(ji)      = 0._wp     
    9997            a_ip_1d(ji)      = 0._wp 
     
    139137      !!                    if no lids:   Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)                      --- from Holland et al 2012 --- 
    140138      !! 
    141       !!              - Overflow:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 --- 
     139      !!              - Flushing:         w = -perm/visc * rho_oce * grav * Hp / Hi                 --- from Flocco et al 2007 --- 
    142140      !!                                     perm = permability of sea-ice 
    143141      !!                                     visc = water viscosity 
     
    145143      !!                                     Hi   = ice thickness thru which there is flushing 
    146144      !! 
    147       !! 
    148145      !!              - Corrections:      remove melt ponds when lid thickness is 10 times the pond thickness 
    149       !! 
    150       !!              - effective pond area: to be used for albedo  
    151146      !! 
    152147      !!              - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip: 
     
    161156      !!                   Holland et al      (J. Clim, 2012) 
    162157      !!------------------------------------------------------------------- 
    163       REAL(wp), DIMENSION(nlay_i) ::   zperm          ! Permeability of sea ice 
     158      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
    164159      !! 
    165160      REAL(wp), PARAMETER ::   zaspect =  0.8_wp      ! pond aspect ratio 
    166161      REAL(wp), PARAMETER ::   zTp     = -2._wp       ! reference temperature 
    167162      REAL(wp), PARAMETER ::   zvisc   =  1.79e-3_wp  ! water viscosity 
    168       REAL(wp), PARAMETER ::   zhl_max =  0.015_wp    ! pond lid thickness above which the ponds disappear from the albedo calculation 
    169       REAL(wp), PARAMETER ::   zhl_min =  0.005_wp    ! pond lid thickness below which the full pond area is used in the albedo calculation 
    170163      !! 
    171164      REAL(wp) ::   zfr_mlt, zdv_mlt                  ! fraction and volume of available meltwater retained for melt ponding 
    172165      REAL(wp) ::   zdv_frz, zdv_flush                ! Amount of melt pond that freezes, flushes 
    173166      REAL(wp) ::   zhp                               ! heigh of top of pond lid wrt ssh 
    174       REAL(wp) ::   v_ip_max                          ! max pond volume allowed 
     167      REAL(wp) ::   zv_ip_max                         ! max pond volume allowed 
    175168      REAL(wp) ::   zdT                               ! zTp-t_su 
    176169      REAL(wp) ::   zsbr                              ! Brine salinity 
     170      REAL(wp) ::   zperm                             ! permeability of sea ice 
    177171      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
    178172      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse 
     
    214208            !    a_ip_max = zfr_mlt * a_i 
    215209            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    216             v_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    217             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) ) 
     210            zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
     211            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    218212 
    219213            ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    220214            !    h_ip_max = 0.5 * h_i 
    221215            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    222             v_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
    223             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) ) 
     216            zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
     217            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    224218             
    225219            !--- Pond growing ---! 
     
    266260            ! v_ip     = h_ip * a_ip 
    267261            ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 
    268             a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) 
     262            a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    269263            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    270264            h_ip_1d(ji)      = zaspect * a_ip_frac_1d(ji) 
     
    283277                     &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
    284278                     &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 ! clem: error here the factor was 0.01878 instead of 0.0178 (cf Flocco 2010) 
    285                   zperm(jk) = MAX( 0._wp, 3.e-08_wp * (sz_i_1d(ji,jk) / zsbr)**3 ) 
     279                  ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
    286280               END DO 
     281               zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
    287282 
    288283               ! Do the drainage using Darcy's law 
    289                zdv_flush   = -MINVAL(zperm(:)) * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 
     284               zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 
    290285               zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
    291286               v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
    292287                
    293288               !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    294                a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) 
     289               a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    295290               a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
    296291               h_ip_1d(ji)      = zaspect * a_ip_frac_1d(ji) 
     
    300295            !--- Corrections and lid thickness ---! 
    301296            IF( ln_pnd_lids ) THEN 
     297               !--- retrieve lid thickness from volume ---! 
     298               IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
     299               ELSE                              ;   h_il_1d(ji) = 0._wp 
     300               ENDIF 
    302301               !--- remove ponds if lids are much larger than ponds ---! 
    303                IF ( v_il_1d(ji) > v_ip_1d(ji) * 10._wp ) THEN 
     302               IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
    304303                  a_ip_1d(ji)      = 0._wp 
    305304                  a_ip_frac_1d(ji) = 0._wp 
    306305                  h_ip_1d(ji)      = 0._wp 
    307                   v_il_1d(ji)      = 0._wp 
    308                ENDIF 
    309                !--- retrieve lid thickness from volume ---! 
    310                IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
    311                ELSE                              ;   h_il_1d(ji) = 0._wp ; 
     306                  h_il_1d(ji)      = 0._wp 
    312307               ENDIF 
    313308            ENDIF 
     
    316311          
    317312      END DO 
    318  
    319       !-------------------------------------------------!             
    320       ! How much melt pond is exposed to the atmosphere ! 
    321       !-------------------------------------------------!             
    322       ! Calculate the melt pond effective area (used for albedo) 
    323       WHERE    ( h_il_1d(1:npti) <= zhl_min )   ;   a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti)       ! lid is very thin.  Expose all the pond 
    324       ELSEWHERE( h_il_1d(1:npti) >= zhl_max )   ;   a_ip_eff_1d(1:npti) = 0._wp                      ! lid is very thick. Cover all the pond up with ice and snow 
    325       ELSEWHERE                                 ;   a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti) * &   ! lid is in between. Expose part of the pond 
    326          &                                                                ( h_il_1d(1:npti) - zhl_min ) / ( zhl_max - zhl_min ) 
    327       END WHERE 
    328313      ! 
    329314   END SUBROUTINE pnd_H12 
Note: See TracChangeset for help on using the changeset viewer.