Changeset 12197


Ignore:
Timestamp:
2019-12-12T00:46:12+01:00 (8 months ago)
Author:
clem
Message:

correct ticket #2349

Location:
NEMO/trunk/src/ICE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_adv.F90

    r11536 r12197  
    8888      CASE( np_advPRA )                ! PRATHER scheme        ! 
    8989         !                             !-----------------------! 
    90          CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, & 
     90         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, h_i, h_s, h_ip, & 
    9191            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
    9292      END SELECT 
  • NEMO/trunk/src/ICE/icedyn_adv_pra.F90

    r11812 r12197  
    5454CONTAINS 
    5555 
    56    SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  & 
     56   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    5757      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    5858      !!---------------------------------------------------------------------- 
     
    7070      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity 
    7171      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity 
     72      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_i       ! ice thickness 
     73      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_s       ! snw thickness 
     74      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   ph_ip      ! ice pond thickness 
    7275      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area 
    7376      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume 
     
    8790      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
    8891      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
     92      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max 
    8993      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
    9094      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
     
    9599      ! 
    96100      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
     101      ! 
     102      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     103      DO jl = 1, jpl 
     104         DO jj = 2, jpjm1 
     105            DO ji = fs_2, fs_jpim1 
     106               zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
     107                  &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
     108                  &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
     109                  &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
     110               zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
     111                  &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
     112                  &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
     113                  &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
     114               zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
     115                  &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
     116                  &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
     117                  &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     118            END DO 
     119         END DO 
     120      END DO 
     121      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 
    97122      ! 
    98123      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     
    239264         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    240265         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     266         ! 
     267         ! --- Make sure ice thickness is not too big --- ! 
     268         !     (because ice thickness can be too large where ice concentration is very small) 
     269         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
    241270         ! 
    242271         ! --- Ensure snow load is not too big --- ! 
     
    588617      ! 
    589618   END SUBROUTINE adv_y 
     619 
     620 
     621   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     622      !!------------------------------------------------------------------- 
     623      !!                  ***  ROUTINE Hbig  *** 
     624      !! 
     625      !! ** Purpose : Thickness correction in case advection scheme creates 
     626      !!              abnormally tick ice or snow 
     627      !! 
     628      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     629      !!                 (before advection) and reduce it by adapting ice concentration 
     630      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     631      !!                 (before advection) and reduce it by sending the excess in the ocean 
     632      !! 
     633      !! ** input   : Max thickness of the surrounding 9-points 
     634      !!------------------------------------------------------------------- 
     635      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
     636      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
     637      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     638      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     639      ! 
     640      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
     641      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     642      !!------------------------------------------------------------------- 
     643      ! 
     644      z1_dt = 1._wp / pdt 
     645      ! 
     646      DO jl = 1, jpl 
     647 
     648         DO jj = 1, jpj 
     649            DO ji = 1, jpi 
     650               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     651                  ! 
     652                  !                               ! -- check h_ip -- ! 
     653                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     654                  IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     655                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     656                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     657                        pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     658                     ENDIF 
     659                  ENDIF 
     660                  ! 
     661                  !                               ! -- check h_i -- ! 
     662                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     663                  zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     664                  IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     665                     pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     666                  ENDIF 
     667                  ! 
     668                  !                               ! -- check h_s -- ! 
     669                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     670                  zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     671                  IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     672                     zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
     673                     ! 
     674                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 
     675                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     676                     ! 
     677                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     678                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     679                  ENDIF            
     680                  !                   
     681               ENDIF 
     682            END DO 
     683         END DO 
     684      END DO  
     685      ! 
     686   END SUBROUTINE Hbig 
    590687 
    591688 
  • NEMO/trunk/src/ICE/icedyn_adv_umx.F90

    r11627 r12197  
    352352         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    353353         ! 
    354          ! Make sure ice thickness is not too big 
    355          !    (because ice thickness can be too large where ice concentration is very small) 
    356          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
    357  
     354         ! --- Make sure ice thickness is not too big --- ! 
     355         !     (because ice thickness can be too large where ice concentration is very small) 
     356         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     357         ! 
     358         ! --- Ensure snow load is not too big --- ! 
     359         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
     360         ! 
    358361      END DO 
    359362      ! 
     
    15141517 
    15151518 
    1516    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     1519   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
    15171520      !!------------------------------------------------------------------- 
    15181521      !!                  ***  ROUTINE Hbig  *** 
     
    15251528      !!              2- check whether snow thickness is larger than the surrounding 9-points 
    15261529      !!                 (before advection) and reduce it by sending the excess in the ocean 
    1527       !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
    1528       !!                 and reduce it by sending the excess in the ocean 
    1529       !!              4- correct pond concentration to avoid a_ip > a_i 
    15301530      !! 
    15311531      !! ** input   : Max thickness of the surrounding 9-points 
     
    15331533      REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    15341534      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    1535       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 
     1535      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
    15361536      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    1537       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
    1538       ! 
    1539       INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
    1540       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zvs_excess, zfra 
    1541       REAL(wp), DIMENSION(jpi,jpj) ::   zswitch 
     1537      ! 
     1538      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
     1539      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
    15421540      !!------------------------------------------------------------------- 
    15431541      ! 
     
    15781576                     pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
    15791577                  ENDIF            
     1578                  !                   
     1579               ENDIF 
     1580            END DO 
     1581         END DO 
     1582      END DO  
     1583      ! 
     1584   END SUBROUTINE Hbig 
     1585 
     1586 
     1587   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s ) 
     1588      !!------------------------------------------------------------------- 
     1589      !!                  ***  ROUTINE Hsnow  *** 
     1590      !! 
     1591      !! ** Purpose : 1- Check snow load after advection 
     1592      !!              2- Correct pond concentration to avoid a_ip > a_i 
     1593      !! 
     1594      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface 
     1595      !!              then put the snow excess in the ocean 
     1596      !! 
     1597      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards 
     1598      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially 
     1599      !!              make the snow very thick (if concentration decreases drastically) 
     1600      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather 
     1601      !!------------------------------------------------------------------- 
     1602      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step 
     1603      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip 
     1604      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
     1605      ! 
     1606      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     1607      REAL(wp) ::   z1_dt, zvs_excess, zfra 
     1608      !!------------------------------------------------------------------- 
     1609      ! 
     1610      z1_dt = 1._wp / pdt 
     1611      ! 
     1612      ! -- check snow load -- ! 
     1613      DO jl = 1, jpl 
     1614         DO jj = 1, jpj 
     1615            DO ji = 1, jpi 
     1616               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
    15801617                  ! 
    1581                   !                               ! -- check snow load -- ! 
    1582                   ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
    1583                   !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
    1584                   !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
    15851618                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
    1586                   IF( zvs_excess > 0._wp ) THEN 
     1619                  ! 
     1620                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface 
     1621                     ! put snow excess in the ocean 
    15871622                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
    15881623                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
    15891624                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
    1590                      ! 
     1625                     ! correct snow volume and heat content 
    15911626                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
    15921627                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
    15931628                  ENDIF 
    1594                    
     1629                  ! 
    15951630               ENDIF 
    15961631            END DO 
    15971632         END DO 
    1598       END DO  
    1599       !                                           !-- correct pond concentration to avoid a_ip > a_i 
     1633      END DO 
     1634      ! 
     1635      !-- correct pond concentration to avoid a_ip > a_i -- ! 
    16001636      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:) 
    16011637      ! 
    1602       ! 
    1603    END SUBROUTINE Hbig 
    1604     
     1638   END SUBROUTINE Hsnow 
     1639 
     1640 
    16051641#else 
    16061642   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.