Changeset 13629


Ignore:
Timestamp:
2020-10-19T08:46:02+02:00 (6 months ago)
Author:
vancop
Message:

Pregurvanic review version, mpp bug fix test

File:
1 edited

Legend:

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

    r13593 r13629  
    777777         ll_u_iterate = .TRUE. 
    778778         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!!! 
    781782 
    782783         DO i_inn = 1, nn_ninn_vp ! inner loop iterations 
     
    789790            ! l_full_nf_update = jter == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    790791 
    791             zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 
    792             zv_b(:,:) = v_ice(:,:) 
    793  
    794             zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp 
    795             zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp 
     792            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 
    796797 
    797798            IF ( ll_u_iterate .OR. ll_v_iterate )   THEN 
     
    800801               IF ( ll_u_iterate ) THEN    ! --- Solve for u-velocity --- ! 
    801802                                           ! ---------------------------- ! 
    802        
     803 
    803804                  ! What follows could be subroutinized... 
    804805       
    805806                  ! Thomas Algorithm  for tridiagonal solver 
    806807                  ! 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 
    807810                   
    808811                  DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) 
     
    826829                           ! see Zhang and Hibler, 1997, Appendix B 
    827830                           zAA3 = 0._wp 
    828                            IF ( ji == 2 )          zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 
    829                            IF ( ji == jpi - 1 )    zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj) 
     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) 
    830833 
    831834                           ! right hand side 
     
    854857                        ELSE ! nn_thomas == 2      ! mitGCM algorithm 
    855858 
    856                            ! ji = 2 
     859                           ! --- ji = 2 
    857860                           zw                   =   zBU(2,jj) 
    858861                           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 ) 
    860863                           zCU_prime(2,jj)      =   zw * zCU(2,jj) 
    861864                           zFU_prime(2,jj)      =   zw * zFU(2,jj) 
    862865 
     866                           ! --- ji = 3 --> jpi-1 
    863867                           DO ji = 3, jpi - 1 
    864868 
    865869                              zw                =   zBU(ji,jj) - zAU(ji,jj) * zCU_prime(ji-1,jj) 
    866870                              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 ) 
    868872                              zCU_prime(ji,jj)  =   zw * zCU(ji,jj)  
    869873                              zFU_prime(ji,jj)  =   zw * ( zFU(ji,jj) - zAU(ji,jj) * zFU_prime(ji-1,jj) ) 
     
    882886                           ! last row  
    883887                           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)     = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) &  
     888                           u_ice(jpi-1,jj)    = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) &  
    885889                                  &            * umask(jpi-1,jj,1) 
    886890 
     
    891895                              ! ji2 = jpi - 2 -> ji = 2 
    892896                              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) ) & ! this 
     897                              u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1)  & ! this 
    894898                                   ! line is guilty 
    895899                                              &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 )  
     
    898902                        ELSE ! nn_thomas == 2 ! mitGCM version 
    899903 
    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) 
    901905 
    902906                           DO ji2 = 2 , jpi-2 ! all other rows       !   MV OPT: could be turned into forward loop (by substituting ji) 
    903907                              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) 
    905911                           END DO 
    906912                      
     
    915921                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        &   ! masking 
    916922                              &         * (           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) 
    918924                            
    919925                        END DO 
     
    935941                  !!! ZH97 explain it is critical for convergence speed 
    936942 
     943                  zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp 
     944 
    937945                  DO jn = 1, nn_zebra_vp ! "zebra" loop 
    938946       
     
    946954                      
    947955                        !------------------------ 
    948                         ! Independent term (zFU) 
     956                        ! Independent term (zFV) 
    949957                        !------------------------ 
    950958                        DO jj = 2, jpj - 1 
     
    970978                        IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 
    971979                      
    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                                                                                         
     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 
    979987                        ELSE ! nn_thomas == 2      ! mitGCM algorithm 
    980988                      
    981989                           ! jj = 2 
    982                            zw                   =   zBV(2,jj) 
     990                           zw                   =   zBV(ji,2) 
    983991                           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) 
    987995 
    988996                           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) 
    991998                              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 ) 
    9931000                              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) ) 
    9961002                           END DO 
    997                         
     1003 
    9981004                        ENDIF 
    9991005                      
     
    10151021                              jj =  jpj - jj2    ! index swap 
    10161022                              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) & 
    10181024                                  &                             / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )        
    10191025                           END DO             
     
    10221028 
    10231029                           ! 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) 
    10251031                            
    10261032                           ! other rows 
    1027                            DO jj2 = 2 , jpi-2 ! all other rows 
     1033                           DO jj2 = 2 , jpj-2 ! all other rows 
    10281034                              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) 
    10301036                           END DO 
    1031                          
     1037 
    10321038                        ENDIF ! nn_thomas 
    10331039                         
     
    11381144      IF ( lwp ) WRITE(numout,*) ' We are out of outer loop ' 
    11391145 
     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 
    11401151      CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG 
    11411152      CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG 
     1153      CALL iom_put( 'zCU_prime'     , zCU_prime      ) ! MV DEBUG 
    11421154      CALL iom_put( 'zFU_prime'     , zFU_prime      ) ! MV DEBUG 
    11431155 
    11441156      CALL iom_put( 'zFV'           , zFV            ) ! MV DEBUG 
    11451157      CALL iom_put( 'zBV_prime'     , zBV_prime      ) ! MV DEBUG 
     1158      CALL iom_put( 'zCV_prime'     , zCV_prime      ) ! MV DEBUG 
    11461159      CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG 
    11471160 
     
    11831196 
    11841197      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
     1198 
     1199      IF ( lwp ) WRITE(numout,*) ' Velocity replaced ' 
    11851200 
    11861201      ! END DEBUG 
     
    12381253         END DO 
    12391254      END DO 
     1255 
     1256      IF ( lwp ) WRITE(numout,*) ' Deformation recalculated ' 
    12401257       
    12411258      CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
     
    12881305          
    12891306      ENDIF 
     1307 
     1308      IF ( lwp ) WRITE(numout,*) ' zs12f recalculated ' 
     1309 
    12901310      ! 
    12911311      !----------------------- 
     
    13411361      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    13421362 
     1363      IF ( lwp ) WRITE(numout,*) 'Some terms recalculated ' 
     1364 
    13431365      ! --- Stress tensor invariants (SIMIP diags) --- ! 
    13441366      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     
    13681390          
    13691391      ENDIF 
     1392 
     1393      IF ( lwp ) WRITE(numout,*) 'SIMIP work done' 
    13701394 
    13711395      ! --- Normalized stress tensor principal components --- ! 
     
    13871411               zfac             =   zp_deltastar_t(ji,jj) 
    13881412               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 
     1413               zsig1            = 0._wp !!! FUCKING DEBUG TEST !!! 
    13891414               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
    13901415               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     
    14001425            END DO 
    14011426         END DO 
     1427         IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done' 
    14021428         ! 
    14031429         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 
    14041430         !       
     1431         IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll ' 
    14051432         ! 
    14061433         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     
    14081435 
    14091436         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 
     1437 
     1438         IF ( lwp ) WRITE(numout,*) ' So what ??? ' 
    14101439          
    14111440      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.