Changeset 13574


Ignore:
Timestamp:
2020-10-07T18:02:40+02:00 (3 months ago)
Author:
stefryn
Message:

corrected bug in shear calculation, updated testcase (example output figs for EVP and EAP included)

Location:
NEMO/branches/2019/dev_r11842_SI3-10_EAP
Files:
2 added
7 edited

Legend:

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

    r13155 r13574  
    760760      !                              !--------------------------------------------------! 
    761761         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 
    762          ismooth = 1 
     762!         ismooth = 1    ! original code 
     763         ismooth = 0    ! try for EAP stability 
    763764         !                           !--------------------------------------------------! 
    764765      ELSE                           ! Zero strength                                    ! 
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90

    r13155 r13574  
    391391               zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
    392392                  &   + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    393                   &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     393!                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)       .25 was an error 
     394                  &   ) * r1_e1e2t(ji,jj) 
    394395               !WRITE(numout,*) 'shear output', ji, jj, zdsT 
    395396                
     
    972973      ! Factor to maintain the same stress as in EVP (see Section 3) 
    973974      ! Can be set to 1 otherwise 
    974       zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
     975!      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
     976      zinvstressconviso = 1._wp 
    975977  
    976978      zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso  
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/make_INITICE.py

    r12263 r13574  
    8484 
    8585# --- add noise in initial thickness to help nucleation of deformation --- 
    86 for y in np.arange(1,1000,1) : 
    87     for x in np.arange(1,1000,1) : 
    88         cdfhti[:,y,x] = rn_hti_ini+0.02*rn_hti_ini*np.random.uniform(-1,1) 
     86#for y in np.arange(1,1000,1) : 
     87#    for x in np.arange(1,1000,1) : 
     88#        cdfhti[:,y,x] = rn_hti_ini+0.02*rn_hti_ini*np.random.uniform(-1,1) 
     89cdfhti[:,:,:]=rn_hti_ini 
    8990 
    9091nc.close() 
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_cfg

    r12263 r13574  
    3232   nn_it000    =       1      !  first time step 
    3333   nn_itend    =     720     !  last  time step (std 5475) 
     34!   nn_itend    =     1440     !  last  time step (std 5475) 
    3435   nn_istate   =       1      !  output the initial state (1) or not (0) 
    3536/ 
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/EXPREF/namelist_ice_cfg

    r12263 r13574  
    4848&namdyn_rhg     !   Ice rheology 
    4949!------------------------------------------------------------------------------ 
    50    ln_rhg_EVP       = .false.          !  EVP rheology 
    51       ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
     50   ln_rhg_EVP       = .true.          !  EVP rheology 
     51      ln_aEVP       = .false.          !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
    5252      rn_creepl     =   2.0e-9        !     creep limit [1/s] 
    5353      rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve           
    54       nn_nevp       = 120             !     number of EVP subcycles                              
     54      nn_nevp       =   200           !     number of EVP subcycles                              
    5555      rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast  
    5656                                      !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 
    57    ln_rhg_EAP       = .true.          !  EVP rheology 
     57   ln_rhg_EAP       = .false.          !  EVP rheology 
    5858/ 
    5959!------------------------------------------------------------------------------ 
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90

    r12263 r13574  
    4949   PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90 
    5050 
    51    REAL(wp), PARAMETER           ::   phi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg) 
     51   REAL(wp), PARAMETER           ::   pphi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg) 
    5252 
    5353   ! look-up table for calculating structure tensor 
     
    6767CONTAINS 
    6868 
    69    SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12 ) 
     69   SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12, prdg_conv ) 
    7070      !!------------------------------------------------------------------- 
    7171      !!                 ***  SUBROUTINE ice_dyn_rhg_eap  *** 
     
    123123      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    124124      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components 
     125      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   prdg_conv                             ! for ridging 
    125126      !! 
    126       INTEGER ::   ji, jj, ij0, ij1, ii0, ii1     ! dummy loop indices 
     127      INTEGER ::   ji, jj       ! dummy loop indices 
    127128      INTEGER ::   jter         ! local integers 
    128129      ! 
     
    141142      REAL(wp) ::   zfac_x, zfac_y 
    142143      REAL(wp) ::   zshear, zdum1, zdum2 
    143       REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components 
    144       REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution 
     144      REAL(wp) ::   zstressptmp, zstressmtmp, zstress12tmpF                ! anisotropic stress tensor components 
     145      REAL(wp) ::   zalphar, zalphas                                      ! for mechanical redistribution 
     146      REAL(wp) ::   zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A        ! for structure tensor evolution 
    145147      REAL(wp) ::   invw                                                ! for test case 
    146148      ! 
    147       REAL(wp), DIMENSION(jpi,jpj) ::   stress12tmp                     ! anisotropic stress tensor component for regridding 
     149      REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                     ! anisotropic stress tensor component for regridding 
     150      REAL(wp), DIMENSION(jpi,jpj) ::   zyield11, zyield22, zyield12       ! yield surface tensor for history 
    148151      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    149152      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
     
    254257      zs2 (:,:) = pstress2_i (:,:) 
    255258      zs12(:,:) = pstress12_i(:,:) 
     259 
     260      ! constants for structure tensor 
     261      z1_dtevp_A = z1_dtevp/10.0_wp 
     262      z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp) 
     263      zp5kth = 0.5_wp * 0.00002_wp 
    256264 
    257265      ! Ice strength 
     
    408416                
    409417               ! --- anisotropic stress calculation --- ! 
    410                CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, & 
     418               CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, & 
    411419                                       zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 
    412                                        stressptmp, stressmtmp, & 
    413                                        stress12tmp(ji,jj), strength(ji,jj), & 
    414                                        alphar, alphas) 
     420                                       zstressptmp, zstressmtmp, & 
     421                                       zstress12tmp(ji,jj), strength(ji,jj), & 
     422                                       zalphar, zalphas) 
     423 
     424               ! structure tensor update 
     425               IF (mod(jter,10) == 0) THEN 
     426                  CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 
     427                                  paniso_11(ji,jj), paniso_12(ji,jj),           & 
     428                                  zmresult11,  zmresult12) 
     429 
     430                  paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A  + zp5kth - zmresult11) * z1dtevpkth ! implicit 
     431                  paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  - zmresult12) * z1dtevpkth ! implicit 
     432               ENDIF 
     433 
     434 
    415435               IF (jter == nn_nevp) THEN 
    416                   rdg_conv(ji,jj) = -min( alphar, 0._wp)     
     436                  zyield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp) 
     437                  zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp) 
     438                  zyield12(ji,jj) = zstress12tmp(ji,jj) 
     439                  prdg_conv(ji,jj) = -min( zalphar, 0._wp)     
    417440               ENDIF 
    418441 
     
    426449               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    427450               !   alpha = beta = sqrt(4*gamma) 
    428                IF( ln_aEVP ) THEN 
    429                   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
    430                   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    431                ENDIF 
     451               !IF( ln_aEVP ) THEN 
     452               !   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     453               !   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     454               !ENDIF 
    432455                
    433456               ! stress at T points (zkt/=0 if landfast) 
    434457               !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
    435458               !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    436                zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1  ) * z1_alph1 
    437                zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1  ) * z1_alph1 
     459               zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1  ) * z1_alph1 
     460               zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1  ) * z1_alph1 
    438461               !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1  ) * z1_alph1 
    439462               !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1  ) * z1_alph1 
     
    442465         END DO 
    443466         !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 
    444          CALL lbc_lnk( 'icedyn_rhg_eap', stress12tmp, 'T', 1. ) 
     467         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 
    445468 
    446469         DO jj = 1, jpjm1 
    447470            DO ji = 1, jpim1 
    448471               ! stress12tmp at F points  
    449                stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
    450                   &   + stress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + stress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  & 
     472               zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
     473                  &   + zstress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + zstress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  & 
    451474                  &   ) * 0.25_wp * r1_e1e2f(ji,jj) 
    452475 
     
    459482               ! stress at F points (zkt/=0 if landfast) 
    460483               !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
    461                zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1  ) * z1_alph1 
     484               zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1  ) * z1_alph1 
    462485               !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1 
    463486 
    464487            END DO 
    465488         END DO 
     489      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    466490 
    467491         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
     
    507531                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    508532                  zTauB = ztauy_base(ji,jj) / zvel 
    509                   !                 !--- OceanBottom-to-Ice stress 
     533                  !                 !--- OceanBottom-to-Ice stress  
    510534                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    511535                  ! 
     
    539563!extra code for test case boundary conditions 
    540564                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    541                      v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj)) 
     565                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
    542566                  END IF 
    543  
    544567               END DO 
    545568            END DO 
     
    557580                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    558581                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    559                   !                 !--- Ocean-to-Ice stress 
     582                  !                 !--- Ocean-to-Ice stress  
    560583                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    561584                  ! 
     
    595618!extra code for test case boundary conditions 
    596619                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    597                      u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj)) 
     620                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
    598621                  END IF 
    599  
    600622               END DO 
    601623            END DO 
     
    615637                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    616638                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    617                   !                 !--- Ocean-to-Ice stress 
     639                  !                 !--- Ocean-to-Ice stress  
    618640                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    619641                  ! 
     
    653675!extra code for test case boundary conditions 
    654676                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    655                      u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj)) 
     677                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
    656678                  END IF 
    657  
    658679               END DO 
    659680            END DO 
     
    709730!extra code for test case boundary conditions 
    710731                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    711                      v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj)) 
     732                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
    712733                  END IF 
    713  
    714734               END DO 
    715735            END DO 
     
    732752!!$         ENDIF 
    733753         ! 
     754 
    734755         !                                                ! ==================== ! 
    735756      END DO                                              !  end loop over jter  ! 
    736757      !                                                   ! ==================== ! 
     758      ! 
     759      CALL lbc_lnk( 'icedyn_rhg_eap', prdg_conv, 'T', 1. )  ! only need this in ridging module after rheology completed 
    737760      ! 
    738761      !------------------------------------------------------------------------------! 
     
    838861!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    839862               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    840                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
     863               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress  
    841864               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    842865            END DO 
     
    850873         ! Stress tensor invariants (normal and shear stress N/m) 
    851874         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    852          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
     875         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress  
    853876 
    854877         DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
    855878      ENDIF 
    856879 
     880      ! --- yieldcurve --- ! 
     881      IF( iom_use('yield11') .OR. iom_use('yield12') .OR. iom_use('yield22')) THEN 
     882 
     883         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1., zyield22, 'T', 1., zyield12, 'T', 1. ) 
     884 
     885         CALL iom_put( 'yield11', zyield11 * zmsk00 ) 
     886         CALL iom_put( 'yield22', zyield22 * zmsk00 ) 
     887         CALL iom_put( 'yield12', zyield12 * zmsk00 ) 
     888      ENDIF 
     889 
    857890      ! --- anisotropy tensor --- ! 
    858       IF( iom_use('aniso') ) THEN 
    859          CALL iom_put( 'aniso' , aniso_11 * zmsk00 ) 
     891      IF( iom_use('aniso') ) THEN        
     892         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1. )  
     893         CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) 
    860894      ENDIF 
    861895       
     
    917951   END SUBROUTINE ice_dyn_rhg_eap 
    918952 
    919    SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu, tension, & 
    920                                  shear, a11, a12, & 
    921                                  stressp,  stressm, & 
    922                                  stress12, strength, & 
    923                                  alphar, alphas) 
     953   SUBROUTINE update_stress_rdg (ksub, kndte, pdivu, ptension, & 
     954                                 pshear, pa11, pa12, & 
     955                                 pstressp,  pstressm, & 
     956                                 pstress12, strength, & 
     957                                 palphar, palphas) 
    924958      !!--------------------------------------------------------------------- 
    925       !!                   ***  SUBROUTINE rhg_eap_rst  *** 
    926       !! Updates the stress depending on values of strain rate and structure 
    927       !! tensor and for ksub=ndte it computes closing and sliding rate 
     959      !!                   ***  SUBROUTINE update_stress_rdg  *** 
     960      !!                      
     961      !! ** Purpose :   Updates the stress depending on values of strain rate and structure 
     962      !!                tensor and for the last subcycle step it computes closing and sliding rate 
    928963      !!--------------------------------------------------------------------- 
    929       INTEGER,  INTENT(in   ) ::   ksub, ndte 
    930       REAL(wp), INTENT(in   ) ::   phi, strength 
    931       REAL(wp), INTENT(in   ) ::   divu, tension, shear 
    932       REAL(wp), INTENT(in   ) ::   a11, a12 
    933       REAL(wp), INTENT(  out) ::   stressp, stressm, stress12 
    934       REAL(wp), INTENT(  out) ::   alphar, alphas 
     964      INTEGER,  INTENT(in   ) ::   ksub, kndte 
     965      REAL(wp), INTENT(in   ) ::   strength 
     966      REAL(wp), INTENT(in   ) ::   pdivu, ptension, pshear 
     967      REAL(wp), INTENT(in   ) ::   pa11, pa12 
     968      REAL(wp), INTENT(  out) ::   pstressp, pstressm, pstress12 
     969      REAL(wp), INTENT(  out) ::   palphar, palphas 
    935970 
    936971      INTEGER  ::   kx ,ky, ka 
    937972 
    938       REAL(wp) ::   stemp11r, stemp12r, stemp22r    
    939       REAL(wp) ::   stemp11s, stemp12s, stemp22s  
    940       REAL(wp) ::   a22, Qd11Qd11, Qd11Qd12, Qd12Qd12  
    941       REAL(wp) ::   Q11Q11, Q11Q12, Q12Q12  
    942       REAL(wp) ::   dtemp11, dtemp12, dtemp22  
    943       REAL(wp) ::   rotstemp11r, rotstemp12r, rotstemp22r    
    944       REAL(wp) ::   rotstemp11s, rotstemp12s, rotstemp22s 
    945       REAL(wp) ::   sig11, sig12, sig22  
    946       REAL(wp) ::   sgprm11, sgprm12, sgprm22  
    947       REAL(wp) ::   invstressconviso  
    948       REAL(wp) ::   Angle_denom_gamma,  Angle_denom_alpha  
    949       REAL(wp) ::   Tany_1, Tany_2  
    950       REAL(wp) ::   x, y, dx, dy, da, kxw, kyw, kaw 
    951       REAL(wp) ::   invdx, invdy, invda    
    952       REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin 
     973      REAL(wp) ::   zstemp11r, zstemp12r, zstemp22r    
     974      REAL(wp) ::   zstemp11s, zstemp12s, zstemp22s  
     975      REAL(wp) ::   za22, zQd11Qd11, zQd11Qd12, zQd12Qd12  
     976      REAL(wp) ::   zQ11Q11, zQ11Q12, zQ12Q12  
     977      REAL(wp) ::   zdtemp11, zdtemp12, zdtemp22  
     978      REAL(wp) ::   zrotstemp11r, zrotstemp12r, zrotstemp22r    
     979      REAL(wp) ::   zrotstemp11s, zrotstemp12s, zrotstemp22s 
     980      REAL(wp) ::   zsig11, zsig12, zsig22  
     981      REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22  
     982      REAL(wp) ::   zinvstressconviso  
     983      REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha  
     984      REAL(wp) ::   zTany_1, zTany_2  
     985      REAL(wp) ::   zx, zy, zdx, zdy, zda, zkxw, kyw, kaw 
     986      REAL(wp) ::   zinvdx, zinvdy, zinvda    
     987      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime, zinvsin 
    953988 
    954989      REAL(wp), PARAMETER ::   kfriction = 0.45_wp 
     
    956991      ! Factor to maintain the same stress as in EVP (see Section 3) 
    957992      ! Can be set to 1 otherwise 
    958       invstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
     993      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
    959994  
    960       invsin = 1._wp/sin(2._wp*phi) * invstressconviso  
     995      zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso  
    961996      !now uses phi as set in higher code 
    962997        
     
    9651000 
    9661001      ! 1) structure tensor 
    967       a22 = 1._wp-a11 
    968       Q11Q11 = 1._wp 
    969       Q12Q12 = rsmall 
    970       Q11Q12 = rsmall 
     1002      za22 = 1._wp-pa11 
     1003      zQ11Q11 = 1._wp 
     1004      zQ12Q12 = rsmall 
     1005      zQ11Q12 = rsmall 
    9711006  
    9721007      ! gamma: angle between general coordiantes and principal axis of A  
    9731008      ! here Tan2gamma = 2 a12 / (a11 - a22)  
    974       if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then       
    975          Angle_denom_gamma = 1._wp/sqrt( ( a11 - a22 )*( a11 - a22) + & 
    976                              4._wp*a12*a12 ) 
     1009      IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN       
     1010         zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + & 
     1011                             4._wp*pa12*pa12 ) 
    9771012    
    978          Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Cos^2  
    979          Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Sin^2 
    980          Q11Q12 = a12*Angle_denom_gamma                     !CosSin  
    981       endif 
     1013         zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2  
     1014         zQ12Q12 = 0.5_wp - ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2 
     1015         zQ11Q12 = pa12*zAngle_denom_gamma                     !CosSin  
     1016      ENDIF 
    9821017  
    9831018      ! rotation Q*atemp*Q^T 
    984       atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22  
     1019      zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22  
    9851020       
    9861021      ! make first principal value the largest 
    987       atempprime = max(atempprime, 1.0_wp - atempprime) 
     1022      zatempprime = max(zatempprime, 1.0_wp - zatempprime) 
    9881023  
    9891024      ! 2) strain rate 
    990       dtemp11 = 0.5_wp*(divu + tension) 
    991       dtemp12 = shear*0.5_wp 
    992       dtemp22 = 0.5_wp*(divu - tension) 
    993  
    994       !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear 
     1025      zdtemp11 = 0.5_wp*(pdivu + ptension) 
     1026      zdtemp12 = pshear*0.5_wp 
     1027      zdtemp22 = 0.5_wp*(pdivu - ptension) 
    9951028 
    9961029      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)  
    9971030 
    998       Qd11Qd11 = 1.0_wp 
    999       Qd12Qd12 = rsmall 
    1000       Qd11Qd12 = rsmall 
    1001  
    1002       if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then 
    1003  
    1004          Angle_denom_alpha = 1.0_wp/sqrt( ( dtemp11 - dtemp22 )* & 
    1005                              ( dtemp11 - dtemp22 ) + 4.0_wp*dtemp12*dtemp12) 
    1006  
    1007          Qd11Qd11 = 0.5_wp + ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Cos^2  
    1008          Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^2 
    1009          Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin 
    1010       endif 
    1011  
    1012       dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp22 
    1013       dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp22 
     1031      zQd11Qd11 = 1.0_wp 
     1032      zQd12Qd12 = rsmall 
     1033      zQd11Qd12 = rsmall 
     1034 
     1035      IF((ABS( zdtemp11 - zdtemp22) > rsmall).OR. (ABS(zdtemp12) > rsmall)) THEN 
     1036 
     1037         zAngle_denom_alpha = 1.0_wp/sqrt( ( zdtemp11 - zdtemp22 )* & 
     1038                             ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12) 
     1039 
     1040         zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2  
     1041         zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2 
     1042         zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin 
     1043      ENDIF 
     1044 
     1045      zdtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22 
     1046      zdtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22 
    10141047      ! In cos and sin values 
    1015       x = 0._wp 
    1016       if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then 
    1017          x = atan2(dtemp2,dtemp1) 
    1018       endif       
     1048      zx = 0._wp 
     1049      IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN 
     1050         zx = atan2(zdtemp2,zdtemp1) 
     1051      ENDIF       
    10191052                
    10201053      ! to ensure the angle lies between pi/4 and 9 pi/4  
    1021       if (x < rpi*0.25_wp) x = x + rpi*2.0_wp 
     1054      IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp 
    10221055       
    10231056      ! y: angle between major principal axis of strain rate and structure 
     
    10271060      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 
    10281061       
    1029       Tany_1 = Q11Q12 - Qd11Qd12 
    1030       Tany_2 = Q11Q11 - Qd12Qd12 
     1062      zTany_1 = zQ11Q12 - zQd11Qd12 
     1063      zTany_2 = zQ11Q11 - zQd12Qd12 
    10311064       
    1032       y = 0._wp 
     1065      zy = 0._wp 
    10331066       
    1034       if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then 
    1035          y = atan2(Tany_1,Tany_2) 
    1036       endif 
     1067      IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN 
     1068         zy = atan2(zTany_1,zTany_2) 
     1069      ENDIF 
    10371070       
    10381071      ! to make sure y is between 0 and pi 
    1039       if (y > rpi) y = y - rpi 
    1040       if (y < 0)  y = y + rpi 
     1072      IF (zy > rpi) zy = zy - rpi 
     1073      IF (zy < 0)  zy = zy + rpi 
    10411074 
    10421075      ! 3) update anisotropic stress tensor  
    1043       dx   = rpi/real(nx_yield-1,kind=wp) 
    1044       dy   = rpi/real(ny_yield-1,kind=wp) 
    1045       da   = 0.5_wp/real(na_yield-1,kind=wp) 
    1046       invdx = 1._wp/dx 
    1047       invdy = 1._wp/dy 
    1048       invda = 1._wp/da 
     1076      zdx   = rpi/real(nx_yield-1,kind=wp) 
     1077      zdy   = rpi/real(ny_yield-1,kind=wp) 
     1078      zda   = 0.5_wp/real(na_yield-1,kind=wp) 
     1079      zinvdx = 1._wp/zdx 
     1080      zinvdy = 1._wp/zdy 
     1081      zinvda = 1._wp/zda 
    10491082 
    10501083      ! % need 8 coords and 8 weights 
    10511084      ! % range in kx 
    1052       kx  = int((x-rpi*0.25_wp-rpi)*invdx) + 1 
    1053       kxw = kx - (x-rpi*0.25_wp-rpi)*invdx  
     1085      kx  = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 
     1086      zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx  
    10541087       
    1055       ky  = int(y*invdy) + 1 
    1056       kyw = ky - y*invdy  
     1088      ky  = int(zy*zinvdy) + 1 
     1089      kyw = ky - zy*zinvdy  
    10571090       
    1058       ka  = int((atempprime-0.5_wp)*invda) + 1 
    1059       kaw = ka - (atempprime-0.5_wp)*invda 
     1091      ka  = int((zatempprime-0.5_wp)*zinvda) + 1 
     1092      kaw = ka - (zatempprime-0.5_wp)*zinvda 
    10601093       
    10611094      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 
    1062       stemp11r =  kxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) & 
    1063         & + (1._wp-kxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) & 
    1064         & + kxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) & 
    1065         & + kxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) & 
    1066         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) & 
    1067         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) & 
    1068         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) & 
    1069         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 
    1070       stemp12r =  kxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) & 
    1071         & + (1._wp-kxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) & 
    1072         & + kxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) & 
    1073         & + kxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) & 
    1074         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) & 
    1075         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) & 
    1076         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) & 
    1077         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 
    1078       stemp22r = kxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) & 
    1079         & + (1._wp-kxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) & 
    1080         & + kxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) & 
    1081         & + kxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) & 
    1082         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) & 
    1083         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) & 
    1084         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) & 
    1085         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 
     1095      zstemp11r =  zkxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) & 
     1096        & + (1._wp-zkxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) & 
     1097        & + zkxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) & 
     1098        & + zkxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) & 
     1099        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) & 
     1100        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) & 
     1101        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) & 
     1102        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 
     1103      zstemp12r =  zkxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) & 
     1104        & + (1._wp-zkxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) & 
     1105        & + zkxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) & 
     1106        & + zkxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) & 
     1107        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) & 
     1108        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) & 
     1109        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) & 
     1110        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 
     1111      zstemp22r = zkxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) & 
     1112        & + (1._wp-zkxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) & 
     1113        & + zkxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) & 
     1114        & + zkxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) & 
     1115        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) & 
     1116        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) & 
     1117        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) & 
     1118        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 
    10861119       
    1087       stemp11s =  kxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) & 
    1088         & + (1._wp-kxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) & 
    1089         & + kxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) & 
    1090         & + kxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) & 
    1091         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) & 
    1092         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) & 
    1093         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) & 
    1094         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 
    1095       stemp12s =  kxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) & 
    1096         & + (1._wp-kxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) & 
    1097         & + kxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) & 
    1098         & + kxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) & 
    1099         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) & 
    1100         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) & 
    1101         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) & 
    1102         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 
    1103       stemp22s =  kxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) & 
    1104         & + (1._wp-kxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) & 
    1105         & + kxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) & 
    1106         & + kxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) & 
    1107         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) & 
    1108         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) & 
    1109         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) & 
    1110         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 
     1120      zstemp11s =  zkxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) & 
     1121        & + (1._wp-zkxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) & 
     1122        & + zkxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) & 
     1123        & + zkxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) & 
     1124        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) & 
     1125        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) & 
     1126        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) & 
     1127        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 
     1128      zstemp12s =  zkxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) & 
     1129        & + (1._wp-zkxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) & 
     1130        & + zkxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) & 
     1131        & + zkxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) & 
     1132        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) & 
     1133        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) & 
     1134        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) & 
     1135        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 
     1136      zstemp22s =  zkxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) & 
     1137        & + (1._wp-zkxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) & 
     1138        & + zkxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) & 
     1139        & + zkxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) & 
     1140        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) & 
     1141        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) & 
     1142        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) & 
     1143        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 
    11111144                    
    11121145      ! Calculate mean ice stress over a collection of floes (Equation 3 in 
    11131146      ! Tsamados 2013) 
    11141147 
    1115       sig11  = strength*(stemp11r + kfriction*stemp11s) * invsin 
    1116       sig12  = strength*(stemp12r + kfriction*stemp12s) * invsin 
    1117       sig22  = strength*(stemp22r + kfriction*stemp22s) * invsin 
    1118  
    1119       !WRITE(numout,*) 'strength inside loop', strength 
     1148      zsig11  = strength*(zstemp11r + kfriction*zstemp11s) * zinvsin 
     1149      zsig12  = strength*(zstemp12r + kfriction*zstemp12s) * zinvsin 
     1150      zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin 
     1151 
    11201152      !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22 
    11211153 
     
    11231155 
    11241156      ! Update stress 
    1125       sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 -        2._wp*Q11Q12 *sig12 
    1126       sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12 
    1127       sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 +        2._wp*Q11Q12 *sig12 
    1128  
    1129       stressp  = sgprm11 + sgprm22 
    1130       stress12 = sgprm12 
    1131       stressm  = sgprm11 - sgprm22 
    1132  
    1133       !WRITE(numout,*) 'stress output inside loop', ksub, stressp 
     1157      zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 -        2._wp*zQ11Q12 *zsig12 
     1158      zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12 
     1159      zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 +        2._wp*zQ11Q12 *zsig12 
     1160 
     1161      pstressp  = zsgprm11 + zsgprm22 
     1162      pstress12 = zsgprm12 
     1163      pstressm  = zsgprm11 - zsgprm22 
    11341164 
    11351165      ! Compute ridging and sliding functions in general coordinates  
    11361166      ! (Equation 11 in Tsamados 2013) 
    1137       if (ksub == ndte) then 
    1138          rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12* stemp12r & 
    1139                      + Q12Q12*stemp22r 
    1140          rotstemp12r = Q11Q11*stemp12r +    Q11Q12*(stemp11r-stemp22r) & 
    1141                      - Q12Q12*stemp12r 
    1142          rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12* stemp12r &  
    1143                      + Q11Q11*stemp22r 
    1144  
    1145          rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12* stemp12s & 
    1146                      + Q12Q12*stemp22s 
    1147          rotstemp12s = Q11Q11*stemp12s +    Q11Q12*(stemp11s-stemp22s) & 
    1148                      - Q12Q12*stemp12s 
    1149          rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12* stemp12s &  
    1150                      + Q11Q11*stemp22s 
    1151  
    1152          alphar =  rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 & 
    1153                  + rotstemp22r*dtemp22 
    1154          alphas =  rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 & 
    1155                  + rotstemp22s*dtemp22 
    1156       endif 
     1167      IF (ksub == kndte) THEN 
     1168         zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r & 
     1169                     + zQ12Q12*zstemp22r 
     1170         zrotstemp12r = zQ11Q11*zstemp12r +    zQ11Q12*(zstemp11r-zstemp22r) & 
     1171                     - zQ12Q12*zstemp12r 
     1172         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r &  
     1173                     + zQ11Q11*zstemp22r 
     1174 
     1175         zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s & 
     1176                     + zQ12Q12*zstemp22s 
     1177         zrotstemp12s = zQ11Q11*zstemp12s +    zQ11Q12*(zstemp11s-zstemp22s) & 
     1178                     - zQ12Q12*zstemp12s 
     1179         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s &  
     1180                     + zQ11Q11*zstemp22s 
     1181 
     1182         palphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 
     1183                 + zrotstemp22r*zdtemp22 
     1184         palphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 
     1185                 + zrotstemp22s*zdtemp22 
     1186      ENDIF 
    11571187   END SUBROUTINE update_stress_rdg  
     1188 
     1189!======================================================================= 
     1190 
     1191 
     1192   SUBROUTINE calc_ffrac (pstressp, pstressm, pstress12, & 
     1193                          pa11, pa12,                   & 
     1194                          pmresult11, pmresult12) 
     1195      !!--------------------------------------------------------------------- 
     1196      !!                   ***  ROUTINE calc_ffrac  *** 
     1197      !!                      
     1198      !! ** Purpose :   Computes term in evolution equation for structure tensor 
     1199      !!                which determines the ice floe re-orientation due to fracture 
     1200      !! ** Method :    Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2 
     1201      !!--------------------------------------------------------------------- 
     1202      REAL (wp), INTENT(in)  ::   pstressp, pstressm, pstress12, pa11, pa12 
     1203      REAL (wp), INTENT(out) ::   pmresult11, pmresult12 
     1204 
     1205      ! local variables 
     1206 
     1207      REAL (wp) ::   zsigma11, zsigma12, zsigma22  ! stress tensor elements 
     1208      REAL (wp) ::   zAngle_denom        ! angle with principal component axis 
     1209      REAL (wp) ::   zsigma_1, zsigma_2   ! principal components of stress 
     1210      REAL (wp) ::   zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 
     1211 
     1212      REAL (wp), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation  
     1213      REAL (wp), PARAMETER ::   threshold = 0.3_wp  ! critical confinement ratio  
     1214      !!--------------------------------------------------------------- 
     1215      ! 
     1216      zsigma11 = 0.5_wp*(pstressp+pstressm)  
     1217      zsigma12 = pstress12  
     1218      zsigma22 = 0.5_wp*(pstressp-pstressm)  
     1219 
     1220      ! Here's the change - no longer calculate gamma, 
     1221      ! use trig formulation, angle signs are kept correct, don't worry 
     1222    
     1223      ! rotate tensor to get into sigma principal axis 
     1224    
     1225      ! here Tan2gamma = 2 sig12 / (sig11 - sig12)  
     1226      ! add rsmall to the denominator to stop 1/0 errors, makes very little  
     1227      ! error to the calculated angles < rsmall 
     1228 
     1229      zQ11Q11 = 0.1_wp 
     1230      zQ12Q12 = rsmall 
     1231      zQ11Q12 = rsmall 
     1232 
     1233      IF((ABS( zsigma11 - zsigma22) > rsmall).OR.(ABS(zsigma12) > rsmall)) THEN 
     1234 
     1235         zAngle_denom = 1.0_wp/sqrt( ( zsigma11 - zsigma22 )*( zsigma11 - & 
     1236                       zsigma22 ) + 4.0_wp*zsigma12*zsigma12) 
     1237 
     1238         zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2  
     1239         zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Sin^2 
     1240         zQ11Q12 = zsigma12*zAngle_denom                      !CosSin  
     1241      ENDIF 
     1242 
     1243      zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12  & 
     1244              + zQ12Q12*zsigma22 ! S(1,1) 
     1245      zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12  & 
     1246              + zQ11Q11*zsigma22 ! S(2,2) 
     1247 
     1248      ! Pure divergence 
     1249      IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp))  THEN 
     1250         pmresult11 = 0.0_wp 
     1251         pmresult12 = 0.0_wp 
     1252 
     1253      ! Unconfined compression: cracking of blocks not along the axial splitting 
     1254      ! direction 
     1255      ! which leads to the loss of their shape, so we again model it through diffusion 
     1256      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN 
     1257         pmresult11 = kfrac * (pa11 - zQ12Q12) 
     1258         pmresult12 = kfrac * (pa12 + zQ11Q12) 
     1259 
     1260      ! Shear faulting 
     1261      ELSEIF (zsigma_2 == 0.0_wp) THEN 
     1262         pmresult11 = 0.0_wp 
     1263         pmresult12 = 0.0_wp 
     1264      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 
     1265         pmresult11 = kfrac * (pa11 - zQ12Q12) 
     1266         pmresult12 = kfrac * (pa12 + zQ11Q12) 
     1267 
     1268      ! Horizontal spalling 
     1269      ELSE  
     1270         pmresult11 = 0.0_wp 
     1271         pmresult12 = 0.0_wp 
     1272      ENDIF 
     1273 
     1274   END SUBROUTINE calc_ffrac 
     1275 
    11581276 
    11591277   SUBROUTINE rhg_eap_rst( cdrw, kt ) 
     
    12341352         s12s(ix,iy,ia) = 0._wp 
    12351353         s22s(ix,iy,ia) = 0._wp 
    1236          IF (ia <= na_yield-1) then 
     1354         IF (ia <= na_yield-1) THEN 
    12371355          DO iz=1,nz 
    12381356          s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    12391357           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1240            s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1358           s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    12411359          s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    12421360           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1243            s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1361           s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    12441362          s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    12451363           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1246            s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)  
     1364           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)  
    12471365          s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    12481366           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1249            s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1367           s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    12501368          s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    12511369           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1252            s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1370           s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    12531371          s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    12541372           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1255            s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1373           s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    12561374          ENDDO 
    12571375          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 
     
    12621380          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 
    12631381         ELSE 
    1264           s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1265           s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1266           s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1267           s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1268           s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1269           s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     1382          s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1383          s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1384          s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1385          s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1386          s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1387          s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
    12701388          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 
    12711389          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 
     
    12951413 
    12961414 
    1297    FUNCTION w1(a) 
     1415   FUNCTION w1(pa) 
    12981416      !!------------------------------------------------------------------- 
    12991417      !! Function : w1 (see Gaussian function psi in Tsamados et al 2013) 
    13001418      !!------------------------------------------------------------------- 
    1301       REAL(wp), INTENT(in   ) ::   a 
     1419      REAL(wp), INTENT(in   ) ::   pa 
    13021420      REAL(wp) ::   w1 
    13031421      !!------------------------------------------------------------------- 
    13041422    
    13051423      w1 = - 223.87569446_wp & 
    1306            + 2361.2198663_wp*a & 
    1307            - 10606.56079975_wp*a*a & 
    1308            + 26315.50025642_wp*a*a*a & 
    1309            - 38948.30444297_wp*a*a*a*a & 
    1310            + 34397.72407466_wp*a*a*a*a*a & 
    1311            - 16789.98003081_wp*a*a*a*a*a*a & 
    1312            + 3495.82839237_wp*a*a*a*a*a*a*a 
     1424       &   + 2361.2198663_wp*pa & 
     1425       &   - 10606.56079975_wp*pa*pa & 
     1426       &   + 26315.50025642_wp*pa*pa*pa & 
     1427       &   - 38948.30444297_wp*pa*pa*pa*pa & 
     1428       &   + 34397.72407466_wp*pa*pa*pa*pa*pa & 
     1429       &   - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 
     1430       &   + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 
    13131431    
    13141432   END FUNCTION w1 
    13151433 
    13161434 
    1317    FUNCTION w2(a) 
     1435   FUNCTION w2(pa) 
    13181436      !!------------------------------------------------------------------- 
    13191437      !! Function : w2 (see Gaussian function psi in Tsamados et al 2013) 
    13201438      !!------------------------------------------------------------------- 
    1321       REAL(wp), INTENT(in   ) ::   a 
     1439      REAL(wp), INTENT(in   ) ::   pa 
    13221440      REAL(wp) ::   w2 
    13231441      !!------------------------------------------------------------------- 
    13241442    
    13251443      w2 = - 6670.68911883_wp & 
    1326            + 70222.33061536_wp*a & 
    1327            - 314871.71525448_wp*a*a & 
    1328            + 779570.02793492_wp*a*a*a & 
    1329            - 1151098.82436864_wp*a*a*a*a & 
    1330            + 1013896.59464498_wp*a*a*a*a*a & 
    1331            - 493379.44906738_wp*a*a*a*a*a*a & 
    1332            + 102356.551518_wp*a*a*a*a*a*a*a 
     1444       &   + 70222.33061536_wp*pa & 
     1445       &   - 314871.71525448_wp*pa*pa & 
     1446       &   + 779570.02793492_wp*pa*pa*pa & 
     1447       &   - 1151098.82436864_wp*pa*pa*pa*pa & 
     1448       &   + 1013896.59464498_wp*pa*pa*pa*pa*pa & 
     1449       &   - 493379.44906738_wp*pa*pa*pa*pa*pa*pa & 
     1450       &   + 102356.551518_wp*pa*pa*pa*pa*pa*pa*pa 
    13331451    
    13341452   END FUNCTION w2 
    13351453 
    1336    FUNCTION s11kr(x,y,z,phi)  
     1454   FUNCTION s11kr(px,py,pz)  
    13371455      !!------------------------------------------------------------------- 
    13381456      !! Function : s11kr 
    13391457      !!------------------------------------------------------------------- 
    1340       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
     1458      REAL(wp), INTENT(in   ) ::   px,py,pz 
    13411459    
    1342       REAL(wp) ::   s11kr, p, pih 
     1460      REAL(wp) ::   s11kr, zpih 
    13431461    
    1344       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1345       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1346       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1347       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1348       REAL(wp) ::   d11, d12, d22 
    1349       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1350       REAL(wp) ::   Hen1t2, Hen2t1 
     1462      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1463      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1464      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1465      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1466      REAL(wp) ::   zd11, zd12, zd22 
     1467      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1468      REAL(wp) ::   zHen1t2, zHen2t1 
    13511469      !!------------------------------------------------------------------- 
    13521470    
    1353       pih = 0.5_wp*rpi 
    1354       p = phi 
     1471      zpih = 0.5_wp*rpi 
    13551472    
    1356       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1357       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1358       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1359       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1360       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1361       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1362       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1363       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1364       t1t2i11 = cos(z-p) * cos(z+p) 
    1365       t1t2i12 = cos(z-p) * sin(z+p) 
    1366       t1t2i21 = sin(z-p) * cos(z+p) 
    1367       t1t2i22 = sin(z-p) * sin(z+p) 
    1368       t2t1i11 = cos(z+p) * cos(z-p) 
    1369       t2t1i12 = cos(z+p) * sin(z-p) 
    1370       t2t1i21 = sin(z+p) * cos(z-p) 
    1371       t2t1i22 = sin(z+p) * sin(z-p) 
     1473      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1474      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1475      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1476      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1477      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1478      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1479      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1480      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1481      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1482      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1483      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1484      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1485      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1486      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1487      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1488      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    13721489   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 
    13731490   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 
    13741491   ! x=x-pi/2 
    1375       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1376       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1377       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1378       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1379       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1380       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1492      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1493      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1494      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1495      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1496      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1497      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    13811498    
    1382       IF (-IIn1t2>=rsmall) then 
    1383       Hen1t2 = 1._wp 
     1499      IF (-zIIn1t2>=rsmall) THEN 
     1500      zHen1t2 = 1._wp 
    13841501      ELSE 
    1385       Hen1t2 = 0._wp 
     1502      zHen1t2 = 0._wp 
    13861503      ENDIF 
    13871504    
    1388       IF (-IIn2t1>=rsmall) then 
    1389       Hen2t1 = 1._wp 
     1505      IF (-zIIn2t1>=rsmall) THEN 
     1506      zHen2t1 = 1._wp 
    13901507      ELSE 
    1391       Hen2t1 = 0._wp 
     1508      zHen2t1 = 0._wp 
    13921509      ENDIF 
    13931510    
    1394       s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 * n2t1i11) 
     1511      s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 
    13951512 
    13961513   END FUNCTION s11kr 
    13971514 
    1398    FUNCTION s12kr(x,y,z,phi) 
     1515   FUNCTION s12kr(px,py,pz) 
    13991516      !!------------------------------------------------------------------- 
    14001517      !! Function : s12kr 
    14011518      !!------------------------------------------------------------------- 
    1402       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1403  
    1404       REAL(wp) ::   s12kr, s12r0, s21r0, p, pih 
    1405  
    1406       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1407       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1408       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1409       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1410       REAL(wp) ::   d11, d12, d22 
    1411       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1412       REAL(wp) ::   Hen1t2, Hen2t1 
    1413       !!------------------------------------------------------------------- 
    1414       pih = 0.5_wp*rpi 
    1415       p = phi 
     1519      REAL(wp), INTENT(in   ) ::   px,py,pz 
     1520 
     1521      REAL(wp) ::   s12kr, zs12r0, zs21r0, zpih 
     1522 
     1523      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1524      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1525      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1526      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1527      REAL(wp) ::   zd11, zd12, zd22 
     1528      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1529      REAL(wp) ::   zHen1t2, zHen2t1 
     1530      !!------------------------------------------------------------------- 
     1531      zpih = 0.5_wp*rpi 
    14161532    
    1417       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1418       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1419       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1420       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1421       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1422       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1423       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1424       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1425       t1t2i11 = cos(z-p) * cos(z+p) 
    1426       t1t2i12 = cos(z-p) * sin(z+p) 
    1427       t1t2i21 = sin(z-p) * cos(z+p) 
    1428       t1t2i22 = sin(z-p) * sin(z+p) 
    1429       t2t1i11 = cos(z+p) * cos(z-p) 
    1430       t2t1i12 = cos(z+p) * sin(z-p) 
    1431       t2t1i21 = sin(z+p) * cos(z-p) 
    1432       t2t1i22 = sin(z+p) * sin(z-p) 
    1433       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1434       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1435       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1436       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1437       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1438       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1533      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1534      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1535      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1536      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1537      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1538      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1539      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1540      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1541      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1542      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1543      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1544      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1545      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1546      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1547      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1548      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
     1549      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1550      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1551      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1552      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1553      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1554      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
    14391555    
    1440       IF (-IIn1t2>=rsmall) then 
    1441       Hen1t2 = 1._wp 
     1556      IF (-zIIn1t2>=rsmall) THEN 
     1557      zHen1t2 = 1._wp 
    14421558      ELSE 
    1443       Hen1t2 = 0._wp 
     1559      zHen1t2 = 0._wp 
    14441560      ENDIF 
    14451561    
    1446       IF (-IIn2t1>=rsmall) then 
    1447       Hen2t1 = 1._wp 
     1562      IF (-zIIn2t1>=rsmall) THEN 
     1563      zHen2t1 = 1._wp 
    14481564      ELSE 
    1449       Hen2t1 = 0._wp 
     1565      zHen2t1 = 0._wp 
    14501566      ENDIF 
    14511567    
    1452       s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 * n2t1i12) 
    1453       s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 * n2t1i21) 
    1454       s12kr=0.5_wp*(s12r0+s21r0) 
     1568      zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 
     1569      zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21) 
     1570      s12kr=0.5_wp*(zs12r0+zs21r0) 
    14551571    
    14561572   END FUNCTION s12kr 
    14571573 
    1458    FUNCTION s22kr(x,y,z,phi) 
     1574   FUNCTION s22kr(px,py,pz) 
    14591575      !!------------------------------------------------------------------- 
    14601576      !! Function : s22kr 
    14611577      !!------------------------------------------------------------------- 
    1462       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1463  
    1464       REAL(wp) ::   s22kr, p, pih 
    1465  
    1466       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1467       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1468       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1469       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1470       REAL(wp) ::   d11, d12, d22 
    1471       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1472       REAL(wp) ::   Hen1t2, Hen2t1 
    1473       !!------------------------------------------------------------------- 
    1474       pih = 0.5_wp*rpi 
    1475       p = phi 
    1476  
    1477       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1478       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1479       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1480       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1481       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1482       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1483       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1484       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1485       t1t2i11 = cos(z-p) * cos(z+p) 
    1486       t1t2i12 = cos(z-p) * sin(z+p) 
    1487       t1t2i21 = sin(z-p) * cos(z+p) 
    1488       t1t2i22 = sin(z-p) * sin(z+p) 
    1489       t2t1i11 = cos(z+p) * cos(z-p) 
    1490       t2t1i12 = cos(z+p) * sin(z-p) 
    1491       t2t1i21 = sin(z+p) * cos(z-p) 
    1492       t2t1i22 = sin(z+p) * sin(z-p) 
    1493       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1494       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1495       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1496       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1497       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1498       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1499  
    1500       IF (-IIn1t2>=rsmall) then 
    1501       Hen1t2 = 1._wp 
     1578      REAL(wp), INTENT(in   ) ::   px,py,pz 
     1579 
     1580      REAL(wp) ::   s22kr, zpih 
     1581 
     1582      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1583      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1584      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1585      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1586      REAL(wp) ::   zd11, zd12, zd22 
     1587      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1588      REAL(wp) ::   zHen1t2, zHen2t1 
     1589      !!------------------------------------------------------------------- 
     1590      zpih = 0.5_wp*rpi 
     1591 
     1592      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1593      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1594      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1595      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1596      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1597      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1598      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1599      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1600      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1601      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1602      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1603      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1604      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1605      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1606      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1607      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
     1608      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1609      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1610      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1611      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1612      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1613      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1614 
     1615      IF (-zIIn1t2>=rsmall) THEN 
     1616      zHen1t2 = 1._wp 
    15021617      ELSE 
    1503       Hen1t2 = 0._wp 
    1504       ENDIF 
    1505  
    1506       IF (-IIn2t1>=rsmall) then 
    1507       Hen2t1 = 1._wp 
     1618      zHen1t2 = 0._wp 
     1619      ENDIF 
     1620 
     1621      IF (-zIIn2t1>=rsmall) THEN 
     1622      zHen2t1 = 1._wp 
    15081623      ELSE 
    1509       Hen2t1 = 0._wp 
    1510       ENDIF 
    1511  
    1512       s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 * n2t1i22) 
     1624      zHen2t1 = 0._wp 
     1625      ENDIF 
     1626 
     1627      s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 
    15131628 
    15141629   END FUNCTION s22kr 
    15151630 
    1516    FUNCTION s11ks(x,y,z,phi) 
     1631   FUNCTION s11ks(px,py,pz) 
    15171632      !!------------------------------------------------------------------- 
    15181633      !! Function : s11ks 
    15191634      !!------------------------------------------------------------------- 
    1520       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1521  
    1522       REAL(wp) ::   s11ks, p, pih 
    1523  
    1524       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1525       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1526       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1527       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1528       REAL(wp) ::   d11, d12, d22 
    1529       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1530       REAL(wp) ::   Hen1t2, Hen2t1 
    1531       !!------------------------------------------------------------------- 
    1532       pih = 0.5_wp*rpi 
    1533       p = phi 
    1534  
    1535       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1536       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1537       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1538       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1539       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1540       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1541       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1542       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1543       t1t2i11 = cos(z-p) * cos(z+p) 
    1544       t1t2i12 = cos(z-p) * sin(z+p) 
    1545       t1t2i21 = sin(z-p) * cos(z+p) 
    1546       t1t2i22 = sin(z-p) * sin(z+p) 
    1547       t2t1i11 = cos(z+p) * cos(z-p) 
    1548       t2t1i12 = cos(z+p) * sin(z-p) 
    1549       t2t1i21 = sin(z+p) * cos(z-p) 
    1550       t2t1i22 = sin(z+p) * sin(z-p) 
    1551       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1552       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1553       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1554       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1555       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1556       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1557  
    1558       IF (-IIn1t2>=rsmall) then 
    1559       Hen1t2 = 1._wp 
     1635      REAL(wp), INTENT(in   ) ::   px,py,pz 
     1636 
     1637      REAL(wp) ::   s11ks, zpih 
     1638 
     1639      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1640      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1641      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1642      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1643      REAL(wp) ::   zd11, zd12, zd22 
     1644      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1645      REAL(wp) ::   zHen1t2, zHen2t1 
     1646      !!------------------------------------------------------------------- 
     1647      zpih = 0.5_wp*rpi 
     1648 
     1649      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1650      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1651      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1652      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1653      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1654      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1655      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1656      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1657      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1658      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1659      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1660      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1661      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1662      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1663      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1664      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
     1665      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1666      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1667      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1668      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1669      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1670      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1671 
     1672      IF (-zIIn1t2>=rsmall) THEN 
     1673      zHen1t2 = 1._wp 
    15601674      ELSE 
    1561       Hen1t2 = 0._wp 
    1562       ENDIF 
    1563  
    1564       IF (-IIn2t1>=rsmall) then 
    1565       Hen2t1 = 1._wp 
     1675      zHen1t2 = 0._wp 
     1676      ENDIF 
     1677 
     1678      IF (-zIIn2t1>=rsmall) THEN 
     1679      zHen2t1 = 1._wp 
    15661680      ELSE 
    1567       Hen2t1 = 0._wp 
    1568       ENDIF 
    1569  
    1570       s11ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i11 + Hen2t1 * t2t1i11) 
     1681      zHen2t1 = 0._wp 
     1682      ENDIF 
     1683 
     1684      s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 
    15711685 
    15721686   END FUNCTION s11ks 
    15731687 
    1574    FUNCTION s12ks(x,y,z,phi) 
     1688   FUNCTION s12ks(px,py,pz) 
    15751689      !!------------------------------------------------------------------- 
    15761690      !! Function : s12ks 
    15771691      !!------------------------------------------------------------------- 
    1578       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1579  
    1580       REAL(wp) ::   s12ks, s12s0, s21s0, p, pih 
    1581  
    1582       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1583       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1584       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1585       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1586       REAL(wp) ::   d11, d12, d22 
    1587       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1588       REAL(wp) ::   Hen1t2, Hen2t1 
    1589       !!------------------------------------------------------------------- 
    1590       pih = 0.5_wp*rpi 
    1591       p =phi 
    1592  
    1593       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1594       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1595       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1596       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1597       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1598       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1599       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1600       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1601       t1t2i11 = cos(z-p) * cos(z+p) 
    1602       t1t2i12 = cos(z-p) * sin(z+p) 
    1603       t1t2i21 = sin(z-p) * cos(z+p) 
    1604       t1t2i22 = sin(z-p) * sin(z+p) 
    1605       t2t1i11 = cos(z+p) * cos(z-p) 
    1606       t2t1i12 = cos(z+p) * sin(z-p) 
    1607       t2t1i21 = sin(z+p) * cos(z-p) 
    1608       t2t1i22 = sin(z+p) * sin(z-p) 
    1609       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1610       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1611       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1612       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1613       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1614       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1615  
    1616       IF (-IIn1t2>=rsmall) then 
    1617       Hen1t2 = 1._wp 
     1692      REAL(wp), INTENT(in   ) ::   px,py,pz 
     1693 
     1694      REAL(wp) ::   s12ks, zs12s0, zs21s0, zpih 
     1695 
     1696      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1697      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1698      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1699      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1700      REAL(wp) ::   zd11, zd12, zd22 
     1701      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1702      REAL(wp) ::   zHen1t2, zHen2t1 
     1703      !!------------------------------------------------------------------- 
     1704      zpih = 0.5_wp*rpi 
     1705 
     1706      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1707      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1708      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1709      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1710      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1711      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1712      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1713      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1714      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1715      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1716      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1717      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1718      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1719      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1720      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1721      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
     1722      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1723      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1724      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1725      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1726      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1727      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1728 
     1729      IF (-zIIn1t2>=rsmall) THEN 
     1730      zHen1t2 = 1._wp 
    16181731      ELSE 
    1619       Hen1t2 = 0._wp 
    1620       ENDIF 
    1621  
    1622       IF (-IIn2t1>=rsmall) then 
    1623       Hen2t1 = 1._wp 
     1732      zHen1t2 = 0._wp 
     1733      ENDIF 
     1734 
     1735      IF (-zIIn2t1>=rsmall) THEN 
     1736      zHen2t1 = 1._wp 
    16241737      ELSE 
    1625       Hen2t1 = 0._wp 
    1626       ENDIF 
    1627  
    1628       s12s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i12 + Hen2t1 * t2t1i12) 
    1629       s21s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i21 + Hen2t1 * t2t1i21) 
    1630       s12ks=0.5_wp*(s12s0+s21s0) 
     1738      zHen2t1 = 0._wp 
     1739      ENDIF 
     1740 
     1741      zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 
     1742      zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21) 
     1743      s12ks=0.5_wp*(zs12s0+zs21s0) 
    16311744 
    16321745   END FUNCTION s12ks 
    16331746 
    1634    FUNCTION s22ks(x,y,z,phi)  
     1747   FUNCTION s22ks(px,py,pz)  
    16351748      !!------------------------------------------------------------------- 
    16361749      !! Function : s22ks 
    16371750      !!------------------------------------------------------------------- 
    1638       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1639  
    1640       REAL(wp) ::   s22ks, p, pih 
    1641  
    1642       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1643       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1644       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1645       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1646       REAL(wp) ::   d11, d12, d22 
    1647       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1648       REAL(wp) ::   Hen1t2, Hen2t1 
    1649       !!------------------------------------------------------------------- 
    1650       pih = 0.5_wp*rpi 
    1651       p =phi 
    1652  
    1653       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1654       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1655       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1656       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1657       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1658       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1659       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1660       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1661       t1t2i11 = cos(z-p) * cos(z+p) 
    1662       t1t2i12 = cos(z-p) * sin(z+p) 
    1663       t1t2i21 = sin(z-p) * cos(z+p) 
    1664       t1t2i22 = sin(z-p) * sin(z+p) 
    1665       t2t1i11 = cos(z+p) * cos(z-p) 
    1666       t2t1i12 = cos(z+p) * sin(z-p) 
    1667       t2t1i21 = sin(z+p) * cos(z-p) 
    1668       t2t1i22 = sin(z+p) * sin(z-p) 
    1669       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1670       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1671       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1672       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1673       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1674       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1675  
    1676       IF (-IIn1t2>=rsmall) THEN 
    1677       Hen1t2 = 1._wp 
     1751      REAL(wp), INTENT(in   ) ::   px,py,pz 
     1752 
     1753      REAL(wp) ::   s22ks, zpih 
     1754 
     1755      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1756      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1757      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1758      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1759      REAL(wp) ::   zd11, zd12, zd22 
     1760      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1761      REAL(wp) ::   zHen1t2, zHen2t1 
     1762      !!------------------------------------------------------------------- 
     1763      zpih = 0.5_wp*rpi 
     1764 
     1765      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1766      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1767      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1768      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1769      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1770      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1771      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1772      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1773      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1774      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1775      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1776      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1777      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1778      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1779      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1780      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
     1781      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1782      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1783      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1784      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1785      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1786      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1787 
     1788      IF (-zIIn1t2>=rsmall) THEN 
     1789      zHen1t2 = 1._wp 
    16781790      ELSE 
    1679       Hen1t2 = 0._wp 
    1680       ENDIF 
    1681  
    1682       IF (-IIn2t1>=rsmall) THEN 
    1683       Hen2t1 = 1._wp 
     1791      zHen1t2 = 0._wp 
     1792      ENDIF 
     1793 
     1794      IF (-zIIn2t1>=rsmall) THEN 
     1795      zHen2t1 = 1._wp 
    16841796      ELSE 
    1685       Hen2t1 = 0._wp 
    1686       ENDIF 
    1687  
    1688       s22ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22) 
     1797      zHen2t1 = 0._wp 
     1798      ENDIF 
     1799 
     1800      s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 
    16891801 
    16901802   END FUNCTION s22ks 
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90

    r12263 r13574  
    112112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    113113      !! 
    114       INTEGER ::   ji, jj, ij0, ij1, ii0, ii1        ! dummy loop indices 
     114      INTEGER ::   ji, jj       ! dummy loop indices 
    115115      INTEGER ::   jter         ! local integers 
    116116      ! 
     
    505505!extra code for test case boundary conditions 
    506506                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    507                      v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj)) 
     507                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
    508508                  END IF 
    509509               END DO 
     
    560560!extra code for test case boundary conditions 
    561561                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    562                      u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj)) 
     562                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
    563563                  END IF 
    564564               END DO 
     
    617617!extra code for test case boundary conditions 
    618618                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    619                      u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj)) 
     619                     u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
    620620                  END IF 
    621621               END DO 
     
    672672!extra code for test case boundary conditions 
    673673                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    674                      v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj)) 
     674                     v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
    675675                  END IF 
    676676               END DO 
Note: See TracChangeset for help on using the changeset viewer.