- Timestamp:
- 2019-12-11T12:02:38+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/ICE/icevar.F90
r10994 r12178 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 … … 46 47 !! ice_var_zapneg : remove negative ice fields 47 48 !! ice_var_roundoff : remove negative values arising from roundoff erros 48 !! ice_var_itd : convert 1-cat to jpl-cat49 !! ice_var_itd2 : convert N-cat to jpl-cat50 49 !! ice_var_bv : brine volume 51 50 !! ice_var_enthalpy : compute ice and snow enthalpies from temperature 52 51 !! ice_var_sshdyn : compute equivalent ssh in lead 52 !! ice_var_itd : convert N-cat to M-cat 53 53 !!---------------------------------------------------------------------- 54 54 USE dom_oce ! ocean space and time domain … … 73 73 PUBLIC ice_var_zapneg 74 74 PUBLIC ice_var_roundoff 75 PUBLIC ice_var_itd76 PUBLIC ice_var_itd277 75 PUBLIC ice_var_bv 78 76 PUBLIC ice_var_enthalpy 79 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 80 83 81 84 !!---------------------------------------------------------------------- … … 101 104 ! 102 105 ! ! integrated values 103 vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 ) 104 vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 ) 105 at_i(:,:) = SUM( a_i(:,:,:) , dim=3 ) 106 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 107 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 ) 108 112 ! 109 113 at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds … … 135 139 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 136 140 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 137 sm_i (:,:) = SUM( sv_i(:,:,:) , dim=3 )* z1_vt_i(:,:)141 sm_i (:,:) = st_i(:,:) * z1_vt_i(:,:) 138 142 ! 139 143 tm_i(:,:) = 0._wp … … 155 159 tm_s (:,:) = rt0 156 160 END WHERE 157 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 ! 158 167 DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 168 ! 159 169 ENDIF 160 170 ! … … 260 270 ! 261 271 ! integrated values 262 vt_i (:,:) = SUM( v_i , dim=3 )263 vt_s (:,:) = SUM( v_s , dim=3 )264 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 ) 265 275 ! 266 276 END SUBROUTINE ice_var_glo2eqv … … 530 540 531 541 ! to be sure that at_i is the sum of a_i(jl) 532 at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 533 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 534 550 535 551 ! open water = 1 if at_i=0 … … 649 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 650 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 651 IF 667 IF( ln_pnd_H12 ) THEN 652 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 653 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 … … 656 672 END SUBROUTINE ice_var_roundoff 657 673 658 659 SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )660 !!-------------------------------------------------------------------661 !! *** ROUTINE ice_var_itd ***662 !!663 !! ** Purpose : converting 1-cat ice to multiple ice categories664 !!665 !! ice thickness distribution follows a gaussian law666 !! around the concentration of the most likely ice thickness667 !! (similar as iceistate.F90)668 !!669 !! ** Method: Iterative procedure670 !!671 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian672 !!673 !! 2) Check whether the distribution conserves area and volume, positivity and674 !! category boundaries675 !!676 !! 3) If not (input ice is too thin), the last category is empty and677 !! the number of categories is reduced (jpl-1)678 !!679 !! 4) Iterate until ok (SUM(itest(:) = 4)680 !!681 !! ** Arguments : zhti: 1-cat ice thickness682 !! zhts: 1-cat snow depth683 !! zati: 1-cat ice concentration684 !!685 !! ** Output : jpl-cat686 !!687 !! (Example of application: BDY forcings when input are cell averaged)688 !!-------------------------------------------------------------------689 INTEGER :: ji, jk, jl ! dummy loop indices690 INTEGER :: idim, i_fill, jl0691 REAL(wp) :: zarg, zV, zconv, zdh, zdv692 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables693 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables694 INTEGER , DIMENSION(4) :: itest695 !!-------------------------------------------------------------------696 !697 ! ----------------------------------------698 ! distribution over the jpl ice categories699 ! ----------------------------------------700 ! a gaussian distribution for ice concentration is used701 ! then we check whether the distribution fullfills702 ! volume and area conservation, positivity and ice categories bounds703 idim = SIZE( zhti , 1 )704 zh_i(1:idim,1:jpl) = 0._wp705 zh_s(1:idim,1:jpl) = 0._wp706 za_i(1:idim,1:jpl) = 0._wp707 !708 DO ji = 1, idim709 !710 IF( zhti(ji) > 0._wp ) THEN711 !712 ! find which category (jl0) the input ice thickness falls into713 jl0 = jpl714 DO jl = 1, jpl715 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN716 jl0 = jl717 CYCLE718 ENDIF719 END DO720 !721 itest(:) = 0722 i_fill = jpl + 1 !------------------------------------723 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories724 ! !------------------------------------725 i_fill = i_fill - 1726 !727 zh_i(ji,1:jpl) = 0._wp728 za_i(ji,1:jpl) = 0._wp729 itest(:) = 0730 !731 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1732 zh_i(ji,1) = zhti(ji)733 za_i (ji,1) = zati (ji)734 ELSE !-- case ice is thicker: fill categories >1735 ! thickness736 DO jl = 1, i_fill - 1737 zh_i(ji,jl) = hi_mean(jl)738 END DO739 !740 ! concentration741 za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))742 DO jl = 1, i_fill - 1743 IF ( jl /= jl0 ) THEN744 zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )745 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)746 ENDIF747 END DO748 !749 ! last category750 za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )751 zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )752 zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )753 !754 ! correction if concentration of upper cat is greater than lower cat755 ! (it should be a gaussian around jl0 but sometimes it is not)756 IF ( jl0 /= jpl ) THEN757 DO jl = jpl, jl0+1, -1758 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN759 zdv = zh_i(ji,jl) * za_i(ji,jl)760 zh_i(ji,jl ) = 0._wp761 za_i (ji,jl ) = 0._wp762 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )763 END IF764 END DO765 ENDIF766 !767 ENDIF768 !769 ! Compatibility tests770 zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )771 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation772 !773 zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )774 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation775 !776 IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?777 !778 itest(4) = 1779 DO jl = 1, i_fill780 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations781 END DO782 ! !----------------------------783 END DO ! end iteration on categories784 ! !----------------------------785 ENDIF786 END DO787 788 ! Add Snow in each category where za_i is not 0789 DO jl = 1, jpl790 DO ji = 1, idim791 IF( za_i(ji,jl) > 0._wp ) THEN792 zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )793 ! In case snow load is in excess that would lead to transformation from snow to ice794 ! Then, transfer the snow excess into the ice (different from icethd_dh)795 zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )796 ! recompute h_i, h_s avoiding out of bounds values797 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )798 zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )799 ENDIF800 END DO801 END DO802 !803 END SUBROUTINE ice_var_itd804 805 806 SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )807 !!-------------------------------------------------------------------808 !! *** ROUTINE ice_var_itd2 ***809 !!810 !! ** Purpose : converting N-cat ice to jpl ice categories811 !!812 !! ice thickness distribution follows a gaussian law813 !! around the concentration of the most likely ice thickness814 !! (similar as iceistate.F90)815 !!816 !! ** Method: Iterative procedure817 !!818 !! 1) Fill ice cat that correspond to input thicknesses819 !! Find the lowest(jlmin) and highest(jlmax) cat that are filled820 !!821 !! 2) Expand the filling to the cat jlmin-1 and jlmax+1822 !! by removing 25% ice area from jlmin and jlmax (resp.)823 !!824 !! 3) Expand the filling to the empty cat between jlmin and jlmax825 !! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)826 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)827 !!828 !! ** Arguments : zhti: N-cat ice thickness829 !! zhts: N-cat snow depth830 !! zati: N-cat ice concentration831 !!832 !! ** Output : jpl-cat833 !!834 !! (Example of application: BDY forcings when inputs have N-cat /= jpl)835 !!-------------------------------------------------------------------836 INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices837 INTEGER :: idim, icat838 REAL(wp), PARAMETER :: ztrans = 0.25_wp839 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables840 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables841 INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2842 INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin843 !!-------------------------------------------------------------------844 !845 idim = SIZE( zhti, 1 )846 icat = SIZE( zhti, 2 )847 !848 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays849 ALLOCATE( jlmin(idim), jlmax(idim) )850 851 ! --- initialize output fields to 0 --- !852 zh_i(1:idim,1:jpl) = 0._wp853 zh_s(1:idim,1:jpl) = 0._wp854 za_i(1:idim,1:jpl) = 0._wp855 !856 ! --- fill the categories --- !857 ! find where cat-input = cat-output and fill cat-output fields858 jlmax(:) = 0859 jlmin(:) = 999860 jlfil(:,:) = 0861 DO jl1 = 1, jpl862 DO jl2 = 1, icat863 DO ji = 1, idim864 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN865 ! fill the right category866 zh_i(ji,jl1) = zhti(ji,jl2)867 zh_s(ji,jl1) = zhts(ji,jl2)868 za_i(ji,jl1) = zati(ji,jl2)869 ! record categories that are filled870 jlmax(ji) = MAX( jlmax(ji), jl1 )871 jlmin(ji) = MIN( jlmin(ji), jl1 )872 jlfil(ji,jl1) = jl1873 ENDIF874 END DO875 END DO876 END DO877 !878 ! --- fill the gaps between categories --- !879 ! transfer from categories filled at the previous step to the empty ones in between880 DO ji = 1, idim881 jl1 = jlmin(ji)882 jl2 = jlmax(ji)883 IF( jl1 > 1 ) THEN884 ! fill the lower cat (jl1-1)885 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)886 zh_i(ji,jl1-1) = hi_mean(jl1-1)887 ! remove from cat jl1888 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)889 ENDIF890 IF( jl2 < jpl ) THEN891 ! fill the upper cat (jl2+1)892 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)893 zh_i(ji,jl2+1) = hi_mean(jl2+1)894 ! remove from cat jl2895 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)896 ENDIF897 END DO898 !899 jlfil2(:,:) = jlfil(:,:)900 ! fill categories from low to high901 DO jl = 2, jpl-1902 DO ji = 1, idim903 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN904 ! fill high905 za_i(ji,jl) = ztrans * za_i(ji,jl-1)906 zh_i(ji,jl) = hi_mean(jl)907 jlfil(ji,jl) = jl908 ! remove low909 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)910 ENDIF911 END DO912 END DO913 !914 ! fill categories from high to low915 DO jl = jpl-1, 2, -1916 DO ji = 1, idim917 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN918 ! fill low919 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)920 zh_i(ji,jl) = hi_mean(jl)921 jlfil2(ji,jl) = jl922 ! remove high923 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)924 ENDIF925 END DO926 END DO927 !928 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays929 DEALLOCATE( jlmin, jlmax )930 !931 END SUBROUTINE ice_var_itd2932 933 674 934 675 SUBROUTINE ice_var_bv … … 992 733 END SUBROUTINE ice_var_enthalpy 993 734 735 994 736 FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 995 737 !!--------------------------------------------------------------------- … … 1038 780 END FUNCTION ice_var_sshdyn 1039 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 1040 1225 1041 1226 #else
Note: See TracChangeset
for help on using the changeset viewer.