Ignore:
Timestamp:
2017-09-08T18:19:17+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part5 - very nearly finished

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90

    r8512 r8515  
    5252CONTAINS 
    5353 
    54    SUBROUTINE ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
     54   SUBROUTINE ice_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i ) 
    5555      !!------------------------------------------------------------------- 
    5656      !!                 ***  SUBROUTINE ice_rhg_evp  *** 
     
    108108      !!------------------------------------------------------------------- 
    109109      INTEGER, INTENT(in) ::   kt      ! time step 
    110       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   stress1_i, stress2_i, stress12_i 
    111       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   u_ice, v_ice, shear_i, divu_i, delta_i  
     110      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i 
     111      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i  
    112112      !! 
    113113      INTEGER ::   ji, jj       ! dummy loop indices 
     
    247247 
    248248      ! Initialise stress tensor  
    249       zs1 (:,:) = stress1_i (:,:)  
    250       zs2 (:,:) = stress2_i (:,:) 
    251       zs12(:,:) = stress12_i(:,:) 
     249      zs1 (:,:) = pstress1_i (:,:)  
     250      zs2 (:,:) = pstress2_i (:,:) 
     251      zs12(:,:) = pstress12_i(:,:) 
    252252 
    253253      ! Ice strength 
     
    340340         IF(ln_ctl) THEN   ! Convergence test 
    341341            DO jj = 1, jpjm1 
    342                zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
    343                zv_ice(:,jj) = v_ice(:,jj) 
     342               zu_ice(:,jj) = pu_ice(:,jj) ! velocity at previous time step 
     343               zv_ice(:,jj) = pv_ice(:,jj) 
    344344            END DO 
    345345         ENDIF 
     
    350350 
    351351               ! shear at F points 
    352                zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    353                   &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     352               zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     353                  &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    354354                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
    355355 
     
    367367               
    368368               ! divergence at T points 
    369                zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    370                   &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     369               zdiv  = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   & 
     370                  &    + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   & 
    371371                  &    ) * r1_e1e2t(ji,jj) 
    372372               zdiv2 = zdiv * zdiv 
    373373                
    374374               ! tension at T points 
    375                zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    376                   &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     375               zdt  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     376                  &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    377377                  &   ) * r1_e1e2t(ji,jj) 
    378378               zdt2 = zdt * zdt 
     
    426426                  ! 
    427427                  !                !--- u_ice at V point 
    428                u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    429                   &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     428               u_iceV(ji,jj) = 0.5_wp * ( ( pu_ice(ji,jj  ) + pu_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     429                  &                     + ( pu_ice(ji,jj+1) + pu_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    430430                  ! 
    431431                  !                !--- v_ice at U point 
    432                v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    433                   &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     432               v_iceU(ji,jj) = 0.5_wp * ( ( pv_ice(ji  ,jj) + pv_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     433                  &                     + ( pv_ice(ji+1,jj) + pv_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    434434               ! 
    435435            END DO 
     
    444444               DO ji = fs_2, fs_jpim1 
    445445                     !                 !--- tau_io/(v_oce - v_ice) 
    446                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     446                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) )  & 
    447447                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    448448                     !                 !--- Ocean-to-Ice stress 
    449                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     449                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 
    450450                     ! 
    451451                     !                 !--- tau_bottom/v_ice 
    452                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     452                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    453453                  zTauB = - tau_icebfr(ji,jj) / zvel 
    454454                     ! 
    455455                     !                 !--- Coriolis at V-points (energy conserving formulation) 
    456456                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    457                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    458                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     457                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  & 
     458                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 
    459459                     ! 
    460460                     !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    465465                     ! 
    466466                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    467                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
    468                      &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    469                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    470                      &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    471                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
    472                      &           ) * zmaskV(ji,jj) 
     467                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity 
     468                     &                                      + zTauE + zTauO * pv_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     469                     &                                      ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     470                     &              + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     471                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin 
     472                     &            ) * zmaskV(ji,jj) 
    473473                     ! 
    474474                  ! Bouillon 2013 
    475                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     475                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) )                  & 
    476476                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    477477                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
     
    479479               END DO 
    480480            END DO 
    481             CALL lbc_lnk( v_ice, 'V', -1. ) 
     481            CALL lbc_lnk( pv_ice, 'V', -1. ) 
    482482            ! 
    483483#if defined key_agrif 
     
    491491                                
    492492                  ! tau_io/(u_oce - u_ice) 
    493                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     493                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) )  & 
    494494                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    495495 
    496496                  ! Ocean-to-Ice stress 
    497                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     497                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 
    498498 
    499499                  ! tau_bottom/u_ice 
    500                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     500                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 
    501501                  zTauB = - tau_icebfr(ji,jj) / zvel 
    502502 
    503503                  ! Coriolis at U-points (energy conserving formulation) 
    504504                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    505                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    506                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     505                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  & 
     506                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 
    507507                   
    508508                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    513513 
    514514                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    515                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
    516                      &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    517                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    518                      &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    519                      &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin  
    520                      &           ) * zmaskU(ji,jj) 
     515                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity 
     516                     &                                      + zTauE + zTauO * pu_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     517                     &                                      ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     518                     &              + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     519                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin  
     520                     &            ) * zmaskU(ji,jj) 
    521521                  ! Bouillon 2013 
    522                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     522                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) )                  & 
    523523                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    524524                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    525525               END DO 
    526526            END DO 
    527             CALL lbc_lnk( u_ice, 'U', -1. ) 
     527            CALL lbc_lnk( pu_ice, 'U', -1. ) 
    528528            ! 
    529529#if defined key_agrif 
     
    539539 
    540540                  ! tau_io/(u_oce - u_ice) 
    541                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     541                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( pu_ice(ji,jj) - u_oce (ji,jj) ) * ( pu_ice(ji,jj) - u_oce (ji,jj) )  & 
    542542                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    543543 
    544544                  ! Ocean-to-Ice stress 
    545                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     545                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - pu_ice(ji,jj) ) 
    546546 
    547547                  ! tau_bottom/u_ice 
    548                   zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     548                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + pu_ice(ji,jj) * pu_ice(ji,jj) ) ) 
    549549                  zTauB = - tau_icebfr(ji,jj) / zvel 
    550550 
    551551                  ! Coriolis at U-points (energy conserving formulation) 
    552552                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    553                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    554                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     553                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * pv_ice(ji  ,jj) + e1v(ji  ,jj-1) * pv_ice(ji  ,jj-1) )  & 
     554                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * pv_ice(ji+1,jj) + e1v(ji+1,jj-1) * pv_ice(ji+1,jj-1) ) ) 
    555555 
    556556                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    561561 
    562562                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    563                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
    564                      &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    565                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    566                      &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    567                      &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
    568                      &           ) * zmaskU(ji,jj) 
     563                  pu_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * pu_ice(ji,jj)                            &  ! previous velocity 
     564                     &                                      + zTauE + zTauO * pu_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     565                     &                                      ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     566                     &              + ( 1._wp - rswitch ) * pu_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     567                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin 
     568                     &            ) * zmaskU(ji,jj) 
    569569                  ! Bouillon 2013 
    570                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     570                  !!pu_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * pu_ice(ji,jj) + pu_ice_b(ji,jj) )                  & 
    571571                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    572572                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
    573573               END DO 
    574574            END DO 
    575             CALL lbc_lnk( u_ice, 'U', -1. ) 
     575            CALL lbc_lnk( pu_ice, 'U', -1. ) 
    576576            ! 
    577577#if defined key_agrif 
     
    585585          
    586586                  ! tau_io/(v_oce - v_ice) 
    587                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     587                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( pv_ice(ji,jj) - v_oce (ji,jj) ) * ( pv_ice(ji,jj) - v_oce (ji,jj) )  & 
    588588                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    589589 
    590590                  ! Ocean-to-Ice stress 
    591                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     591                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - pv_ice(ji,jj) ) 
    592592 
    593593                  ! tau_bottom/v_ice 
    594                   zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     594                  zvel  = MAX( zepsi, SQRT( pv_ice(ji,jj) * pv_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    595595                  ztauB = - tau_icebfr(ji,jj) / zvel 
    596596                   
    597597                  ! Coriolis at V-points (energy conserving formulation) 
    598598                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    599                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    600                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     599                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * pu_ice(ji,jj  ) + e2u(ji-1,jj  ) * pu_ice(ji-1,jj  ) )  & 
     600                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * pu_ice(ji,jj+1) + e2u(ji-1,jj+1) * pu_ice(ji-1,jj+1) ) ) 
    601601 
    602602                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     
    607607                   
    608608                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
    609                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
    610                      &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    611                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
    612                      &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    613                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
    614                      &           ) * zmaskV(ji,jj) 
     609                  pv_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * pv_ice(ji,jj)                            &  ! previous velocity 
     610                     &                                      + zTauE + zTauO * pv_ice(ji,jj)                            &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     611                     &                                      ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
     612                     &              + ( 1._wp - rswitch ) * pv_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     613                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                   &  ! v_ice = v_oce if mass < zmmin 
     614                     &            ) * zmaskV(ji,jj) 
    615615                  ! Bouillon 2013 
    616                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     616                  !!pv_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * pv_ice(ji,jj) + pv_ice_b(ji,jj) )                  & 
    617617                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    618618                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    619619               END DO 
    620620            END DO 
    621             CALL lbc_lnk( v_ice, 'V', -1. ) 
     621            CALL lbc_lnk( pv_ice, 'V', -1. ) 
    622622            ! 
    623623#if defined key_agrif 
     
    631631         IF(ln_ctl) THEN   ! Convergence test 
    632632            DO jj = 2 , jpjm1 
    633                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
     633               zresr(:,jj) = MAX( ABS( pu_ice(:,jj) - zu_ice(:,jj) ), ABS( pv_ice(:,jj) - zv_ice(:,jj) ) ) 
    634634            END DO 
    635635            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 
     
    648648 
    649649            ! shear at F points 
    650             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    651                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     650            zds(ji,jj) = ( ( pu_ice(ji,jj+1) * r1_e1u(ji,jj+1) - pu_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     651               &         + ( pv_ice(ji+1,jj) * r1_e2v(ji+1,jj) - pv_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    652652               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
    653653 
     
    660660             
    661661            ! tension**2 at T points 
    662             zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    663                &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     662            zdt  = ( ( pu_ice(ji,jj) * r1_e2u(ji,jj) - pu_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     663               &   - ( pv_ice(ji,jj) * r1_e1v(ji,jj) - pv_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    664664               &   ) * r1_e1e2t(ji,jj) 
    665665            zdt2 = zdt * zdt 
     
    671671             
    672672            ! shear at T points 
    673             shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     673            pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
    674674 
    675675            ! divergence at T points 
    676             divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    677                &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    678                &            ) * r1_e1e2t(ji,jj) 
     676            pdivu_i(ji,jj) = ( e2u(ji,jj) * pu_ice(ji,jj) - e2u(ji-1,jj) * pu_ice(ji-1,jj)   & 
     677               &             + e1v(ji,jj) * pv_ice(ji,jj) - e1v(ji,jj-1) * pv_ice(ji,jj-1)   & 
     678               &             ) * r1_e1e2t(ji,jj) 
    679679             
    680680            ! delta at T points 
    681             zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
     681            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    682682            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    683             delta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     683            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
    684684 
    685685         END DO 
    686686      END DO 
    687       CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 
     687      CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
    688688       
    689689      ! --- Store the stress tensor for the next time step --- ! 
    690       stress1_i (:,:) = zs1 (:,:) 
    691       stress2_i (:,:) = zs2 (:,:) 
    692       stress12_i(:,:) = zs12(:,:) 
     690      pstress1_i (:,:) = zs1 (:,:) 
     691      pstress2_i (:,:) = zs2 (:,:) 
     692      pstress12_i(:,:) = zs12(:,:) 
    693693      ! 
    694694 
     
    703703 
    704704      ! --- divergence, shear and strength --- ! 
    705       IF( iom_use('idive')  )   CALL iom_put( "idive"  , divu_i(:,:)   * zswi(:,:) )   ! divergence 
    706       IF( iom_use('ishear') )   CALL iom_put( "ishear" , shear_i(:,:)  * zswi(:,:) )   ! shear 
     705      IF( iom_use('idive')  )   CALL iom_put( "idive"  , pdivu_i(:,:)   * zswi(:,:) )   ! divergence 
     706      IF( iom_use('ishear') )   CALL iom_put( "ishear" , pshear_i(:,:)  * zswi(:,:) )   ! shear 
    707707      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength 
    708708 
     
    714714         DO jj = 2, jpjm1 
    715715            DO ji = 2, jpim1 
    716                zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    717                   &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  & 
    718                   &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
    719  
    720                zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
     716               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
     717                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
     718                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
     719 
     720               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    721721 
    722722               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    723723 
    724 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    725 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    726 !!               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
     724!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
     725!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
     726!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    727727!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    728                zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
     728               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    729729               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress 
    730                zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
     730               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    731731            END DO 
    732732         END DO 
     
    773773                
    774774               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    775                zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
    776                zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
     775               zfac_x = 0.5 * pu_ice(ji,jj) * e2u(ji,jj) * rswitch 
     776               zfac_y = 0.5 * pv_ice(ji,jj) * e1v(ji,jj) * rswitch 
    777777                
    778778               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
Note: See TracChangeset for help on using the changeset viewer.