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 13911 for NEMO – NEMO

Changeset 13911 for NEMO


Ignore:
Timestamp:
2020-11-30T10:49:09+01:00 (3 years ago)
Author:
vancop
Message:

Level-ice pond routine and two compilation bugfixes

Location:
NEMO/branches/2020/SI3_martin_ponds/src/ICE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icedyn_adv_umx.F90

    r13908 r13911  
    361361            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
    362362               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 
    363          ELSEIF( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 
     363         ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 
    364364            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
    365365               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 
  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_pnd.F90

    r13909 r13911  
    7373      diag_dvpn_lid_1d,  &  !! exchange with lid / refreezing   [m/day] 
    7474      diag_dvpn_rnf_1d      !! meltwater pond lost to runoff    [m/day] 
    75    ! 
    76  
     75 
     76   !! * Substitutions 
     77#  include "do_loop_substitute.h90" 
    7778   !!---------------------------------------------------------------------- 
    7879   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    213214      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect 
    214215      !! 
    215       !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 
     216      !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 
    216217      !!  
    217       !! ** Note       :   mostly stolen from CICE 
    218       !! 
    219       !! ** References :   Flocco and Feltham (JGR, 2007) 
    220       !!                   Flocco et al       (JGR, 2010) 
    221       !!                   Holland et al      (J. Clim, 2012) 
    222       !!------------------------------------------------------------------- 
     218      !! ** Note       :   Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds.  
     219      !! 
     220      !! ** References :   Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 
     221      !!                    
     222      !!------------------------------------------------------------------- 
     223       
    223224      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
    224225      !! 
     
    236237      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
    237238      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse 
    238       !! 
    239       INTEGER  ::   ji, jk                            ! loop indices 
    240       !!------------------------------------------------------------------- 
     239      REAL(wp) ::   zvold                             ! 
     240      !! 
     241      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     242       
     243      !!------------------------------------------------------------------- 
     244       
    241245      z1_rhow   = 1._wp / rhow  
    242246      z1_aspect = 1._wp / zaspect 
    243247      z1_Tp     = 1._wp / zTp  
    244  
    245       DO ji = 1, npti 
    246          !                                                            !----------------------------------------------------! 
    247          IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
    248             !                                                         !----------------------------------------------------! 
    249             !--- Remove ponds on thin ice or tiny ice fractions 
    250             a_ip_1d(ji)      = 0._wp 
    251             h_ip_1d(ji)      = 0._wp 
    252             h_il_1d(ji)      = 0._wp 
    253             !                                                         !--------------------------------! 
    254          ELSE                                                         ! Case ice thickness >= rn_himin ! 
    255             !                                                         !--------------------------------! 
    256             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
    257             v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
    258             ! 
    259             !------------------! 
    260             ! case ice melting ! 
    261             !------------------! 
    262             ! 
    263             !--- available meltwater for melt ponding ---! 
    264             zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    265             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 
    266             zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
    267             ! 
    268             !--- overflow ---! 
    269             ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
    270             !    a_ip_max = zfr_mlt * a_i 
    271             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    272             zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    273             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    274  
    275             ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    276             !    h_ip_max = 0.5 * h_i 
    277             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    278             zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
    279             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    280              
    281             !--- Pond growing ---! 
    282             v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
    283             ! 
    284             !--- Lid melting ---! 
    285             IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
    286             ! 
    287             !--- mass flux ---! 
    288             IF( zdv_mlt > 0._wp ) THEN 
    289                zfac = zdv_mlt * rhow * r1_Dt_ice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
    290                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    291                ! 
    292                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
    293                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    294                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    295             ENDIF 
    296  
    297             !-------------------! 
    298             ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
    299             !-------------------! 
    300             ! 
    301             zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
    302             ! 
    303             !--- Pond contraction (due to refreezing) ---! 
    304             IF( ln_pnd_lids ) THEN 
    305                ! 
    306                !--- Lid growing and subsequent pond shrinking ---!  
    307                zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
    308                   &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
    309                 
    310                ! Lid growing 
    311                v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
    312                 
    313                ! Pond shrinking 
    314                v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
    315  
    316             ELSE 
    317                ! Pond shrinking 
    318                v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
    319             ENDIF 
    320             ! 
    321             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    322             ! v_ip     = h_ip * a_ip 
    323             ! 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) 
    324             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 
    325             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    326  
    327             !---------------!             
    328             ! Pond flushing ! 
    329             !---------------! 
    330             ! height of top of the pond above sea-level 
    331             zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 
    332              
    333             ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
     248       
     249      !----------------------------------------------------------------------------------------------- 
     250      !  Identify grid cells with ice 
     251      !----------------------------------------------------------------------------------------------- 
     252      at_i(:,:) = SUM( a_i, dim=3 ) 
     253      ! 
     254      npti = 0   ;   nptidx(:) = 0 
     255      DO_2D( 1, 1, 1, 1 ) 
     256         IF ( at_i(ji,jj) > epsi10 ) THEN 
     257            npti = npti + 1 
     258            nptidx( npti ) = (jj - 1) * jpi + ji 
     259         ENDIF 
     260      END_2D 
     261       
     262      !----------------------------------------------------------------------------------------------- 
     263      ! Prepare 1D arrays 
     264      !----------------------------------------------------------------------------------------------- 
     265 
     266      IF( npti > 0 ) THEN 
     267       
     268         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum      ) 
     269         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti)   , wfx_sum          ) 
     270         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti)   , wfx_pnd          ) 
     271          
     272         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 
     273         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 
     274         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 
     275         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 
     276       
     277         DO jl = 1, jpl 
     278          
     279            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d      (1:npti), a_i      (:,:,jl) ) 
     280            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d      (1:npti), h_i      (:,:,jl) ) 
     281            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d     (1:npti), t_su     (:,:,jl) ) 
     282            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
     283            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
     284            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
     285 
     286            CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum    (1:npti), dh_i_sum_2d     (:,:,jl) ) 
     287            CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt    (1:npti), dh_s_mlt_2d     (:,:,jl) ) 
     288             
    334289            DO jk = 1, nlay_i 
    335                zsbr = - 1.2_wp                                  & 
    336                   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
    337                   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
    338                   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
    339                ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     290               CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  ) 
    340291            END DO 
    341             zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
    342              
    343             ! Do the drainage using Darcy's law 
    344             zdv_flush   = -zperm * rho0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 
    345             zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
    346             v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
    347              
    348             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    349             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 
    350             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    351  
    352             !--- Corrections and lid thickness ---! 
    353             IF( ln_pnd_lids ) THEN 
    354                !--- retrieve lid thickness from volume ---! 
    355                IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
    356                ELSE                              ;   h_il_1d(ji) = 0._wp 
    357                ENDIF 
    358                !--- remove ponds if lids are much larger than ponds ---! 
    359                IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     292             
     293            !----------------------------------------------------------------------------------------------- 
     294            ! Go for ponds 
     295            !----------------------------------------------------------------------------------------------- 
     296 
     297            DO ji = 1, npti 
     298               !                                                            !----------------------------------------------------! 
     299               IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
     300                  !                                                         !----------------------------------------------------! 
     301                  !--- Remove ponds on thin ice or tiny ice fractions 
    360302                  a_ip_1d(ji)      = 0._wp 
    361303                  h_ip_1d(ji)      = 0._wp 
    362304                  h_il_1d(ji)      = 0._wp 
     305                  !                                                         !--------------------------------! 
     306               ELSE                                                         ! Case ice thickness >= rn_himin ! 
     307                  !                                                         !--------------------------------! 
     308                  v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
     309                  v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     310                  ! 
     311                  !------------------! 
     312                  ! case ice melting ! 
     313                  !------------------! 
     314                  ! 
     315                  !--- available meltwater for melt ponding ---! 
     316                  zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
     317                  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 
     318                  zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
     319                   
     320                  diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_Dt_ice                      ! surface melt input diag 
     321                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_Dt_ice   ! runoff diag 
     322                  ! 
     323                  !--- overflow ---! 
     324                  ! 
     325                  ! 1) area driven overflow 
     326                  ! 
     327                  ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 
     328                  !    a_ip_max = zfr_mlt * a_i 
     329                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     330                  zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
     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                  ! 
     337                  ! If pond depth exceeds half the ice thickness then reduce the pond volume 
     338                  !    h_ip_max = 0.5 * h_i 
     339                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     340                  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                  zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     342                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_Dt_ice       ! runoff diag - overflow contribution 
     343             
     344                  !--- Pond growing ---! 
     345                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     346                  ! 
     347                  !--- Lid melting ---! 
     348                  IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
     349                  ! 
     350                  !--- mass flux ---! 
     351                  ! MV I would recommend to remove that 
     352                  ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 
     353                   
     354                  IF( zdv_mlt > 0._wp ) THEN 
     355                     zfac = zdv_mlt * rhow * r1_Dt_ice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
     356                     wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
     357                     ! 
     358                     zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
     359                     wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
     360                     wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
     361                  ENDIF 
     362 
     363                  !-------------------! 
     364                  ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
     365                  !-------------------! 
     366                  ! 
     367                  zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
     368                  ! 
     369                  !--- Pond contraction (due to refreezing) ---! 
     370                  zvold       = v_ip_1d(ji) ! for diag 
     371 
     372                  IF( ln_pnd_lids ) THEN 
     373                     ! 
     374                     !--- Lid growing and subsequent pond shrinking ---!  
     375                     zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
     376                        &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
     377                
     378                     ! Lid growing 
     379                     v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
     380                
     381                     ! Pond shrinking 
     382                     v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
     383                      
     384                  ELSE 
     385                   
     386                     ! Pond shrinking 
     387                     v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
     388                      
     389                  ENDIF 
     390 
     391                  diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_Dt_ice   ! shrinking counted as lid diagnostic 
     392 
     393                  ! 
     394                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     395                  ! v_ip     = h_ip * a_ip 
     396                  ! 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) 
     397                  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 
     398                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     399 
     400                  !------------------------------------------------!             
     401                  ! Pond drainage through brine network (flushing) ! 
     402                  !------------------------------------------------! 
     403                  ! height of top of the pond above sea-level 
     404                  zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 
     405             
     406                  ! Calculate the permeability of the ice  
     407                  ! mv- note here that a linear expression for permeability is fine,  
     408                  ! simpler and more consistent with the rest of SI3 code 
     409                  DO jk = 1, nlay_i 
     410                     ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     411                  END DO 
     412                  zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
     413             
     414                  ! Do the drainage using Darcy's law 
     415                  zdv_flush   = - zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 
     416                  zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
     417                  ! zdv_flush   = 0._wp ! MV remove pond drainage for now 
     418                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
     419                   
     420                  diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_Dt_ice   ! shrinking counted as lid diagnostic 
     421             
     422                  ! MV --- why pond drainage does not give back water into freshwater flux ?             
     423                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     424                  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 
     425                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     426 
     427                  !--- Corrections and lid thickness ---! 
     428                  IF( ln_pnd_lids ) THEN 
     429                     !--- retrieve lid thickness from volume ---! 
     430                     IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
     431                     ELSE                              ;   h_il_1d(ji) = 0._wp 
     432                     ENDIF 
     433                     !--- remove ponds if lids are much larger than ponds ---! 
     434                     IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     435                        a_ip_1d(ji)      = 0._wp 
     436                        h_ip_1d(ji)      = 0._wp 
     437                        h_il_1d(ji)      = 0._wp 
     438                     ENDIF 
     439                  ENDIF 
     440                  ! 
    363441               ENDIF 
    364             ENDIF 
    365             ! 
    366          ENDIF 
     442                
     443            END DO ! ji 
     444 
     445            !----------------------------------------------------------------------------------------------- 
     446            ! Retrieve 2D arrays 
     447            !----------------------------------------------------------------------------------------------- 
     448             
     449            v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
     450            v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 
     451             
     452            CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
     453            CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
     454            CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
     455            CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d     (1:npti), v_ip     (:,:,jl) ) 
     456            CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d     (1:npti), v_il     (:,:,jl) ) 
     457            DO jk = 1, nlay_i 
     458               CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  ) 
     459            END DO 
     460    
     461         END DO ! ji 
     462                   
     463         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 
     464         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        ) 
     465         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    367466          
    368       END DO 
    369       ! 
     467         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        ) 
     468         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        ) 
     469         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        ) 
     470         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        )                   
     471                            
     472      ! 
     473      ENDIF 
     474       
    370475   END SUBROUTINE pnd_LEV 
    371  
    372476 
    373477   SUBROUTINE ice_thd_pnd_init  
  • NEMO/branches/2020/SI3_martin_ponds/src/ICE/icevar.F90

    r13908 r13911  
    690690      WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    691691      WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    692       IF( ln_pnd_LEV .OR ln_pnd_TOPO ) THEN 
     692      IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    693693         WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    694694         WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 
Note: See TracChangeset for help on using the changeset viewer.