- Timestamp:
- 2020-10-19T08:46:02+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
r13593 r13629 777 777 ll_u_iterate = .TRUE. 778 778 ll_v_iterate = .TRUE. 779 780 nn_thomas = 2 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm 779 780 ! nn_thomas = 2 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm ! ---> not working for now 781 nn_thomas = 1 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm ! ---> BUG MPP!!! 781 782 782 783 DO i_inn = 1, nn_ninn_vp ! inner loop iterations … … 789 790 ! l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 790 791 791 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate792 zv_b(:,:) = v_ice(:,:)793 794 zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp795 zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp792 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 793 zv_b(:,:) = v_ice(:,:) 794 795 ! zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 796 ! zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp 796 797 797 798 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN … … 800 801 IF ( ll_u_iterate ) THEN ! --- Solve for u-velocity --- ! 801 802 ! ---------------------------- ! 802 803 803 804 ! What follows could be subroutinized... 804 805 805 806 ! Thomas Algorithm for tridiagonal solver 806 807 ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 808 809 zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; zCU_prime(:,:) = 0._wp 807 810 808 811 DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) … … 826 829 ! see Zhang and Hibler, 1997, Appendix B 827 830 zAA3 = 0._wp 828 IF ( ji == 2 ) 829 IF ( ji == jpi - 1 ) 831 IF ( ji == 2 ) zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 832 IF ( ji == jpi - 1 ) zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj) 830 833 831 834 ! right hand side … … 854 857 ELSE ! nn_thomas == 2 ! mitGCM algorithm 855 858 856 ! ji = 2859 ! --- ji = 2 857 860 zw = zBU(2,jj) 858 861 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 859 zw = zfac * 1.0_wp/ MAX ( ABS( zw ) , epsi20 )862 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 860 863 zCU_prime(2,jj) = zw * zCU(2,jj) 861 864 zFU_prime(2,jj) = zw * zFU(2,jj) 862 865 866 ! --- ji = 3 --> jpi-1 863 867 DO ji = 3, jpi - 1 864 868 865 869 zw = zBU(ji,jj) - zAU(ji,jj) * zCU_prime(ji-1,jj) 866 870 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 867 zw = zfac * 1.0_wp/ MAX ( ABS( zw ) , epsi20 )871 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 868 872 zCU_prime(ji,jj) = zw * zCU(ji,jj) 869 873 zFU_prime(ji,jj) = zw * ( zFU(ji,jj) - zAU(ji,jj) * zFU_prime(ji-1,jj) ) … … 882 886 ! last row 883 887 zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 884 u_ice(jpi-1,jj) 888 u_ice(jpi-1,jj) = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 885 889 & * umask(jpi-1,jj,1) 886 890 … … 891 895 ! ji2 = jpi - 2 -> ji = 2 892 896 zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 893 u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) * umask(ji,jj,1) )& ! this897 u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) & ! this 894 898 ! line is guilty 895 899 & / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) … … 898 902 ELSE ! nn_thomas == 2 ! mitGCM version 899 903 900 u_ice(jpi-1,jj) = zFU_prime(jpi-1,jj)904 u_ice(jpi-1,jj) = zFU_prime(jpi-1,jj) * umask(jpi-1,jj,1) 901 905 902 906 DO ji2 = 2 , jpi-2 ! all other rows ! MV OPT: could be turned into forward loop (by substituting ji) 903 907 ji = jpi - ji2 904 u_ice(ji,jj) = zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj) 908 ! ji2 = 2, ji = jpi - 2 909 ! ji2 = jpi- 2, ji = 2 910 u_ice(ji,jj) = ( zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) 905 911 END DO 906 912 … … 915 921 u_ice(ji,jj) = zmsk00x(ji,jj) & ! masking 916 922 & * ( zmsk01x(ji,jj) * u_ice(ji,jj) & 917 & + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) * umask(ji,jj,1)923 & + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) * umask(ji,jj,1) 918 924 919 925 END DO … … 935 941 !!! ZH97 explain it is critical for convergence speed 936 942 943 zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp 944 937 945 DO jn = 1, nn_zebra_vp ! "zebra" loop 938 946 … … 946 954 947 955 !------------------------ 948 ! Independent term (zF U)956 ! Independent term (zFV) 949 957 !------------------------ 950 958 DO jj = 2, jpj - 1 … … 970 978 IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 971 979 972 973 974 975 976 977 978 980 DO jj = 3, jpj - 1 981 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 982 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 ) 983 zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 984 zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 985 END DO 986 979 987 ELSE ! nn_thomas == 2 ! mitGCM algorithm 980 988 981 989 ! jj = 2 982 zw = zBV( 2,jj)990 zw = zBV(ji,2) 983 991 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)992 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 993 zCV_prime(ji,2) = zw * zCV(ji,2) 994 zFV_prime(ji,2) = zw * zFV(ji,2) 987 995 988 996 DO jj = 3, jpj - 1 989 990 zw = zBV(ji,jj) - zAV(ji,jj) * zCV_prime(ji-1,jj) 997 zw = zBV(ji,jj) - zAV(ji,jj) * zCV_prime(ji,jj-1) 991 998 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 992 zw = zfac * 1.0_wp/ MAX ( ABS( zw ) , epsi20 )999 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 993 1000 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 1001 zFV_prime(ji,jj) = zw * ( zFV(ji,jj) - zAV(ji,jj) * zFV_prime(ji,jj-1) ) 996 1002 END DO 997 1003 998 1004 ENDIF 999 1005 … … 1015 1021 jj = jpj - jj2 ! index swap 1016 1022 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)) &1023 v_ice(ji,jj) = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) & 1018 1024 & / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 ) 1019 1025 END DO … … 1022 1028 1023 1029 ! last row 1024 v_ice(ji,jpj-1) = zFV_prime(ji,jpj )1030 v_ice(ji,jpj-1) = zFV_prime(ji,jpj-1) * vmask(ji,jpj-1,1) 1025 1031 1026 1032 ! other rows 1027 DO jj2 = 2 , jp i-2 ! all other rows1033 DO jj2 = 2 , jpj-2 ! all other rows 1028 1034 jj = jpj - jj2 1029 v_ice(ji,jj) = zFV_prime(ji,jj) - zCV_prime(ji,jj) * v_ice(ji,jj+1)1035 v_ice(ji,jj) = ( zFV_prime(ji,jj) - zCV_prime(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) 1030 1036 END DO 1031 1037 1032 1038 ENDIF ! nn_thomas 1033 1039 … … 1138 1144 IF ( lwp ) WRITE(numout,*) ' We are out of outer loop ' 1139 1145 1146 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU , 'U', 1., zFV , 'V', 1. ) 1147 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zBV_prime , 'V', 1. ) 1148 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime , 'U', 1., zFV_prime , 'V', 1. ) 1149 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU_prime , 'U', 1., zCV_prime , 'V', 1. ) 1150 1140 1151 CALL iom_put( 'zFU' , zFU ) ! MV DEBUG 1141 1152 CALL iom_put( 'zBU_prime' , zBU_prime ) ! MV DEBUG 1153 CALL iom_put( 'zCU_prime' , zCU_prime ) ! MV DEBUG 1142 1154 CALL iom_put( 'zFU_prime' , zFU_prime ) ! MV DEBUG 1143 1155 1144 1156 CALL iom_put( 'zFV' , zFV ) ! MV DEBUG 1145 1157 CALL iom_put( 'zBV_prime' , zBV_prime ) ! MV DEBUG 1158 CALL iom_put( 'zCV_prime' , zCV_prime ) ! MV DEBUG 1146 1159 CALL iom_put( 'zFV_prime' , zFV_prime ) ! MV DEBUG 1147 1160 … … 1183 1196 1184 1197 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1198 1199 IF ( lwp ) WRITE(numout,*) ' Velocity replaced ' 1185 1200 1186 1201 ! END DEBUG … … 1238 1253 END DO 1239 1254 END DO 1255 1256 IF ( lwp ) WRITE(numout,*) ' Deformation recalculated ' 1240 1257 1241 1258 CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) … … 1288 1305 1289 1306 ENDIF 1307 1308 IF ( lwp ) WRITE(numout,*) ' zs12f recalculated ' 1309 1290 1310 ! 1291 1311 !----------------------- … … 1341 1361 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 1342 1362 1363 IF ( lwp ) WRITE(numout,*) 'Some terms recalculated ' 1364 1343 1365 ! --- Stress tensor invariants (SIMIP diags) --- ! 1344 1366 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN … … 1368 1390 1369 1391 ENDIF 1392 1393 IF ( lwp ) WRITE(numout,*) 'SIMIP work done' 1370 1394 1371 1395 ! --- Normalized stress tensor principal components --- ! … … 1387 1411 zfac = zp_deltastar_t(ji,jj) 1388 1412 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 1413 zsig1 = 0._wp !!! FUCKING DEBUG TEST !!! 1389 1414 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 1390 1415 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) … … 1400 1425 END DO 1401 1426 END DO 1427 IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done' 1402 1428 ! 1403 1429 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 1404 1430 ! 1431 IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll ' 1405 1432 ! 1406 1433 CALL iom_put( 'sig1_pnorm' , zsig1_p ) … … 1408 1435 1409 1436 DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 1437 1438 IF ( lwp ) WRITE(numout,*) ' So what ??? ' 1410 1439 1411 1440 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.