Changeset 11332
- Timestamp:
- 2019-07-23T16:26:43+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90
r11223 r11332 773 773 !! ** Purpose : converting N-cat ice to jpl ice categories 774 774 !!------------------------------------------------------------------- 775 SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 775 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, & 776 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 776 777 !!------------------------------------------------------------------- 777 778 !! ** Purpose : converting 1-cat ice to 1 ice category 778 779 !!------------------------------------------------------------------- 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(:) 780 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 781 REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 782 REAL(wp), DIMENSION(:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 783 REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 784 !!------------------------------------------------------------------- 785 ! == thickness and concentration == ! 786 ph_i(:) = phti(:) 787 ph_s(:) = phts(:) 788 pa_i(:) = pati(:) 789 ! 790 ! == temperature and salinity == ! 791 IF( PRESENT( pt_i ) ) pt_i (:) = ptmi (:) 792 IF( PRESENT( pt_s ) ) pt_s (:) = ptms (:) 793 IF( PRESENT( pt_su ) ) pt_su(:) = ptmsu(:) 794 IF( PRESENT( ps_i ) ) ps_i (:) = psmi (:) 795 785 796 END SUBROUTINE ice_var_itd_1c1c 786 797 787 SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i ) 798 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, & 799 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 788 800 !!------------------------------------------------------------------- 789 801 !! ** Purpose : converting N-cat ice to 1 ice category 790 802 !!------------------------------------------------------------------- 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 803 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 804 REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 805 REAL(wp), DIMENSION(:,:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 806 REAL(wp), DIMENSION(:) , INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 807 ! 808 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs 809 ! 810 INTEGER :: idim 811 !!------------------------------------------------------------------- 812 ! 813 idim = SIZE( phti, 1 ) 814 ! 815 ! == thickness and concentration == ! 816 ALLOCATE( z1_ai(idim) ) 817 ! 818 pa_i(:) = SUM( pati(:,:), dim=2 ) 819 820 WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:) 821 ELSEWHERE ; z1_ai(:) = 0._wp 803 822 END WHERE 823 824 ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 825 ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 826 ! 827 ! == temperature and salinity == ! 828 IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 829 ! 830 ALLOCATE( z1_vi(idim), z1_vs(idim) ) 831 ! 832 WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 833 ELSEWHERE ; z1_vi(:) = 0._wp 834 END WHERE 835 WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 836 ELSEWHERE ; z1_vs(:) = 0._wp 837 END WHERE 838 ! 839 IF( PRESENT( pt_i ) ) pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 840 IF( PRESENT( pt_s ) ) pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 841 IF( PRESENT( pt_su ) ) pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 842 IF( PRESENT( ps_i ) ) ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 843 ! 844 DEALLOCATE( z1_vi, z1_vs ) 845 ! 846 ENDIF 847 ! 848 DEALLOCATE( z1_ai ) 804 849 ! 805 850 END SUBROUTINE ice_var_itd_Nc1c 806 851 807 SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 852 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 853 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 808 854 !!------------------------------------------------------------------- 809 855 !! … … 826 872 !! 4) Iterate until ok (SUM(itest(:) = 4) 827 873 !! 828 !! ** Arguments : zhti: 1-cat ice thickness829 !! zhts: 1-cat snow depth830 !! zati: 1-cat ice concentration874 !! ** Arguments : phti: 1-cat ice thickness 875 !! phts: 1-cat snow depth 876 !! pati: 1-cat ice concentration 831 877 !! 832 878 !! ** Output : jpl-cat … … 834 880 !! (Example of application: BDY forcings when input are cell averaged) 835 881 !!------------------------------------------------------------------- 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 ! ---------------------------------------- 882 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 883 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 884 REAL(wp), DIMENSION(:) , INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 885 REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 886 ! 887 INTEGER , DIMENSION(4) :: itest 888 INTEGER :: ji, jk, jl 889 INTEGER :: idim, i_fill, jl0 890 REAL(wp) :: zarg, zV, zconv, zdh, zdv 891 !!------------------------------------------------------------------- 892 ! 893 ! == thickness and concentration == ! 845 894 ! 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 ! 855 IF( jpl == 1 ) THEN 856 CALL ice_var_itd_1c1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 857 RETURN 858 ENDIF 895 ! a gaussian distribution for ice concentration is used 896 ! then we check whether the distribution fullfills 897 ! volume and area conservation, positivity and ice categories bounds 898 idim = SIZE( phti , 1 ) 899 ! 900 ph_i(1:idim,1:jpl) = 0._wp 901 ph_s(1:idim,1:jpl) = 0._wp 902 pa_i(1:idim,1:jpl) = 0._wp 859 903 ! 860 904 DO ji = 1, idim 861 905 ! 862 IF( zhti(ji) > 0._wp ) THEN906 IF( phti(ji) > 0._wp ) THEN 863 907 ! 864 908 ! find which category (jl0) the input ice thickness falls into 865 909 jl0 = jpl 866 910 DO jl = 1, jpl 867 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN911 IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN 868 912 jl0 = jl 869 913 CYCLE … … 877 921 i_fill = i_fill - 1 878 922 ! 879 zh_i(ji,1:jpl) = 0._wp880 za_i(ji,1:jpl) = 0._wp923 ph_i(ji,1:jpl) = 0._wp 924 pa_i(ji,1:jpl) = 0._wp 881 925 itest(:) = 0 882 926 ! 883 927 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1 884 zh_i(ji,1) = zhti(ji)885 za_i (ji,1) = zati (ji)928 ph_i(ji,1) = phti(ji) 929 pa_i (ji,1) = pati (ji) 886 930 ELSE !-- case ice is thicker: fill categories >1 887 931 ! thickness 888 932 DO jl = 1, i_fill - 1 889 zh_i(ji,jl) = hi_mean(jl)933 ph_i(ji,jl) = hi_mean(jl) 890 934 END DO 891 935 ! 892 936 ! concentration 893 za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))937 pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl)) 894 938 DO jl = 1, i_fill - 1 895 939 IF ( jl /= jl0 ) THEN 896 zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )897 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)940 zarg = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp ) 941 pa_i(ji,jl) = pa_i (ji,jl0) * EXP(-zarg**2) 898 942 ENDIF 899 943 END DO 900 944 ! 901 945 ! last category 902 za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )903 zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )904 zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )946 pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) ) 947 zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) ) 948 ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 ) 905 949 ! 906 950 ! correction if concentration of upper cat is greater than lower cat … … 908 952 IF ( jl0 /= jpl ) THEN 909 953 DO jl = jpl, jl0+1, -1 910 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN911 zdv = zh_i(ji,jl) * za_i(ji,jl)912 zh_i(ji,jl ) = 0._wp913 za_i (ji,jl ) = 0._wp914 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )954 IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN 955 zdv = ph_i(ji,jl) * pa_i(ji,jl) 956 ph_i(ji,jl ) = 0._wp 957 pa_i (ji,jl ) = 0._wp 958 pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 ) 915 959 END IF 916 960 END DO … … 920 964 ! 921 965 ! Compatibility tests 922 zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )966 zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) ) 923 967 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation 924 968 ! 925 zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )969 zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) ) 926 970 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation 927 971 ! 928 IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?972 IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ? 929 973 ! 930 974 itest(4) = 1 931 975 DO jl = 1, i_fill 932 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations976 IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations 933 977 END DO 934 978 ! !---------------------------- … … 938 982 END DO 939 983 940 ! Add Snow in each category where za_i is not 0984 ! Add Snow in each category where pa_i is not 0 941 985 DO jl = 1, jpl 942 986 DO ji = 1, idim 943 IF( za_i(ji,jl) > 0._wp ) THEN944 zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )987 IF( pa_i(ji,jl) > 0._wp ) THEN 988 ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) ) 945 989 ! In case snow load is in excess that would lead to transformation from snow to ice 946 990 ! Then, transfer the snow excess into the ice (different from icethd_dh) 947 zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )991 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 948 992 ! recompute h_i, h_s avoiding out of bounds values 949 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )950 zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )993 ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) 994 ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos ) 951 995 ENDIF 952 996 END DO 953 997 END DO 954 998 ! 999 ! == temperature and salinity == ! 1000 IF( PRESENT( pt_i ) ) THEN 1001 DO jl = 1, jpl 1002 pt_i(:,jl) = ptmi (:) 1003 END DO 1004 ENDIF 1005 IF( PRESENT( pt_s ) ) THEN 1006 DO jl = 1, jpl 1007 pt_s (:,jl) = ptms (:) 1008 END DO 1009 ENDIF 1010 IF( PRESENT( pt_su ) ) THEN 1011 DO jl = 1, jpl 1012 pt_su(:,jl) = ptmsu(:) 1013 END DO 1014 ENDIF 1015 IF( PRESENT( ps_i ) ) THEN 1016 DO jl = 1, jpl 1017 ps_i (:,jl) = psmi (:) 1018 END DO 1019 ENDIF 1020 ! 955 1021 END SUBROUTINE ice_var_itd_1cMc 956 1022 957 SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i ) 1023 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, & 1024 & ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i ) 958 1025 !!------------------------------------------------------------------- 959 1026 !! … … 976 1043 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 977 1044 !! 978 !! ** Arguments : zhti: N-cat ice thickness979 !! zhts: N-cat snow depth980 !! zati: N-cat ice concentration1045 !! ** Arguments : phti: N-cat ice thickness 1046 !! phts: N-cat snow depth 1047 !! pati: N-cat ice concentration 981 1048 !! 982 1049 !! ** Output : jpl-cat … … 984 1051 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 985 1052 !!------------------------------------------------------------------- 986 INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices 1053 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 1054 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 1055 REAL(wp), DIMENSION(:,:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal 1056 REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal 1057 ! 1058 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2 1059 INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin 1060 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp 1061 ! 1062 REAL(wp), PARAMETER :: ztrans = 0.25_wp 1063 INTEGER :: ji, jl, jl1, jl2 987 1064 INTEGER :: idim, icat 988 REAL(wp), PARAMETER :: ztrans = 0.25_wp 989 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables 990 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables 991 INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2 992 INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin 993 !!------------------------------------------------------------------- 994 ! 995 idim = SIZE( zhti, 1 ) 996 icat = SIZE( zhti, 2 ) 1065 !!------------------------------------------------------------------- 1066 ! 1067 idim = SIZE( phti, 1 ) 1068 icat = SIZE( phti, 2 ) 1069 ! 1070 ! == thickness and concentration == ! 997 1071 ! ! ---------------------- ! 998 1072 IF( icat == jpl ) THEN ! input cat = output cat ! 999 1073 ! ! ---------------------- ! 1000 zh_i(:,:) = zhti(:,:) 1001 zh_s(:,:) = zhts(:,:) 1002 za_i(:,:) = zati(:,:) 1074 ph_i(:,:) = phti(:,:) 1075 ph_s(:,:) = phts(:,:) 1076 pa_i(:,:) = pati(:,:) 1077 ! 1078 ! == temperature and salinity == ! 1079 IF( PRESENT( pt_i ) ) pt_i (:,:) = ptmi (:,:) 1080 IF( PRESENT( pt_s ) ) pt_s (:,:) = ptms (:,:) 1081 IF( PRESENT( pt_su ) ) pt_su(:,:) = ptmsu(:,:) 1082 IF( PRESENT( ps_i ) ) ps_i (:,:) = psmi (:,:) 1003 1083 ! ! ---------------------- ! 1004 ELSEIF( icat == 1 ) THEN ! specific case if N = 1!1084 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! 1005 1085 ! ! ---------------------- ! 1006 ! 1007 CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i, zh_s, za_i ) 1008 ! 1086 CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), ph_i(:,:), ph_s(:,:), pa_i (:,:), & 1087 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) ) 1009 1088 ! ! ---------------------- ! 1010 ELSEIF( jpl == 1 ) THEN ! specific case if M = 1!1089 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! 1011 1090 ! ! ---------------------- ! 1012 ! 1013 CALL ice_var_itd_Nc1c( zhti, zhts, zati, zh_i(:,1), zh_s(:,1), za_i(:,1) ) 1014 ! 1091 CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), ph_i(:,1), ph_s(:,1), pa_i (:,1), & 1092 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) ) 1015 1093 ! ! ----------------------- ! 1016 1094 ELSE ! input cat /= output cat ! … … 1021 1099 1022 1100 ! --- initialize output fields to 0 --- ! 1023 zh_i(1:idim,1:jpl) = 0._wp1024 zh_s(1:idim,1:jpl) = 0._wp1025 za_i(1:idim,1:jpl) = 0._wp1101 ph_i(1:idim,1:jpl) = 0._wp 1102 ph_s(1:idim,1:jpl) = 0._wp 1103 pa_i(1:idim,1:jpl) = 0._wp 1026 1104 ! 1027 1105 ! --- fill the categories --- ! … … 1033 1111 DO jl2 = 1, icat 1034 1112 DO ji = 1, idim 1035 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN1113 IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 1036 1114 ! fill the right category 1037 zh_i(ji,jl1) = zhti(ji,jl2)1038 zh_s(ji,jl1) = zhts(ji,jl2)1039 za_i(ji,jl1) = zati(ji,jl2)1115 ph_i(ji,jl1) = phti(ji,jl2) 1116 ph_s(ji,jl1) = phts(ji,jl2) 1117 pa_i(ji,jl1) = pati(ji,jl2) 1040 1118 ! record categories that are filled 1041 1119 jlmax(ji) = MAX( jlmax(ji), jl1 ) … … 1054 1132 IF( jl1 > 1 ) THEN 1055 1133 ! fill the lower cat (jl1-1) 1056 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)1057 zh_i(ji,jl1-1) = hi_mean(jl1-1)1134 pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 1135 ph_i(ji,jl1-1) = hi_mean(jl1-1) 1058 1136 ! remove from cat jl1 1059 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)1137 pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 1060 1138 ENDIF 1061 1139 IF( jl2 < jpl ) THEN 1062 1140 ! fill the upper cat (jl2+1) 1063 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)1064 zh_i(ji,jl2+1) = hi_mean(jl2+1)1141 pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 1142 ph_i(ji,jl2+1) = hi_mean(jl2+1) 1065 1143 ! remove from cat jl2 1066 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)1144 pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 1067 1145 ENDIF 1068 1146 END DO … … 1074 1152 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 1075 1153 ! fill high 1076 za_i(ji,jl) = ztrans * za_i(ji,jl-1)1077 zh_i(ji,jl) = hi_mean(jl)1154 pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 1155 ph_i(ji,jl) = hi_mean(jl) 1078 1156 jlfil(ji,jl) = jl 1079 1157 ! remove low 1080 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) 1081 1159 ENDIF 1082 1160 END DO … … 1088 1166 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 1089 1167 ! fill low 1090 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)1091 zh_i(ji,jl) = hi_mean(jl)1168 pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 1169 ph_i(ji,jl) = hi_mean(jl) 1092 1170 jlfil2(ji,jl) = jl 1093 1171 ! remove high 1094 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)1172 pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 1095 1173 ENDIF 1096 1174 END DO … … 1099 1177 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 1100 1178 DEALLOCATE( jlmin, jlmax ) 1179 ! 1180 ! == temperature and salinity == ! 1181 ! 1182 IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN 1183 ! 1184 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 1185 ! 1186 WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 1187 ELSEWHERE ; z1_ai(:) = 0._wp 1188 END WHERE 1189 WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 1190 ELSEWHERE ; z1_vi(:) = 0._wp 1191 END WHERE 1192 WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 1193 ELSEWHERE ; z1_vs(:) = 0._wp 1194 END WHERE 1195 ! 1196 ! fill all the categories with the same value 1197 IF( PRESENT( pt_i ) ) THEN 1198 ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1199 DO jl = 1, jpl 1200 pt_i (:,jl) = ztmp(:) 1201 END DO 1202 ENDIF 1203 IF( PRESENT( pt_s ) ) THEN 1204 ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 1205 DO jl = 1, jpl 1206 pt_s (:,jl) = ztmp(:) 1207 END DO 1208 ENDIF 1209 IF( PRESENT( pt_su ) ) THEN 1210 ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 1211 DO jl = 1, jpl 1212 pt_su(:,jl) = ztmp(:) 1213 END DO 1214 ENDIF 1215 IF( PRESENT( ps_i ) ) THEN 1216 ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1217 DO jl = 1, jpl 1218 ps_i (:,jl) = ztmp(:) 1219 END DO 1220 ENDIF 1221 ! 1222 DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 1223 ! 1224 ENDIF 1101 1225 ! 1102 1226 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
