- Timestamp:
- 2019-09-18T16:11:52+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/ICE/icevar.F90
r10945 r11564 32 32 !! - vt_s(jpi,jpj) 33 33 !! - at_i(jpi,jpj) 34 !! - st_i(jpi,jpj) 34 35 !! - et_s(jpi,jpj) total snow heat content 35 36 !! - et_i(jpi,jpj) total ice thermal content … … 44 45 !! ice_var_salprof1d : salinity profile in the ice 1D 45 46 !! ice_var_zapsmall : remove very small area and volume 46 !! ice_var_zapneg : remove negative ice fields (to debug the advection scheme UM3-5) 47 !! ice_var_itd : convert 1-cat to jpl-cat 48 !! ice_var_itd2 : convert N-cat to jpl-cat 47 !! ice_var_zapneg : remove negative ice fields 48 !! ice_var_roundoff : remove negative values arising from roundoff erros 49 49 !! ice_var_bv : brine volume 50 50 !! ice_var_enthalpy : compute ice and snow enthalpies from temperature 51 51 !! ice_var_sshdyn : compute equivalent ssh in lead 52 !! ice_var_itd : convert N-cat to M-cat 52 53 !!---------------------------------------------------------------------- 53 54 USE dom_oce ! ocean space and time domain … … 71 72 PUBLIC ice_var_zapsmall 72 73 PUBLIC ice_var_zapneg 73 PUBLIC ice_var_itd 74 PUBLIC ice_var_itd2 74 PUBLIC ice_var_roundoff 75 75 PUBLIC ice_var_bv 76 76 PUBLIC ice_var_enthalpy 77 77 PUBLIC ice_var_sshdyn 78 PUBLIC ice_var_itd 79 80 INTERFACE ice_var_itd 81 MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc 82 END INTERFACE 78 83 79 84 !!---------------------------------------------------------------------- … … 99 104 ! 100 105 ! ! integrated values 101 vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 ) 102 vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 ) 103 at_i(:,:) = SUM( a_i(:,:,:) , dim=3 ) 104 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 105 et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 106 vt_i(:,:) = SUM( v_i (:,:,:) , dim=3 ) 107 vt_s(:,:) = SUM( v_s (:,:,:) , dim=3 ) 108 st_i(:,:) = SUM( sv_i(:,:,:) , dim=3 ) 109 at_i(:,:) = SUM( a_i (:,:,:) , dim=3 ) 110 et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 111 et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 106 112 ! 107 113 at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds … … 133 139 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 134 140 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 135 sm_i (:,:) = SUM( sv_i(:,:,:) , dim=3 )* z1_vt_i(:,:)141 sm_i (:,:) = st_i(:,:) * z1_vt_i(:,:) 136 142 ! 137 143 tm_i(:,:) = 0._wp … … 153 159 tm_s (:,:) = rt0 154 160 END WHERE 155 161 ! 162 ! ! mean melt pond depth 163 WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 164 ELSEWHERE ; hm_ip(:,:) = 0._wp 165 END WHERE 166 ! 156 167 DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 168 ! 157 169 ENDIF 158 170 ! … … 258 270 ! 259 271 ! integrated values 260 vt_i (:,:) = SUM( v_i , dim=3 )261 vt_s (:,:) = SUM( v_s , dim=3 )262 at_i (:,:) = SUM( a_i , dim=3 )272 vt_i (:,:) = SUM( v_i , dim=3 ) 273 vt_s (:,:) = SUM( v_s , dim=3 ) 274 at_i (:,:) = SUM( a_i , dim=3 ) 263 275 ! 264 276 END SUBROUTINE ice_var_glo2eqv … … 528 540 529 541 ! to be sure that at_i is the sum of a_i(jl) 530 at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 531 vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) 542 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) 543 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) 544 !!clem add? 545 ! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 ) 546 ! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 ) 547 ! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 548 ! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 549 !!clem 532 550 533 551 ! open water = 1 if at_i=0 … … 622 640 END SUBROUTINE ice_var_zapneg 623 641 642 643 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 644 !!------------------------------------------------------------------- 645 !! *** ROUTINE ice_var_roundoff *** 646 !! 647 !! ** Purpose : Remove negative sea ice values arising from roundoff errors 648 !!------------------------------------------------------------------- 649 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration 650 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume 651 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume 652 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content 653 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content 654 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 655 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume 656 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content 657 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content 658 !!------------------------------------------------------------------- 659 ! 660 WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0 661 WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0 662 WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0 663 WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0 664 WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0 665 WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0 666 WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0 667 IF( ln_pnd_H12 ) THEN 668 WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 669 WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 670 ENDIF 671 ! 672 END SUBROUTINE ice_var_roundoff 624 673 625 SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )626 !!-------------------------------------------------------------------627 !! *** ROUTINE ice_var_itd ***628 !!629 !! ** Purpose : converting 1-cat ice to multiple ice categories630 !!631 !! ice thickness distribution follows a gaussian law632 !! around the concentration of the most likely ice thickness633 !! (similar as iceistate.F90)634 !!635 !! ** Method: Iterative procedure636 !!637 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian638 !!639 !! 2) Check whether the distribution conserves area and volume, positivity and640 !! category boundaries641 !!642 !! 3) If not (input ice is too thin), the last category is empty and643 !! the number of categories is reduced (jpl-1)644 !!645 !! 4) Iterate until ok (SUM(itest(:) = 4)646 !!647 !! ** Arguments : zhti: 1-cat ice thickness648 !! zhts: 1-cat snow depth649 !! zati: 1-cat ice concentration650 !!651 !! ** Output : jpl-cat652 !!653 !! (Example of application: BDY forcings when input are cell averaged)654 !!-------------------------------------------------------------------655 INTEGER :: ji, jk, jl ! dummy loop indices656 INTEGER :: idim, i_fill, jl0657 REAL(wp) :: zarg, zV, zconv, zdh, zdv658 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables659 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables660 INTEGER , DIMENSION(4) :: itest661 !!-------------------------------------------------------------------662 !663 ! ----------------------------------------664 ! distribution over the jpl ice categories665 ! ----------------------------------------666 ! a gaussian distribution for ice concentration is used667 ! then we check whether the distribution fullfills668 ! volume and area conservation, positivity and ice categories bounds669 idim = SIZE( zhti , 1 )670 zh_i(1:idim,1:jpl) = 0._wp671 zh_s(1:idim,1:jpl) = 0._wp672 za_i(1:idim,1:jpl) = 0._wp673 !674 DO ji = 1, idim675 !676 IF( zhti(ji) > 0._wp ) THEN677 !678 ! find which category (jl0) the input ice thickness falls into679 jl0 = jpl680 DO jl = 1, jpl681 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN682 jl0 = jl683 CYCLE684 ENDIF685 END DO686 !687 itest(:) = 0688 i_fill = jpl + 1 !------------------------------------689 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories690 ! !------------------------------------691 i_fill = i_fill - 1692 !693 zh_i(ji,1:jpl) = 0._wp694 za_i(ji,1:jpl) = 0._wp695 itest(:) = 0696 !697 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1698 zh_i(ji,1) = zhti(ji)699 za_i (ji,1) = zati (ji)700 ELSE !-- case ice is thicker: fill categories >1701 ! thickness702 DO jl = 1, i_fill - 1703 zh_i(ji,jl) = hi_mean(jl)704 END DO705 !706 ! concentration707 za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))708 DO jl = 1, i_fill - 1709 IF ( jl /= jl0 ) THEN710 zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )711 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)712 ENDIF713 END DO714 !715 ! last category716 za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )717 zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )718 zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )719 !720 ! correction if concentration of upper cat is greater than lower cat721 ! (it should be a gaussian around jl0 but sometimes it is not)722 IF ( jl0 /= jpl ) THEN723 DO jl = jpl, jl0+1, -1724 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN725 zdv = zh_i(ji,jl) * za_i(ji,jl)726 zh_i(ji,jl ) = 0._wp727 za_i (ji,jl ) = 0._wp728 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )729 END IF730 END DO731 ENDIF732 !733 ENDIF734 !735 ! Compatibility tests736 zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )737 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation738 !739 zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )740 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation741 !742 IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?743 !744 itest(4) = 1745 DO jl = 1, i_fill746 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations747 END DO748 ! !----------------------------749 END DO ! end iteration on categories750 ! !----------------------------751 ENDIF752 END DO753 754 ! Add Snow in each category where za_i is not 0755 DO jl = 1, jpl756 DO ji = 1, idim757 IF( za_i(ji,jl) > 0._wp ) THEN758 zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )759 ! In case snow load is in excess that would lead to transformation from snow to ice760 ! Then, transfer the snow excess into the ice (different from icethd_dh)761 zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )762 ! recompute h_i, h_s avoiding out of bounds values763 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )764 zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )765 ENDIF766 END DO767 END DO768 !769 END SUBROUTINE ice_var_itd770 771 772 SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )773 !!-------------------------------------------------------------------774 !! *** ROUTINE ice_var_itd2 ***775 !!776 !! ** Purpose : converting N-cat ice to jpl ice categories777 !!778 !! ice thickness distribution follows a gaussian law779 !! around the concentration of the most likely ice thickness780 !! (similar as iceistate.F90)781 !!782 !! ** Method: Iterative procedure783 !!784 !! 1) Fill ice cat that correspond to input thicknesses785 !! Find the lowest(jlmin) and highest(jlmax) cat that are filled786 !!787 !! 2) Expand the filling to the cat jlmin-1 and jlmax+1788 !! by removing 25% ice area from jlmin and jlmax (resp.)789 !!790 !! 3) Expand the filling to the empty cat between jlmin and jlmax791 !! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)792 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)793 !!794 !! ** Arguments : zhti: N-cat ice thickness795 !! zhts: N-cat snow depth796 !! zati: N-cat ice concentration797 !!798 !! ** Output : jpl-cat799 !!800 !! (Example of application: BDY forcings when inputs have N-cat /= jpl)801 !!-------------------------------------------------------------------802 INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices803 INTEGER :: idim, icat804 REAL(wp), PARAMETER :: ztrans = 0.25_wp805 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables806 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables807 INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2808 INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin809 !!-------------------------------------------------------------------810 !811 idim = SIZE( zhti, 1 )812 icat = SIZE( zhti, 2 )813 !814 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays815 ALLOCATE( jlmin(idim), jlmax(idim) )816 817 ! --- initialize output fields to 0 --- !818 zh_i(1:idim,1:jpl) = 0._wp819 zh_s(1:idim,1:jpl) = 0._wp820 za_i(1:idim,1:jpl) = 0._wp821 !822 ! --- fill the categories --- !823 ! find where cat-input = cat-output and fill cat-output fields824 jlmax(:) = 0825 jlmin(:) = 999826 jlfil(:,:) = 0827 DO jl1 = 1, jpl828 DO jl2 = 1, icat829 DO ji = 1, idim830 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN831 ! fill the right category832 zh_i(ji,jl1) = zhti(ji,jl2)833 zh_s(ji,jl1) = zhts(ji,jl2)834 za_i(ji,jl1) = zati(ji,jl2)835 ! record categories that are filled836 jlmax(ji) = MAX( jlmax(ji), jl1 )837 jlmin(ji) = MIN( jlmin(ji), jl1 )838 jlfil(ji,jl1) = jl1839 ENDIF840 END DO841 END DO842 END DO843 !844 ! --- fill the gaps between categories --- !845 ! transfer from categories filled at the previous step to the empty ones in between846 DO ji = 1, idim847 jl1 = jlmin(ji)848 jl2 = jlmax(ji)849 IF( jl1 > 1 ) THEN850 ! fill the lower cat (jl1-1)851 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)852 zh_i(ji,jl1-1) = hi_mean(jl1-1)853 ! remove from cat jl1854 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)855 ENDIF856 IF( jl2 < jpl ) THEN857 ! fill the upper cat (jl2+1)858 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)859 zh_i(ji,jl2+1) = hi_mean(jl2+1)860 ! remove from cat jl2861 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)862 ENDIF863 END DO864 !865 jlfil2(:,:) = jlfil(:,:)866 ! fill categories from low to high867 DO jl = 2, jpl-1868 DO ji = 1, idim869 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN870 ! fill high871 za_i(ji,jl) = ztrans * za_i(ji,jl-1)872 zh_i(ji,jl) = hi_mean(jl)873 jlfil(ji,jl) = jl874 ! remove low875 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)876 ENDIF877 END DO878 END DO879 !880 ! fill categories from high to low881 DO jl = jpl-1, 2, -1882 DO ji = 1, idim883 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN884 ! fill low885 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)886 zh_i(ji,jl) = hi_mean(jl)887 jlfil2(ji,jl) = jl888 ! remove high889 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)890 ENDIF891 END DO892 END DO893 !894 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays895 DEALLOCATE( jlmin, jlmax )896 !897 END SUBROUTINE ice_var_itd2898 899 674 900 675 SUBROUTINE ice_var_bv … … 958 733 END SUBROUTINE ice_var_enthalpy 959 734 735 960 736 FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 961 737 !!--------------------------------------------------------------------- … … 1004 780 END FUNCTION ice_var_sshdyn 1005 781 782 783 !!------------------------------------------------------------------- 784 !! *** INTERFACE ice_var_itd *** 785 !! 786 !! ** Purpose : converting N-cat ice to jpl ice categories 787 !!------------------------------------------------------------------- 788 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, & 789 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 790 !!------------------------------------------------------------------- 791 !! ** Purpose : converting 1-cat ice to 1 ice category 792 !!------------------------------------------------------------------- 793 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 794 REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 795 REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 796 REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 797 !!------------------------------------------------------------------- 798 ! == thickness and concentration == ! 799 ph_i(:) = phti(:) 800 ph_s(:) = phts(:) 801 pa_i(:) = pati(:) 802 ! 803 ! == temperature and salinity and ponds == ! 804 pt_i (:) = ptmi (:) 805 pt_s (:) = ptms (:) 806 pt_su(:) = ptmsu(:) 807 ps_i (:) = psmi (:) 808 pa_ip(:) = patip(:) 809 ph_ip(:) = phtip(:) 810 811 END SUBROUTINE ice_var_itd_1c1c 812 813 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, & 814 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 815 !!------------------------------------------------------------------- 816 !! ** Purpose : converting N-cat ice to 1 ice category 817 !!------------------------------------------------------------------- 818 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 819 REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 820 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 821 REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 822 ! 823 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs 824 ! 825 INTEGER :: idim 826 !!------------------------------------------------------------------- 827 ! 828 idim = SIZE( phti, 1 ) 829 ! 830 ! == thickness and concentration == ! 831 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) ) 832 ! 833 pa_i(:) = SUM( pati(:,:), dim=2 ) 834 835 WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:) 836 ELSEWHERE ; z1_ai(:) = 0._wp 837 END WHERE 838 839 ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 840 ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 841 ! 842 ! == temperature and salinity == ! 843 WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 844 ELSEWHERE ; z1_vi(:) = 0._wp 845 END WHERE 846 WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 847 ELSEWHERE ; z1_vs(:) = 0._wp 848 END WHERE 849 pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 850 pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 851 pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 852 ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 853 854 ! == ponds == ! 855 pa_ip(:) = SUM( patip(:,:), dim=2 ) 856 WHERE( pa_ip(:) /= 0._wp ) ; ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 857 ELSEWHERE ; ph_ip(:) = 0._wp 858 END WHERE 859 ! 860 DEALLOCATE( z1_ai, z1_vi, z1_vs ) 861 ! 862 END SUBROUTINE ice_var_itd_Nc1c 863 864 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 865 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 866 !!------------------------------------------------------------------- 867 !! 868 !! ** Purpose : converting 1-cat ice to jpl ice categories 869 !! 870 !! 871 !! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015) 872 !! it has the property of conserving total concentration and volume 873 !! 874 !! 875 !! ** Arguments : phti: 1-cat ice thickness 876 !! phts: 1-cat snow depth 877 !! pati: 1-cat ice concentration 878 !! 879 !! ** Output : jpl-cat 880 !! 881 !! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015. 882 !! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice. 883 !! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614 884 !!------------------------------------------------------------------- 885 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 886 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 887 REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 888 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 889 ! 890 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti 891 INTEGER :: ji, jk, jl 892 INTEGER :: idim 893 REAL(wp) :: zv, zdh 894 !!------------------------------------------------------------------- 895 ! 896 idim = SIZE( phti , 1 ) 897 ! 898 ph_i(1:idim,1:jpl) = 0._wp 899 ph_s(1:idim,1:jpl) = 0._wp 900 pa_i(1:idim,1:jpl) = 0._wp 901 ! 902 ALLOCATE( z1_hti(idim) ) 903 WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:) 904 ELSEWHERE ; z1_hti(:) = 0._wp 905 END WHERE 906 ! 907 ! == thickness and concentration == ! 908 ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl) 909 DO jl = 1, jpl-1 910 DO ji = 1, idim 911 ! 912 IF( phti(ji) > 0._wp ) THEN 913 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 914 pa_i(ji,jl) = pati(ji) * z1_hti(ji) * ( ( phti(ji) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 915 & - ( phti(ji) + 2.*hi_max(jl ) ) * EXP( -2.*hi_max(jl )*z1_hti(ji) ) ) 916 ! 917 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 918 zv = pati(ji) * z1_hti(ji) * ( ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl-1) + 2.*hi_max(jl-1)*hi_max(jl-1) ) & 919 & * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 920 & - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) & 921 & * EXP(-2.*hi_max(jl)*z1_hti(ji)) ) 922 ! thickness 923 IF( pa_i(ji,jl) > epsi06 ) THEN 924 ph_i(ji,jl) = zv / pa_i(ji,jl) 925 ELSE 926 ph_i(ji,jl) = 0. 927 pa_i(ji,jl) = 0. 928 ENDIF 929 ENDIF 930 ! 931 ENDDO 932 ENDDO 933 ! 934 ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity 935 DO ji = 1, idim 936 ! 937 IF( phti(ji) > 0._wp ) THEN 938 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity 939 pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 940 941 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity 942 zv = pati(ji) * z1_hti(ji) * ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jpl-1) + 2.*hi_max(jpl-1)*hi_max(jpl-1) ) & 943 & * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 944 ! thickness 945 IF( pa_i(ji,jpl) > epsi06 ) THEN 946 ph_i(ji,jpl) = zv / pa_i(ji,jpl) 947 else 948 ph_i(ji,jpl) = 0. 949 pa_i(ji,jpl) = 0. 950 ENDIF 951 ENDIF 952 ! 953 ENDDO 954 ! 955 ! Add Snow in each category where pa_i is not 0 956 DO jl = 1, jpl 957 DO ji = 1, idim 958 IF( pa_i(ji,jl) > 0._wp ) THEN 959 ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 960 ! In case snow load is in excess that would lead to transformation from snow to ice 961 ! Then, transfer the snow excess into the ice (different from icethd_dh) 962 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 963 ! recompute h_i, h_s avoiding out of bounds values 964 ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) 965 ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos ) 966 ENDIF 967 END DO 968 END DO 969 ! 970 DEALLOCATE( z1_hti ) 971 ! 972 ! == temperature and salinity == ! 973 DO jl = 1, jpl 974 pt_i (:,jl) = ptmi (:) 975 pt_s (:,jl) = ptms (:) 976 pt_su(:,jl) = ptmsu(:) 977 ps_i (:,jl) = psmi (:) 978 ps_i (:,jl) = psmi (:) 979 END DO 980 ! 981 ! == ponds == ! 982 ALLOCATE( zfra(idim) ) 983 ! keep the same pond fraction atip/ati for each category 984 WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:) 985 ELSEWHERE ; zfra(:) = 0._wp 986 END WHERE 987 DO jl = 1, jpl 988 pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 989 END DO 990 ! keep the same v_ip/v_i ratio for each category 991 WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) ) 992 ELSEWHERE ; zfra(:) = 0._wp 993 END WHERE 994 DO jl = 1, jpl 995 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 996 ELSEWHERE ; ph_ip(:,jl) = 0._wp 997 END WHERE 998 END DO 999 DEALLOCATE( zfra ) 1000 ! 1001 END SUBROUTINE ice_var_itd_1cMc 1002 1003 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, & 1004 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 1005 !!------------------------------------------------------------------- 1006 !! 1007 !! ** Purpose : converting N-cat ice to jpl ice categories 1008 !! 1009 !! ice thickness distribution follows a gaussian law 1010 !! around the concentration of the most likely ice thickness 1011 !! (similar as iceistate.F90) 1012 !! 1013 !! ** Method: Iterative procedure 1014 !! 1015 !! 1) Fill ice cat that correspond to input thicknesses 1016 !! Find the lowest(jlmin) and highest(jlmax) cat that are filled 1017 !! 1018 !! 2) Expand the filling to the cat jlmin-1 and jlmax+1 1019 !! by removing 25% ice area from jlmin and jlmax (resp.) 1020 !! 1021 !! 3) Expand the filling to the empty cat between jlmin and jlmax 1022 !! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax) 1023 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 1024 !! 1025 !! ** Arguments : phti: N-cat ice thickness 1026 !! phts: N-cat snow depth 1027 !! pati: N-cat ice concentration 1028 !! 1029 !! ** Output : jpl-cat 1030 !! 1031 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 1032 !!------------------------------------------------------------------- 1033 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 1034 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 1035 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 1036 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 1037 ! 1038 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2 1039 INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin 1040 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra 1041 ! 1042 REAL(wp), PARAMETER :: ztrans = 0.25_wp 1043 INTEGER :: ji, jl, jl1, jl2 1044 INTEGER :: idim, icat 1045 !!------------------------------------------------------------------- 1046 ! 1047 idim = SIZE( phti, 1 ) 1048 icat = SIZE( phti, 2 ) 1049 ! 1050 ! == thickness and concentration == ! 1051 ! ! ---------------------- ! 1052 IF( icat == jpl ) THEN ! input cat = output cat ! 1053 ! ! ---------------------- ! 1054 ph_i(:,:) = phti(:,:) 1055 ph_s(:,:) = phts(:,:) 1056 pa_i(:,:) = pati(:,:) 1057 ! 1058 ! == temperature and salinity and ponds == ! 1059 pt_i (:,:) = ptmi (:,:) 1060 pt_s (:,:) = ptms (:,:) 1061 pt_su(:,:) = ptmsu(:,:) 1062 ps_i (:,:) = psmi (:,:) 1063 pa_ip(:,:) = patip(:,:) 1064 ph_ip(:,:) = phtip(:,:) 1065 ! ! ---------------------- ! 1066 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! 1067 ! ! ---------------------- ! 1068 CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 1069 & ph_i(:,:), ph_s(:,:), pa_i (:,:), & 1070 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 1071 & pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:) ) 1072 ! ! ---------------------- ! 1073 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! 1074 ! ! ---------------------- ! 1075 CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 1076 & ph_i(:,1), ph_s(:,1), pa_i (:,1), & 1077 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 1078 & pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1) ) 1079 ! ! ----------------------- ! 1080 ELSE ! input cat /= output cat ! 1081 ! ! ----------------------- ! 1082 1083 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays 1084 ALLOCATE( jlmin(idim), jlmax(idim) ) 1085 1086 ! --- initialize output fields to 0 --- ! 1087 ph_i(1:idim,1:jpl) = 0._wp 1088 ph_s(1:idim,1:jpl) = 0._wp 1089 pa_i(1:idim,1:jpl) = 0._wp 1090 ! 1091 ! --- fill the categories --- ! 1092 ! find where cat-input = cat-output and fill cat-output fields 1093 jlmax(:) = 0 1094 jlmin(:) = 999 1095 jlfil(:,:) = 0 1096 DO jl1 = 1, jpl 1097 DO jl2 = 1, icat 1098 DO ji = 1, idim 1099 IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 1100 ! fill the right category 1101 ph_i(ji,jl1) = phti(ji,jl2) 1102 ph_s(ji,jl1) = phts(ji,jl2) 1103 pa_i(ji,jl1) = pati(ji,jl2) 1104 ! record categories that are filled 1105 jlmax(ji) = MAX( jlmax(ji), jl1 ) 1106 jlmin(ji) = MIN( jlmin(ji), jl1 ) 1107 jlfil(ji,jl1) = jl1 1108 ENDIF 1109 END DO 1110 END DO 1111 END DO 1112 ! 1113 ! --- fill the gaps between categories --- ! 1114 ! transfer from categories filled at the previous step to the empty ones in between 1115 DO ji = 1, idim 1116 jl1 = jlmin(ji) 1117 jl2 = jlmax(ji) 1118 IF( jl1 > 1 ) THEN 1119 ! fill the lower cat (jl1-1) 1120 pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 1121 ph_i(ji,jl1-1) = hi_mean(jl1-1) 1122 ! remove from cat jl1 1123 pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 1124 ENDIF 1125 IF( jl2 < jpl ) THEN 1126 ! fill the upper cat (jl2+1) 1127 pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 1128 ph_i(ji,jl2+1) = hi_mean(jl2+1) 1129 ! remove from cat jl2 1130 pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 1131 ENDIF 1132 END DO 1133 ! 1134 jlfil2(:,:) = jlfil(:,:) 1135 ! fill categories from low to high 1136 DO jl = 2, jpl-1 1137 DO ji = 1, idim 1138 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 1139 ! fill high 1140 pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 1141 ph_i(ji,jl) = hi_mean(jl) 1142 jlfil(ji,jl) = jl 1143 ! remove low 1144 pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 1145 ENDIF 1146 END DO 1147 END DO 1148 ! 1149 ! fill categories from high to low 1150 DO jl = jpl-1, 2, -1 1151 DO ji = 1, idim 1152 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 1153 ! fill low 1154 pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 1155 ph_i(ji,jl) = hi_mean(jl) 1156 jlfil2(ji,jl) = jl 1157 ! remove high 1158 pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 1159 ENDIF 1160 END DO 1161 END DO 1162 ! 1163 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 1164 DEALLOCATE( jlmin, jlmax ) 1165 ! 1166 ! == temperature and salinity == ! 1167 ! 1168 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 1169 ! 1170 WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 1171 ELSEWHERE ; z1_ai(:) = 0._wp 1172 END WHERE 1173 WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 1174 ELSEWHERE ; z1_vi(:) = 0._wp 1175 END WHERE 1176 WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 1177 ELSEWHERE ; z1_vs(:) = 0._wp 1178 END WHERE 1179 ! 1180 ! fill all the categories with the same value 1181 ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1182 DO jl = 1, jpl 1183 pt_i (:,jl) = ztmp(:) 1184 END DO 1185 ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 1186 DO jl = 1, jpl 1187 pt_s (:,jl) = ztmp(:) 1188 END DO 1189 ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 1190 DO jl = 1, jpl 1191 pt_su(:,jl) = ztmp(:) 1192 END DO 1193 ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1194 DO jl = 1, jpl 1195 ps_i (:,jl) = ztmp(:) 1196 END DO 1197 ! 1198 DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 1199 ! 1200 ! == ponds == ! 1201 ALLOCATE( zfra(idim) ) 1202 ! keep the same pond fraction atip/ati for each category 1203 WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 ) 1204 ELSEWHERE ; zfra(:) = 0._wp 1205 END WHERE 1206 DO jl = 1, jpl 1207 pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 1208 END DO 1209 ! keep the same v_ip/v_i ratio for each category 1210 WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 1211 zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 1212 ELSEWHERE 1213 zfra(:) = 0._wp 1214 END WHERE 1215 DO jl = 1, jpl 1216 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 1217 ELSEWHERE ; ph_ip(:,jl) = 0._wp 1218 END WHERE 1219 END DO 1220 DEALLOCATE( zfra ) 1221 ! 1222 ENDIF 1223 ! 1224 END SUBROUTINE ice_var_itd_NcMc 1006 1225 1007 1226 #else
Note: See TracChangeset
for help on using the changeset viewer.
