Changeset 13128


Ignore:
Timestamp:
2020-06-18T16:55:01+02:00 (7 weeks ago)
Author:
stefryn
Message:

variable name updates, to be continued

File:
1 edited

Legend:

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

    r13113 r13128  
    142142      REAL(wp) ::   zfac_x, zfac_y 
    143143      REAL(wp) ::   zshear, zdum1, zdum2 
    144       REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components 
    145       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 
    146146      REAL(wp) ::   mresult11, mresult12, z1dtevpkth, p5kth, z1_dtevp_A        ! for structure tensor evolution 
    147147      ! 
    148       REAL(wp), DIMENSION(jpi,jpj) ::   stress12tmp                     ! anisotropic stress tensor component for regridding 
     148      REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                     ! anisotropic stress tensor component for regridding 
    149149      REAL(wp), DIMENSION(jpi,jpj) ::   yield11, yield22, yield12       ! yield surface tensor for history 
    150150      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
     
    415415               CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, & 
    416416                                       zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 
    417                                        stressptmp, stressmtmp, & 
    418                                        stress12tmp(ji,jj), strength(ji,jj), & 
    419                                        alphar, alphas) 
     417                                       zstressptmp, zstressmtmp, & 
     418                                       zstress12tmp(ji,jj), strength(ji,jj), & 
     419                                       zalphar, zalphas) 
    420420 
    421421               ! structure tensor update 
    422422               IF (mod(jter,10) == 0) THEN 
    423                   CALL calc_ffrac(stressptmp, stressmtmp, stress12tmp(ji,jj), & 
     423                  CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 
    424424                                  paniso_11(ji,jj), paniso_12(ji,jj),           & 
    425425                                  mresult11,  mresult12) 
     
    431431 
    432432               IF (jter == nn_nevp) THEN 
    433                   yield11(ji,jj) = 0.5_wp * (stressptmp + stressmtmp) 
    434                   yield22(ji,jj) = 0.5_wp * (stressptmp - stressmtmp) 
    435                   yield12(ji,jj) = stress12tmp(ji,jj) 
    436                   prdg_conv(ji,jj) = -min( alphar, 0._wp)     
     433                  yield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp) 
     434                  yield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp) 
     435                  yield12(ji,jj) = zstress12tmp(ji,jj) 
     436                  prdg_conv(ji,jj) = -min( zalphar, 0._wp)     
    437437               ENDIF 
    438438 
     
    454454               !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
    455455               !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    456                zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1  ) * z1_alph1 
    457                zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1  ) * z1_alph1 
     456               zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1  ) * z1_alph1 
     457               zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1  ) * z1_alph1 
    458458               !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1  ) * z1_alph1 
    459459               !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1  ) * z1_alph1 
     
    462462         END DO 
    463463         !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 
    464          CALL lbc_lnk_multi( 'icedyn_rhg_eap', stress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 
     464         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 
    465465 
    466466         DO jj = 1, jpjm1 
    467467            DO ji = 1, jpim1 
    468468               ! stress12tmp at F points  
    469                stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
    470                   &   + stress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + stress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  & 
     469               zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
     470                  &   + zstress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + zstress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  & 
    471471                  &   ) * 0.25_wp * r1_e1e2f(ji,jj) 
    472472 
     
    479479               ! stress at F points (zkt/=0 if landfast) 
    480480               !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
    481                zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1  ) * z1_alph1 
     481               zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1  ) * z1_alph1 
    482482               !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1 
    483483 
     
    527527                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    528528                  zTauB = ztauy_base(ji,jj) / zvel 
    529                   !                 !--- OceanBottom-to-Ice stress 
     529                  !                 !--- OceanBottom-to-Ice stress  
    530530                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    531531                  ! 
     
    572572                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    573573                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    574                   !                 !--- Ocean-to-Ice stress 
     574                  !                 !--- Ocean-to-Ice stress  
    575575                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    576576                  ! 
     
    625625                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    626626                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    627                   !                 !--- Ocean-to-Ice stress 
     627                  !                 !--- Ocean-to-Ice stress  
    628628                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    629629                  ! 
     
    841841!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    842842               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    843                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
     843               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress  
    844844               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    845845            END DO 
     
    853853         ! Stress tensor invariants (normal and shear stress N/m) 
    854854         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    855          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
     855         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress  
    856856 
    857857         DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
     
    935935                                 stressp,  stressm, & 
    936936                                 stress12, strength, & 
    937                                  alphar, alphas) 
     937                                 palphar, palphas) 
    938938      !!--------------------------------------------------------------------- 
    939939      !!                   ***  SUBROUTINE rhg_eap_rst  *** 
    940940      !!                      
    941       !! ** Purpose :   Updates the stress depending on values of strain rate and structure 
     941      !! ** Purpose :   Updates the depending on values of strain rate and structure 
    942942      !!                tensor and for ksub=ndte it computes closing and sliding rate 
    943943      !!--------------------------------------------------------------------- 
     
    947947      REAL(wp), INTENT(in   ) ::   a11, a12 
    948948      REAL(wp), INTENT(  out) ::   stressp, stressm, stress12 
    949       REAL(wp), INTENT(  out) ::   alphar, alphas 
     949      REAL(wp), INTENT(  out) ::   palphar, palphas 
    950950 
    951951      INTEGER  ::   kx ,ky, ka 
     
    964964      REAL(wp) ::   zTany_1, zTany_2  
    965965      REAL(wp) ::   zx, zy, zdx, zdy, zda, zkxw, kyw, kaw 
    966       REAL(wp) ::   invdx, invdy, invda    
    967       REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin 
     966      REAL(wp) ::   zinvdx, zinvdy, zinvda    
     967      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime, zinvsin 
    968968 
    969969      REAL(wp), PARAMETER ::   kfriction = 0.45_wp 
     
    973973      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
    974974  
    975       invsin = 1._wp/sin(2._wp*phi) * zinvstressconviso  
     975      zinvsin = 1._wp/sin(2._wp*phi) * zinvstressconviso  
    976976      !now uses phi as set in higher code 
    977977        
     
    997997  
    998998      ! rotation Q*atemp*Q^T 
    999       atempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22  
     999      zatempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22  
    10001000       
    10011001      ! make first principal value the largest 
    1002       atempprime = max(atempprime, 1.0_wp - atempprime) 
     1002      zatempprime = max(zatempprime, 1.0_wp - zatempprime) 
    10031003  
    10041004      ! 2) strain rate 
     
    10251025      ENDIF 
    10261026 
    1027       dtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22 
    1028       dtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22 
     1027      zdtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22 
     1028      zdtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22 
    10291029      ! In cos and sin values 
    10301030      zx = 0._wp 
    1031       IF ((ABS(dtemp1) > rsmall).OR.(ABS(dtemp2) > rsmall)) THEN 
    1032          zx = atan2(dtemp2,dtemp1) 
     1031      IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN 
     1032         zx = atan2(zdtemp2,zdtemp1) 
    10331033      ENDIF       
    10341034                
     
    10591059      zdy   = rpi/real(ny_yield-1,kind=wp) 
    10601060      zda   = 0.5_wp/real(na_yield-1,kind=wp) 
    1061       invdx = 1._wp/zdx 
    1062       invdy = 1._wp/zdy 
    1063       invda = 1._wp/zda 
     1061      zinvdx = 1._wp/zdx 
     1062      zinvdy = 1._wp/zdy 
     1063      zinvda = 1._wp/zda 
    10641064 
    10651065      ! % need 8 coords and 8 weights 
    10661066      ! % range in kx 
    1067       kx  = int((zx-rpi*0.25_wp-rpi)*invdx) + 1 
    1068       zkxw = kx - (zx-rpi*0.25_wp-rpi)*invdx  
     1067      kx  = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 
     1068      zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx  
    10691069       
    1070       ky  = int(zy*invdy) + 1 
    1071       kyw = ky - zy*invdy  
     1070      ky  = int(zy*zinvdy) + 1 
     1071      kyw = ky - zy*zinvdy  
    10721072       
    1073       ka  = int((atempprime-0.5_wp)*invda) + 1 
    1074       kaw = ka - (atempprime-0.5_wp)*invda 
     1073      ka  = int((zatempprime-0.5_wp)*zinvda) + 1 
     1074      kaw = ka - (zatempprime-0.5_wp)*zinvda 
    10751075       
    10761076      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 
     
    11281128      ! Tsamados 2013) 
    11291129 
    1130       zsig11  = strength*(zstemp11r + kfriction*zstemp11s) * invsin 
    1131       zsig12  = strength*(zstemp12r + kfriction*zstemp12s) * invsin 
    1132       zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * invsin 
     1130      zsig11  = strength*(zstemp11r + kfriction*zstemp11s) * zinvsin 
     1131      zsig12  = strength*(zstemp12r + kfriction*zstemp12s) * zinvsin 
     1132      zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin 
    11331133 
    11341134      !WRITE(numout,*) 'strength inside loop', strength 
     
    11461146      stressm  = zsgprm11 - zsgprm22 
    11471147 
    1148       !WRITE(numout,*) 'stress output inside loop', ksub, stressp 
     1148      !WRITE(numout,*) ' output inside loop', ksub, stressp 
    11491149 
    11501150      ! Compute ridging and sliding functions in general coordinates  
     
    11651165                     + zQ11Q11*zstemp22s 
    11661166 
    1167          alphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 
     1167         palphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 
    11681168                 + zrotstemp22r*zdtemp22 
    1169          alphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 
     1169         palphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 
    11701170                 + zrotstemp22s*zdtemp22 
    11711171      ENDIF 
     
    11901190      ! local variables 
    11911191 
    1192       REAL (wp) ::   sigma11, sigma12, sigma22  ! stress tensor elements 
    1193       REAL (wp) ::   Angle_denom        ! angle with principal component axis 
    1194       REAL (wp) ::   sigma_1, sigma_2   ! principal components of stress 
    1195       REAL (wp) ::   Q11, Q12, Q11Q11, Q11Q12, Q12Q12 
     1192      REAL (wp) ::   zsigma11, zsigma12, zsigma22  ! stress tensor elements 
     1193      REAL (wp) ::   zAngle_denom        ! angle with principal component axis 
     1194      REAL (wp) ::   zsigma_1, zsigma_2   ! principal components of stress 
     1195      REAL (wp) ::   zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 
    11961196 
    11971197      REAL (wp), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation  
     
    11991199      !!--------------------------------------------------------------- 
    12001200      ! 
    1201       sigma11 = 0.5_wp*(stressp+stressm)  
    1202       sigma12 = stress12  
    1203       sigma22 = 0.5_wp*(stressp-stressm)  
     1201      zsigma11 = 0.5_wp*(stressp+stressm)  
     1202      zsigma12 = stress12  
     1203      zsigma22 = 0.5_wp*(stressp-stressm)  
    12041204 
    12051205      ! Here's the change - no longer calculate gamma, 
     
    12121212      ! error to the calculated angles < rsmall 
    12131213 
    1214       Q11Q11 = 0.1_wp 
    1215       Q12Q12 = rsmall 
    1216       Q11Q12 = rsmall 
    1217  
    1218       IF((ABS( sigma11 - sigma22) > rsmall).OR.(ABS(sigma12) > rsmall)) THEN 
    1219  
    1220          Angle_denom = 1.0_wp/sqrt( ( sigma11 - sigma22 )*( sigma11 - & 
    1221                        sigma22 ) + 4.0_wp*sigma12*sigma12) 
    1222  
    1223          Q11Q11 = 0.5_wp + ( sigma11 - sigma22 )*0.5_wp*Angle_denom   !Cos^2  
    1224          Q12Q12 = 0.5_wp - ( sigma11 - sigma22 )*0.5_wp*Angle_denom   !Sin^2 
    1225          Q11Q12 = sigma12*Angle_denom                      !CosSin  
    1226       ENDIF 
    1227  
    1228       sigma_1 = Q11Q11*sigma11 + 2.0_wp*Q11Q12*sigma12  & 
    1229               + Q12Q12*sigma22 ! S(1,1) 
    1230       sigma_2 = Q12Q12*sigma11 - 2.0_wp*Q11Q12*sigma12  & 
    1231               + Q11Q11*sigma22 ! S(2,2) 
     1214      zQ11Q11 = 0.1_wp 
     1215      zQ12Q12 = rsmall 
     1216      zQ11Q12 = rsmall 
     1217 
     1218      IF((ABS( zsigma11 - zsigma22) > rsmall).OR.(ABS(zsigma12) > rsmall)) THEN 
     1219 
     1220         zAngle_denom = 1.0_wp/sqrt( ( zsigma11 - zsigma22 )*( zsigma11 - & 
     1221                       zsigma22 ) + 4.0_wp*zsigma12*zsigma12) 
     1222 
     1223         zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2  
     1224         zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Sin^2 
     1225         zQ11Q12 = zsigma12*zAngle_denom                      !CosSin  
     1226      ENDIF 
     1227 
     1228      zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12  & 
     1229              + zQ12Q12*zsigma22 ! S(1,1) 
     1230      zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12  & 
     1231              + zQ11Q11*zsigma22 ! S(2,2) 
    12321232 
    12331233      ! Pure divergence 
    1234       IF ((sigma_1 >= 0.0_wp).AND.(sigma_2 >= 0.0_wp))  THEN 
     1234      IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp))  THEN 
    12351235         mresult11 = 0.0_wp 
    12361236         mresult12 = 0.0_wp 
     
    12391239      ! direction 
    12401240      ! which leads to the loss of their shape, so we again model it through diffusion 
    1241       ELSEIF ((sigma_1 >= 0.0_wp).AND.(sigma_2 < 0.0_wp))  THEN 
    1242          mresult11 = kfrac * (a11 - Q12Q12) 
    1243          mresult12 = kfrac * (a12 + Q11Q12) 
     1241      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN 
     1242         mresult11 = kfrac * (a11 - zQ12Q12) 
     1243         mresult12 = kfrac * (a12 + zQ11Q12) 
    12441244 
    12451245      ! Shear faulting 
    1246       ELSEIF (sigma_2 == 0.0_wp) THEN 
     1246      ELSEIF (zsigma_2 == 0.0_wp) THEN 
    12471247         mresult11 = 0.0_wp 
    12481248         mresult12 = 0.0_wp 
    1249       ELSEIF ((sigma_1 <= 0.0_wp).AND.(sigma_1/sigma_2 <= threshold)) THEN 
    1250          mresult11 = kfrac * (a11 - Q12Q12) 
    1251          mresult12 = kfrac * (a12 + Q11Q12) 
     1249      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 
     1250         mresult11 = kfrac * (a11 - zQ12Q12) 
     1251         mresult12 = kfrac * (a12 + zQ11Q12) 
    12521252 
    12531253      ! Horizontal spalling 
     
    13981398 
    13991399 
    1400    FUNCTION w1(a) 
     1400   FUNCTION w1(pa) 
    14011401      !!------------------------------------------------------------------- 
    14021402      !! Function : w1 (see Gaussian function psi in Tsamados et al 2013) 
    14031403      !!------------------------------------------------------------------- 
    1404       REAL(wp), INTENT(in   ) ::   a 
     1404      REAL(wp), INTENT(in   ) ::   pa 
    14051405      REAL(wp) ::   w1 
    14061406      !!------------------------------------------------------------------- 
    14071407    
    14081408      w1 = - 223.87569446_wp & 
    1409        &   + 2361.2198663_wp*a & 
    1410        &   - 10606.56079975_wp*a*a & 
    1411        &   + 26315.50025642_wp*a*a*a & 
    1412        &   - 38948.30444297_wp*a*a*a*a & 
    1413        &   + 34397.72407466_wp*a*a*a*a*a & 
    1414        &   - 16789.98003081_wp*a*a*a*a*a*a & 
    1415        &   + 3495.82839237_wp*a*a*a*a*a*a*a 
     1409       &   + 2361.2198663_wp*pa & 
     1410       &   - 10606.56079975_wp*pa*pa & 
     1411       &   + 26315.50025642_wp*pa*pa*pa & 
     1412       &   - 38948.30444297_wp*pa*pa*pa*pa & 
     1413       &   + 34397.72407466_wp*pa*pa*pa*pa*pa & 
     1414       &   - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 
     1415       &   + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 
    14161416    
    14171417   END FUNCTION w1 
    14181418 
    14191419 
    1420    FUNCTION w2(a) 
     1420   FUNCTION w2(pa) 
    14211421      !!------------------------------------------------------------------- 
    14221422      !! Function : w2 (see Gaussian function psi in Tsamados et al 2013) 
    14231423      !!------------------------------------------------------------------- 
    1424       REAL(wp), INTENT(in   ) ::   a 
     1424      REAL(wp), INTENT(in   ) ::   pa 
    14251425      REAL(wp) ::   w2 
    14261426      !!------------------------------------------------------------------- 
    14271427    
    14281428      w2 = - 6670.68911883_wp & 
    1429        &   + 70222.33061536_wp*a & 
    1430        &   - 314871.71525448_wp*a*a & 
    1431        &   + 779570.02793492_wp*a*a*a & 
    1432        &   - 1151098.82436864_wp*a*a*a*a & 
    1433        &   + 1013896.59464498_wp*a*a*a*a*a & 
    1434        &   - 493379.44906738_wp*a*a*a*a*a*a & 
    1435        &   + 102356.551518_wp*a*a*a*a*a*a*a 
     1429       &   + 70222.33061536_wp*pa & 
     1430       &   - 314871.71525448_wp*pa*pa & 
     1431       &   + 779570.02793492_wp*pa*pa*pa & 
     1432       &   - 1151098.82436864_wp*pa*pa*pa*pa & 
     1433       &   + 1013896.59464498_wp*pa*pa*pa*pa*pa & 
     1434       &   - 493379.44906738_wp*pa*pa*pa*pa*pa*pa & 
     1435       &   + 102356.551518_wp*pa*pa*pa*pa*pa*pa*pa 
    14361436    
    14371437   END FUNCTION w2 
    14381438 
    1439    FUNCTION s11kr(x,y,z,phi)  
     1439   FUNCTION s11kr(px,py,pz,phi)  
    14401440      !!------------------------------------------------------------------- 
    14411441      !! Function : s11kr 
    14421442      !!------------------------------------------------------------------- 
    1443       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1444     
    1445       REAL(wp) ::   s11kr, p, pih 
    1446     
     1443      REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1444    
     1445      REAL(wp) ::   s11kr, zpih 
     1446    
     1447      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1448      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1449      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1450      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1451      REAL(wp) ::   zd11, zd12, zd22 
     1452      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1453      REAL(wp) ::   zHen1t2, zHen2t1 
     1454      !!------------------------------------------------------------------- 
     1455    
     1456      zpih = 0.5_wp*rpi 
     1457    
     1458      zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 
     1459      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 
     1460      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 
     1461      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 
     1462      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 
     1463      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 
     1464      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 
     1465      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 
     1466      zt1t2i11 = cos(pz-phi) * cos(pz+phi) 
     1467      zt1t2i12 = cos(pz-phi) * sin(pz+phi) 
     1468      zt1t2i21 = sin(pz-phi) * cos(pz+phi) 
     1469      zt1t2i22 = sin(pz-phi) * sin(pz+phi) 
     1470      zt2t1i11 = cos(pz+phi) * cos(pz-phi) 
     1471      zt2t1i12 = cos(pz+phi) * sin(pz-phi) 
     1472      zt2t1i21 = sin(pz+phi) * cos(pz-phi) 
     1473      zt2t1i22 = sin(pz+phi) * sin(pz-phi) 
     1474   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 
     1475   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 
     1476   ! x=x-pi/2 
     1477      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1478      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1479      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1480      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1481      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1482      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1483    
     1484      IF (-zIIn1t2>=rsmall) THEN 
     1485      zHen1t2 = 1._wp 
     1486      ELSE 
     1487      zHen1t2 = 0._wp 
     1488      ENDIF 
     1489    
     1490      IF (-zIIn2t1>=rsmall) THEN 
     1491      zHen2t1 = 1._wp 
     1492      ELSE 
     1493      zHen2t1 = 0._wp 
     1494      ENDIF 
     1495    
     1496      s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 
     1497 
     1498   END FUNCTION s11kr 
     1499 
     1500   FUNCTION s12kr(px,py,pz,phi) 
     1501      !!------------------------------------------------------------------- 
     1502      !! Function : s12kr 
     1503      !!------------------------------------------------------------------- 
     1504      REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1505 
     1506      REAL(wp) ::   s12kr, zs12r0, zs21r0, zpih 
     1507 
     1508      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1509      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1510      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1511      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1512      REAL(wp) ::   zd11, zd12, zd22 
     1513      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1514      REAL(wp) ::   zHen1t2, zHen2t1 
     1515      !!------------------------------------------------------------------- 
     1516      zpih = 0.5_wp*rpi 
     1517    
     1518      zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 
     1519      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 
     1520      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 
     1521      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 
     1522      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 
     1523      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 
     1524      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 
     1525      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 
     1526      zt1t2i11 = cos(pz-phi) * cos(pz+phi) 
     1527      zt1t2i12 = cos(pz-phi) * sin(pz+phi) 
     1528      zt1t2i21 = sin(pz-phi) * cos(pz+phi) 
     1529      zt1t2i22 = sin(pz-phi) * sin(pz+phi) 
     1530      zt2t1i11 = cos(pz+phi) * cos(pz-phi) 
     1531      zt2t1i12 = cos(pz+phi) * sin(pz-phi) 
     1532      zt2t1i21 = sin(pz+phi) * cos(pz-phi) 
     1533      zt2t1i22 = sin(pz+phi) * sin(pz-phi) 
     1534      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1535      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1536      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1537      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1538      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1539      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1540    
     1541      IF (-zIIn1t2>=rsmall) THEN 
     1542      zHen1t2 = 1._wp 
     1543      ELSE 
     1544      zHen1t2 = 0._wp 
     1545      ENDIF 
     1546    
     1547      IF (-zIIn2t1>=rsmall) THEN 
     1548      zHen2t1 = 1._wp 
     1549      ELSE 
     1550      zHen2t1 = 0._wp 
     1551      ENDIF 
     1552    
     1553      zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 
     1554      zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21) 
     1555      s12kr=0.5_wp*(zs12r0+zs21r0) 
     1556    
     1557   END FUNCTION s12kr 
     1558 
     1559   FUNCTION s22kr(px,py,pz,phi) 
     1560      !!------------------------------------------------------------------- 
     1561      !! Function : s22kr 
     1562      !!------------------------------------------------------------------- 
     1563      REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1564 
     1565      REAL(wp) ::   s22kr, zpih 
     1566 
     1567      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1568      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1569      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1570      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1571      REAL(wp) ::   zd11, zd12, zd22 
     1572      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1573      REAL(wp) ::   zHen1t2, zHen2t1 
     1574      !!------------------------------------------------------------------- 
     1575      zpih = 0.5_wp*rpi 
     1576 
     1577      zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 
     1578      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 
     1579      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 
     1580      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 
     1581      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 
     1582      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 
     1583      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 
     1584      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 
     1585      zt1t2i11 = cos(pz-phi) * cos(pz+phi) 
     1586      zt1t2i12 = cos(pz-phi) * sin(pz+phi) 
     1587      zt1t2i21 = sin(pz-phi) * cos(pz+phi) 
     1588      zt1t2i22 = sin(pz-phi) * sin(pz+phi) 
     1589      zt2t1i11 = cos(pz+phi) * cos(pz-phi) 
     1590      zt2t1i12 = cos(pz+phi) * sin(pz-phi) 
     1591      zt2t1i21 = sin(pz+phi) * cos(pz-phi) 
     1592      zt2t1i22 = sin(pz+phi) * sin(pz-phi) 
     1593      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1594      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1595      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1596      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1597      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1598      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1599 
     1600      IF (-zIIn1t2>=rsmall) THEN 
     1601      zHen1t2 = 1._wp 
     1602      ELSE 
     1603      zHen1t2 = 0._wp 
     1604      ENDIF 
     1605 
     1606      IF (-zIIn2t1>=rsmall) THEN 
     1607      zHen2t1 = 1._wp 
     1608      ELSE 
     1609      zHen2t1 = 0._wp 
     1610      ENDIF 
     1611 
     1612      s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 
     1613 
     1614   END FUNCTION s22kr 
     1615 
     1616   FUNCTION s11ks(px,py,pz,phi) 
     1617      !!------------------------------------------------------------------- 
     1618      !! Function : s11ks 
     1619      !!------------------------------------------------------------------- 
     1620      REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1621 
     1622      REAL(wp) ::   s11ks, zpih 
     1623 
     1624      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1625      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1626      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1627      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1628      REAL(wp) ::   zd11, zd12, zd22 
     1629      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1630      REAL(wp) ::   zHen1t2, zHen2t1 
     1631      !!------------------------------------------------------------------- 
     1632      zpih = 0.5_wp*rpi 
     1633 
     1634      zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 
     1635      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 
     1636      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 
     1637      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 
     1638      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 
     1639      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 
     1640      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 
     1641      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 
     1642      zt1t2i11 = cos(pz-phi) * cos(pz+phi) 
     1643      zt1t2i12 = cos(pz-phi) * sin(pz+phi) 
     1644      zt1t2i21 = sin(pz-phi) * cos(pz+phi) 
     1645      zt1t2i22 = sin(pz-phi) * sin(pz+phi) 
     1646      zt2t1i11 = cos(pz+phi) * cos(pz-phi) 
     1647      zt2t1i12 = cos(pz+phi) * sin(pz-phi) 
     1648      zt2t1i21 = sin(pz+phi) * cos(pz-phi) 
     1649      zt2t1i22 = sin(pz+phi) * sin(pz-phi) 
     1650      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1651      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1652      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1653      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1654      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1655      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1656 
     1657      IF (-zIIn1t2>=rsmall) THEN 
     1658      zHen1t2 = 1._wp 
     1659      ELSE 
     1660      zHen1t2 = 0._wp 
     1661      ENDIF 
     1662 
     1663      IF (-zIIn2t1>=rsmall) THEN 
     1664      zHen2t1 = 1._wp 
     1665      ELSE 
     1666      zHen2t1 = 0._wp 
     1667      ENDIF 
     1668 
     1669      s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 
     1670 
     1671   END FUNCTION s11ks 
     1672 
     1673   FUNCTION s12ks(px,py,pz,phi) 
     1674      !!------------------------------------------------------------------- 
     1675      !! Function : s12ks 
     1676      !!------------------------------------------------------------------- 
     1677      REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1678 
     1679      REAL(wp) ::   s12ks, zs12s0, zs21s0, zpih 
     1680 
     1681      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1682      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1683      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1684      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1685      REAL(wp) ::   zd11, zd12, zd22 
     1686      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1687      REAL(wp) ::   zHen1t2, zHen2t1 
     1688      !!------------------------------------------------------------------- 
     1689      zpih = 0.5_wp*rpi 
     1690 
     1691      zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 
     1692      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 
     1693      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 
     1694      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 
     1695      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 
     1696      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 
     1697      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 
     1698      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 
     1699      zt1t2i11 = cos(pz-phi) * cos(pz+phi) 
     1700      zt1t2i12 = cos(pz-phi) * sin(pz+phi) 
     1701      zt1t2i21 = sin(pz-phi) * cos(pz+phi) 
     1702      zt1t2i22 = sin(pz-phi) * sin(pz+phi) 
     1703      zt2t1i11 = cos(pz+phi) * cos(pz-phi) 
     1704      zt2t1i12 = cos(pz+phi) * sin(pz-phi) 
     1705      zt2t1i21 = sin(pz+phi) * cos(pz-phi) 
     1706      zt2t1i22 = sin(pz+phi) * sin(pz-phi) 
     1707      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1708      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1709      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1710      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1711      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1712      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1713 
     1714      IF (-zIIn1t2>=rsmall) THEN 
     1715      zHen1t2 = 1._wp 
     1716      ELSE 
     1717      zHen1t2 = 0._wp 
     1718      ENDIF 
     1719 
     1720      IF (-zIIn2t1>=rsmall) THEN 
     1721      zHen2t1 = 1._wp 
     1722      ELSE 
     1723      zHen2t1 = 0._wp 
     1724      ENDIF 
     1725 
     1726      zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 
     1727      zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21) 
     1728      s12ks=0.5_wp*(zs12s0+zs21s0) 
     1729 
     1730   END FUNCTION s12ks 
     1731 
     1732   FUNCTION s22ks(px,py,pz,phi)  
     1733      !!------------------------------------------------------------------- 
     1734      !! Function : s22ks 
     1735      !!------------------------------------------------------------------- 
     1736      REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1737 
     1738      REAL(wp) ::   s22ks, zpih 
     1739 
    14471740      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    14481741      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
     
    14531746      REAL(wp) ::   Hen1t2, Hen2t1 
    14541747      !!------------------------------------------------------------------- 
    1455     
    1456       pih = 0.5_wp*rpi 
    1457       p = phi 
    1458     
    1459       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1460       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1461       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1462       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1463       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1464       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1465       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1466       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1467       t1t2i11 = cos(z-p) * cos(z+p) 
    1468       t1t2i12 = cos(z-p) * sin(z+p) 
    1469       t1t2i21 = sin(z-p) * cos(z+p) 
    1470       t1t2i22 = sin(z-p) * sin(z+p) 
    1471       t2t1i11 = cos(z+p) * cos(z-p) 
    1472       t2t1i12 = cos(z+p) * sin(z-p) 
    1473       t2t1i21 = sin(z+p) * cos(z-p) 
    1474       t2t1i22 = sin(z+p) * sin(z-p) 
    1475    ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 
    1476    ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 
    1477    ! x=x-pi/2 
    1478       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1479       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1480       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
     1748      zpih = 0.5_wp*rpi 
     1749 
     1750      n1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 
     1751      n1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 
     1752      n1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 
     1753      n1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 
     1754      n2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 
     1755      n2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 
     1756      n2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 
     1757      n2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 
     1758      t1t2i11 = cos(pz-phi) * cos(pz+phi) 
     1759      t1t2i12 = cos(pz-phi) * sin(pz+phi) 
     1760      t1t2i21 = sin(pz-phi) * cos(pz+phi) 
     1761      t1t2i22 = sin(pz-phi) * sin(pz+phi) 
     1762      t2t1i11 = cos(pz+phi) * cos(pz-phi) 
     1763      t2t1i12 = cos(pz+phi) * sin(pz-phi) 
     1764      t2t1i21 = sin(pz+phi) * cos(pz-phi) 
     1765      t2t1i22 = sin(pz+phi) * sin(pz-phi) 
     1766      d11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1767      d12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1768      d22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
    14811769      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    14821770      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    14831771      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1484     
     1772 
    14851773      IF (-IIn1t2>=rsmall) THEN 
    14861774      Hen1t2 = 1._wp 
     
    14881776      Hen1t2 = 0._wp 
    14891777      ENDIF 
    1490     
     1778 
    14911779      IF (-IIn2t1>=rsmall) THEN 
    14921780      Hen2t1 = 1._wp 
     
    14941782      Hen2t1 = 0._wp 
    14951783      ENDIF 
    1496     
    1497       s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 * n2t1i11) 
    1498  
    1499    END FUNCTION s11kr 
    1500  
    1501    FUNCTION s12kr(x,y,z,phi) 
    1502       !!------------------------------------------------------------------- 
    1503       !! Function : s12kr 
    1504       !!------------------------------------------------------------------- 
    1505       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1506  
    1507       REAL(wp) ::   s12kr, s12r0, s21r0, p, pih 
    1508  
    1509       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1510       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1511       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1512       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1513       REAL(wp) ::   d11, d12, d22 
    1514       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1515       REAL(wp) ::   Hen1t2, Hen2t1 
    1516       !!------------------------------------------------------------------- 
    1517       pih = 0.5_wp*rpi 
    1518       p = phi 
    1519     
    1520       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1521       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1522       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1523       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1524       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1525       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1526       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1527       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1528       t1t2i11 = cos(z-p) * cos(z+p) 
    1529       t1t2i12 = cos(z-p) * sin(z+p) 
    1530       t1t2i21 = sin(z-p) * cos(z+p) 
    1531       t1t2i22 = sin(z-p) * sin(z+p) 
    1532       t2t1i11 = cos(z+p) * cos(z-p) 
    1533       t2t1i12 = cos(z+p) * sin(z-p) 
    1534       t2t1i21 = sin(z+p) * cos(z-p) 
    1535       t2t1i22 = sin(z+p) * sin(z-p) 
    1536       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1537       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1538       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1539       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1540       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1541       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1542     
    1543       IF (-IIn1t2>=rsmall) THEN 
    1544       Hen1t2 = 1._wp 
    1545       ELSE 
    1546       Hen1t2 = 0._wp 
    1547       ENDIF 
    1548     
    1549       IF (-IIn2t1>=rsmall) THEN 
    1550       Hen2t1 = 1._wp 
    1551       ELSE 
    1552       Hen2t1 = 0._wp 
    1553       ENDIF 
    1554     
    1555       s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 * n2t1i12) 
    1556       s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 * n2t1i21) 
    1557       s12kr=0.5_wp*(s12r0+s21r0) 
    1558     
    1559    END FUNCTION s12kr 
    1560  
    1561    FUNCTION s22kr(x,y,z,phi) 
    1562       !!------------------------------------------------------------------- 
    1563       !! Function : s22kr 
    1564       !!------------------------------------------------------------------- 
    1565       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1566  
    1567       REAL(wp) ::   s22kr, p, pih 
    1568  
    1569       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1570       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1571       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1572       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1573       REAL(wp) ::   d11, d12, d22 
    1574       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1575       REAL(wp) ::   Hen1t2, Hen2t1 
    1576       !!------------------------------------------------------------------- 
    1577       pih = 0.5_wp*rpi 
    1578       p = phi 
    1579  
    1580       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1581       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1582       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1583       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1584       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1585       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1586       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1587       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1588       t1t2i11 = cos(z-p) * cos(z+p) 
    1589       t1t2i12 = cos(z-p) * sin(z+p) 
    1590       t1t2i21 = sin(z-p) * cos(z+p) 
    1591       t1t2i22 = sin(z-p) * sin(z+p) 
    1592       t2t1i11 = cos(z+p) * cos(z-p) 
    1593       t2t1i12 = cos(z+p) * sin(z-p) 
    1594       t2t1i21 = sin(z+p) * cos(z-p) 
    1595       t2t1i22 = sin(z+p) * sin(z-p) 
    1596       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1597       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1598       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1599       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1600       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1601       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1602  
    1603       IF (-IIn1t2>=rsmall) THEN 
    1604       Hen1t2 = 1._wp 
    1605       ELSE 
    1606       Hen1t2 = 0._wp 
    1607       ENDIF 
    1608  
    1609       IF (-IIn2t1>=rsmall) THEN 
    1610       Hen2t1 = 1._wp 
    1611       ELSE 
    1612       Hen2t1 = 0._wp 
    1613       ENDIF 
    1614  
    1615       s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 * n2t1i22) 
    1616  
    1617    END FUNCTION s22kr 
    1618  
    1619    FUNCTION s11ks(x,y,z,phi) 
    1620       !!------------------------------------------------------------------- 
    1621       !! Function : s11ks 
    1622       !!------------------------------------------------------------------- 
    1623       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1624  
    1625       REAL(wp) ::   s11ks, p, pih 
    1626  
    1627       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1628       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1629       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1630       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1631       REAL(wp) ::   d11, d12, d22 
    1632       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1633       REAL(wp) ::   Hen1t2, Hen2t1 
    1634       !!------------------------------------------------------------------- 
    1635       pih = 0.5_wp*rpi 
    1636       p = phi 
    1637  
    1638       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1639       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1640       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1641       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1642       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1643       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1644       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1645       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1646       t1t2i11 = cos(z-p) * cos(z+p) 
    1647       t1t2i12 = cos(z-p) * sin(z+p) 
    1648       t1t2i21 = sin(z-p) * cos(z+p) 
    1649       t1t2i22 = sin(z-p) * sin(z+p) 
    1650       t2t1i11 = cos(z+p) * cos(z-p) 
    1651       t2t1i12 = cos(z+p) * sin(z-p) 
    1652       t2t1i21 = sin(z+p) * cos(z-p) 
    1653       t2t1i22 = sin(z+p) * sin(z-p) 
    1654       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1655       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1656       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1657       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1658       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1659       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1660  
    1661       IF (-IIn1t2>=rsmall) THEN 
    1662       Hen1t2 = 1._wp 
    1663       ELSE 
    1664       Hen1t2 = 0._wp 
    1665       ENDIF 
    1666  
    1667       IF (-IIn2t1>=rsmall) THEN 
    1668       Hen2t1 = 1._wp 
    1669       ELSE 
    1670       Hen2t1 = 0._wp 
    1671       ENDIF 
    1672  
    1673       s11ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i11 + Hen2t1 * t2t1i11) 
    1674  
    1675    END FUNCTION s11ks 
    1676  
    1677    FUNCTION s12ks(x,y,z,phi) 
    1678       !!------------------------------------------------------------------- 
    1679       !! Function : s12ks 
    1680       !!------------------------------------------------------------------- 
    1681       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1682  
    1683       REAL(wp) ::   s12ks, s12s0, s21s0, p, pih 
    1684  
    1685       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1686       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1687       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1688       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1689       REAL(wp) ::   d11, d12, d22 
    1690       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1691       REAL(wp) ::   Hen1t2, Hen2t1 
    1692       !!------------------------------------------------------------------- 
    1693       pih = 0.5_wp*rpi 
    1694       p =phi 
    1695  
    1696       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1697       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1698       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1699       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1700       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1701       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1702       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1703       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1704       t1t2i11 = cos(z-p) * cos(z+p) 
    1705       t1t2i12 = cos(z-p) * sin(z+p) 
    1706       t1t2i21 = sin(z-p) * cos(z+p) 
    1707       t1t2i22 = sin(z-p) * sin(z+p) 
    1708       t2t1i11 = cos(z+p) * cos(z-p) 
    1709       t2t1i12 = cos(z+p) * sin(z-p) 
    1710       t2t1i21 = sin(z+p) * cos(z-p) 
    1711       t2t1i22 = sin(z+p) * sin(z-p) 
    1712       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1713       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1714       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1715       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1716       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1717       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1718  
    1719       IF (-IIn1t2>=rsmall) THEN 
    1720       Hen1t2 = 1._wp 
    1721       ELSE 
    1722       Hen1t2 = 0._wp 
    1723       ENDIF 
    1724  
    1725       IF (-IIn2t1>=rsmall) THEN 
    1726       Hen2t1 = 1._wp 
    1727       ELSE 
    1728       Hen2t1 = 0._wp 
    1729       ENDIF 
    1730  
    1731       s12s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i12 + Hen2t1 * t2t1i12) 
    1732       s21s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i21 + Hen2t1 * t2t1i21) 
    1733       s12ks=0.5_wp*(s12s0+s21s0) 
    1734  
    1735    END FUNCTION s12ks 
    1736  
    1737    FUNCTION s22ks(x,y,z,phi)  
    1738       !!------------------------------------------------------------------- 
    1739       !! Function : s22ks 
    1740       !!------------------------------------------------------------------- 
    1741       REAL(wp), INTENT(in   ) ::   x,y,z,phi 
    1742  
    1743       REAL(wp) ::   s22ks, p, pih 
    1744  
    1745       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1746       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1747       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1748       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1749       REAL(wp) ::   d11, d12, d22 
    1750       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1751       REAL(wp) ::   Hen1t2, Hen2t1 
    1752       !!------------------------------------------------------------------- 
    1753       pih = 0.5_wp*rpi 
    1754       p =phi 
    1755  
    1756       n1t2i11 = cos(z+pih-p) * cos(z+p) 
    1757       n1t2i12 = cos(z+pih-p) * sin(z+p) 
    1758       n1t2i21 = sin(z+pih-p) * cos(z+p) 
    1759       n1t2i22 = sin(z+pih-p) * sin(z+p) 
    1760       n2t1i11 = cos(z-pih+p) * cos(z-p) 
    1761       n2t1i12 = cos(z-pih+p) * sin(z-p) 
    1762       n2t1i21 = sin(z-pih+p) * cos(z-p) 
    1763       n2t1i22 = sin(z-pih+p) * sin(z-p) 
    1764       t1t2i11 = cos(z-p) * cos(z+p) 
    1765       t1t2i12 = cos(z-p) * sin(z+p) 
    1766       t1t2i21 = sin(z-p) * cos(z+p) 
    1767       t1t2i22 = sin(z-p) * sin(z+p) 
    1768       t2t1i11 = cos(z+p) * cos(z-p) 
    1769       t2t1i12 = cos(z+p) * sin(z-p) 
    1770       t2t1i21 = sin(z+p) * cos(z-p) 
    1771       t2t1i22 = sin(z+p) * sin(z-p) 
    1772       d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
    1773       d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
    1774       d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
    1775       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1776       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1777       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1778  
    1779       IF (-IIn1t2>=rsmall) THEN 
    1780       Hen1t2 = 1._wp 
    1781       ELSE 
    1782       Hen1t2 = 0._wp 
    1783       ENDIF 
    1784  
    1785       IF (-IIn2t1>=rsmall) THEN 
    1786       Hen2t1 = 1._wp 
    1787       ELSE 
    1788       Hen2t1 = 0._wp 
    1789       ENDIF 
    17901784 
    17911785      s22ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22) 
Note: See TracChangeset for help on using the changeset viewer.