New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13155 – NEMO

Changeset 13155


Ignore:
Timestamp:
2020-06-24T17:18:26+02:00 (4 years ago)
Author:
stefryn
Message:

final variable updates

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

Legend:

Unmodified
Added
Removed

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

    r12261 r13155  
    139139      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not 
    140140      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i 
    141       REAL(wp), DIMENSION(jpij) ::   conv          ! 1D rdg_conv (if EAP rheology) 
     141      REAL(wp), DIMENSION(jpij) ::   zconv         ! 1D rdg_conv (if EAP rheology) 
    142142      ! 
    143143      INTEGER, PARAMETER ::   jp_itermax = 20     
    … …  
    177177         ! just needed here 
    178178         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i  ) 
    179          CALL tab_2d_1d( npti, nptidx(1:npti), conv    (1:npti)      , rdg_conv ) 
     179         CALL tab_2d_1d( npti, nptidx(1:npti), zconv   (1:npti)      , rdg_conv ) 
    180180         ! needed here and in the iteration loop 
    181181         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i) 
    … …  
    189189            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    190190            IF( ln_rhg_EVP )  closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    191             IF( ln_rhg_EAP )  closing_net(ji) = conv(ji) 
     191            IF( ln_rhg_EAP )  closing_net(ji) = zconv(ji) 
    192192            ! 
    193193            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough 

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

    r13128 r13155  
    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 
    … …  
    144144      REAL(wp) ::   zstressptmp, zstressmtmp, zstress12tmpF                ! anisotropic stress tensor components 
    145145      REAL(wp) ::   zalphar, zalphas                                      ! for mechanical redistribution 
    146       REAL(wp) ::   mresult11, mresult12, z1dtevpkth, p5kth, z1_dtevp_A        ! for structure tensor evolution 
     146      REAL(wp) ::   zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A        ! for structure tensor evolution 
    147147      ! 
    148148      REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                     ! anisotropic stress tensor component for regridding 
    149       REAL(wp), DIMENSION(jpi,jpj) ::   yield11, yield22, yield12       ! yield surface tensor for history 
     149      REAL(wp), DIMENSION(jpi,jpj) ::   zyield11, zyield22, zyield12       ! yield surface tensor for history 
    150150      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    151151      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    … …  
    258258      z1_dtevp_A = z1_dtevp/10.0_wp 
    259259      z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp) 
    260       p5kth = 0.5_wp * 0.00002_wp 
     260      zp5kth = 0.5_wp * 0.00002_wp 
    261261 
    262262      ! Ice strength 
    … …  
    413413                
    414414               ! --- anisotropic stress calculation --- ! 
    415                CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, & 
     415               CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, & 
    416416                                       zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 
    417417                                       zstressptmp, zstressmtmp, & 
    … …  
    423423                  CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 
    424424                                  paniso_11(ji,jj), paniso_12(ji,jj),           & 
    425                                   mresult11,  mresult12) 
    426  
    427                   paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A  + p5kth - mresult11) * z1dtevpkth ! implicit 
    428                   paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  - mresult12) * z1dtevpkth ! implicit 
     425                                  zmresult11,  zmresult12) 
     426 
     427                  paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A  + zp5kth - zmresult11) * z1dtevpkth ! implicit 
     428                  paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  - zmresult12) * z1dtevpkth ! implicit 
    429429               ENDIF 
    430430 
    431431 
    432432               IF (jter == nn_nevp) THEN 
    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) 
     433                  zyield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp) 
     434                  zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp) 
     435                  zyield12(ji,jj) = zstress12tmp(ji,jj) 
    436436                  prdg_conv(ji,jj) = -min( zalphar, 0._wp)     
    437437               ENDIF 
    … …  
    484484            END DO 
    485485         END DO 
     486      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    486487 
    487488         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
    … …  
    861862      IF( iom_use('yield11') .OR. iom_use('yield12') .OR. iom_use('yield22')) THEN 
    862863 
    863          CALL lbc_lnk_multi( 'icedyn_rhg_eap', yield11, 'T', 1., yield22, 'T', 1., yield12, 'T', 1. ) 
    864  
    865          CALL iom_put( 'yield11', yield11 * zmsk00 ) 
    866          CALL iom_put( 'yield22', yield22 * zmsk00 ) 
    867          CALL iom_put( 'yield12', yield12 * zmsk00 ) 
     864         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1., zyield22, 'T', 1., zyield12, 'T', 1. ) 
     865 
     866         CALL iom_put( 'yield11', zyield11 * zmsk00 ) 
     867         CALL iom_put( 'yield22', zyield22 * zmsk00 ) 
     868         CALL iom_put( 'yield12', zyield12 * zmsk00 ) 
    868869      ENDIF 
    869870 
    … …  
    931932   END SUBROUTINE ice_dyn_rhg_eap 
    932933 
    933    SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu, tension, & 
    934                                  shear, a11, a12, & 
    935                                  stressp,  stressm, & 
    936                                  stress12, strength, & 
     934   SUBROUTINE update_stress_rdg (ksub, kndte, pdivu, ptension, & 
     935                                 pshear, pa11, pa12, & 
     936                                 pstressp,  pstressm, & 
     937                                 pstress12, strength, & 
    937938                                 palphar, palphas) 
    938939      !!--------------------------------------------------------------------- 
    939       !!                   ***  SUBROUTINE rhg_eap_rst  *** 
     940      !!                   ***  SUBROUTINE update_stress_rdg  *** 
    940941      !!                      
    941       !! ** Purpose :   Updates the depending on values of strain rate and structure 
    942       !!                tensor and for ksub=ndte it computes closing and sliding rate 
     942      !! ** Purpose :   Updates the stress depending on values of strain rate and structure 
     943      !!                tensor and for the last subcycle step it computes closing and sliding rate 
    943944      !!--------------------------------------------------------------------- 
    944       INTEGER,  INTENT(in   ) ::   ksub, ndte 
    945       REAL(wp), INTENT(in   ) ::   phi, strength 
    946       REAL(wp), INTENT(in   ) ::   divu, tension, shear 
    947       REAL(wp), INTENT(in   ) ::   a11, a12 
    948       REAL(wp), INTENT(  out) ::   stressp, stressm, stress12 
     945      INTEGER,  INTENT(in   ) ::   ksub, kndte 
     946      REAL(wp), INTENT(in   ) ::   strength 
     947      REAL(wp), INTENT(in   ) ::   pdivu, ptension, pshear 
     948      REAL(wp), INTENT(in   ) ::   pa11, pa12 
     949      REAL(wp), INTENT(  out) ::   pstressp, pstressm, pstress12 
    949950      REAL(wp), INTENT(  out) ::   palphar, palphas 
    950951 
    … …  
    973974      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
    974975  
    975       zinvsin = 1._wp/sin(2._wp*phi) * zinvstressconviso  
     976      zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso  
    976977      !now uses phi as set in higher code 
    977978        
    … …  
    980981 
    981982      ! 1) structure tensor 
    982       za22 = 1._wp-a11 
     983      za22 = 1._wp-pa11 
    983984      zQ11Q11 = 1._wp 
    984985      zQ12Q12 = rsmall 
    … …  
    987988      ! gamma: angle between general coordiantes and principal axis of A  
    988989      ! here Tan2gamma = 2 a12 / (a11 - a22)  
    989       IF((ABS(a11 - za22) > rsmall).OR.(ABS(a12) > rsmall)) THEN       
    990          zAngle_denom_gamma = 1._wp/sqrt( ( a11 - za22 )*( a11 - za22) + & 
    991                              4._wp*a12*a12 ) 
    992     
    993          zQ11Q11 = 0.5_wp + ( a11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2  
    994          zQ12Q12 = 0.5_wp - ( a11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2 
    995          zQ11Q12 = a12*zAngle_denom_gamma                     !CosSin  
     990      IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN       
     991         zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + & 
     992                             4._wp*pa12*pa12 ) 
     993    
     994         zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2  
     995         zQ12Q12 = 0.5_wp - ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2 
     996         zQ11Q12 = pa12*zAngle_denom_gamma                     !CosSin  
    996997      ENDIF 
    997998  
    998999      ! rotation Q*atemp*Q^T 
    999       zatempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22  
     1000      zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22  
    10001001       
    10011002      ! make first principal value the largest 
    … …  
    10031004  
    10041005      ! 2) strain rate 
    1005       zdtemp11 = 0.5_wp*(divu + tension) 
    1006       zdtemp12 = shear*0.5_wp 
    1007       zdtemp22 = 0.5_wp*(divu - tension) 
    1008  
    1009       !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear 
     1006      zdtemp11 = 0.5_wp*(pdivu + ptension) 
     1007      zdtemp12 = pshear*0.5_wp 
     1008      zdtemp22 = 0.5_wp*(pdivu - ptension) 
    10101009 
    10111010      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)  
    … …  
    11321131      zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin 
    11331132 
    1134       !WRITE(numout,*) 'strength inside loop', strength 
    11351133      !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22 
    11361134 
    … …  
    11421140      zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 +        2._wp*zQ11Q12 *zsig12 
    11431141 
    1144       stressp  = zsgprm11 + zsgprm22 
    1145       stress12 = zsgprm12 
    1146       stressm  = zsgprm11 - zsgprm22 
    1147  
    1148       !WRITE(numout,*) ' output inside loop', ksub, stressp 
     1142      pstressp  = zsgprm11 + zsgprm22 
     1143      pstress12 = zsgprm12 
     1144      pstressm  = zsgprm11 - zsgprm22 
    11491145 
    11501146      ! Compute ridging and sliding functions in general coordinates  
    11511147      ! (Equation 11 in Tsamados 2013) 
    1152       IF (ksub == ndte) THEN 
     1148      IF (ksub == kndte) THEN 
    11531149         zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r & 
    11541150                     + zQ12Q12*zstemp22r 
    … …  
    11751171 
    11761172 
    1177    SUBROUTINE calc_ffrac (stressp, stressm, stress12, & 
    1178                           a11, a12,                   & 
    1179                           mresult11, mresult12) 
     1173   SUBROUTINE calc_ffrac (pstressp, pstressm, pstress12, & 
     1174                          pa11, pa12,                   & 
     1175                          pmresult11, pmresult12) 
    11801176      !!--------------------------------------------------------------------- 
    11811177      !!                   ***  ROUTINE calc_ffrac  *** 
    … …  
    11851181      !! ** Method :    Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2 
    11861182      !!--------------------------------------------------------------------- 
    1187       REAL (wp), INTENT(in)  ::   stressp, stressm, stress12, a11, a12 
    1188       REAL (wp), INTENT(out) ::   mresult11, mresult12 
     1183      REAL (wp), INTENT(in)  ::   pstressp, pstressm, pstress12, pa11, pa12 
     1184      REAL (wp), INTENT(out) ::   pmresult11, pmresult12 
    11891185 
    11901186      ! local variables 
    … …  
    11991195      !!--------------------------------------------------------------- 
    12001196      ! 
    1201       zsigma11 = 0.5_wp*(stressp+stressm)  
    1202       zsigma12 = stress12  
    1203       zsigma22 = 0.5_wp*(stressp-stressm)  
     1197      zsigma11 = 0.5_wp*(pstressp+pstressm)  
     1198      zsigma12 = pstress12  
     1199      zsigma22 = 0.5_wp*(pstressp-pstressm)  
    12041200 
    12051201      ! Here's the change - no longer calculate gamma, 
    … …  
    12331229      ! Pure divergence 
    12341230      IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp))  THEN 
    1235          mresult11 = 0.0_wp 
    1236          mresult12 = 0.0_wp 
     1231         pmresult11 = 0.0_wp 
     1232         pmresult12 = 0.0_wp 
    12371233 
    12381234      ! Unconfined compression: cracking of blocks not along the axial splitting 
    … …  
    12401236      ! which leads to the loss of their shape, so we again model it through diffusion 
    12411237      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN 
    1242          mresult11 = kfrac * (a11 - zQ12Q12) 
    1243          mresult12 = kfrac * (a12 + zQ11Q12) 
     1238         pmresult11 = kfrac * (pa11 - zQ12Q12) 
     1239         pmresult12 = kfrac * (pa12 + zQ11Q12) 
    12441240 
    12451241      ! Shear faulting 
    12461242      ELSEIF (zsigma_2 == 0.0_wp) THEN 
    1247          mresult11 = 0.0_wp 
    1248          mresult12 = 0.0_wp 
     1243         pmresult11 = 0.0_wp 
     1244         pmresult12 = 0.0_wp 
    12491245      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 
    1250          mresult11 = kfrac * (a11 - zQ12Q12) 
    1251          mresult12 = kfrac * (a12 + zQ11Q12) 
     1246         pmresult11 = kfrac * (pa11 - zQ12Q12) 
     1247         pmresult12 = kfrac * (pa12 + zQ11Q12) 
    12521248 
    12531249      ! Horizontal spalling 
    12541250      ELSE  
    1255          mresult11 = 0.0_wp 
    1256          mresult12 = 0.0_wp 
     1251         pmresult11 = 0.0_wp 
     1252         pmresult12 = 0.0_wp 
    12571253      ENDIF 
    12581254 
    … …  
    13411337          s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    13421338           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1343            s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1339           s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    13441340          s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    13451341           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1346            s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1342           s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    13471343          s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    13481344           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1349            s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)  
     1345           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)  
    13501346          s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    13511347           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1352            s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1348           s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    13531349          s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    13541350           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1355            s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1351           s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    13561352          s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
    13571353           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
    1358            s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     1354           s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
    13591355          ENDDO 
    13601356          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 
    … …  
    13651361          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 
    13661362         ELSE 
    1367           s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1368           s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1369           s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1370           s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1371           s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
    1372           s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     1363          s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1364          s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1365          s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1366          s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1367          s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
     1368          s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 
    13731369          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 
    13741370          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 
    … …  
    14371433   END FUNCTION w2 
    14381434 
    1439    FUNCTION s11kr(px,py,pz,phi)  
     1435   FUNCTION s11kr(px,py,pz)  
    14401436      !!------------------------------------------------------------------- 
    14411437      !! Function : s11kr 
    14421438      !!------------------------------------------------------------------- 
    1443       REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1439      REAL(wp), INTENT(in   ) ::   px,py,pz 
    14441440    
    14451441      REAL(wp) ::   s11kr, zpih 
    … …  
    14561452      zpih = 0.5_wp*rpi 
    14571453    
    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) 
     1454      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1455      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1456      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1457      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1458      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1459      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1460      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1461      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1462      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1463      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1464      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1465      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1466      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1467      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1468      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1469      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    14741470   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 
    14751471   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 
    … …  
    14981494   END FUNCTION s11kr 
    14991495 
    1500    FUNCTION s12kr(px,py,pz,phi) 
     1496   FUNCTION s12kr(px,py,pz) 
    15011497      !!------------------------------------------------------------------- 
    15021498      !! Function : s12kr 
    15031499      !!------------------------------------------------------------------- 
    1504       REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1500      REAL(wp), INTENT(in   ) ::   px,py,pz 
    15051501 
    15061502      REAL(wp) ::   s12kr, zs12r0, zs21r0, zpih 
    … …  
    15161512      zpih = 0.5_wp*rpi 
    15171513    
    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) 
     1514      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1515      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1516      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1517      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1518      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1519      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1520      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1521      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1522      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1523      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1524      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1525      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1526      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1527      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1528      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1529      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    15341530      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    15351531      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    … …  
    15571553   END FUNCTION s12kr 
    15581554 
    1559    FUNCTION s22kr(px,py,pz,phi) 
     1555   FUNCTION s22kr(px,py,pz) 
    15601556      !!------------------------------------------------------------------- 
    15611557      !! Function : s22kr 
    15621558      !!------------------------------------------------------------------- 
    1563       REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1559      REAL(wp), INTENT(in   ) ::   px,py,pz 
    15641560 
    15651561      REAL(wp) ::   s22kr, zpih 
    … …  
    15751571      zpih = 0.5_wp*rpi 
    15761572 
    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) 
     1573      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1574      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1575      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1576      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1577      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1578      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1579      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1580      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1581      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1582      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1583      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1584      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1585      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1586      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1587      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1588      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    15931589      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    15941590      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    … …  
    16141610   END FUNCTION s22kr 
    16151611 
    1616    FUNCTION s11ks(px,py,pz,phi) 
     1612   FUNCTION s11ks(px,py,pz) 
    16171613      !!------------------------------------------------------------------- 
    16181614      !! Function : s11ks 
    16191615      !!------------------------------------------------------------------- 
    1620       REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1616      REAL(wp), INTENT(in   ) ::   px,py,pz 
    16211617 
    16221618      REAL(wp) ::   s11ks, zpih 
    … …  
    16321628      zpih = 0.5_wp*rpi 
    16331629 
    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) 
     1630      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1631      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1632      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1633      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1634      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1635      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1636      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1637      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1638      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1639      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1640      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1641      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1642      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1643      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1644      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1645      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    16501646      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    16511647      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    … …  
    16711667   END FUNCTION s11ks 
    16721668 
    1673    FUNCTION s12ks(px,py,pz,phi) 
     1669   FUNCTION s12ks(px,py,pz) 
    16741670      !!------------------------------------------------------------------- 
    16751671      !! Function : s12ks 
    16761672      !!------------------------------------------------------------------- 
    1677       REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1673      REAL(wp), INTENT(in   ) ::   px,py,pz 
    16781674 
    16791675      REAL(wp) ::   s12ks, zs12s0, zs21s0, zpih 
    … …  
    16891685      zpih = 0.5_wp*rpi 
    16901686 
    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) 
     1687      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1688      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1689      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1690      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1691      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1692      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1693      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1694      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1695      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1696      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1697      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1698      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1699      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1700      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1701      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1702      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
    17071703      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
    17081704      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
    … …  
    17301726   END FUNCTION s12ks 
    17311727 
    1732    FUNCTION s22ks(px,py,pz,phi)  
     1728   FUNCTION s22ks(px,py,pz)  
    17331729      !!------------------------------------------------------------------- 
    17341730      !! Function : s22ks 
    17351731      !!------------------------------------------------------------------- 
    1736       REAL(wp), INTENT(in   ) ::   px,py,pz,phi 
     1732      REAL(wp), INTENT(in   ) ::   px,py,pz 
    17371733 
    17381734      REAL(wp) ::   s22ks, zpih 
    17391735 
    1740       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
    1741       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
    1742       REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
    1743       REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
    1744       REAL(wp) ::   d11, d12, d22 
    1745       REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
    1746       REAL(wp) ::   Hen1t2, Hen2t1 
     1736      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 
     1737      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 
     1738      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 
     1739      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 
     1740      REAL(wp) ::   zd11, zd12, zd22 
     1741      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2 
     1742      REAL(wp) ::   zHen1t2, zHen2t1 
    17471743      !!------------------------------------------------------------------- 
    17481744      zpih = 0.5_wp*rpi 
    17491745 
    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)) 
    1769       IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
    1770       IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
    1771       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    1772  
    1773       IF (-IIn1t2>=rsmall) THEN 
    1774       Hen1t2 = 1._wp 
     1746      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 
     1747      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 
     1748      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 
     1749      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 
     1750      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 
     1751      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 
     1752      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 
     1753      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 
     1754      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 
     1755      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 
     1756      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 
     1757      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 
     1758      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 
     1759      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 
     1760      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 
     1761      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 
     1762      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 
     1763      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 
     1764      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 
     1765      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 
     1766      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 
     1767      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 
     1768 
     1769      IF (-zIIn1t2>=rsmall) THEN 
     1770      zHen1t2 = 1._wp 
    17751771      ELSE 
    1776       Hen1t2 = 0._wp 
    1777       ENDIF 
    1778  
    1779       IF (-IIn2t1>=rsmall) THEN 
    1780       Hen2t1 = 1._wp 
     1772      zHen1t2 = 0._wp 
     1773      ENDIF 
     1774 
     1775      IF (-zIIn2t1>=rsmall) THEN 
     1776      zHen2t1 = 1._wp 
    17811777      ELSE 
    1782       Hen2t1 = 0._wp 
    1783       ENDIF 
    1784  
    1785       s22ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22) 
     1778      zHen2t1 = 0._wp 
     1779      ENDIF 
     1780 
     1781      s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 
    17861782 
    17871783   END FUNCTION s22ks 
Note: See TracChangeset for help on using the changeset viewer.