Changeset 13681


Ignore:
Timestamp:
2020-10-27T13:54:43+01:00 (6 months ago)
Author:
vancop
Message:

Post-gurvanian subroutine volume 2 including several bugfixes for Clement check

Location:
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE
Files:
2 edited

Legend:

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

    r13629 r13681  
    4343   !! for convergence tests 
    4444   INTEGER ::   ncvgid        ! netcdf file id 
    45    INTEGER ::   nvarid_ucvg   ! netcdf variable id 
    4645   INTEGER ::   nvarid_ures 
    4746   INTEGER ::   nvarid_vres 
     
    4948   INTEGER ::   nvarid_udif 
    5049   INTEGER ::   nvarid_vdif 
     50   INTEGER ::   nvarid_veldif 
    5151   INTEGER ::   nvarid_mke 
    52    INTEGER ::   nvarid_usub, nvarid_vsub 
     52   INTEGER ::   nvarid_ures_xy, nvarid_vres_xy 
    5353 
    5454   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
     
    136136      INTEGER ::   nn_zebra_vp         ! number of zebra steps 
    137137 
    138       INTEGER ::   nn_thomas           ! tridiagonal system version (1=mitgcm, 2=martin) 
    139  
    140138      INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)       
    141139      ! 
     
    144142      REAL(wp) ::   zglob_area                                          ! global ice area for diagnostics 
    145143      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    146       REAL(wp) ::   zm2, zm3, zmassU, zmassV                            ! ice/snow mass and volume 
     144      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass and volume 
    147145      REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars 
    148146      REAL(wp) ::   zp_deltastar_f                                      !  
     
    294292         END DO 
    295293      END DO 
    296       CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp ) 
    297294 
    298295      ! Lateral boundary conditions on velocity (modify zfmask) 
     
    335332      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
    336333      zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 
     334 
     335      zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:)       ! Snow and ice mass at T-point 
     336      zmf(:,:) = zmt(:,:) * ff_t(:,:)                      ! Coriolis factor at T points (m*f) 
    337337       
    338338      DO jj = 2, jpj - 1 
     
    344344 
    345345            ! Snow and ice mass at U-V points 
    346             zmt(ji,jj)      = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )  
    347             zm2             = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) ) 
    348             zm3             = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) ) 
     346            zm1             = zmt(ji,jj) 
     347            zm2             = zmt(ji+1,jj) 
     348            zm3             = zmt(ji,jj+1) 
    349349             
    350             zmassU          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    351             zmassV          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     350            zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     351            zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    352352                           
    353353            ! Mass per unit area divided by time step 
     
    362362            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 
    363363            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 
    364  
    365             ! Coriolis factor at T points (m*f) 
    366             zmf(ji,jj)      = zmt(ji,jj) * ff_t(ji,jj) 
    367364             
    368365            ! Wind stress 
     
    379376 
    380377            ! Mask for lots of ice (1) or little ice (0) 
    381             IF( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
    382             ELSE                                                     ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
    383             IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
    384             ELSE                                                     ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
     378            IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
     379            ELSE                                                      ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     380            IF ( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
     381            ELSE                                                      ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
     382 
     383! MV TEST DEBUG 
     384            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji+1,jj)  <= zmmin )     .AND.  & 
     385               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji+1,jj) <= zamin ) )    THEN   ;   zmsk01x(ji,jj) = 0._wp 
     386            ELSE                                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     387 
     388            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji,jj+1)  <= zmmin )     .AND.  & 
     389               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji,jj+1) <= zamin ) )    THEN   ;   zmsk01y(ji,jj) = 0._wp 
     390            ELSE                                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
     391! END MV TEST DEBUG 
    385392 
    386393         END DO 
     
    391398      CALL iom_put( 'zmsk01x'    , zmsk01x  )   ! MV DEBUG 
    392399      CALL iom_put( 'zmsk01y'    , zmsk01y  )   ! MV DEBUG 
    393       CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG 
    394400      CALL iom_put( 'ztaux_ai'   , ztaux_ai )   ! MV DEBUG 
    395401      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG 
     
    449455         IF( lwp )   WRITE(numout,*) ' outer loop  1a i_out : ', i_out 
    450456 
    451          DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    452             DO ji = 2, jpi - 1 !  
     457         !DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
     458         !   DO ji = 2, jpi - 1 !  
     459 
     460! MV DEBUG 
     461         DO jj = 2, jpj        ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
     462            DO ji = 2, jpi     !  
     463! END MV DEBUG 
    453464 
    454465               ! shear**2 at T points (doc eq. A16) 
     
    500511                
    501512               ! Temporary zef factor at F-point 
    502                zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 
     513               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) 
    503514 
    504515            END DO 
     
    545556 
    546557         IF( lwp )   WRITE(numout,*) ' outer loop  1d i_out : ', i_out 
    547          CALL lbc_lnk( 'icedyn_rhg_vp', zCwU, 'U', 1. ) 
    548          CALL lbc_lnk( 'icedyn_rhg_vp', zCwV, 'V', 1. ) 
     558          
     559         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCwU ,  'U', -1., zCwV, 'V', -1. ) 
     560         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCorU,  'U', -1., zCorV, 'V', -1. ) 
    549561 
    550562         CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG 
    551563         CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG 
    552          IF( lwp )   WRITE(numout,*) ' outer loop  1e i_out : ', i_out 
    553564         CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG 
    554565         CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG 
     566 
    555567         IF( lwp )   WRITE(numout,*) ' outer loop  1f i_out : ', i_out 
    556568          
     
    652664         END DO 
    653665          
    654          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', 1., zrhsv, 'V',  1. ) 
    655  
     666         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V',  -1.) 
     667         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V',  -1.) 
     668         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V',  -1.) 
     669 
     670         CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG 
     671         CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG 
     672         CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG 
     673         CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG 
    656674         CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG 
    657675         CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG 
     
    778796         ll_v_iterate = .TRUE. 
    779797 
    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!!! 
    782  
    783798         DO i_inn = 1, nn_ninn_vp ! inner loop iterations 
    784799 
     
    840855                        END DO 
    841856 
    842                         !--------------- 
    843                         ! Forward sweep 
    844                         !--------------- 
    845           
    846                         IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 
    847  
    848                            DO ji = 3, jpi - 1 
    849  
    850                               zfac             = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) ) 
    851                               zw               = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 )  
    852                               zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 
    853                               zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 
    854  
    855                            END DO 
    856                          
    857                         ELSE ! nn_thomas == 2      ! mitGCM algorithm 
    858  
    859                            ! --- ji = 2 
    860                            zw                   =   zBU(2,jj) 
    861                            zfac                 =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 
    862                            zw                   =   zfac /  MAX ( ABS( zw ) , epsi20 ) 
    863                            zCU_prime(2,jj)      =   zw * zCU(2,jj) 
    864                            zFU_prime(2,jj)      =   zw * zFU(2,jj) 
    865  
    866                            ! --- ji = 3 --> jpi-1 
    867                            DO ji = 3, jpi - 1 
    868  
    869                               zw                =   zBU(ji,jj) - zAU(ji,jj) * zCU_prime(ji-1,jj) 
    870                               zfac              =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 
    871                               zw                =   zfac /  MAX ( ABS( zw ) , epsi20 ) 
    872                               zCU_prime(ji,jj)  =   zw * zCU(ji,jj)  
    873                               zFU_prime(ji,jj)  =   zw * ( zFU(ji,jj) - zAU(ji,jj) * zFU_prime(ji-1,jj) ) 
    874  
    875                            END DO 
    876  
    877                         ENDIF 
     857                     END DO 
     858                      
     859                     CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U',  1. ) 
     860                      
     861                     !--------------- 
     862                     ! Forward sweep 
     863                     !---------------    
     864                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
     865       
     866                        DO ji = 3, jpi - 1 
     867 
     868                           zfac             = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) ) 
     869                           zw               = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 )  
     870                           zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 
     871                           zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 
     872 
     873                        END DO 
     874 
     875                     END DO 
     876 
     877                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime, 'U',  1., zBU_prime, 'U', 1. ) 
    878878  
    879                         !----------------------------- 
    880                         ! Backward sweep & relaxation 
    881                         !----------------------------- 
     879                     !----------------------------- 
     880                     ! Backward sweep & relaxation 
     881                     !----------------------------- 
     882 
     883                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
    882884                     
    883885                        ! --- Backward sweep  
    884                         IF ( nn_thomas == 1 ) THEN ! MV version 
    885  
    886                            ! last row  
    887                            zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , 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 ) &  
    889                                   &            * umask(jpi-1,jj,1) 
    890  
    891                            DO ji2 = 2 , jpi-2 ! all other rows        
    892                            !DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop 
    893                               ji =  jpi - ji2 
    894                               ! ji2 = 2       -> ji = jpi - 2 
    895                               ! ji2 = jpi - 2 -> ji = 2 
    896                               zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 
    897                               u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1)   & ! this 
    898                                    ! line is guilty 
    899                                               &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 )  
    900                            END DO             
    901  
    902                         ELSE ! nn_thomas == 2 ! mitGCM version 
    903  
    904                            u_ice(jpi-1,jj)    =   zFU_prime(jpi-1,jj) * umask(jpi-1,jj,1) 
    905  
    906                            DO ji2 = 2 , jpi-2 ! all other rows       !   MV OPT: could be turned into forward loop (by substituting ji) 
    907                               ji =  jpi - ji2 
    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) 
    911                            END DO 
    912                       
    913                         ENDIF 
     886                        ! last row  
     887                        zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , 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 ) &  
     889                                           &            * umask(jpi-1,jj,1) 
     890                        DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop 
     891                           zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 
     892                           u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1)   &  
     893                                           &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 )  
     894                        END DO 
    914895 
    915896                        !--- Relaxation 
     
    972953                        END DO 
    973954                         
    974                         !--------------- 
    975                         ! Forward sweep 
    976                         !--------------- 
     955                     END DO 
     956 
     957                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV, 'V',  1.) 
    977958                      
    978                         IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 
     959                     !--------------- 
     960                     ! Forward sweep 
     961                     !--------------- 
     962                     DO ji = ji_min, jpi - 1, nn_zebra_vp  
    979963                      
    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  
    987                         ELSE ! nn_thomas == 2      ! mitGCM algorithm 
     964                        DO jj = 3, jpj - 1 
     965 
     966                           zfac             = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 
     967                           zw               = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 ) 
     968                           zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 
     969                           zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1)  
     970 
     971                        END DO 
     972 
     973                     END DO 
     974 
     975                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV_prime, 'V',  1., zBV_prime, 'V', 1. ) 
    988976                      
    989                            ! jj = 2 
    990                            zw                   =   zBV(ji,2) 
    991                            zfac                 =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 
    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) 
    995  
    996                            DO jj = 3, jpj - 1 
    997                               zw                =   zBV(ji,jj) - zAV(ji,jj) * zCV_prime(ji,jj-1) 
    998                               zfac              =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 
    999                               zw                =   zfac /  MAX ( ABS( zw ) , epsi20 ) 
    1000                               zCV_prime(ji,jj)  =   zw * zCV(ji,jj)  
    1001                               zFV_prime(ji,jj)  =   zw * ( zFV(ji,jj) - zAV(ji,jj) * zFV_prime(ji,jj-1) ) 
    1002                            END DO 
    1003  
    1004                         ENDIF 
    1005                       
    1006                         !----------------------------- 
    1007                         ! Backward sweep & relaxation 
    1008                         !----------------------------- 
     977                     !----------------------------- 
     978                     ! Backward sweep & relaxation 
     979                     !----------------------------- 
     980                     DO ji = ji_min, jpi - 1, nn_zebra_vp  
    1009981                     
    1010982                        ! --- Backward sweep  
    1011                         IF ( nn_thomas == 1 ) THEN ! nn_thomas = 1    ! --- MV version seemingly more than mitGCM algorithm 
    1012               
    1013                            ! last row 
    1014                            zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) ) 
    1015                            v_ice(ji,jpj-1)  = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) &  
    1016                                   &         * vmask(ji,jpj-1,1)  ! last row 
    1017  
    1018                            ! other rows 
    1019                            !DO jj = jpj-2, 2, -1 ! original back loop 
    1020                            DO jj2 = 2 , jpj-2    ! modified forward loop 
    1021                               jj =  jpj - jj2    ! index swap 
    1022                               zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) ) 
    1023                               v_ice(ji,jj)   = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) & 
    1024                                   &                             / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )        
    1025                            END DO             
    1026                             
    1027                         ELSE                       ! nn_thomas == 2   ! --- mitGCM algorithm 
    1028  
    1029                            ! last row                         
    1030                            v_ice(ji,jpj-1)   =   zFV_prime(ji,jpj-1) * vmask(ji,jpj-1,1) 
    1031                             
    1032                            ! other rows 
    1033                            DO jj2 = 2 , jpj-2 ! all other rows 
    1034                               jj =  jpj - jj2 
    1035                               v_ice(ji,jj)   = ( zFV_prime(ji,jj) - zCV_prime(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) 
    1036                            END DO 
    1037  
    1038                         ENDIF ! nn_thomas 
    1039                          
     983                        ! last row 
     984                        zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) ) 
     985                        v_ice(ji,jpj-1)  = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) &  
     986                                         &         * vmask(ji,jpj-1,1)  ! last row 
     987 
     988                        ! other rows 
     989                        DO jj = jpj-2, 2, -1 ! original back loop 
     990                           zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) ) 
     991                           v_ice(ji,jj)   = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) & 
     992                                          &  / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )        
     993                        END DO             
     994                                                    
    1040995                        ! --- Relaxation  & masking (should it be now or later) 
    1041996                        DO jj = 2, jpj - 1 
     
    11251080                   
    11261081               ENDIF ! ll_v_iterate 
    1127                 
    1128                !--------------------------------------------------------------------------------------- 
    1129                ! 
    1130                ! --- Calculate extra convergence diagnostics and save them 
    1131                ! 
    1132                !--------------------------------------------------------------------------------------- 
    1133                IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
    1134                           &                         zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 
    1135  
    1136                IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 
    1137                 
     1082 
    11381083            ENDIF ! ---    end ll_u_iterate or ll_v_iterate 
     1084                
     1085            !--------------------------------------------------------------------------------------- 
     1086            ! 
     1087            ! --- Calculate extra convergence diagnostics and save them 
     1088            ! 
     1089            !--------------------------------------------------------------------------------------- 
     1090 
     1091            IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
     1092                      &                         zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 
     1093 
     1094            IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 
    11391095 
    11401096         END DO ! i_inn, end of inner loop 
     
    15961552            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim ) 
    15971553 
    1598             istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid_ucvg ) ! the name uice_cvg sucks, no ? 
    15991554            ! i suggest vel_dif instead 
    1600             istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE   , (/ idtime /), nvarid_ures ) 
    1601             istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE   , (/ idtime /), nvarid_vres ) 
    1602             istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE   , (/ idtime /), nvarid_velres ) 
    1603             istatus = NF90_DEF_VAR( ncvgid, 'u_dif', NF90_DOUBLE   , (/ idtime /), nvarid_udif ) 
    1604             istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE   , (/ idtime /), nvarid_vdif ) 
    1605             istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE   , (/ idtime /), nvarid_mke ) 
    1606  
    1607             istatus = NF90_DEF_VAR( ncvgid, 'usub'    , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_usub) ! --> u-velocity in sub-iterations 
    1608             istatus = NF90_DEF_VAR( ncvgid, 'vsub'    , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_vsub) ! --> v-velocity in sub-iterations 
     1555            istatus = NF90_DEF_VAR( ncvgid, 'u_res'  , NF90_DOUBLE  , (/ idtime /), nvarid_ures ) 
     1556            istatus = NF90_DEF_VAR( ncvgid, 'v_res'  , NF90_DOUBLE  , (/ idtime /), nvarid_vres ) 
     1557            istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE  , (/ idtime /), nvarid_velres ) 
     1558            istatus = NF90_DEF_VAR( ncvgid, 'u_dif'  , NF90_DOUBLE  , (/ idtime /), nvarid_udif ) 
     1559            istatus = NF90_DEF_VAR( ncvgid, 'v_dif'  , NF90_DOUBLE  , (/ idtime /), nvarid_vdif ) 
     1560            istatus = NF90_DEF_VAR( ncvgid, 'vel_dif', NF90_DOUBLE  , (/ idtime /), nvarid_veldif ) 
     1561            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE  , (/ idtime /), nvarid_mke ) 
     1562 
     1563            istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy) 
     1564            istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy) 
    16091565 
    16101566            istatus = NF90_ENDDEF(ncvgid) 
     
    16201576                                            ! if puerrmask and pverrmax are masked at 15% (TEST) 
    16211577 
    1622       IF ( lwp ) WRITE(numout,*) ' zfeldif : ', zveldif 
    1623        
    16241578      ! ---  Mean residual and kinetic energy 
    16251579      IF ( kiter == 1 ) THEN 
    16261580 
    1627               zu_res_mean = 0._wp 
    1628               zv_res_mean = 0._wp 
    1629               zvelres     = 0._wp 
    1630               zmke        = 0._wp 
     1581         zu_res_mean = 0._wp 
     1582         zv_res_mean = 0._wp 
     1583         zvelres     = 0._wp 
     1584         zmke        = 0._wp 
    16311585 
    16321586      ELSE 
     
    16411595 
    16421596      z1_pglob_area = 1._wp / pglob_area 
     1597 
     1598      zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 
    16431599 
    16441600      DO jj = 2, jpj - 1 
     
    16531609            zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 
    16541610 
    1655           END DO 
    1656        END DO                   
    1657  
    1658        IF ( lwp ) WRITE(numout,*) ' TEST 2 '  
    1659        zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 
    1660        zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 
    1661        IF ( lwp ) WRITE(numout,*) ' TEST 3 '  
    1662        zvelres     = 0.5_wp * ( zu_res_mean + zv_res_mean ) 
    1663  
    1664        IF ( lwp ) WRITE(numout,*) ' TEST 4 '  
     1611         END DO 
     1612      END DO                   
     1613 
     1614      IF ( lwp ) WRITE(numout,*) ' TEST 2 '  
     1615      zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 
     1616      zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 
     1617      IF ( lwp ) WRITE(numout,*) ' TEST 3 '  
     1618      zvelres     = 0.5_wp * ( zu_res_mean + zv_res_mean ) 
     1619 
     1620      IF ( lwp ) WRITE(numout,*) ' TEST 4 '  
    16651621                                           
    1666        ! -- Global mean kinetic energy per unit area (J/m2)   
    1667        DO jj = 2, jpj - 1 
    1668           DO ji = 2, jpi - 1                    
    1669              zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 
    1670              zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 
    1671              zvel2(:,:)  = zu*zu + zv*zv              ! square of ice velocity at T-point   
    1672           END DO 
    1673        END DO 
     1622      ! -- Global mean kinetic energy per unit area (J/m2) 
     1623      zvel2(:,:) = 0._wp 
     1624      DO jj = 2, jpj - 1 
     1625         DO ji = 2, jpi - 1                    
     1626            zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 
     1627            zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 
     1628            zvel2(ji,jj)  = zu*zu + zv*zv              ! square of ice velocity at T-point   
     1629         END DO 
     1630      END DO 
    16741631        
    1675          IF ( lwp ) WRITE(numout,*) ' TEST 5 '  
    1676        zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 
    1677          IF ( lwp ) WRITE(numout,*) ' TEST 6 '  
    1678  
    1679       ENDIF 
     1632      IF ( lwp ) WRITE(numout,*) ' TEST 5 '  
     1633 
     1634      zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 
     1635 
     1636      IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 
     1637 
     1638      ENDIF ! kiter 
    16801639                   
    16811640         !                                                ! ==================== ! 
     
    16871646      IF( lwm ) THEN 
    16881647         ! write variables 
    1689             istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) )     ! max U or V velocity diff between subiterations 
    1690             istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 
    1691             istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 
    1692             istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) )   ! average of u- and v- residuals 
    1693             istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) )   ! max U velocity difference, inner iterations 
    1694             istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) )   ! max V velocity difference, inner iterations 
    1695             istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) )         ! mean kinetic energy 
    1696  
    1697             istatus = NF90_PUT_VAR( ncvgid, nvarid_usub, (/pu/), (/it/) )          ! u-velocity 
    1698             istatus = NF90_PUT_VAR( ncvgid, nvarid_vsub, (/pv/), (/it/) )          ! v-velocity 
     1648         istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 
     1649         istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 
     1650         istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) )   ! average of u- and v- residuals 
     1651         istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) )   ! max U velocity difference, inner iterations 
     1652         istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) )   ! max V velocity difference, inner iterations 
     1653         istatus = NF90_PUT_VAR( ncvgid, nvarid_veldif, (/zveldif/), (/it/), (/1/) )   ! max U or V velocity diff between subiterations 
     1654         istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) )         ! mean kinetic energy 
     1655 
     1656         ! 
     1657         IF ( kiter == kitermax ) THEN 
     1658            WRITE(numout,*) ' Should plot the spatially dependent residual ' 
     1659            istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) )          ! U-residual, spatially dependent 
     1660            istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) )          ! V-residual, spatially dependent 
     1661         ENDIF 
    16991662 
    17001663         ! close file 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/iceistate.F90

    r12988 r13681  
    165165      u_ice (:,:) = 0._wp 
    166166      v_ice (:,:) = 0._wp 
     167 
     168      CALL lbc_lnk_multi( 'ice_istate' , u_ice, 'U', -1., v_ice, 'V', -1. ) 
    167169      ! 
    168170      !------------------------------------------------------------------------ 
Note: See TracChangeset for help on using the changeset viewer.