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 13466 for NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90 – NEMO

Ignore:
Timestamp:
2020-09-15T09:27:47+02:00 (4 years ago)
Author:
smasson
Message:

r4_trunk: merge r4 13280:13310, see #2523

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90

    r13271 r13466  
    4141   USE prtctl         ! Print control 
    4242 
     43   USE netcdf         ! NetCDF library for convergence test 
    4344   IMPLICIT NONE 
    4445   PRIVATE 
     
    4950   !! * Substitutions 
    5051#  include "vectopt_loop_substitute.h90" 
     52 
     53   !! for convergence tests 
     54   INTEGER ::   ncvgid   ! netcdf file id 
     55   INTEGER ::   nvarid   ! netcdf variable id 
     56   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
    5157   !!---------------------------------------------------------------------- 
    5258   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    119125      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
    120126      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017 
     127      REAl(wp) ::   zbetau, zbetav 
    121128      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    122129      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
     
    125132      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
    126133      ! 
    127       REAL(wp) ::   zresm                                               ! Maximal error on ice velocity 
    128134      REAL(wp) ::   zintb, zintn                                        ! dummy argument 
    129135      REAL(wp) ::   zfac_x, zfac_y 
     
    141147      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
    142148      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    143 !!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence 
    144149      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    145150      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
     
    160165      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small 
    161166      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
     167      !! --- check convergence 
     168      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    162169      !! --- diags 
    163       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00 
    164170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
    165171      !! --- SIMIP diags 
     
    174180      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 
    175181      ! 
    176 !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
     182      ! for diagnostics and convergence tests 
     183      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 
     184      DO jj = 1, jpj 
     185         DO ji = 1, jpi 
     186            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
     187            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     188         END DO 
     189      END DO 
     190      ! 
     191      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
    177192      !------------------------------------------------------------------------------! 
    178193      ! 0) mask at F points for the ice 
     
    222237      z1_ecc2 = 1._wp / ecc2 
    223238 
    224       ! Time step for subcycling 
    225       zdtevp   = rdt_ice / REAL( nn_nevp ) 
    226       z1_dtevp = 1._wp / zdtevp 
    227  
    228239      ! alpha parameters (Bouillon 2009) 
    229240      IF( .NOT. ln_aEVP ) THEN 
    230          zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     241         zdtevp   = rdt_ice / REAL( nn_nevp ) 
     242         zalph1 =   2._wp * rn_relast * REAL( nn_nevp ) 
    231243         zalph2 = zalph1 * z1_ecc2 
    232244 
    233245         z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    234246         z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     247      ELSE 
     248         zdtevp   = rdt_ice 
     249         ! zalpha parameters set later on adaptatively 
    235250      ENDIF 
     251      z1_dtevp = 1._wp / zdtevp 
    236252          
    237253      ! Initialise stress tensor  
     
    244260 
    245261      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    246       IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     262      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile 
    247263      ELSE                         ;   zkt = 0._wp 
    248264      ENDIF 
     
    315331               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    316332               ! ice-bottom stress at U points 
    317                zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 
    318                ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     333               zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 
     334               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    319335               ! ice-bottom stress at V points 
    320                zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 
    321                ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     336               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 
     337               ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    322338               ! ice_bottom stress at T points 
    323                zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 
    324                tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     339               zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 
     340               tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    325341            END DO 
    326342         END DO 
     
    345361         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    346362         ! 
    347 !!$         IF(ln_ctl) THEN   ! Convergence test 
    348 !!$            DO jj = 1, jpjm1 
    349 !!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
    350 !!$               zv_ice(:,jj) = v_ice(:,jj) 
    351 !!$            END DO 
    352 !!$         ENDIF 
     363         ! convergence test 
     364         IF( ln_rhg_chkcvg ) THEN 
     365            DO jj = 1, jpj 
     366               DO ji = 1, jpi 
     367                  zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 
     368                  zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 
     369               END DO 
     370            END DO 
     371         ENDIF 
    353372 
    354373         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     
    391410               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    392411 
    393                ! alpha & beta for aEVP 
     412               ! alpha for aEVP 
    394413               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    395414               !   alpha = beta = sqrt(4*gamma) 
     
    399418                  zalph2   = zalph1 
    400419                  z1_alph2 = z1_alph1 
     420                  ! explicit: 
     421                  ! z1_alph1 = 1._wp / zalph1 
     422                  ! z1_alph2 = 1._wp / zalph1 
     423                  ! zalph1 = zalph1 - 1._wp 
     424                  ! zalph2 = zalph1 
    401425               ENDIF 
    402426                
     
    409433         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 
    410434 
     435         ! Save beta at T-points for further computations 
     436         IF( ln_aEVP ) THEN 
     437            DO jj = 1, jpj 
     438               DO ji = 1, jpi 
     439                  zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     440               END DO 
     441            END DO 
     442         ENDIF 
     443          
    411444         DO jj = 1, jpjm1 
    412445            DO ji = 1, jpim1 
    413446 
    414                ! alpha & beta for aEVP 
     447               ! alpha for aEVP 
    415448               IF( ln_aEVP ) THEN 
    416                   zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     449                  zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 
    417450                  z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    418                   zbeta(ji,jj) = zalph2 
     451                  ! explicit: 
     452                  ! z1_alph2 = 1._wp / zalph2 
     453                  ! zalph2 = zalph2 - 1._wp 
    419454               ENDIF 
    420455                
     
    486521                  ! 
    487522                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    488                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    489                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    490                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    491                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    492                         &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     523                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 
     524                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     525                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     526                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     527                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &  
     528                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     529                        &                                    ) / ( zbetav + 1._wp )                                              & 
     530                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    493531                        &           )   * zmsk00y(ji,jj) 
    494532                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    495                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
    496                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    497                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    498                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    499                         &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    500                         &            )   * zmsk00y(ji,jj) 
     533                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
     534                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     535                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     536                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     537                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     538                        &            )  * zmsk00y(ji,jj) 
    501539                  ENDIF 
    502540               END DO 
     
    537575                  ! 
    538576                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    539                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    540                         &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    541                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    542                         &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    543                         &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     577                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 
     578                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     579                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     580                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     581                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   & 
     582                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     583                        &                                    ) / ( zbetau + 1._wp )                                              & 
     584                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    544585                        &           )   * zmsk00x(ji,jj) 
    545586                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    546                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
    547                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    548                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    549                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    550                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    551                         &            )   * zmsk00x(ji,jj) 
     587                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
     588                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     589                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     590                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     591                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     592                        &           )   * zmsk00x(ji,jj) 
    552593                  ENDIF 
    553594               END DO 
     
    590631                  ! 
    591632                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    592                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    593                         &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    594                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    595                         &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    596                         &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     633                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 
     634                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     635                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     636                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     637                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   & 
     638                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     639                        &                                    ) / ( zbetau + 1._wp )                                              & 
     640                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    597641                        &           )   * zmsk00x(ji,jj) 
    598642                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    599                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
    600                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    601                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    602                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    603                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    604                         &            )   * zmsk00x(ji,jj) 
     643                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
     644                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     645                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     646                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     647                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     648                        &           )   * zmsk00x(ji,jj) 
    605649                  ENDIF 
    606650               END DO 
     
    641685                  ! 
    642686                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    643                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    644                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    645                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    646                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    647                         &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     687                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 
     688                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     689                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     690                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     691                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
     692                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     693                        &                                    ) / ( zbetav + 1._wp )                                              &  
     694                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    648695                        &           )   * zmsk00y(ji,jj) 
    649696                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    650                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
    651                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    652                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    653                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    654                         &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    655                         &            )   * zmsk00y(ji,jj) 
     697                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
     698                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     699                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     700                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     701                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     702                        &           )   * zmsk00y(ji,jj) 
    656703                  ENDIF 
    657704               END DO 
     
    667714         ENDIF 
    668715 
    669 !!$         IF(ln_ctl) THEN   ! Convergence test 
    670 !!$            DO jj = 2 , jpjm1 
    671 !!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    672 !!$            END DO 
    673 !!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 
    674 !!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
    675 !!$         ENDIF 
     716         ! convergence test 
     717         IF( ln_rhg_chkcvg )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 
    676718         ! 
    677719         !                                                ! ==================== ! 
    678720      END DO                                              !  end loop over jter  ! 
    679721      !                                                   ! ==================== ! 
     722      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta ) 
    680723      ! 
    681724      !------------------------------------------------------------------------------! 
     
    734777      ! 5) diagnostics 
    735778      !------------------------------------------------------------------------------! 
    736       DO jj = 1, jpj 
    737          DO ji = 1, jpi 
    738             zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    739          END DO 
    740       END DO 
    741  
    742779      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
    743780      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
     
    796833         DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
    797834      ENDIF 
    798        
     835 
    799836      ! --- SIMIP --- ! 
    800837      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
     
    852889      ENDIF 
    853890      ! 
     891      ! --- convergence tests --- ! 
     892      IF( ln_rhg_chkcvg ) THEN 
     893         IF( iom_use('uice_cvg') ) THEN 
     894            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
     895               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 
     896                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     897            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
     898               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 
     899                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     900            ENDIF 
     901         ENDIF 
     902      ENDIF       
     903      ! 
     904      DEALLOCATE( zmsk00, zmsk15 ) 
     905      ! 
    854906   END SUBROUTINE ice_dyn_rhg_evp 
     907 
     908 
     909   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 
     910      !!---------------------------------------------------------------------- 
     911      !!                    ***  ROUTINE rhg_cvg  *** 
     912      !!                      
     913      !! ** Purpose :   check convergence of oce rheology 
     914      !! 
     915      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity 
     916      !!                during the sub timestepping of rheology so as: 
     917      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 
     918      !!                This routine is called every sub-iteration, so it is cpu expensive 
     919      !! 
     920      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     921      !!---------------------------------------------------------------------- 
     922      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index 
     923      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities 
     924      !! 
     925      INTEGER           ::   it, idtime, istatus 
     926      INTEGER           ::   ji, jj          ! dummy loop indices 
     927      REAL(wp)          ::   zresm           ! local real  
     928      CHARACTER(len=20) ::   clname 
     929      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence 
     930      !!---------------------------------------------------------------------- 
     931 
     932      ! create file 
     933      IF( kt == nit000 .AND. kiter == 1 ) THEN 
     934         ! 
     935         IF( lwp ) THEN 
     936            WRITE(numout,*) 
     937            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 
     938            WRITE(numout,*) '~~~~~~~' 
     939         ENDIF 
     940         ! 
     941         IF( lwm ) THEN 
     942            clname = 'ice_cvg.nc' 
     943            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
     944            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 
     945            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime ) 
     946            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid ) 
     947            istatus = NF90_ENDDEF(ncvgid) 
     948         ENDIF 
     949         ! 
     950      ENDIF 
     951 
     952      ! time 
     953      it = ( kt - 1 ) * kitermax + kiter 
     954       
     955      ! convergence 
     956      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 
     957         zresm = 0._wp 
     958      ELSE 
     959         DO jj = 1, jpj 
     960            DO ji = 1, jpi 
     961               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     962                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 
     963            END DO 
     964         END DO 
     965         zresm = MAXVAL( zres ) 
     966         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
     967      ENDIF 
     968 
     969      IF( lwm ) THEN 
     970         ! write variables 
     971         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
     972         ! close file 
     973         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid) 
     974      ENDIF 
     975       
     976   END SUBROUTINE rhg_cvg 
    855977 
    856978 
     
    9101032   END SUBROUTINE rhg_evp_rst 
    9111033 
     1034    
    9121035#else 
    9131036   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.