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 13891 – NEMO

Changeset 13891


Ignore:
Timestamp:
2020-11-26T17:02:26+01:00 (4 years ago)
Author:
vancop
Message:

fix a few bugs, add pond volume budget diagnostics

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90

    r13879 r13891  
    4545      zTd     = 0.15_wp       , & ! temperature difference for freeze-up (C) 
    4646      zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
     47   !  
     48   ! Pond volume per area budget diags 
     49   !   
     50   ! The idea of diags is the volume of ponds per grid cell area is 
     51   ! 
     52   ! dV/dt = mlt + drn + lid + rnf 
     53   ! mlt   = input from surface melting 
     54   ! drn   = drainage through brine network 
     55   ! lid   = lid growth & melt 
     56   ! rnf   = runoff (water directly removed out of surface melting + overflow) 
     57   ! 
     58   ! In topo mode, the pond water lost because it is in the snow is not included in the budget 
     59   ! 
     60   ! In level mode 
     61 
     62   REAL(wp), DIMENSION(jpi,jpj) ::   ! pond volume budget diagnostics 
     63      diag_dvpn_mlt,  &     !! meltwater pond volume input      [m/day] 
     64      diag_dvpn_drn,  &     !! pond volume lost by drainage     [m/day] 
     65      diag_dvpn_lid,  &     !! pond volume lost by refreezing   [m/day] 
     66      diag_dvpn_rnf         !! meltwater pond lost to runoff    [m/day] 
     67       
     68   REAL(wp), DIMENSION(jpij) ::   ! pond volume budget diagnostics (1d) 
     69      diag_dvpn_mlt_1d,  &  !! meltwater pond volume input      [m/day] 
     70      diag_dvpn_drn_1d,  &  !! pond volume lost by drainage     [m/day] 
     71      diag_dvpn_lid_1d,  &  !! pond volume lost by refreezing   [m/day] 
     72      diag_dvpn_rnf_1d      !! meltwater pond lost to runoff    [m/day] 
     73 
    4774 
    4875   !! * Substitutions 
     
    76103      !!------------------------------------------------------------------- 
    77104      ! 
     105      diag_dvpn_mlt(:,:) = 0._wp      ;   diag_dvpn_drn(:,:)    = 0._wp 
     106      diag_dvpn_lid(:,:) = 0._wp      ;   diag_dvpn_rnf(:,:)    = 0._wp 
     107      diag_dvpn_mlt_1d(:,:) = 0._wp   ;   diag_dvpn_drn_1d(:,:) = 0._wp 
     108      diag_dvpn_lid_1d(:,:) = 0._wp   ;   diag_dvpn_rnf_1d(:,:) = 0._wp 
     109 
    78110      SELECT CASE ( nice_pnd ) 
    79111      ! 
     
    82114      CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==! 
    83115         ! 
    84       CASE (np_pndTOPO)  ;   CALL pnd_TOPO                  !&         !==  Topographic melt ponds  ==! 
     116      CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==! 
    85117         ! 
    86118      END SELECT 
    87119      ! 
     120 
     121      IF( iom_use('dvpn_mlt'  ) )   CALL iom_put( 'dvpn_mlt', dvpn_mlt * 100._wp ) ! input from melting 
     122      IF( iom_use('dvpn_lid'  ) )   CALL iom_put( 'dvpn_lid', dvpn_lid * 100._wp ) ! exchanges with lid 
     123      IF( iom_use('dvpn_drn'  ) )   CALL iom_put( 'dvpn_drn', dvpn_drn * 100._wp ) ! vertical drainage 
     124      IF( iom_use('dvpn_rnf'  ) )   CALL iom_put( 'dvpn_rnf', dvpn_rnf * 100._wp ) ! runoff + overflow 
     125       
    88126   END SUBROUTINE ice_thd_pnd  
     127 
    89128 
    90129 
     
    191230      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
    192231      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse 
     232      REAL(wp) ::   z1_rdtice                         ! inverse time step 
     233      REAL(wp) ::   zvold                             ! 
    193234      !! 
    194235      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     
    199240      z1_aspect = 1._wp / zaspect 
    200241      z1_Tp     = 1._wp / zTp  
     242      z1_rdtice = 1._wp / rdt_ice 
    201243       
    202244      !----------------------------------------------------------------------------------------------- 
     
    221263      IF( npti > 0 ) THEN 
    222264       
    223          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum   ) 
    224          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum          ) 
    225          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          ) 
     265         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum      ) 
     266         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti)   , wfx_sum          ) 
     267         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti)   , wfx_pnd          ) 
     268          
     269         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 
     270         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 
     271         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 
     272         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 
    226273       
    227274         DO jl = 1, jpl 
     
    268315                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed 
    269316                  zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
     317                   
     318                  diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * z1_rdtice                      ! surface melt input diag 
     319                   
     320                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * z1_rdtice   ! runoff diag 
     321                   
    270322                  ! 
    271323                  !--- overflow ---! 
    272                   ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
     324                  ! 
     325                  ! 1) area driven overflow 
     326                  ! 
     327                  ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 
    273328                  !    a_ip_max = zfr_mlt * a_i 
    274329                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    275330                  zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    276                   zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    277  
     331                  zvold     = zdv_mlt 
     332                  zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )                      
     333                   
     334                  !  
     335                  ! 2) depth driven overflow 
     336                  ! 
    278337                  ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    279338                  !    h_ip_max = 0.5 * h_i 
    280339                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    281340                  zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear 
     341                                    
    282342                  zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     343                   
     344                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * z1_rdtice       ! runoff diag - overflow contribution 
    283345             
    284346                  !--- Pond growing ---! 
     
    308370                  ! 
    309371                  !--- Pond contraction (due to refreezing) ---! 
     372                  zvold       = v_ip_1d(ji) ! for diag 
     373 
    310374                  IF( ln_pnd_lids ) THEN 
    311375                     ! 
     
    319383                     ! Pond shrinking 
    320384                     v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
    321  
     385                      
    322386                  ELSE 
     387                   
    323388                     ! Pond shrinking 
    324389                     v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
     390                      
    325391                  ENDIF 
     392 
     393                  diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * z1_rdtice   ! shrinking counted as lid diagnostic 
     394 
    326395                  ! 
    327396                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     
    331400                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    332401 
    333                   !---------------!             
    334                   ! Pond flushing ! 
    335                   !---------------! 
     402                  !------------------------------------------------!             
     403                  ! Pond drainage through brine network (flushing) ! 
     404                  !------------------------------------------------! 
    336405                  ! height of top of the pond above sea-level 
    337406                  zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 
     
    352421                  zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 
    353422                  zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
    354                   zdv_flush   = 0._wp ! MV remove pond drainage for now 
     423                  ! zdv_flush   = 0._wp ! MV remove pond drainage for now 
    355424                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
     425                   
     426                  diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * z1_rdtice   ! shrinking counted as lid diagnostic 
    356427             
    357                   ! MV --- why pond drainage does not give back water into freshwater flux ? 
    358              
     428                  ! MV --- why pond drainage does not give back water into freshwater flux ?             
    359429                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    360430                  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 
     
    385455            v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
    386456            v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 
     457             
    387458            CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
    388459            CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
     
    399470         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        ) 
    400471         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    401                    
     472          
     473         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        ) 
     474         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        ) 
     475         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        ) 
     476         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        )                   
     477                            
    402478      ! 
    403479      ENDIF 
     
    438514      ! local variables 
    439515      REAL(wp) :: & 
    440          zdHui,   &   ! change in thickness of ice lid (m) 
    441          zomega,  &   ! conduction 
    442          zdTice,  &   ! temperature difference across ice lid (C) 
    443          zdvice,  &   ! change in ice volume (m) 
    444          zTavg,   &   ! mean surface temperature across categories (C) 
    445          zfsurf,  &   ! net heat flux, excluding conduction and transmitted radiation (W/m2) 
    446          zTp,     &   ! pond freezing temperature (C) 
    447          zrhoi_L, &   ! volumetric latent heat of sea ice (J/m^3) 
    448          zfr_mlt, &   ! fraction and volume of available meltwater retained for melt ponding 
    449          z1_rhow, &   ! inverse water density 
    450          zpond  , &   ! dummy variable 
    451          zdum         ! dummy variable 
    452  
    453       REAL(wp), DIMENSION(jpi,jpj) ::   zvolp, & !! total melt pond water available before redistribution and drainage 
    454                                         zvolp_res 
    455  
     516         zdHui,   &      ! change in thickness of ice lid (m) 
     517         zomega,  &      ! conduction 
     518         zdTice,  &      ! temperature difference across ice lid (C) 
     519         zdvice,  &      ! change in ice volume (m) 
     520         zTavg,   &      ! mean surface temperature across categories (C) 
     521         zfsurf,  &      ! net heat flux, excluding conduction and transmitted radiation (W/m2) 
     522         zTp,     &      ! pond freezing temperature (C) 
     523         zrhoi_L, &      ! volumetric latent heat of sea ice (J/m^3) 
     524         zfr_mlt, &      ! fraction and volume of available meltwater retained for melt ponding 
     525         z1_rhow, &      ! inverse water density 
     526         z1_rdtice, &    ! inverse time step 
     527         zv_pnd  , &     ! volume of meltwater contributing to ponds 
     528         zv_mlt          ! total amount of meltwater produced 
     529 
     530      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp, &     !! total melt pond water available before redistribution and drainage 
     531                                        zvolp_res    !! remaining melt pond water available after drainage 
     532                                         
    456533      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i 
    457534 
     
    473550      zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline) 
    474551      z1_rhow   = 1._wp / rhow  
     552      z1_rdtice = 1._wp / rdt_ice 
    475553 
    476554      ! Set required ice variables (hard-coded here for now) 
     
    491569      ! Change 2D to 1D 
    492570      !--------------------------------------------------------------- 
     571       
     572      ! use what we have in iceitd.F90 (incremental remapping) 
    493573 
    494574      !-------------------------------------------------------------- 
    495       ! Collect total available pond water 
     575      ! Collect total available pond water volume 
    496576      !-------------------------------------------------------------- 
    497  
     577      ! Assuming that meltwater (+rain in principle) runsoff the surface 
     578      ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 
     579      ! I cite her words, they are very talkative 
     580      ! "grid cells with very little ice cover (and hence more open water area)  
     581      ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 
     582      ! "This results in the same runoff fraction r for each ice category within a grid cell" 
     583       
    498584      zvolp(:,:) = 0._wp 
    499585 
     
    505591             
    506592                  !--- Available meltwater for melt ponding ---! 
    507                   zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * a_i(ji,jj,jl) !  = ( 1 - r ) = fraction of melt water that is not flushed 
    508                   zdum    = -( dh_i_sum_2d(ji,jj,jl)*rhoi + dh_s_mlt_2d(ji,jj,jl)*rhos ) * z1_rhow * a_i(ji,jj,jl) 
    509                   zpond   = zfr_mlt * zdum ! check 
     593                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds 
     594 
     595                  zv_mlt  = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) &        ! total volume of surface melt water 
     596                     &    * z1_rhow * a_i(ji,jj,jl) 
     597                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! volume of meltwater contributing to ponds 
     598                   
     599                  diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * z1_rdtice                     ! diagnostics 
     600                   
     601                  diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * z1_rdtice    
    510602 
    511603                  !--- Create possible new ponds 
     
    513605                  IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 
    514606                     h_ip(ji,jj,jl)       = 0._wp 
    515                      a_ip_frac(ji,jj,jl) = 1.0_wp    ! pond fraction of sea ice (apnd for CICE) 
     607                     a_ip_frac(ji,jj,jl)  = 1.0_wp    ! pond fraction of sea ice (apnd for CICE) 
     608                     a_ip(ji,jj,jl)       = a_i(ji,jj,jl) 
    516609                  ENDIF 
    517610                   
    518                   !--- Deepen existing ponds before redistribution and drainage(later on) 
    519                   h_ip(ji,jj,jl) = ( zpond + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) ) / a_ip_frac(ji,jj,jl) 
    520                   zvolp(ji,jj)   = zvolp(ji,jj) + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)  
    521 !                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) !useless 
     611                  !--- Deepen existing ponds before redistribution and drainage 
     612                  ! without pond fraction 
     613                  v_ip(ji,jj,jl) = v_i_p(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness 
     614                      
     615                  h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     616                   
     617                  zvolp(ji,jj)   = zvolp(ji,jj) + v_ip(ji,jj,jl) 
     618                   
     619!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 
    522620                    
    523621               ENDIF ! a_i 
     
    526624         END DO ! jj 
    527625      END DO ! ji 
    528           
    529       h_ip(:,:,:)  = 0._wp ! pond thickness 
    530       a_ip(:,:,:)  = 0._wp ! pond fraction per grid area   
    531           
     626                   
    532627      !-------------------------------------------------------------- 
    533628      ! Redistribute and drain water from ponds 
     
    564659 
    565660                        zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 
    566                         IF ( zdvice > epsi10 ) then 
     661                         
     662                        IF ( zdvice > epsi10 ) THEN 
     663                         
    567664                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice 
    568                            v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice 
     665                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice -  
     666                                                                        ! as it is already counted in surface melt 
    569667!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar) 
    570668!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw) 
     
    577675                           ENDIF 
    578676                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 
     677                            
     678                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag 
     679                            
    579680                        ENDIF 
    580                    
     681                         
    581682                     !---------------------------------------------------------------- 
    582683                     ! Freeze pre-existing lid  
     
    598699!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice 
    599700                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     701                            
     702                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag 
     703                            
    600704                        ENDIF 
    601705                   
     
    624728                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
    625729                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     730                         
     731                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag 
    626732!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice 
    627733!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice 
     
    849955             DO ns = 1, jl 
    850956                cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 
    851                    - rhos/rhow * &   ! free air fraction that can be filled by water 
     957                   - rhos/rhow * &     ! free air fraction that can be filled by water 
    852958                     asnon(ns)  * &    ! effective areal fraction of snow in that category 
    853959                     max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 
     
    870976          drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 
    871977          zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 
     978 
     979          diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice) 
     980           
    872981          zdvolp(ji,jj) = drain         ! this is the drained water 
    873982          IF (zvolp(ji,jj) < epsi10) THEN 
     
    878987 
    879988       ! height and area corresponding to the remaining volume 
     989       ! routine leaves zvolp unchanged 
    880990       CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
    881991 
     
    9171027 
    9181028                perm = 0._wp ! MV ugly dummy patch 
    919                 ! CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof  
     1029                CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof  
    9201030                IF (perm > 0._wp) permflag = 1 
    9211031 
     
    9241034                zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 
    9251035                zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 
     1036 
     1037                diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 
     1038                 
    9261039                IF (zvolp(ji,jj) < epsi10) THEN 
    9271040                   zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     
    9501063       ! snow in melt ponds is not melted 
    9511064       !------------------------------------------------------------------------ 
     1065        
     1066       ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 
     1067       ! how much, so I did not diagnose it 
     1068       ! so if there is a problem here, nobody is going to see it... 
     1069        
    9521070 
    9531071       ! Calculate pond volume for lower categories 
    9541072       DO jl = 1,m_index-1 
    9551073          v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 
    956                      - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 
     1074                         - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 
    9571075       END DO 
    9581076 
Note: See TracChangeset for help on using the changeset viewer.