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 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 – NEMO

Ignore:
Timestamp:
2015-04-13T15:08:59+02:00 (9 years ago)
Author:
davestorkey
Message:

Merge in changes from trunk up to 5021.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r4688 r5208  
    5050   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2) 
    5151 
    52    REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    53        
    5452   !! * Substitutions 
    5553#  include "vectopt_loop_substitute.h90" 
     
    119117      CHARACTER (len=50) ::   charout 
    120118      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         ! 
    121       REAL(wp) ::   za, zstms, zsang, zmask   ! local scalars 
     119      REAL(wp) ::   za, zstms, zmask   ! local scalars 
     120      REAL(wp) ::   zc1, zc2, zc3             ! ice mass 
    122121 
    123122      REAL(wp) ::   dtevp              ! time step for subcycling 
     
    125124      REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars 
    126125      REAL(wp) ::   zu_ice2, zv_ice1   ! 
    127       REAL(wp) ::   zddc, zdtc, zzdst   ! delta on corners and on centre 
     126      REAL(wp) ::   zddc, zdtc         ! delta on corners and on centre 
     127      REAL(wp) ::   zdst               ! shear at the center of the grid point 
    128128      REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface 
    129129      REAL(wp) ::   sigma1, sigma2     ! internal ice stress 
    130130 
    131131      REAL(wp) ::   zresm         ! Maximal error on ice velocity 
    132       REAL(wp) ::   zindb         ! ice (1) or not (0)       
    133132      REAL(wp) ::   zdummy        ! dummy argument 
    134133      REAL(wp) ::   zintb, zintn  ! dummy argument 
     
    140139      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points 
    141140      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   zc1              ! ice mass 
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   zusw             ! temporary weight for ice strength computation 
    144141      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                            
    145142      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points 
     
    147144      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses 
    148145       
    149       REAL(wp), POINTER, DIMENSION(:,:) ::   zdd   , zdt      ! Divergence and tension at centre of grid cells 
     146      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells 
    150147      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells 
    151       REAL(wp), POINTER, DIMENSION(:,:) ::   zdst             ! Shear on centre of grid cells 
    152       REAL(wp), POINTER, DIMENSION(:,:) ::   deltat, deltac   ! Delta at centre and corners of grid cells 
    153148      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2  
    154149      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12 
     
    160155 
    161156      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    162       CALL wrk_alloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                ) 
    163       CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  ) 
    164       CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
     157      CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                ) 
     158      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
     159      CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
    165160 
    166161#if  defined key_lim2 && ! defined key_lim2_vp 
     
    179174      ! 
    180175      !------------------------------------------------------------------------------! 
    181       ! 1) Ice-Snow mass (zc1), ice strength (zpresh)                                ! 
     176      ! 1) Ice strength (zpresh)                                ! 
    182177      !------------------------------------------------------------------------------! 
    183178      ! 
    184179      ! Put every vector to 0 
    185       zpresh (:,:) = 0._wp   ;   zc1   (:,:) = 0._wp 
     180      delta_i(:,:) = 0._wp   ; 
     181      zpresh (:,:) = 0._wp   ;   
    186182      zpreshc(:,:) = 0._wp 
    187183      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp 
    188       zdd    (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp 
     184      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp 
     185      shear_i(:,:) = 0._wp 
    189186 
    190187#if defined key_lim3 
     
    196193!CDIR NOVERRCHK 
    197194         DO ji = 1 , jpi 
    198             zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) 
    199195#if defined key_lim3 
    200196            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) 
     
    218214               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
    219215               &              tms(ji,jj)     * wght(ji+1,jj+1,1,1) 
    220             zusw(ji,jj)    = 1.0 / MAX( zstms, epsd ) 
    221216            zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 
    222217               &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + & 
    223218               &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &  
    224219               &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   & 
    225                &             ) * zusw(ji,jj) 
     220               &             ) / MAX( zstms, epsd ) 
    226221         END DO 
    227222      END DO 
     
    265260         DO ji = fs_2, fs_jpim1 
    266261 
     262            zc1 = tms(ji  ,jj  ) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
     263            zc2 = tms(ji+1,jj  ) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
     264            zc3 = tms(ji  ,jj+1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
     265 
    267266            zt11 = tms(ji  ,jj) * e1t(ji  ,jj) 
    268267            zt12 = tms(ji+1,jj) * e1t(ji+1,jj) 
     
    275274 
    276275            ! Mass, coriolis coeff. and currents 
    277             zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd) 
    278             zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd) 
     276            zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+epsd) 
     277            zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+epsd) 
    279278            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   & 
    280279               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd ) 
     
    344343               !   
    345344               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 
    346                !- zdd(:,:), zdt(:,:): divergence and tension at centre of grid cells 
     345               !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 
    347346               !- zds(:,:): shear on northeast corner of grid cells 
    348347               ! 
     
    353352               !                      bugs (Martin, for Miguel). 
    354353               ! 
    355                !- ALSO: arrays zdd, zdt, zds and delta could  
     354               !- ALSO: arrays zdt, zds and delta could  
    356355               !  be removed in the future to minimise memory demand. 
    357356               ! 
     
    361360               ! 
    362361               ! 
    363                zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    364                   &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
    365                   &          +e1v(ji,jj)*v_ice(ji,jj)                      & 
    366                   &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
    367                   &          )                                             & 
    368                   &         / area(ji,jj) 
     362               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
     363                  &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
     364                  &             +e1v(ji,jj)*v_ice(ji,jj)                      & 
     365                  &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
     366                  &             )                                             & 
     367                  &            / area(ji,jj) 
    369368 
    370369               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    & 
     
    408407 
    409408               !- Calculate Delta at centre of grid cells 
    410                zzdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          & 
     409               zdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          & 
    411410                  &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          & 
    412411                  &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          & 
     
    415414                  &         / area(ji,jj) 
    416415 
    417                delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 )   
    418                ! MV rewriting 
    419                ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 
    420                !!gm faster to replace the line above with simply: 
    421                !!                deltat(ji,jj) = MAX( delta, creepl ) 
    422                !!gm end   
    423                deltat(ji,jj) = delta + creepl 
    424                ! END MV 
     416               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )   
     417               delta_i(ji,jj) = delta + creepl 
    425418               !-Calculate stress tensor components zs1 and zs2  
    426419               !-at centre of grid cells (see section 3.5 of CICE user's guide). 
    427                !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   & 
    428                !   &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &        
    429                !   &          / ( 1._wp + alphaevp * dtotel ) 
    430  
    431                !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   & 
    432                !              zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   & 
    433                !   &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
    434  
    435                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
    436                zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  & 
     420               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) )  & 
    437421                  &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
    438                zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  & 
     422               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) )  & 
    439423                  &         / ( 1._wp + dtotel ) 
    440424 
     
    468452                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) 
    469453 
    470                deltac(ji,jj) = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 
     454               zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 
    471455 
    472456               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 
    473                !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  & 
    474                !   &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
    475                !   &          / ( 1._wp + alphaevp * ecc2 * dtotel )  
    476  
    477                ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
    478457               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
    479                   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
     458                  &          ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) )  & 
    480459                  &          / ( 1.0 + dtotel )  
    481460 
     
    513492               DO ji = fs_2, fs_jpim1 
    514493                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    515                   zsang        = SIGN ( 1.0 , fcor(ji,jj) ) * sangvg 
    516494                  z0           = zmass1(ji,jj)/dtevp 
    517495 
     
    523501                     (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 
    524502                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 
    525                      za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 
    526                   zcca         = z0+za*cangvg 
    527                   zccb         = zcorl1(ji,jj)+za*zsang 
     503                     za*(u_oce1(ji,jj)) 
     504                  zcca         = z0+za 
     505                  zccb         = zcorl1(ji,jj) 
    528506                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask  
    529507 
     
    536514#endif 
    537515#if defined key_bdy 
    538          ! clem: change u_ice and v_ice at the boundary for each iteration 
    539516         CALL bdy_ice_lim_dyn( 'U' ) 
    540517#endif          
     
    546523 
    547524                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 
    548                   zsang        = SIGN(1.0,fcor(ji,jj))*sangvg 
    549525                  z0           = zmass2(ji,jj)/dtevp 
    550526                  ! SB modif because ocean has no slip boundary condition 
     
    555531                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 
    556532                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + & 
    557                      za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 
    558                   zcca         = z0+za*cangvg 
    559                   zccb         = zcorl2(ji,jj)+za*zsang 
     533                     za2ct(ji,jj) + za*(v_oce2(ji,jj)) 
     534                  zcca         = z0+za 
     535                  zccb         = zcorl2(ji,jj) 
    560536                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 
    561537 
     
    568544#endif 
    569545#if defined key_bdy 
    570          ! clem: change u_ice and v_ice at the boundary for each iteration 
    571546         CALL bdy_ice_lim_dyn( 'V' ) 
    572547#endif          
     
    578553               DO ji = fs_2, fs_jpim1 
    579554                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 
    580                   zsang        = SIGN(1.0,fcor(ji,jj))*sangvg 
    581555                  z0           = zmass2(ji,jj)/dtevp 
    582556                  ! SB modif because ocean has no slip boundary condition 
     
    588562                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 
    589563                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + & 
    590                      za2ct(ji,jj) + za*(cangvg*v_oce2(ji,jj)+zsang*u_oce2(ji,jj)) 
    591                   zcca         = z0+za*cangvg 
    592                   zccb         = zcorl2(ji,jj)+za*zsang 
     564                     za2ct(ji,jj) + za*(v_oce2(ji,jj)) 
     565                  zcca         = z0+za 
     566                  zccb         = zcorl2(ji,jj) 
    593567                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 
    594568 
     
    601575#endif 
    602576#if defined key_bdy 
    603          ! clem: change u_ice and v_ice at the boundary for each iteration 
    604577         CALL bdy_ice_lim_dyn( 'V' ) 
    605578#endif          
     
    610583               DO ji = fs_2, fs_jpim1 
    611584                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    612                   zsang        = SIGN(1.0,fcor(ji,jj))*sangvg 
    613585                  z0           = zmass1(ji,jj)/dtevp 
    614                   ! SB modif because ocean has no slip boundary condition 
    615                   ! GG Bug 
    616                   !                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)      & 
    617                   !                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   & 
    618                   !                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    619586                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      & 
    620587                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   & 
     
    624591                     (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 
    625592                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 
    626                      za*(cangvg*u_oce1(ji,jj)-zsang*v_oce1(ji,jj)) 
    627                   zcca         = z0+za*cangvg 
    628                   zccb         = zcorl1(ji,jj)+za*zsang 
     593                     za*(u_oce1(ji,jj)) 
     594                  zcca         = z0+za 
     595                  zccb         = zcorl1(ji,jj) 
    629596                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask  
    630597               END DO ! ji 
     
    636603#endif 
    637604#if defined key_bdy 
    638          ! clem: change u_ice and v_ice at the boundary for each iteration 
    639605         CALL bdy_ice_lim_dyn( 'U' ) 
    640606#endif          
     
    666632!CDIR NOVERRCHK 
    667633         DO ji = fs_2, fs_jpim1 
    668             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
    669             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    670634            zdummy = vt_i(ji,jj) 
    671635            IF ( zdummy .LE. hminrhg ) THEN 
     
    683647#endif 
    684648#if defined key_bdy 
    685       ! clem: change u_ice and v_ice at the boundary 
    686649      CALL bdy_ice_lim_dyn( 'U' ) 
    687650      CALL bdy_ice_lim_dyn( 'V' ) 
     
    690653      DO jj = k_j1+1, k_jpj-1  
    691654         DO ji = fs_2, fs_jpim1 
    692             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
    693             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    694655            zdummy = vt_i(ji,jj) 
    695656            IF ( zdummy .LE. hminrhg ) THEN 
     
    713674!CDIR NOVERRCHK 
    714675         DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi 
    715             !- zdd(:,:), zdt(:,:): divergence and tension at centre  
     676            !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    716677            !- zds(:,:): shear on northeast corner of grid cells 
    717             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
    718             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    719678            zdummy = vt_i(ji,jj) 
    720679            IF ( zdummy .LE. hminrhg ) THEN 
    721680 
    722                zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    723                   &          -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
    724                   &          +e1v(ji,jj)*v_ice(ji,jj)                      & 
    725                   &          -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
    726                   &         )                                              & 
    727                   &         / area(ji,jj) 
     681               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
     682                  &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
     683                  &             +e1v(ji,jj)*v_ice(ji,jj)                      & 
     684                  &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
     685                  &            )                                              & 
     686                  &            / area(ji,jj) 
    728687 
    729688               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    & 
     
    747706                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
    748707 
    749                zdst(ji,jj) = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    & 
     708               zdst = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    & 
    750709                  &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    & 
    751710                  &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    & 
    752711                  &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj) 
    753712 
    754 !              deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   &  
    755 !                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 &  
    756 !                  &                          ) + creepl 
    757                ! MV rewriting 
    758                delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 )   
    759                deltat(ji,jj) = delta + creepl 
    760                ! END MV 
     713               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )   
     714               delta_i(ji,jj) = delta + creepl 
    761715             
    762716            ENDIF ! zdummy 
     
    773727      DO jj = k_j1+1, k_jpj-1 
    774728         DO ji = fs_2, fs_jpim1 
    775             divu_i (ji,jj) = zdd   (ji,jj) 
    776             delta_i(ji,jj) = deltat(ji,jj) 
    777729            ! begin TECLIM change  
    778             zdst(ji,jj)= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &    
     730            zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &    
    779731               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &    
    780732               &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &    
    781733               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj)  
    782             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) ) 
     734            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    783735            ! end TECLIM change 
    784736         END DO 
     
    834786      ! 
    835787      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    836       CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                ) 
    837       CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  ) 
    838       CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
     788      CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                ) 
     789      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
     790      CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
    839791 
    840792   END SUBROUTINE lim_rhg 
Note: See TracChangeset for help on using the changeset viewer.