Ignore:
Timestamp:
2008-03-14T15:44:32+01:00 (13 years ago)
Author:
rblod
Message:

Correct a bug in ice rheology, see ticket #78

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/LIM_SRC_3/limrhg.F90

    r834 r866  
    176176         zresr                         !: Local error on velocity 
    177177 
    178       INTEGER :: & 
    179          stress_mean_swi = 1 
    180178! 
    181179!------------------------------------------------------------------------------! 
     
    203201            zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) 
    204202            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) / 2. 
    205             ! tmi = 1 where ther is ice or on land 
    206             ! Claude sees a bug, next line : tmi(ji,jj) 
    207             tmi           = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & 
     203            ! tmi = 1 where there is ice or on land 
     204            tmi(ji,jj)    = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & 
    208205                                    epsd ) ) ) * tms(ji,jj) 
    209206         END DO 
     
    213210      CALL lbc_lnk( zpresh(:,:), 'T',  1. ) 
    214211 
    215 ! Claude sees a bug 
    216 !     CALL lbc_lnk( tmi(:,:), 'T',  1. ) 
     212      CALL lbc_lnk( tmi(:,:), 'T',  1. ) 
    217213 
    218214      ! Ice strength on grid cell corners (zpreshc) 
     
    280276 
    281277            ! Ocean has no slip boundary condition 
    282             v_oce1(ji,jj)  = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji+1,jj)    & 
    283                 &                 +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) & 
     278! GG bug 
     279!           v_oce1(ji,jj)  = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji+1,jj)    & 
     280!               &                 +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji,jj)) & 
     281!               &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)   
     282            v_oce1(ji,jj)  = 0.5*( (v_io(ji,jj)+v_io(ji,jj-1))*e1t(ji,jj)    & 
     283                &                 +(v_io(ji+1,jj)+v_io(ji+1,jj-1))*e1t(ji+1,jj)) & 
    284284                &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)   
    285285 
    286             u_oce2(ji,jj)  = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1)     & 
    287                 &                 +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) & 
     286! GG bug 
     287!           u_oce2(ji,jj)  = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj+1)     & 
     288!               &                 +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj)) & 
     289!               &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
     290            u_oce2(ji,jj)  = 0.5*((u_io(ji,jj)+u_io(ji-1,jj))*e2t(ji,jj)     & 
     291                &                 +(u_io(ji,jj+1)+u_io(ji-1,jj+1))*e2t(ji,jj+1)) & 
    288292                &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
    289293 
     
    322326 
    323327      !-Initialise stress tensor  
    324       IF ( stress_mean_swi .NE. 1 ) THEN 
    325          zs1(:,:)  = 0.0 
    326          zs2(:,:)  = 0.0 
    327          zs12(:,:) = 0.0 
    328       ELSE 
    329          zs1(:,:)  = stress1_i(:,:) 
    330          zs2(:,:)  = stress2_i(:,:) 
    331          zs12(:,:) = stress12_i(:,:) 
    332       ENDIF 
     328      zs1(:,:)  = stress1_i(:,:) 
     329      zs2(:,:)  = stress2_i(:,:) 
     330      zs12(:,:) = stress12_i(:,:) 
    333331 
    334332      v_ice1(:,:) = 0.0 
     
    512510 
    513511                    ! SB modif because ocean has no slip boundary condition 
    514                     zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)     & 
    515                       &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   & 
     512! GG bug 
     513!                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)     & 
     514!                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   & 
     515!                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
     516                    zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         & 
     517                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   & 
    516518                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    517519                    za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 
     
    535537                    z0           = zmass2(ji,jj)/dtevp 
    536538                    ! SB modif because ocean has no slip boundary condition 
    537                     zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)     & 
    538                 &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj))   & 
     539! GG bug 
     540!                   zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)     & 
     541!               &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj))   & 
     542!               &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
     543                    zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     & 
     544                &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   & 
    539545                &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
    540546                    za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &  
     
    558564                    z0           = zmass2(ji,jj)/dtevp 
    559565                    ! SB modif because ocean has no slip boundary condition 
    560                     zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)      & 
    561                        &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj))   & 
     566! GG Bug 
     567!                   zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)      & 
     568!                      &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj))   & 
     569!                      &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)    
     570                    zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      & 
     571                       &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   & 
    562572                       &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)    
    563573 
     
    581591                    z0           = zmass1(ji,jj)/dtevp 
    582592                    ! SB modif because ocean has no slip boundary condition 
    583                     zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)      & 
    584                        &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   & 
     593! GG Bug 
     594!                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)      & 
     595!                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj))   & 
     596!                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
     597                    zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      & 
     598                       &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   & 
    585599                       &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    586600     
     
    628642         END DO 
    629643      END DO 
    630 ! 
    631 !------------------------------------------------------------------------------! 
    632 ! 5) Compute stress rain invariants and store strain tensor 
    633 !------------------------------------------------------------------------------! 
    634 ! 
    635       ! Compute delta, shear and div, inputs for mechanical redistribution  
     644 
     645      ! Recompute delta, shear and div, inputs for mechanical redistribution  
    636646      DO jj = k_j1+1, k_jpj-1 
    637647         DO ji = 2, jpim1 
     
    690700     &                          ) + creepl 
    691701 
    692              divu_i(ji,jj)  = zdd(ji,jj)  
    693              delta_i(ji,jj) = deltat(ji,jj)  
    694              shear_i(ji,jj) = zds(ji,jj) 
    695  
    696702             ENDIF ! zdummy 
    697703 
    698704         END DO !jj 
    699705      END DO !ji 
    700  
    701 ! MV This should be removed as it is the same as previous lines 
    702 !     ! dynamical invariants 
    703 !     IF ( stress_mean_swi .EQ. 0 ) THEN 
    704 ! the following should not be necessary 
    705          DO jj = k_j1+1, k_jpj-1 
    706             DO ji = 2, jpim1 
    707                divu_i(ji,jj)  = zdd(ji,jj) 
    708                delta_i(ji,jj) = deltat (ji,jj) 
    709                shear_i(ji,jj) = zds (ji,jj) 
    710             END DO 
     706! 
     707!------------------------------------------------------------------------------! 
     708! 5) Store stress tensor and its invariants 
     709!------------------------------------------------------------------------------! 
     710! 
     711      ! * Invariants of the stress tensor are required for limitd_me 
     712      ! * Store the stress tensor for the next time step 
     713      ! accelerates convergence and improves stability 
     714      DO jj = k_j1+1, k_jpj-1 
     715         DO ji = 2, jpim1 
     716            divu_i(ji,jj)  = zdd(ji,jj) 
     717            delta_i(ji,jj) = deltat (ji,jj) 
     718            shear_i(ji,jj) = zds (ji,jj) 
     719            stress1_i(ji,jj)  = zs1(ji,jj) 
     720            stress2_i(ji,jj)  = zs2(ji,jj) 
     721            stress12_i(ji,jj) = zs12(ji,jj) 
    711722         END DO 
    712 !     ENDIF 
    713  
     723      END DO 
     724 
     725      ! Lateral boundary condition 
    714726      CALL lbc_lnk( divu_i(:,:) , 'T', 1. ) 
    715727      CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 
    716728      CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 
    717  
    718  
    719       ! Store the stress tensor for next time step 
    720       ! this accelerates convergence and improves stability 
    721       IF ( stress_mean_swi .EQ. 1 ) THEN 
    722          DO jj = k_j1+1, k_jpj-1 
    723             DO ji = 2, jpim1 
    724                stress1_i(ji,jj)  = zs1(ji,jj) 
    725                stress2_i(ji,jj)  = zs2(ji,jj) 
    726                stress12_i(ji,jj) = zs12(ji,jj) 
    727             END DO 
    728          END DO 
    729       ENDIF 
    730  
    731       ! Lateral boundary condition 
    732729      CALL lbc_lnk( stress1_i(:,:), 'T', 1. ) 
    733730      CALL lbc_lnk( stress2_i(:,:), 'T', 1. ) 
    734731      CALL lbc_lnk( stress12_i(:,:), 'F', 1. ) 
    735   
     732 
    736733! 
    737734!------------------------------------------------------------------------------! 
Note: See TracChangeset for help on using the changeset viewer.