 Timestamp:
 20201014T17:29:11+02:00 (4 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/branches/2020/SI303_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(i1,j)+B*u(i,j)+C*u(i+1,j) = F 810 807 811 808 DO jn = 1, nn_zebra_vp ! "zebra" loop (! redblacksor!!! ) … … 819 816 IF ( lwp ) WRITE(numout,*) ' Into the Uzebra 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(i1,j)+B*u(i,j)+C*u(i+1,j) = F 831 832 !  Righthand 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(ji1,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 , jpi2 ! all other rows 914 !DO ji = jpi2 , 2, 1 ! all other rows ! > original backward loop888 !DO ji = jpi2 , 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(jpi1,jpj) = zFU_prime(jpi1,jpj) 898 ELSE ! nn_thomas == 2 ! mitGCM version 899 900 u_ice(jpi1,jj) = zFU_prime(jpi1,jj) 931 901 932 902 DO ji2 = 2 , jpi2 ! 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 littleice and noice 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 littleice and noice 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,j1)+B*v(i,j)+C*v(i,j+1) = F 934 ! It is intentional to have a ji then jj loop for Vvelocity 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 Vvelocity 981 !!! ZH97 explain it is critical for convergence speed 982 983 ! 984 !  Tridiagonal system solver for each row 985 ! 986 ! A*v(i,j1)+B*v(i,j)+C*v(i,j+1) = F 987 988 !  Righthand 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,jj1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj1) )  epsi20 ) )1010 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj1) ) , epsi20 )1011 zBV_prime(ji,jj) = zBV(ji,jj)  zw * zCV(ji,jj1)1012 zFV_prime(ji,jj) = zFV(ji,jj)  zw * zFV(ji,jj1)1013 END DO1014 1015 ! Backward sweep1016 zfac = SIGN( 1._wp , zBV_prime(ji,jpj1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj1) )  epsi20 ) )1017 v_ice(ji,jpj1) = zfac * zFV_prime(ji,jpj1) / MAX ( ABS(zBV_prime(ji,jpj1) ) , epsi20 ) &1018 & * vmask(ji,jpj1,1) ! last row1019 1020 DO jj = jpj2, 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,jj1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj1) )  epsi20 ) ) 974 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj1) ) , epsi20 ) 975 zBV_prime(ji,jj) = zBV(ji,jj)  zw * zCV(ji,jj1) 976 zFV_prime(ji,jj) = zFV(ji,jj)  zw * zFV(ji,jj1) 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(ji1,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(ji1,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,jpj1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj1) )  epsi20 ) ) 1009 v_ice(ji,jpj1) = zfac * zFV_prime(ji,jpj1) / MAX ( ABS(zBV_prime(ji,jpj1) ) , epsi20 ) & 1010 & * vmask(ji,jpj1,1) ! last row 1011 1012 ! other rows 1013 !DO jj = jpj2, 2, 1 ! original back loop 1014 DO jj2 = 2 , jpj2 ! 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,jpj1) = zFV_prime(ji,jpj) 1025 1026 ! other rows 1027 DO jj2 = 2 , jpi2 ! 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 noice and littleice 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.