- Timestamp:
- 2019-07-25T14:02:55+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90
r11334 r11348 778 778 !! ** Purpose : converting N-cat ice to jpl ice categories 779 779 !!------------------------------------------------------------------- 780 SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 780 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, & 781 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 781 782 !!------------------------------------------------------------------- 782 783 !! ** Purpose : converting 1-cat ice to 1 ice category 783 784 !!------------------------------------------------------------------- 784 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 785 REAL(wp), DIMENSION(:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 786 !!------------------------------------------------------------------- 787 zh_i(:) = zhti(:) 788 zh_s(:) = zhts(:) 789 za_i(:) = zati(:) 785 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 786 REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 787 REAL(wp), DIMENSION(:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 788 REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 789 !!------------------------------------------------------------------- 790 ! == thickness and concentration == ! 791 ph_i(:) = phti(:) 792 ph_s(:) = phts(:) 793 pa_i(:) = pati(:) 794 ! 795 ! == temperature and salinity == ! 796 IF( PRESENT( pt_i ) ) pt_i (:) = ptmi (:) 797 IF( PRESENT( pt_s ) ) pt_s (:) = ptms (:) 798 IF( PRESENT( pt_su ) ) pt_su(:) = ptmsu(:) 799 IF( PRESENT( ps_i ) ) ps_i (:) = psmi (:) 800 790 801 END SUBROUTINE ice_var_itd_1c1c 791 802 792 SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 803 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, & 804 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 793 805 !!------------------------------------------------------------------- 794 806 !! ** Purpose : converting N-cat ice to 1 ice category 795 807 !!------------------------------------------------------------------- 796 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 797 REAL(wp), DIMENSION(:) , INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 798 !!------------------------------------------------------------------- 799 ! 800 za_i(:) = SUM( zati(:,:), dim=2 ) 801 ! 802 WHERE( za_i(:) /= 0._wp ) 803 zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:) 804 zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:) 805 ELSEWHERE 806 zh_i(:) = 0._wp 807 zh_s(:) = 0._wp 808 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 809 REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 810 REAL(wp), DIMENSION(:,:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 811 REAL(wp), DIMENSION(:) , INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 812 ! 813 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs 814 ! 815 INTEGER :: idim 816 !!------------------------------------------------------------------- 817 ! 818 idim = SIZE( phti, 1 ) 819 ! 820 ! == thickness and concentration == ! 821 ALLOCATE( z1_ai(idim) ) 822 ! 823 pa_i(:) = SUM( pati(:,:), dim=2 ) 824 825 WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:) 826 ELSEWHERE ; z1_ai(:) = 0._wp 808 827 END WHERE 828 829 ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 830 ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 831 ! 832 ! == temperature and salinity == ! 833 IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 834 ! 835 ALLOCATE( z1_vi(idim), z1_vs(idim) ) 836 ! 837 WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 838 ELSEWHERE ; z1_vi(:) = 0._wp 839 END WHERE 840 WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 841 ELSEWHERE ; z1_vs(:) = 0._wp 842 END WHERE 843 ! 844 IF( PRESENT( pt_i ) ) pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 845 IF( PRESENT( pt_s ) ) pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 846 IF( PRESENT( pt_su ) ) pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 847 IF( PRESENT( ps_i ) ) ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 848 ! 849 DEALLOCATE( z1_vi, z1_vs ) 850 ! 851 ENDIF 852 ! 853 DEALLOCATE( z1_ai ) 809 854 ! 810 855 END SUBROUTINE ice_var_itd_Nc1c 811 856 812 SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 857 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 858 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 813 859 !!------------------------------------------------------------------- 814 860 !! … … 831 877 !! 4) Iterate until ok (SUM(itest(:) = 4) 832 878 !! 833 !! ** Arguments : zhti: 1-cat ice thickness834 !! zhts: 1-cat snow depth835 !! zati: 1-cat ice concentration879 !! ** Arguments : phti: 1-cat ice thickness 880 !! phts: 1-cat snow depth 881 !! pati: 1-cat ice concentration 836 882 !! 837 883 !! ** Output : jpl-cat … … 839 885 !! (Example of application: BDY forcings when input are cell averaged) 840 886 !!------------------------------------------------------------------- 841 INTEGER :: ji, jk, jl ! dummy loop indices 842 INTEGER :: idim, i_fill, jl0 843 REAL(wp) :: zarg, zV, zconv, zdh, zdv 844 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 845 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 846 INTEGER , DIMENSION(4) :: itest 847 !!------------------------------------------------------------------- 848 ! 849 ! ---------------------------------------- 887 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 888 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 889 REAL(wp), DIMENSION(:) , INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 890 REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 891 ! 892 INTEGER , DIMENSION(4) :: itest 893 INTEGER :: ji, jk, jl 894 INTEGER :: idim, i_fill, jl0 895 REAL(wp) :: zarg, zV, zconv, zdh, zdv 896 !!------------------------------------------------------------------- 897 ! 898 ! == thickness and concentration == ! 850 899 ! distribution over the jpl ice categories 851 ! ---------------------------------------- 852 ! a gaussian distribution for ice concentration is used 853 ! then we check whether the distribution fullfills 854 ! volume and area conservation, positivity and ice categories bounds 855 idim = SIZE( zhti , 1 ) 856 zh_i(1:idim,1:jpl) = 0._wp 857 zh_s(1:idim,1:jpl) = 0._wp 858 za_i(1:idim,1:jpl) = 0._wp 859 ! 860 IF( jpl == 1 ) THEN 861 CALL ice_var_itd_1c1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 862 RETURN 863 ENDIF 900 ! a gaussian distribution for ice concentration is used 901 ! then we check whether the distribution fullfills 902 ! volume and area conservation, positivity and ice categories bounds 903 idim = SIZE( phti , 1 ) 904 ! 905 ph_i(1:idim,1:jpl) = 0._wp 906 ph_s(1:idim,1:jpl) = 0._wp 907 pa_i(1:idim,1:jpl) = 0._wp 864 908 ! 865 909 DO ji = 1, idim 866 910 ! 867 IF( zhti(ji) > 0._wp ) THEN911 IF( phti(ji) > 0._wp ) THEN 868 912 ! 869 913 ! find which category (jl0) the input ice thickness falls into 870 914 jl0 = jpl 871 915 DO jl = 1, jpl 872 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN916 IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN 873 917 jl0 = jl 874 918 CYCLE … … 882 926 i_fill = i_fill - 1 883 927 ! 884 zh_i(ji,1:jpl) = 0._wp885 za_i(ji,1:jpl) = 0._wp928 ph_i(ji,1:jpl) = 0._wp 929 pa_i(ji,1:jpl) = 0._wp 886 930 itest(:) = 0 887 931 ! 888 932 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1 889 zh_i(ji,1) = zhti(ji)890 za_i (ji,1) = zati (ji)933 ph_i(ji,1) = phti(ji) 934 pa_i (ji,1) = pati (ji) 891 935 ELSE !-- case ice is thicker: fill categories >1 892 936 ! thickness 893 937 DO jl = 1, i_fill - 1 894 zh_i(ji,jl) = hi_mean(jl)938 ph_i(ji,jl) = hi_mean(jl) 895 939 END DO 896 940 ! 897 941 ! concentration 898 za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))942 pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl)) 899 943 DO jl = 1, i_fill - 1 900 944 IF ( jl /= jl0 ) THEN 901 zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )902 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)945 zarg = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp ) 946 pa_i(ji,jl) = pa_i (ji,jl0) * EXP(-zarg**2) 903 947 ENDIF 904 948 END DO 905 949 ! 906 950 ! last category 907 za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )908 zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )909 zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )951 pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) ) 952 zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) ) 953 ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 ) 910 954 ! 911 955 ! correction if concentration of upper cat is greater than lower cat … … 913 957 IF ( jl0 /= jpl ) THEN 914 958 DO jl = jpl, jl0+1, -1 915 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN916 zdv = zh_i(ji,jl) * za_i(ji,jl)917 zh_i(ji,jl ) = 0._wp918 za_i (ji,jl ) = 0._wp919 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )959 IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN 960 zdv = ph_i(ji,jl) * pa_i(ji,jl) 961 ph_i(ji,jl ) = 0._wp 962 pa_i (ji,jl ) = 0._wp 963 pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 ) 920 964 END IF 921 965 END DO … … 925 969 ! 926 970 ! Compatibility tests 927 zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )971 zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) ) 928 972 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation 929 973 ! 930 zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )974 zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) ) 931 975 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation 932 976 ! 933 IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?977 IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ? 934 978 ! 935 979 itest(4) = 1 936 980 DO jl = 1, i_fill 937 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations981 IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations 938 982 END DO 939 983 ! !---------------------------- … … 943 987 END DO 944 988 945 ! Add Snow in each category where za_i is not 0989 ! Add Snow in each category where pa_i is not 0 946 990 DO jl = 1, jpl 947 991 DO ji = 1, idim 948 IF( za_i(ji,jl) > 0._wp ) THEN949 zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )992 IF( pa_i(ji,jl) > 0._wp ) THEN 993 ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) ) 950 994 ! In case snow load is in excess that would lead to transformation from snow to ice 951 995 ! Then, transfer the snow excess into the ice (different from icethd_dh) 952 zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )996 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 953 997 ! recompute h_i, h_s avoiding out of bounds values 954 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )955 zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )998 ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) 999 ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos ) 956 1000 ENDIF 957 1001 END DO 958 1002 END DO 959 1003 ! 1004 ! == temperature and salinity == ! 1005 IF( PRESENT( pt_i ) ) THEN 1006 DO jl = 1, jpl 1007 pt_i(:,jl) = ptmi (:) 1008 END DO 1009 ENDIF 1010 IF( PRESENT( pt_s ) ) THEN 1011 DO jl = 1, jpl 1012 pt_s (:,jl) = ptms (:) 1013 END DO 1014 ENDIF 1015 IF( PRESENT( pt_su ) ) THEN 1016 DO jl = 1, jpl 1017 pt_su(:,jl) = ptmsu(:) 1018 END DO 1019 ENDIF 1020 IF( PRESENT( ps_i ) ) THEN 1021 DO jl = 1, jpl 1022 ps_i (:,jl) = psmi (:) 1023 END DO 1024 ENDIF 1025 ! 960 1026 END SUBROUTINE ice_var_itd_1cMc 961 1027 962 SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 1028 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, & 1029 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 963 1030 !!------------------------------------------------------------------- 964 1031 !! … … 981 1048 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 982 1049 !! 983 !! ** Arguments : zhti: N-cat ice thickness984 !! zhts: N-cat snow depth985 !! zati: N-cat ice concentration1050 !! ** Arguments : phti: N-cat ice thickness 1051 !! phts: N-cat snow depth 1052 !! pati: N-cat ice concentration 986 1053 !! 987 1054 !! ** Output : jpl-cat … … 989 1056 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 990 1057 !!------------------------------------------------------------------- 991 INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices 1058 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 1059 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 1060 REAL(wp), DIMENSION(:,:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 1061 REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 1062 ! 1063 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2 1064 INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin 1065 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp 1066 ! 1067 REAL(wp), PARAMETER :: ztrans = 0.25_wp 1068 INTEGER :: ji, jl, jl1, jl2 992 1069 INTEGER :: idim, icat 993 REAL(wp), PARAMETER :: ztrans = 0.25_wp 994 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 995 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 996 INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2 997 INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin 998 !!------------------------------------------------------------------- 999 ! 1000 idim = SIZE( zhti, 1 ) 1001 icat = SIZE( zhti, 2 ) 1070 !!------------------------------------------------------------------- 1071 ! 1072 idim = SIZE( phti, 1 ) 1073 icat = SIZE( phti, 2 ) 1074 ! 1075 ! == thickness and concentration == ! 1002 1076 ! ! ---------------------- ! 1003 1077 IF( icat == jpl ) THEN ! input cat = output cat ! 1004 1078 ! ! ---------------------- ! 1005 zh_i(:,:) = zhti(:,:) 1006 zh_s(:,:) = zhts(:,:) 1007 za_i(:,:) = zati(:,:) 1079 ph_i(:,:) = phti(:,:) 1080 ph_s(:,:) = phts(:,:) 1081 pa_i(:,:) = pati(:,:) 1082 ! 1083 ! == temperature and salinity == ! 1084 IF( PRESENT( pt_i ) ) pt_i (:,:) = ptmi (:,:) 1085 IF( PRESENT( pt_s ) ) pt_s (:,:) = ptms (:,:) 1086 IF( PRESENT( pt_su ) ) pt_su(:,:) = ptmsu(:,:) 1087 IF( PRESENT( ps_i ) ) ps_i (:,:) = psmi (:,:) 1008 1088 ! ! ---------------------- ! 1009 ELSEIF( icat == 1 ) THEN ! specific case if N = 1!1089 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! 1010 1090 ! ! ---------------------- ! 1011 ! 1012 CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i, zh_s, za_i ) 1013 ! 1091 CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), ph_i(:,:), ph_s(:,:), pa_i (:,:), & 1092 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) ) 1014 1093 ! ! ---------------------- ! 1015 ELSEIF( jpl == 1 ) THEN ! specific case if M = 1!1094 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! 1016 1095 ! ! ---------------------- ! 1017 ! 1018 CALL ice_var_itd_Nc1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 1019 ! 1096 CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), ph_i(:,1), ph_s(:,1), pa_i (:,1), & 1097 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) ) 1020 1098 ! ! ----------------------- ! 1021 1099 ELSE ! input cat /= output cat ! … … 1026 1104 1027 1105 ! --- initialize output fields to 0 --- ! 1028 zh_i(1:idim,1:jpl) = 0._wp1029 zh_s(1:idim,1:jpl) = 0._wp1030 za_i(1:idim,1:jpl) = 0._wp1106 ph_i(1:idim,1:jpl) = 0._wp 1107 ph_s(1:idim,1:jpl) = 0._wp 1108 pa_i(1:idim,1:jpl) = 0._wp 1031 1109 ! 1032 1110 ! --- fill the categories --- ! … … 1038 1116 DO jl2 = 1, icat 1039 1117 DO ji = 1, idim 1040 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN1118 IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 1041 1119 ! fill the right category 1042 zh_i(ji,jl1) = zhti(ji,jl2)1043 zh_s(ji,jl1) = zhts(ji,jl2)1044 za_i(ji,jl1) = zati(ji,jl2)1120 ph_i(ji,jl1) = phti(ji,jl2) 1121 ph_s(ji,jl1) = phts(ji,jl2) 1122 pa_i(ji,jl1) = pati(ji,jl2) 1045 1123 ! record categories that are filled 1046 1124 jlmax(ji) = MAX( jlmax(ji), jl1 ) … … 1059 1137 IF( jl1 > 1 ) THEN 1060 1138 ! fill the lower cat (jl1-1) 1061 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)1062 zh_i(ji,jl1-1) = hi_mean(jl1-1)1139 pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 1140 ph_i(ji,jl1-1) = hi_mean(jl1-1) 1063 1141 ! remove from cat jl1 1064 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)1142 pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 1065 1143 ENDIF 1066 1144 IF( jl2 < jpl ) THEN 1067 1145 ! fill the upper cat (jl2+1) 1068 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)1069 zh_i(ji,jl2+1) = hi_mean(jl2+1)1146 pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 1147 ph_i(ji,jl2+1) = hi_mean(jl2+1) 1070 1148 ! remove from cat jl2 1071 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)1149 pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 1072 1150 ENDIF 1073 1151 END DO … … 1079 1157 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 1080 1158 ! fill high 1081 za_i(ji,jl) = ztrans * za_i(ji,jl-1)1082 zh_i(ji,jl) = hi_mean(jl)1159 pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 1160 ph_i(ji,jl) = hi_mean(jl) 1083 1161 jlfil(ji,jl) = jl 1084 1162 ! remove low 1085 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)1163 pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 1086 1164 ENDIF 1087 1165 END DO … … 1093 1171 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 1094 1172 ! fill low 1095 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)1096 zh_i(ji,jl) = hi_mean(jl)1173 pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 1174 ph_i(ji,jl) = hi_mean(jl) 1097 1175 jlfil2(ji,jl) = jl 1098 1176 ! remove high 1099 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)1177 pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 1100 1178 ENDIF 1101 1179 END DO … … 1104 1182 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 1105 1183 DEALLOCATE( jlmin, jlmax ) 1184 ! 1185 ! == temperature and salinity == ! 1186 ! 1187 IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 1188 ! 1189 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 1190 ! 1191 WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 1192 ELSEWHERE ; z1_ai(:) = 0._wp 1193 END WHERE 1194 WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 1195 ELSEWHERE ; z1_vi(:) = 0._wp 1196 END WHERE 1197 WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 1198 ELSEWHERE ; z1_vs(:) = 0._wp 1199 END WHERE 1200 ! 1201 ! fill all the categories with the same value 1202 IF( PRESENT( pt_i ) ) THEN 1203 ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1204 DO jl = 1, jpl 1205 pt_i (:,jl) = ztmp(:) 1206 END DO 1207 ENDIF 1208 IF( PRESENT( pt_s ) ) THEN 1209 ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 1210 DO jl = 1, jpl 1211 pt_s (:,jl) = ztmp(:) 1212 END DO 1213 ENDIF 1214 IF( PRESENT( pt_su ) ) THEN 1215 ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 1216 DO jl = 1, jpl 1217 pt_su(:,jl) = ztmp(:) 1218 END DO 1219 ENDIF 1220 IF( PRESENT( ps_i ) ) THEN 1221 ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1222 DO jl = 1, jpl 1223 ps_i (:,jl) = ztmp(:) 1224 END DO 1225 ENDIF 1226 ! 1227 DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 1228 ! 1229 ENDIF 1106 1230 ! 1107 1231 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
