Changeset 13593
- Timestamp:
- 2020-10-14T17:29:11+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r13591 r13593 131 131 LOGICAL :: ll_u_iterate, ll_v_iterate ! continue iteration or not 132 132 ! 133 INTEGER :: ji, ji2, jj, j n ! dummy loop indices133 INTEGER :: ji, ji2, jj, jj2, jn ! dummy loop indices 134 134 INTEGER :: jter, i_out, i_inn ! 135 135 INTEGER :: ji_min, jj_min ! … … 777 777 ll_u_iterate = .TRUE. 778 778 ll_v_iterate = .TRUE. 779 780 ! zBU 781 ! it is probably the way zBU is used that is problematic... 782 ! because zBU itself does not have the zipper!!! 783 ! BUG only occurs when zBU is different from zero 784 ! H1 - the jj loop does not work well with jpj = 149 ? 785 ! zBU(:,:) = 0._wp; zBV(:,:) = 0._wp 779 780 nn_thomas = 2 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm 786 781 787 782 DO i_inn = 1, nn_ninn_vp ! inner loop iterations … … 806 801 ! ---------------------------- ! 807 802 803 ! What follows could be subroutinized... 804 808 805 ! Thomas Algorithm for tridiagonal solver 809 ! what follows could be subroutinized...806 ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 810 807 811 808 DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) … … 819 816 IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min 820 817 821 !----------------------------------------------------------822 ! -- Boundary condition (link with neighbouring processor)823 !----------------------------------------------------------824 825 818 DO jj = jj_min, jpj - 1, nn_zebra_vp 826 827 ! 828 ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan 829 ! 830 ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 831 832 ! - Right-hand side of tridiagonal system (zFU) 819 820 !------------------------ 821 ! Independent term (zFU) 822 !------------------------ 833 823 DO ji = 2, jpi - 1 834 824 835 825 ! boundary condition substitution 836 826 ! see Zhang and Hibler, 1997, Appendix B 837 ! MV NOTE possibly not fully appropriate838 827 zAA3 = 0._wp 839 828 IF ( ji == 2 ) zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) … … 848 837 END DO 849 838 850 END DO 851 852 ! CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. ) 853 854 !--------------- 855 ! Forward sweep 856 !--------------- 857 nn_thomas = 2 ! 1 = martin, 2 = mitgcm 858 859 IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 860 861 DO jj = jj_min, jpj - 1, nn_zebra_vp 839 !--------------- 840 ! Forward sweep 841 !--------------- 842 843 IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 862 844 863 845 DO ji = 3, jpi - 1 … … 870 852 END DO 871 853 872 END DO 873 874 ELSE ! nn_thomas == 2 ! mitGCM algorithm 875 876 DO jj = jj_min, jpj - 1, nn_zebra_vp 854 ELSE ! nn_thomas == 2 ! mitGCM algorithm 877 855 878 856 ! ji = 2 … … 893 871 END DO 894 872 895 END DO 896 897 ENDIF 898 899 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zFU_prime, 'U', 1. ) 900 901 !---------------- 902 ! Backward sweep 903 !---------------- 904 IF ( nn_thomas == 1 ) THEN ! MV version 905 906 DO jj = jj_min, jpj - 1, nn_zebra_vp 873 ENDIF 874 875 !----------------------------- 876 ! Backward sweep & relaxation 877 !----------------------------- 878 879 ! --- Backward sweep 880 IF ( nn_thomas == 1 ) THEN ! MV version 907 881 908 882 ! last row … … 912 886 913 887 DO ji2 = 2 , jpi-2 ! all other rows 914 !DO ji = jpi-2 , 2, -1 ! all other rows ! ---> original backward loop888 !DO ji = jpi-2 , 2, -1 ! all other rows ! ---> original backward loop 915 889 ji = jpi - ji2 916 890 ! ji2 = 2 -> ji = jpi - 2 … … 922 896 END DO 923 897 924 END DO 925 926 ELSE ! nn_thomas == 2 ! mitGCM version 927 928 DO jj = jj_min, jpj - 1, nn_zebra_vp 929 930 u_ice(jpi-1,jpj) = zFU_prime(jpi-1,jpj) 898 ELSE ! nn_thomas == 2 ! mitGCM version 899 900 u_ice(jpi-1,jj) = zFU_prime(jpi-1,jj) 931 901 932 902 DO ji2 = 2 , jpi-2 ! all other rows ! MV OPT: could be turned into forward loop (by substituting ji) … … 935 905 END DO 936 906 937 END DO 938 939 ENDIF 940 941 !--------------- 942 ! -- Relaxation 943 !--------------- 944 945 DO jj = jj_min, jpj - 1, nn_zebra_vp 946 907 ENDIF 908 909 !--- Relaxation 910 ! and velocity masking for little-ice and no-ice cases 947 911 DO ji = 2, jpi - 1 948 912 949 u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) 913 u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) ! relaxation 950 914 951 ! velocity masking for little-ice and no-ice cases 952 u_ice(ji,jj) = zmsk00x(ji,jj) & 915 u_ice(ji,jj) = zmsk00x(ji,jj) & ! masking 953 916 & * ( zmsk01x(ji,jj) * u_ice(ji,jj) & 954 917 & + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) * umask(ji,jj,1) … … 966 929 ! ! ---------------------------- ! 967 930 968 ! MV OPT: what follows could be subroutinized... 969 931 ! MV OPT: what follows could be subroutinized... 932 ! Thomas Algorithm for tridiagonal solver 933 ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F 934 ! It is intentional to have a ji then jj loop for V-velocity 935 !!! ZH97 explain it is critical for convergence speed 936 970 937 DO jn = 1, nn_zebra_vp ! "zebra" loop 971 938 … … 977 944 978 945 DO ji = ji_min, jpi - 1, nn_zebra_vp 979 980 !!! It is intentional to have a ji then jj loop for V-velocity 981 !!! ZH97 explain it is critical for convergence speed 982 983 !------------------------------------------- 984 ! -- Tridiagonal system solver for each row 985 !------------------------------------------- 986 ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F 987 988 ! --- Right-hand side of tridiagonal system (zFU) 946 947 !------------------------ 948 ! Independent term (zFU) 949 !------------------------ 989 950 DO jj = 2, jpj - 1 990 951 … … 1002 963 1003 964 END DO 1004 1005 ! --- Thomas Algorithm 1006 ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 965 966 !--------------- 1007 967 ! Forward sweep 1008 DO jj = 3, jpj - 11009 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) )1010 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 )1011 zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1)1012 zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1)1013 END DO1014 1015 ! Backward sweep1016 zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) )1017 v_ice(ji,jpj-1) = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) &1018 & * vmask(ji,jpj-1,1) ! last row1019 1020 DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed1021 zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) )1022 v_ice(ji,jj) = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) ) &1023 & / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )1024 END DO1025 1026 !---------------1027 ! -- Relaxation1028 968 !--------------- 969 970 IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 971 972 DO jj = 3, jpj - 1 973 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 974 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 ) 975 zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 976 zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 977 END DO 978 979 ELSE ! nn_thomas == 2 ! mitGCM algorithm 980 981 ! jj = 2 982 zw = zBV(2,jj) 983 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 984 zw = zfac * 1.0_wp / MAX ( ABS( zw ) , epsi20 ) 985 zCV_prime(2,jj) = zw * zCV(2,jj) 986 zFV_prime(2,jj) = zw * zFV(2,jj) 987 988 DO jj = 3, jpj - 1 989 990 zw = zBV(ji,jj) - zAV(ji,jj) * zCV_prime(ji-1,jj) 991 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 992 zw = zfac * 1.0_wp / MAX ( ABS( zw ) , epsi20 ) 993 zCV_prime(ji,jj) = zw * zCV(ji,jj) 994 zFV_prime(ji,jj) = zw * ( zFV(ji,jj) - zAV(ji,jj) * zFV_prime(ji-1,jj) ) 995 996 END DO 997 998 ENDIF 999 1000 !----------------------------- 1001 ! Backward sweep & relaxation 1002 !----------------------------- 1003 1004 ! --- Backward sweep 1005 IF ( nn_thomas == 1 ) THEN ! nn_thomas = 1 ! --- MV version seemingly more than mitGCM algorithm 1006 1007 ! last row 1008 zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) ) 1009 v_ice(ji,jpj-1) = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) & 1010 & * vmask(ji,jpj-1,1) ! last row 1011 1012 ! other rows 1013 !DO jj = jpj-2, 2, -1 ! original back loop 1014 DO jj2 = 2 , jpj-2 ! modified forward loop 1015 jj = jpj - jj2 ! index swap 1016 zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) ) 1017 v_ice(ji,jj) = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) ) & 1018 & / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 ) 1019 END DO 1020 1021 ELSE ! nn_thomas == 2 ! --- mitGCM algorithm 1022 1023 ! last row 1024 v_ice(ji,jpj-1) = zFV_prime(ji,jpj) 1025 1026 ! other rows 1027 DO jj2 = 2 , jpi-2 ! all other rows 1028 jj = jpj - jj2 1029 v_ice(ji,jj) = zFV_prime(ji,jj) - zCV_prime(ji,jj) * v_ice(ji,jj+1) 1030 END DO 1031 1032 ENDIF ! nn_thomas 1033 1034 ! --- Relaxation & masking (should it be now or later) 1029 1035 DO jj = 2, jpj - 1 1030 v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) )1031 1032 ! mask velocity for no-ice and little-ice cases1033 v_ice(ji,jj) = zmsk00y(ji,jj) &1034 & * ( zmsk01y(ji,jj) * v_ice(ji,jj) &1036 1037 v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) ) ! relaxation 1038 1039 v_ice(ji,jj) = zmsk00y(ji,jj) & ! masking 1040 & * ( zmsk01y(ji,jj) * v_ice(ji,jj) & 1035 1041 & + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp ) * vmask(ji,jj,1) 1036 1042 1037 END DO 1038 1043 END DO ! jj 1044 1039 1045 END DO ! ji 1040 1046 1041 1047 END DO ! zebra loop 1042 1048
Note: See TracChangeset
for help on using the changeset viewer.