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 6763 for branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 – NEMO

Ignore:
Timestamp:
2016-06-30T17:17:35+02:00 (8 years ago)
Author:
clem
Message:

new parameterization from Lupkes et al. (2012) of ice-atm and ocean-atm drags (dependent on ice morphology) => ln_Cd_L12 in namelist_ref

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r6746 r6763  
    116116      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2  ! alpha and beta from Bouillon 2009 and 2013 
    117117      REAL(wp) ::   zm1, zm2, zm3                              ! ice/snow mass 
    118       REAL(wp) ::   delta, zp_delf, zds2 
     118      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 
    119119      REAL(wp) ::   zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel  ! temporary scalars 
    120120 
     
    134134      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2               ! internal stresses 
    135135       
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zdt, zds                  ! tension and shear 
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                       ! shear 
    137137      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12            ! stress tensor components 
    138138      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr     ! check convergence 
     
    149149      CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
    150150      CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
    151       CALL wrk_alloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     151      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    152152      CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2 ) 
    153153 
     
    257257         END DO 
    258258      END DO 
    259  
    260259      ! 
    261260      !------------------------------------------------------------------------------! 
     
    266265      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    267266         !                                            !----------------------!         
    268          ! Convergence test 
    269          IF(ln_ctl) THEN 
     267         IF(ln_ctl) THEN   ! Convergence test 
    270268            DO jj = 1, jpjm1 
    271269               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
     
    275273 
    276274         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
    277          DO jj = 2, jpjm1 
    278             DO ji = fs_2, fs_jpim1 
    279                 
    280                ! divergence at T points 
    281                divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    282                   &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    283                   &            ) * r1_e12t(ji,jj) 
    284  
    285                ! tension at T points 
    286                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    287                   &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    288                   &         ) * r1_e12t(ji,jj) 
     275         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
     276            DO ji = 1, fs_jpim1 
    289277 
    290278               ! shear at F points 
     
    292280                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    293281                  &         ) * r1_e12f(ji,jj) 
     282 
     283            END DO 
     284         END DO 
     285         CALL lbc_lnk( zds, 'F', 1. ) 
     286 
     287         DO jj = 2, jpjm1 
     288            DO ji = 2, jpim1 ! no vector loop 
     289 
     290               ! shear**2 at T points (doc eq. A16) 
     291               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     292                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     293                  &   ) * 0.25_wp * r1_e12t(ji,jj) 
     294               
     295               ! divergence at T points 
     296               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     297                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     298                  &    ) * r1_e12t(ji,jj) 
     299               zdiv2 = zdiv * zdiv 
     300                
     301               ! tension at T points 
     302               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)   & 
     303                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     304                  &   ) * r1_e12t(ji,jj) 
     305               zdt2 = zdt * zdt 
     306                
     307               ! delta at T points 
     308               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )   
     309 
     310               ! P/delta at T points 
     311               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl ) 
     312                
     313               ! stress at T points 
     314               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 
     315               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 
     316              
     317            END DO 
     318         END DO 
     319         CALL lbc_lnk( zp_delt, 'T', 1. ) 
     320 
     321         DO jj = 1, jpjm1 
     322            DO ji = 1, jpim1 
     323 
     324               ! P/delta at F points 
     325               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
     326                
     327               ! stress at F points 
     328               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 
     329 
     330            END DO 
     331         END DO 
     332         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     333  
     334 
     335         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
     336         DO jj = 2, jpjm1 
     337            DO ji = fs_2, fs_jpim1                
     338 
     339               ! U points 
     340               zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     341                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     342                  &                    ) * r1_e2u(ji,jj)                                                                      & 
     343                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     344                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     345                  &                  ) * r1_e12u(ji,jj) 
     346 
     347               ! V points 
     348               zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     349                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     350                  &                    ) * r1_e1v(ji,jj)                                                                      & 
     351                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     352                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     353                  &                  ) * r1_e12v(ji,jj) 
    294354 
    295355               ! u_ice at V point 
     
    303363            END DO 
    304364         END DO 
    305          CALL lbc_lnk( zds, 'F', 1. ) 
    306  
    307          DO jj = 2, jpjm1 
    308             DO ji = 2, jpim1 ! no vector loop 
    309  
    310                ! shear**2 at T points (doc eq. A16) 
    311                zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
    312                   &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
    313                   &   ) * 0.25_wp * r1_e12t(ji,jj) 
    314                
    315                ! delta at T points 
    316                delta          = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zds2 ) * usecc2 )   
    317                delta_i(ji,jj) = delta + rn_creepl 
    318  
    319                ! P/delta at T points 
    320                zp_delt(ji,jj) = zpresh(ji,jj) / delta_i(ji,jj) 
    321             END DO 
    322          END DO 
    323          CALL lbc_lnk( zp_delt, 'T', 1. ) 
    324  
    325          ! --- Stress tensor --- ! 
    326          DO jj = 2, jpjm1 
    327             DO ji = 2, jpim1 ! no vector loop 
    328                 
    329                ! P/delta at F points 
    330                zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
    331                 
    332                ! stress tensor at T points 
    333                zs1(ji,jj) = ( zs1 (ji,jj) * zalph1 + zp_delt(ji,jj) * ( divu_i(ji,jj) - (delta_i(ji,jj)-rn_creepl) ) ) * z1_alph1 
    334                zs2(ji,jj) = ( zs2 (ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt(ji,jj) * z1_ecc2  )                      ) * z1_alph2 
    335                ! F points 
    336                zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf        * ( zds(ji,jj) * z1_ecc2  ) * 0.5_wp             ) * z1_alph2 
    337             END DO 
    338          END DO 
    339          CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    340   
    341          ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
    342          DO jj = 2, jpjm1 
    343             DO ji = fs_2, fs_jpim1                
    344                ! U points 
    345                zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
    346                   &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
    347                   &                    ) * r1_e2u(ji,jj)                                                                      & 
    348                   &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
    349                   &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    350                   &                  ) * r1_e12u(ji,jj) 
    351  
    352                ! V points 
    353                zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    354                   &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
    355                   &                    ) * r1_e1v(ji,jj)                                                                      & 
    356                   &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
    357                   &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    358                   &                  ) * r1_e12v(ji,jj) 
    359             END DO 
    360          END DO 
    361365         ! 
    362366         ! --- Computation of ice velocity --- ! 
     
    368372               DO ji = fs_2, fs_jpim1 
    369373 
    370                   zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    371                      &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
    372  
     374                  zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 
     375                   
    373376                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
    374                   zTauO   =   za2(ji,jj) * rhoco * zvel 
    375                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     377                  zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     378                     &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
     379                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    376380                  zCor    =   zcor2(ji,jj)  * u_ice2(ji,jj) 
    377381                  zmt     =   zmass2(ji,jj) * z1_dtevp 
     
    403407                     &               + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)  
    404408                   
    405                   zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    406                      &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
    407  
     409                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     410                   
    408411                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
    409                   zTauO   =   za1(ji,jj) * rhoco * zvel 
    410                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     412                  zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     413                     &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     414                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    411415                  zCor    =   zcor1(ji,jj)  * zv_ice1 
    412416                  zmt     =   zmass1(ji,jj) * z1_dtevp 
     
    435439               DO ji = fs_2, fs_jpim1 
    436440 
    437                   zvel    = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    438                      &          + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     441                  zvel    = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    439442 
    440443                  zTauA   =   za1(ji,jj) * utau_ice(ji,jj) 
    441                   zTauO   =   za1(ji,jj) * rhoco * zvel 
    442                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     444                  zTauO   =   za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     445                     &                                 + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 
     446                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    443447                  zCor    =   zcor1(ji,jj)  * v_ice1(ji,jj) 
    444448                  zmt     =   zmass1(ji,jj) * z1_dtevp 
     
    469473                     &               + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    470474 
    471                   zvel    = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    472                      &          + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
     475                  zvel    = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 
    473476          
    474477                  zTauA   =   za2(ji,jj) * vtau_ice(ji,jj) 
    475                   zTauO   =   za2(ji,jj) * rhoco * zvel 
    476                   ztauB   = - tau_icebfr(ji,jj) / zvel 
     478                  zTauO   =   za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     479                     &                                 + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 
     480                  ztauB   = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 
    477481                  zCor    =   zcor2(ji,jj)  * zu_ice2 
    478482                  zmt     =   zmass2(ji,jj) * z1_dtevp 
     
    485489                  !!   &           + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj)  & 
    486490                  !!   &           ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 
    487  
    488491               END DO 
    489492            END DO 
     
    499502         ENDIF 
    500503          
    501          ! Convergence test 
    502          IF(ln_ctl) THEN 
     504         IF(ln_ctl) THEN   ! Convergence test 
    503505            DO jj = 2 , jpjm1 
    504506               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
     
    539541      ! 5) Recompute delta, shear and div (inputs for mechanical redistribution)  
    540542      !------------------------------------------------------------------------------! 
    541  
    542       ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     543      DO jj = 1, jpjm1 
     544         DO ji = 1, fs_jpim1 
     545 
     546            ! shear at F points 
     547            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)   & 
     548               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     549               &         ) * r1_e12f(ji,jj) 
     550 
     551         END DO 
     552      END DO            
     553      CALL lbc_lnk( zds, 'F', 1. ) 
     554       
    543555      DO jj = 2, jpjm1 
    544          DO ji = fs_2, fs_jpim1 
    545              
     556         DO ji = 2, jpim1 ! no vector loop 
     557             
     558            ! tension**2 at T points 
     559            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)   & 
     560               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     561               &   ) * r1_e12t(ji,jj) 
     562            zdt2 = zdt * zdt 
     563             
     564            ! shear**2 at T points (doc eq. A16) 
     565            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     566               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     567               &   ) * 0.25_wp * r1_e12t(ji,jj) 
     568             
     569            ! shear at T points 
     570            shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     571 
    546572            ! divergence at T points 
    547573            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     
    549575               &            ) * r1_e12t(ji,jj) 
    550576             
    551             ! tension at T points 
    552             zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    553                &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    554                &         ) * r1_e12t(ji,jj) 
    555              
    556             ! shear at F points 
    557             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)   & 
    558                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    559                &         ) * r1_e12f(ji,jj) 
    560              
     577            ! delta at T points 
     578            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )   
     579            delta_i(ji,jj) = zdelta + rn_creepl 
     580 
    561581         END DO 
    562582      END DO 
    563       CALL lbc_lnk_multi( divu_i, 'T', 1., zds, 'F', 1. ) 
    564        
    565  
    566       DO jj = 2, jpjm1 
    567          DO ji = 2, jpim1 ! no vector loop 
    568              
    569             ! shear**2 at T points (doc eq. A16) 
    570             zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
    571                &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
    572                &   ) * 0.25_wp * r1_e12t(ji,jj) 
    573              
    574             ! delta at T points 
    575             delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zds2 ) * usecc2 )   
    576             delta_i(ji,jj) = delta + rn_creepl 
    577              
    578             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zds2 ) 
    579          END DO 
    580       END DO 
    581       CALL lbc_lnk_multi( delta_i, 'T', 1.,  shear_i, 'T', 1. ) 
     583      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 
    582584       
    583585      ! --- Store the stress tensor for the next time step --- ! 
     
    626628      CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 
    627629      CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 
    628       CALL wrk_dealloc( jpi,jpj, zds, zdt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     630      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
    629631      CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2 ) 
    630632 
Note: See TracChangeset for help on using the changeset viewer.