Changeset 12197
- Timestamp:
- 2019-12-12T00:46:12+01:00 (5 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv.F90
r11536 r12197 88 88 CASE( np_advPRA ) ! PRATHER scheme ! 89 89 ! !-----------------------! 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, & 91 91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 92 92 END SELECT -
NEMO/trunk/src/ICE/icedyn_adv_pra.F90
r11812 r12197 54 54 CONTAINS 55 55 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, & 57 57 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 58 58 !!---------------------------------------------------------------------- … … 70 70 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity 71 71 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 72 75 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 73 76 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 87 90 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 88 91 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx 92 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 89 93 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea 90 94 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi … … 95 99 ! 96 100 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. ) 97 122 ! 98 123 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! … … 239 264 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 240 265 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 ) 241 270 ! 242 271 ! --- Ensure snow load is not too big --- ! … … 588 617 ! 589 618 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 590 687 591 688 -
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r11627 r12197 352 352 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 ) 353 353 ! 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 ! 358 361 END DO 359 362 ! … … 1514 1517 1515 1518 1516 SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, p sv_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 ) 1517 1520 !!------------------------------------------------------------------- 1518 1521 !! *** ROUTINE Hbig *** … … 1525 1528 !! 2- check whether snow thickness is larger than the surrounding 9-points 1526 1529 !! (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 ocean1529 !! 4- correct pond concentration to avoid a_ip > a_i1530 1530 !! 1531 1531 !! ** input : Max thickness of the surrounding 9-points … … 1533 1533 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1534 1534 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, p sv_i, poa_i, pa_i, pa_ip, pv_ip1535 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip 1536 1536 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 1542 1540 !!------------------------------------------------------------------- 1543 1541 ! … … 1578 1576 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1579 1577 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 1580 1617 ! 1581 ! ! -- check snow load -- !1582 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean1583 ! 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)1585 1618 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 1587 1622 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1588 1623 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 1589 1624 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 1591 1626 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1592 1627 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1593 1628 ENDIF 1594 1629 ! 1595 1630 ENDIF 1596 1631 END DO 1597 1632 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 -- ! 1600 1636 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1601 1637 ! 1602 !1603 END SUBROUTINE Hbig 1604 1638 END SUBROUTINE Hsnow 1639 1640 1605 1641 #else 1606 1642 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.