New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15375 – NEMO

Changeset 15375


Ignore:
Timestamp:
2021-10-14T19:52:46+02:00 (3 years ago)
Author:
clem
Message:

change convergence test from the max over the domain to the average. This way it resembles more to what exists in the literature. Also, add alternate directions for the upstream advection scheme used in the rheology. However this scheme seems to slow down convergence, hence it should not be used

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

    r15334 r15375  
    181181      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s) 
    182182      !! -- advect fields at the rheology time step for the calculation of strength 
     183      !!    it seems that convergence is worse when ll_advups=true. So it not really a good idea 
    183184      LOGICAL  ::   ll_advups = .FALSE. 
    184185      REAL(wp) ::   zdt_ups 
    … …  
    704705            ENDIF 
    705706            ! 
    706             CALL rhg_upstream( zdt_ups, u_ice, v_ice, za_i_ups )   ! upstream advection: a_i 
    707             CALL rhg_upstream( zdt_ups, u_ice, v_ice, zv_i_ups )   ! upstream advection: v_i 
     707            CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, za_i_ups )   ! upstream advection: a_i 
     708            CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, zv_i_ups )   ! upstream advection: v_i 
    708709            ! 
    709710            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! strength 
    … …  
    967968      REAL(wp)          ::   zresm           ! local real 
    968969      CHARACTER(len=20) ::   clname 
     970      LOGICAL           ::   ll_maxcvg 
     971      REAL(wp), DIMENSION(jpi,jpj,2) ::   zres 
     972      REAL(wp), DIMENSION(2)         ::   ztmp 
    969973      !!---------------------------------------------------------------------- 
    970  
     974      ll_maxcvg = .FALSE. 
     975      ! 
    971976      ! create file 
    972977      IF( kt == nit000 .AND. kiter == 1 ) THEN 
    … …  
    983988            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 
    984989            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime ) 
    985             istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid ) 
     990            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) 
    986991            istatus = NF90_ENDDEF(ncvgid) 
    987992         ENDIF 
    … …  
    9971002      ELSE 
    9981003         zresm = 0._wp 
    999          DO_2D( 0, 0, 0, 0 ) 
    1000             zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
    1001                &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) 
    1002          END_2D 
    1003          CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
     1004         IF( ll_maxcvg ) THEN   ! error max over the domain 
     1005            DO_2D( 0, 0, 0, 0 ) 
     1006               zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1007                  &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) 
     1008            END_2D 
     1009            CALL mpp_max( 'icedyn_rhg_evp', zresm ) 
     1010         ELSE                   ! error averaged over the domain 
     1011            DO_2D( 0, 0, 0, 0 ) 
     1012               zres(ji,jj,1) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1013                  &                 ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) 
     1014               zres(ji,jj,2) = pmsk15(ji,jj) 
     1015            END_2D 
     1016            ztmp(:) = glob_sum_vec( 'icedyn_rhg_evp', zres ) 
     1017            IF( ztmp(2) /= 0._wp )   zresm = ztmp(1) / ztmp(2) 
     1018         ENDIF 
    10041019      ENDIF 
    10051020 
    … …  
    10691084   END SUBROUTINE rhg_evp_rst 
    10701085 
    1071    SUBROUTINE rhg_upstream( pdt, pu, pv, pt ) 
     1086   SUBROUTINE rhg_upstream( jter, pdt, pu, pv, pt ) 
    10721087      !!--------------------------------------------------------------------- 
    10731088      !!                    ***  ROUTINE rhg_upstream  *** 
    … …  
    10751090      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer 
    10761091      !!---------------------------------------------------------------------- 
     1092      INTEGER                    , INTENT(in   ) ::   jter 
    10771093      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step 
    10781094      REAL(wp), DIMENSION(:,:  ) , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
    … …  
    10811097      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    10821098      REAL(wp) ::   ztra          ! local scalar 
    1083       REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups   ! upstream fluxes 
     1099      LOGICAL  ::   ll_upsxy = .TRUE. 
     1100      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zpt   ! upstream fluxes and tracer guess 
    10841101      !!---------------------------------------------------------------------- 
    10851102      DO jl = 1, jpl 
    1086          DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    1087             zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + & 
    1088                &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj  ,jl) 
    1089             zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + & 
    1090                &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj+1,jl) 
    1091          END_2D 
     1103         IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **! 
     1104            ! 
     1105            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
     1106               zfu_ups(ji,jj) = MAX(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji+1,jj,jl) 
     1107               zfv_ups(ji,jj) = MAX(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj+1,jl) 
     1108            END_2D 
     1109            ! 
     1110         ELSE                              !** alternate directions **! 
     1111            ! 
     1112            IF( MOD(jter,2) == 1 ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     1113               ! 
     1114               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )       !-- flux in x-direction 
     1115                  zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj,jl) + & 
     1116                     &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     1117               END_2D 
     1118               ! 
     1119               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls )     !-- first guess of tracer from u-flux 
     1120                  ztra       = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) 
     1121                  zpt(ji,jj) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1122               END_2D 
     1123               ! 
     1124               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )   !-- flux in y-direction 
     1125                  zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj  ) + & 
     1126                     &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj+1) 
     1127               END_2D 
     1128               ! 
     1129            ELSE                          !==  even ice time step:  adv_y then adv_x  ==! 
     1130               ! 
     1131               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )       !-- flux in y-direction 
     1132                  zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj  ,jl) + & 
     1133                     &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     1134               END_2D 
     1135               ! 
     1136               DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 )     !-- first guess of tracer from v-flux 
     1137                  ztra       = - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 
     1138                  zpt(ji,jj) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1139               END_2D 
     1140               ! 
     1141               DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )   !-- flux in x-direction 
     1142                  zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji  ,jj) + & 
     1143                     &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji+1,jj) 
     1144               END_2D 
     1145               ! 
     1146            ENDIF 
     1147            ! 
     1148         ENDIF 
     1149         ! 
    10921150         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    1093             ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 
    1094             ! 
    1095             pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1151            ztra         = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 
     1152            pt(ji,jj,jl) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    10961153         END_2D 
    10971154      END DO 
Note: See TracChangeset for help on using the changeset viewer.