Ignore:
Timestamp:
2020-10-14T17:29:11+02:00 (4 months ago)
Author:
vancop
Message:

Modified tridiagonal system

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90

    r13591 r13593  
    131131      LOGICAL ::   ll_u_iterate, ll_v_iterate   ! continue iteration or not 
    132132      ! 
    133       INTEGER ::   ji, ji2, jj, jn          ! dummy loop indices 
     133      INTEGER ::   ji, ji2, jj, jj2, jn          ! dummy loop indices 
    134134      INTEGER ::   jter, i_out, i_inn  !  
    135135      INTEGER ::   ji_min, jj_min      ! 
     
    777777         ll_u_iterate = .TRUE. 
    778778         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 
    786781 
    787782         DO i_inn = 1, nn_ninn_vp ! inner loop iterations 
     
    806801                                           ! ---------------------------- ! 
    807802       
     803                  ! What follows could be subroutinized... 
     804       
    808805                  ! 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 
    810807                   
    811808                  DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) 
     
    819816                     IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min 
    820817 
    821                      !---------------------------------------------------------- 
    822                      ! -- Boundary condition (link with neighbouring processor) 
    823                      !---------------------------------------------------------- 
    824  
    825818                     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                        !------------------------ 
    833823                        DO ji = 2, jpi - 1     
    834824 
    835825                           ! boundary condition substitution 
    836826                           ! see Zhang and Hibler, 1997, Appendix B 
    837                            ! MV NOTE possibly not fully appropriate 
    838827                           zAA3 = 0._wp 
    839828                           IF ( ji == 2 )          zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 
     
    848837                        END DO 
    849838 
    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 
    862844 
    863845                           DO ji = 3, jpi - 1 
     
    870852                           END DO 
    871853                         
    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 
    877855 
    878856                           ! ji = 2 
     
    893871                           END DO 
    894872 
    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 
    907881 
    908882                           ! last row  
     
    912886 
    913887                           DO ji2 = 2 , jpi-2 ! all other rows        
    914                           !DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop 
     888                           !DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop 
    915889                              ji =  jpi - ji2 
    916890                              ! ji2 = 2       -> ji = jpi - 2 
     
    922896                           END DO             
    923897 
    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) 
    931901 
    932902                           DO ji2 = 2 , jpi-2 ! all other rows       !   MV OPT: could be turned into forward loop (by substituting ji) 
     
    935905                           END DO 
    936906                      
    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 
    947911                        DO ji = 2, jpi - 1     
    948912                         
    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 
    950914                            
    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 
    953916                              &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             & 
    954917                              &           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) * umask(ji,jj,1) 
     
    966929               !                           ! ---------------------------- ! 
    967930                                           
    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 
    970937                  DO jn = 1, nn_zebra_vp ! "zebra" loop 
    971938       
     
    977944          
    978945                     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                        !------------------------ 
    989950                        DO jj = 2, jpj - 1 
    990951 
     
    1002963 
    1003964                        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                        !--------------- 
    1007967                        ! Forward sweep 
    1008                         DO jj = 3, jpj - 1 
    1009                            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 DO 
    1014              
    1015                         ! Backward sweep 
    1016                         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 row 
    1019  
    1020                         DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed 
    1021                            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 DO             
    1025  
    1026                         !---------------        
    1027                         ! -- Relaxation 
    1028968                        !--------------- 
     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) 
    10291035                        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 cases                           
    1033                            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)              & 
    10351041                              &           + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) * vmask(ji,jj,1) 
    10361042 
    1037                         END DO 
    1038  
     1043                        END DO ! jj 
     1044                         
    10391045                     END DO ! ji 
    1040  
     1046                      
    10411047                  END DO ! zebra loop 
    10421048                                     
Note: See TracChangeset for help on using the changeset viewer.