Changeset 11573 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/icevar.F90
- Timestamp:
- 2019-09-19T11:18:03+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/ICE/icevar.F90
r11229 r11573 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 … … 104 104 ! 105 105 ! ! integrated values 106 vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 ) 107 vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 ) 108 at_i(:,:) = SUM( a_i(:,:,:) , dim=3 ) 109 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 110 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 ) 111 112 ! 112 113 at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds … … 138 139 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 139 140 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 140 sm_i (:,:) = SUM( sv_i(:,:,:) , dim=3 )* z1_vt_i(:,:)141 sm_i (:,:) = st_i(:,:) * z1_vt_i(:,:) 141 142 ! 142 143 tm_i(:,:) = 0._wp … … 158 159 tm_s (:,:) = rt0 159 160 END WHERE 160 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 ! 161 167 DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 168 ! 162 169 ENDIF 163 170 ! … … 263 270 ! 264 271 ! integrated values 265 vt_i (:,:) = SUM( v_i , dim=3 )266 vt_s (:,:) = SUM( v_s , dim=3 )267 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 ) 268 275 ! 269 276 END SUBROUTINE ice_var_glo2eqv … … 533 540 534 541 ! to be sure that at_i is the sum of a_i(jl) 535 at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 536 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 537 550 538 551 ! open water = 1 if at_i=0 … … 652 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 653 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 654 IF 667 IF( ln_pnd_H12 ) THEN 655 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 656 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 … … 773 786 !! ** Purpose : converting N-cat ice to jpl ice categories 774 787 !!------------------------------------------------------------------- 775 SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 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 ) 776 790 !!------------------------------------------------------------------- 777 791 !! ** Purpose : converting 1-cat ice to 1 ice category 778 792 !!------------------------------------------------------------------- 779 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 780 REAL(wp), DIMENSION(:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 781 !!------------------------------------------------------------------- 782 zh_i(:) = zhti(:) 783 zh_s(:) = zhts(:) 784 za_i(:) = zati(:) 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 785 811 END SUBROUTINE ice_var_itd_1c1c 786 812 787 SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 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 ) 788 815 !!------------------------------------------------------------------- 789 816 !! ** Purpose : converting N-cat ice to 1 ice category 790 817 !!------------------------------------------------------------------- 791 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 792 REAL(wp), DIMENSION(:) , INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 793 !!------------------------------------------------------------------- 794 ! 795 za_i(:) = SUM( zati(:,:), dim=2 ) 796 ! 797 WHERE( za_i(:) /= 0._wp ) 798 zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:) 799 zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:) 800 ELSEWHERE 801 zh_i(:) = 0._wp 802 zh_s(:) = 0._wp 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 803 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 ) 804 861 ! 805 862 END SUBROUTINE ice_var_itd_Nc1c 806 863 807 SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 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 ) 808 866 !!------------------------------------------------------------------- 809 867 !! 810 868 !! ** Purpose : converting 1-cat ice to jpl ice categories 811 869 !! 812 !! ice thickness distribution follows a gaussian law 813 !! around the concentration of the most likely ice thickness 814 !! (similar as iceistate.F90) 815 !! 816 !! ** Method: Iterative procedure 817 !! 818 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 819 !! 820 !! 2) Check whether the distribution conserves area and volume, positivity and 821 !! category boundaries 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 822 873 !! 823 !! 3) If not (input ice is too thin), the last category is empty and 824 !! the number of categories is reduced (jpl-1) 825 !! 826 !! 4) Iterate until ok (SUM(itest(:) = 4) 827 !! 828 !! ** Arguments : zhti: 1-cat ice thickness 829 !! zhts: 1-cat snow depth 830 !! zati: 1-cat ice concentration 874 !! 875 !! ** Arguments : phti: 1-cat ice thickness 876 !! phts: 1-cat snow depth 877 !! pati: 1-cat ice concentration 831 878 !! 832 879 !! ** Output : jpl-cat 833 880 !! 834 !! (Example of application: BDY forcings when input are cell averaged) 835 !!------------------------------------------------------------------- 836 INTEGER :: ji, jk, jl ! dummy loop indices 837 INTEGER :: idim, i_fill, jl0 838 REAL(wp) :: zarg, zV, zconv, zdh, zdv 839 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 840 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 841 INTEGER , DIMENSION(4) :: itest 842 !!------------------------------------------------------------------- 843 ! 844 ! ---------------------------------------- 845 ! distribution over the jpl ice categories 846 ! ---------------------------------------- 847 ! a gaussian distribution for ice concentration is used 848 ! then we check whether the distribution fullfills 849 ! volume and area conservation, positivity and ice categories bounds 850 idim = SIZE( zhti , 1 ) 851 zh_i(1:idim,1:jpl) = 0._wp 852 zh_s(1:idim,1:jpl) = 0._wp 853 za_i(1:idim,1:jpl) = 0._wp 854 ! 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 855 935 DO ji = 1, idim 856 936 ! 857 IF( zhti(ji) > 0._wp ) THEN 858 ! 859 ! find which category (jl0) the input ice thickness falls into 860 jl0 = jpl 861 DO jl = 1, jpl 862 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 863 jl0 = jl 864 CYCLE 865 ENDIF 866 END DO 867 ! 868 itest(:) = 0 869 i_fill = jpl + 1 !------------------------------------ 870 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 871 ! !------------------------------------ 872 i_fill = i_fill - 1 873 ! 874 zh_i(ji,1:jpl) = 0._wp 875 za_i(ji,1:jpl) = 0._wp 876 itest(:) = 0 877 ! 878 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1 879 zh_i(ji,1) = zhti(ji) 880 za_i (ji,1) = zati (ji) 881 ELSE !-- case ice is thicker: fill categories >1 882 ! thickness 883 DO jl = 1, i_fill - 1 884 zh_i(ji,jl) = hi_mean(jl) 885 END DO 886 ! 887 ! concentration 888 za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl)) 889 DO jl = 1, i_fill - 1 890 IF ( jl /= jl0 ) THEN 891 zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 892 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 893 ENDIF 894 END DO 895 ! 896 ! last category 897 za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) ) 898 zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) ) 899 zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 900 ! 901 ! correction if concentration of upper cat is greater than lower cat 902 ! (it should be a gaussian around jl0 but sometimes it is not) 903 IF ( jl0 /= jpl ) THEN 904 DO jl = jpl, jl0+1, -1 905 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 906 zdv = zh_i(ji,jl) * za_i(ji,jl) 907 zh_i(ji,jl ) = 0._wp 908 za_i (ji,jl ) = 0._wp 909 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 910 END IF 911 END DO 912 ENDIF 913 ! 914 ENDIF 915 ! 916 ! Compatibility tests 917 zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) ) 918 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation 919 ! 920 zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) ) 921 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation 922 ! 923 IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ? 924 ! 925 itest(4) = 1 926 DO jl = 1, i_fill 927 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations 928 END DO 929 ! !---------------------------- 930 END DO ! end iteration on categories 931 ! !---------------------------- 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 932 951 ENDIF 933 END DO 934 935 ! Add Snow in each category where za_i is not 0 952 ! 953 ENDDO 954 ! 955 ! Add Snow in each category where pa_i is not 0 936 956 DO jl = 1, jpl 937 957 DO ji = 1, idim 938 IF( za_i(ji,jl) > 0._wp ) THEN939 zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji))958 IF( pa_i(ji,jl) > 0._wp ) THEN 959 ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 940 960 ! In case snow load is in excess that would lead to transformation from snow to ice 941 961 ! Then, transfer the snow excess into the ice (different from icethd_dh) 942 zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )962 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 943 963 ! recompute h_i, h_s avoiding out of bounds values 944 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )945 zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )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 ) 946 966 ENDIF 947 967 END DO 948 968 END DO 949 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 ! 950 1001 END SUBROUTINE ice_var_itd_1cMc 951 1002 952 SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 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 ) 953 1005 !!------------------------------------------------------------------- 954 1006 !! … … 971 1023 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 972 1024 !! 973 !! ** Arguments : zhti: N-cat ice thickness974 !! zhts: N-cat snow depth975 !! zati: N-cat ice concentration1025 !! ** Arguments : phti: N-cat ice thickness 1026 !! phts: N-cat snow depth 1027 !! pati: N-cat ice concentration 976 1028 !! 977 1029 !! ** Output : jpl-cat … … 979 1031 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 980 1032 !!------------------------------------------------------------------- 981 INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices 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 982 1044 INTEGER :: idim, icat 983 REAL(wp), PARAMETER :: ztrans = 0.25_wp 984 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 985 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 986 INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2 987 INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin 988 !!------------------------------------------------------------------- 989 ! 990 idim = SIZE( zhti, 1 ) 991 icat = SIZE( zhti, 2 ) 1045 !!------------------------------------------------------------------- 1046 ! 1047 idim = SIZE( phti, 1 ) 1048 icat = SIZE( phti, 2 ) 1049 ! 1050 ! == thickness and concentration == ! 992 1051 ! ! ---------------------- ! 993 1052 IF( icat == jpl ) THEN ! input cat = output cat ! 994 1053 ! ! ---------------------- ! 995 zh_i(:,:) = zhti(:,:) 996 zh_s(:,:) = zhts(:,:) 997 za_i(:,:) = zati(:,:) 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(:,:) 998 1065 ! ! ---------------------- ! 999 1066 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! 1000 1067 ! ! ---------------------- ! 1001 CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i(:,:), zh_s(:,:), za_i(:,:) ) 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(:,:) ) 1002 1072 ! ! ---------------------- ! 1003 1073 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! 1004 1074 ! ! ---------------------- ! 1005 CALL ice_var_itd_Nc1c( zhti(:,:), zhts(:,:), zati(:,:), zh_i(:,1), zh_s(:,1), za_i(:,1) ) 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) ) 1006 1079 ! ! ----------------------- ! 1007 1080 ELSE ! input cat /= output cat ! … … 1012 1085 1013 1086 ! --- initialize output fields to 0 --- ! 1014 zh_i(1:idim,1:jpl) = 0._wp1015 zh_s(1:idim,1:jpl) = 0._wp1016 za_i(1:idim,1:jpl) = 0._wp1087 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 1017 1090 ! 1018 1091 ! --- fill the categories --- ! … … 1024 1097 DO jl2 = 1, icat 1025 1098 DO ji = 1, idim 1026 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN1099 IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 1027 1100 ! fill the right category 1028 zh_i(ji,jl1) = zhti(ji,jl2)1029 zh_s(ji,jl1) = zhts(ji,jl2)1030 za_i(ji,jl1) = zati(ji,jl2)1101 ph_i(ji,jl1) = phti(ji,jl2) 1102 ph_s(ji,jl1) = phts(ji,jl2) 1103 pa_i(ji,jl1) = pati(ji,jl2) 1031 1104 ! record categories that are filled 1032 1105 jlmax(ji) = MAX( jlmax(ji), jl1 ) … … 1045 1118 IF( jl1 > 1 ) THEN 1046 1119 ! fill the lower cat (jl1-1) 1047 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)1048 zh_i(ji,jl1-1) = hi_mean(jl1-1)1120 pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 1121 ph_i(ji,jl1-1) = hi_mean(jl1-1) 1049 1122 ! remove from cat jl1 1050 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)1123 pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 1051 1124 ENDIF 1052 1125 IF( jl2 < jpl ) THEN 1053 1126 ! fill the upper cat (jl2+1) 1054 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)1055 zh_i(ji,jl2+1) = hi_mean(jl2+1)1127 pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 1128 ph_i(ji,jl2+1) = hi_mean(jl2+1) 1056 1129 ! remove from cat jl2 1057 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)1130 pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 1058 1131 ENDIF 1059 1132 END DO … … 1065 1138 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 1066 1139 ! fill high 1067 za_i(ji,jl) = ztrans * za_i(ji,jl-1)1068 zh_i(ji,jl) = hi_mean(jl)1140 pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 1141 ph_i(ji,jl) = hi_mean(jl) 1069 1142 jlfil(ji,jl) = jl 1070 1143 ! remove low 1071 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)1144 pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 1072 1145 ENDIF 1073 1146 END DO … … 1079 1152 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 1080 1153 ! fill low 1081 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)1082 zh_i(ji,jl) = hi_mean(jl)1154 pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 1155 ph_i(ji,jl) = hi_mean(jl) 1083 1156 jlfil2(ji,jl) = jl 1084 1157 ! remove high 1085 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)1158 pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 1086 1159 ENDIF 1087 1160 END DO … … 1090 1163 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 1091 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 ) 1092 1221 ! 1093 1222 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.