Changeset 11951


Ignore:
Timestamp:
2019-11-22T16:52:49+01:00 (2 months ago)
Author:
stefryn
Message:

fixed bug in stress calculation

Location:
NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/ice.F90

    r11877 r11951  
    226226   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field                  [s-1] 
    227227   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aniso_11, aniso_12   !: structure tensor elements 
     228   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdg_conv 
    228229   ! 
    229230   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]      
     
    409410         &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      & 
    410411         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) ,                      & 
    411          &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , STAT=ierr(ii) ) 
     412         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv  (jpi,jpj) , STAT=ierr(ii) ) 
    412413 
    413414      ii = ii + 1 
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90

    r11933 r11951  
    6767CONTAINS 
    6868 
    69    SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, aniso_11, aniso_12 ) 
     69   SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12 ) 
    7070      !!------------------------------------------------------------------- 
    7171      !!                 ***  SUBROUTINE ice_dyn_rhg_eap  *** 
     
    122122      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   ! 
    123123      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    124       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   aniso_11, aniso_12   ! structure tensor components 
     124      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components 
    125125      !! 
    126126      INTEGER ::   ji, jj       ! dummy loop indices 
     
    141141      REAL(wp) ::   zfac_x, zfac_y 
    142142      REAL(wp) ::   zshear, zdum1, zdum2 
    143       REAL(wp) ::   stressptmp, stressmtmp                              ! anisotropic stress tensor components 
     143      REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components 
    144144      REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution 
    145145      ! 
     
    404404               ! --- anisotropic stress calculation --- ! 
    405405               CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, & 
    406                                        zdsT, aniso_11(ji,jj), aniso_12(ji,jj), & 
     406                                       zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 
    407407                                       stressptmp, stressmtmp, & 
    408408                                       stress12tmp(ji,jj), strength(ji,jj), & 
    409409                                       alphar, alphas) 
     410               IF (jter == nn_nevp) THEN 
     411                  rdg_conv(ji,jj) = -min( alphar, 0._wp)     
     412               ENDIF 
    410413 
    411414               ! delta at T points 
    412                !zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     415               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    413416 
    414417               ! P/delta at T points 
    415                !zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
     418               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    416419 
    417420               ! alpha & beta for aEVP 
     
    428431               zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1  ) * z1_alph1 
    429432               zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1  ) * z1_alph1 
    430               
     433               !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1  ) * z1_alph1 
     434               !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1  ) * z1_alph1 
     435               !WRITE(numout,*) 'stress output', ji, jj, zs1(ji,jj) 
    431436            END DO 
    432437         END DO 
     
    437442            DO ji = 1, jpim1 
    438443               ! stress12tmp at F points  
    439                zdsT = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
     444               stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
    440445                  &   + stress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + stress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  & 
    441446                  &   ) * 0.25_wp * r1_e1e2f(ji,jj) 
     
    449454               ! stress at F points (zkt/=0 if landfast) 
    450455               !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
    451                zs12(ji,jj) = ( zs12(ji,jj) + stress12tmp(ji,jj) * zalph1  ) * z1_alph1 
     456               zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1  ) * z1_alph1 
     457               !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1 
    452458 
    453459            END DO 
     
    749755      END DO 
    750756      CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
     757      !CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1. ) 
    751758       
    752759      ! --- Store the stress tensor for the next time step --- ! 
     
    922929      REAL(wp), PARAMETER ::   kfriction = 0.45_wp 
    923930      !!--------------------------------------------------------------------- 
     931      ! Factor to maintain the same stress as in EVP (see Section 3) 
     932      ! Can be set to 1 otherwise 
     933      invstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
     934  
     935      invsin = 1._wp/sin(2._wp*phi) * invstressconviso  
     936      !now uses phi as set in higher code 
     937        
     938      ! compute eigenvalues, eigenvectors and angles for structure tensor, strain 
     939      ! rates 
     940 
    924941      ! 1) structure tensor 
    925942      a22 = 1._wp-a11 
     
    949966      dtemp12 = shear*0.5_wp 
    950967      dtemp22 = 0.5_wp*(divu - tension) 
     968 
     969      !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear 
    951970 
    952971      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)  
     
    10731092      sig22  = strength*(stemp22r + kfriction*stemp22s) * invsin 
    10741093 
     1094      !WRITE(numout,*) 'strength inside loop', strength 
     1095      !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22 
     1096 
    10751097      ! Back - rotation of the stress from principal axes into general coordinates 
    10761098 
     
    10831105      stress12 = sgprm12 
    10841106      stressm  = sgprm11 - sgprm22 
     1107 
     1108      !WRITE(numout,*) 'stress output inside loop', ksub, stressp 
    10851109 
    10861110      ! Compute ridging and sliding functions in general coordinates  
Note: See TracChangeset for help on using the changeset viewer.