Changeset 11920


Ignore:
Timestamp:
2019-11-15T18:38:36+01:00 (9 months ago)
Author:
stefryn
Message:

add calculation of the look-up table for the yield surface

File:
1 edited

Legend:

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

    r11878 r11920  
    4646   PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90 
    4747   PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90 
     48 
     49   ! look-up table for calculating structure tensor 
     50   INTEGER, PARAMETER ::   nx_yield = 41 
     51   INTEGER, PARAMETER ::   ny_yield = 41 
     52   INTEGER, PARAMETER ::   na_yield = 21 
     53 
     54   REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) ::   s11r, s12r, s22r, s11s, s12s, s22s 
    4855 
    4956   !! * Substitutions 
     
    872879      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step 
    873880      ! 
    874       INTEGER  ::   iter            ! local integer 
     881      INTEGER  ::   iter                      ! local integer 
    875882      INTEGER  ::   id1, id2, id3, id4, id5   ! local integers 
     883      INTEGER  ::   ix, iy, ip, iz, n, ia     ! local integers 
     884     
     885      INTEGER, PARAMETER            ::    nz = 100 
     886       
     887      REAL(wp) ::   ainit, xinit, yinit, pinit, zinit 
     888      REAL(wp) ::   da, dx, dy, dp, dz, a1 
     889 
     890      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp 
     891      REAL(wp), PARAMETER           ::   phi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg) 
    876892      !!---------------------------------------------------------------------- 
    877893      ! 
     
    911927         ENDIF 
    912928         ! 
     929 
     930      da = 0.5_wp/real(na_yield-1,kind=wp) 
     931      ainit = 0.5_wp - da 
     932      dx = rpi/real(nx_yield-1,kind=wp) 
     933      xinit = rpi + 0.25_wp*rpi - dx 
     934      dz = rpi/real(nz,kind=wp) 
     935      zinit = -rpi*0.5_wp 
     936      dy = rpi/real(ny_yield-1,kind=wp) 
     937      yinit = -dy 
     938 
     939      DO ia=1,na_yield 
     940       DO ix=1,nx_yield 
     941        DO iy=1,ny_yield 
     942         s11r(ix,iy,ia) = 0._wp 
     943         s12r(ix,iy,ia) = 0._wp 
     944         s22r(ix,iy,ia) = 0._wp 
     945         s11s(ix,iy,ia) = 0._wp 
     946         s12s(ix,iy,ia) = 0._wp 
     947         s22s(ix,iy,ia) = 0._wp 
     948         IF (ia <= na_yield-1) then 
     949          DO iz=1,nz 
     950          s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
     951           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
     952           s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     953          s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
     954           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
     955           s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     956          s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
     957           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
     958           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)  
     959          s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
     960           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
     961           s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     962          s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
     963           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
     964           s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     965          s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
     966           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 
     967           s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi) 
     968          ENDDO 
     969          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 
     970          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 
     971          IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 
     972          IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 
     973          IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 
     974          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 
     975         ELSE 
     976          s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     977          s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     978          s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     979          s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     980          s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     981          s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi) 
     982          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 
     983          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp 
     984          IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp 
     985          IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp 
     986          IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp 
     987          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 
     988         ENDIF 
     989        ENDDO 
     990       ENDDO 
     991      ENDDO 
     992 
    913993      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    914994         !                                   ! ------------------- 
     
    9261006   END SUBROUTINE rhg_eap_rst 
    9271007 
     1008 
     1009   FUNCTION w1(a) 
     1010      !!------------------------------------------------------------------- 
     1011      !! Function : w1 (see Gaussian function psi in Tsamados et al 2013) 
     1012      !!------------------------------------------------------------------- 
     1013      REAL(wp), INTENT(in   ) ::   a 
     1014      REAL(wp) ::   w1 
     1015      !!------------------------------------------------------------------- 
     1016    
     1017      w1 = - 223.87569446_wp & 
     1018           + 2361.2198663_wp*a & 
     1019           - 10606.56079975_wp*a*a & 
     1020           + 26315.50025642_wp*a*a*a & 
     1021           - 38948.30444297_wp*a*a*a*a & 
     1022           + 34397.72407466_wp*a*a*a*a*a & 
     1023           - 16789.98003081_wp*a*a*a*a*a*a & 
     1024           + 3495.82839237_wp*a*a*a*a*a*a*a 
     1025    
     1026   END FUNCTION w1 
     1027 
     1028 
     1029   FUNCTION w2(a) 
     1030      !!------------------------------------------------------------------- 
     1031      !! Function : w2 (see Gaussian function psi in Tsamados et al 2013) 
     1032      !!------------------------------------------------------------------- 
     1033      REAL(wp), INTENT(in   ) ::   a 
     1034      REAL(wp) ::   w2 
     1035      !!------------------------------------------------------------------- 
     1036    
     1037      w2 = - 6670.68911883_wp & 
     1038           + 70222.33061536_wp*a & 
     1039           - 314871.71525448_wp*a*a & 
     1040           + 779570.02793492_wp*a*a*a & 
     1041           - 1151098.82436864_wp*a*a*a*a & 
     1042           + 1013896.59464498_wp*a*a*a*a*a & 
     1043           - 493379.44906738_wp*a*a*a*a*a*a & 
     1044           + 102356.551518_wp*a*a*a*a*a*a*a 
     1045    
     1046   END FUNCTION w2 
     1047 
     1048   FUNCTION s11kr(x,y,z,phi)  
     1049      !!------------------------------------------------------------------- 
     1050      !! Function : s11kr 
     1051      !!------------------------------------------------------------------- 
     1052      REAL(wp), INTENT(in   ) ::   x,y,z,phi 
     1053    
     1054      REAL(wp) ::   s11kr, p, pih 
     1055    
     1056      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
     1057      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
     1058      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
     1059      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
     1060      REAL(wp) ::   d11, d12, d22 
     1061      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
     1062      REAL(wp) ::   Hen1t2, Hen2t1 
     1063      !!------------------------------------------------------------------- 
     1064    
     1065      pih = 0.5_wp*rpi 
     1066      p = phi 
     1067    
     1068      n1t2i11 = cos(z+pih-p) * cos(z+p) 
     1069      n1t2i12 = cos(z+pih-p) * sin(z+p) 
     1070      n1t2i21 = sin(z+pih-p) * cos(z+p) 
     1071      n1t2i22 = sin(z+pih-p) * sin(z+p) 
     1072      n2t1i11 = cos(z-pih+p) * cos(z-p) 
     1073      n2t1i12 = cos(z-pih+p) * sin(z-p) 
     1074      n2t1i21 = sin(z-pih+p) * cos(z-p) 
     1075      n2t1i22 = sin(z-pih+p) * sin(z-p) 
     1076      t1t2i11 = cos(z-p) * cos(z+p) 
     1077      t1t2i12 = cos(z-p) * sin(z+p) 
     1078      t1t2i21 = sin(z-p) * cos(z+p) 
     1079      t1t2i22 = sin(z-p) * sin(z+p) 
     1080      t2t1i11 = cos(z+p) * cos(z-p) 
     1081      t2t1i12 = cos(z+p) * sin(z-p) 
     1082      t2t1i21 = sin(z+p) * cos(z-p) 
     1083      t2t1i22 = sin(z+p) * sin(z-p) 
     1084   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 
     1085   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 
     1086   ! x=x-pi/2 
     1087      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
     1088      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
     1089      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
     1090      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
     1091      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
     1092      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1093    
     1094      IF (-IIn1t2>=rsmall) then 
     1095      Hen1t2 = 1._wp 
     1096      ELSE 
     1097      Hen1t2 = 0._wp 
     1098      ENDIF 
     1099    
     1100      IF (-IIn2t1>=rsmall) then 
     1101      Hen2t1 = 1._wp 
     1102      ELSE 
     1103      Hen2t1 = 0._wp 
     1104      ENDIF 
     1105    
     1106      s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 * n2t1i11) 
     1107 
     1108   END FUNCTION s11kr 
     1109 
     1110   FUNCTION s12kr(x,y,z,phi) 
     1111      !!------------------------------------------------------------------- 
     1112      !! Function : s12kr 
     1113      !!------------------------------------------------------------------- 
     1114      REAL(wp), INTENT(in   ) ::   x,y,z,phi 
     1115 
     1116      REAL(wp) ::   s12kr, s12r0, s21r0, p, pih 
     1117 
     1118      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
     1119      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
     1120      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
     1121      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
     1122      REAL(wp) ::   d11, d12, d22 
     1123      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
     1124      REAL(wp) ::   Hen1t2, Hen2t1 
     1125      !!------------------------------------------------------------------- 
     1126      pih = 0.5_wp*rpi 
     1127      p = phi 
     1128    
     1129      n1t2i11 = cos(z+pih-p) * cos(z+p) 
     1130      n1t2i12 = cos(z+pih-p) * sin(z+p) 
     1131      n1t2i21 = sin(z+pih-p) * cos(z+p) 
     1132      n1t2i22 = sin(z+pih-p) * sin(z+p) 
     1133      n2t1i11 = cos(z-pih+p) * cos(z-p) 
     1134      n2t1i12 = cos(z-pih+p) * sin(z-p) 
     1135      n2t1i21 = sin(z-pih+p) * cos(z-p) 
     1136      n2t1i22 = sin(z-pih+p) * sin(z-p) 
     1137      t1t2i11 = cos(z-p) * cos(z+p) 
     1138      t1t2i12 = cos(z-p) * sin(z+p) 
     1139      t1t2i21 = sin(z-p) * cos(z+p) 
     1140      t1t2i22 = sin(z-p) * sin(z+p) 
     1141      t2t1i11 = cos(z+p) * cos(z-p) 
     1142      t2t1i12 = cos(z+p) * sin(z-p) 
     1143      t2t1i21 = sin(z+p) * cos(z-p) 
     1144      t2t1i22 = sin(z+p) * sin(z-p) 
     1145      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
     1146      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
     1147      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
     1148      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
     1149      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
     1150      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1151    
     1152      IF (-IIn1t2>=rsmall) then 
     1153      Hen1t2 = 1._wp 
     1154      ELSE 
     1155      Hen1t2 = 0._wp 
     1156      ENDIF 
     1157    
     1158      IF (-IIn2t1>=rsmall) then 
     1159      Hen2t1 = 1._wp 
     1160      ELSE 
     1161      Hen2t1 = 0._wp 
     1162      ENDIF 
     1163    
     1164      s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 * n2t1i12) 
     1165      s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 * n2t1i21) 
     1166      s12kr=0.5_wp*(s12r0+s21r0) 
     1167    
     1168   END FUNCTION s12kr 
     1169 
     1170   FUNCTION s22kr(x,y,z,phi) 
     1171      !!------------------------------------------------------------------- 
     1172      !! Function : s22kr 
     1173      !!------------------------------------------------------------------- 
     1174      REAL(wp), INTENT(in   ) ::   x,y,z,phi 
     1175 
     1176      REAL(wp) ::   s22kr, p, pih 
     1177 
     1178      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
     1179      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
     1180      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
     1181      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
     1182      REAL(wp) ::   d11, d12, d22 
     1183      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
     1184      REAL(wp) ::   Hen1t2, Hen2t1 
     1185      !!------------------------------------------------------------------- 
     1186      pih = 0.5_wp*rpi 
     1187      p = phi 
     1188 
     1189      n1t2i11 = cos(z+pih-p) * cos(z+p) 
     1190      n1t2i12 = cos(z+pih-p) * sin(z+p) 
     1191      n1t2i21 = sin(z+pih-p) * cos(z+p) 
     1192      n1t2i22 = sin(z+pih-p) * sin(z+p) 
     1193      n2t1i11 = cos(z-pih+p) * cos(z-p) 
     1194      n2t1i12 = cos(z-pih+p) * sin(z-p) 
     1195      n2t1i21 = sin(z-pih+p) * cos(z-p) 
     1196      n2t1i22 = sin(z-pih+p) * sin(z-p) 
     1197      t1t2i11 = cos(z-p) * cos(z+p) 
     1198      t1t2i12 = cos(z-p) * sin(z+p) 
     1199      t1t2i21 = sin(z-p) * cos(z+p) 
     1200      t1t2i22 = sin(z-p) * sin(z+p) 
     1201      t2t1i11 = cos(z+p) * cos(z-p) 
     1202      t2t1i12 = cos(z+p) * sin(z-p) 
     1203      t2t1i21 = sin(z+p) * cos(z-p) 
     1204      t2t1i22 = sin(z+p) * sin(z-p) 
     1205      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
     1206      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
     1207      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
     1208      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
     1209      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
     1210      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1211 
     1212      IF (-IIn1t2>=rsmall) then 
     1213      Hen1t2 = 1._wp 
     1214      ELSE 
     1215      Hen1t2 = 0._wp 
     1216      ENDIF 
     1217 
     1218      IF (-IIn2t1>=rsmall) then 
     1219      Hen2t1 = 1._wp 
     1220      ELSE 
     1221      Hen2t1 = 0._wp 
     1222      ENDIF 
     1223 
     1224      s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 * n2t1i22) 
     1225 
     1226   END FUNCTION s22kr 
     1227 
     1228   FUNCTION s11ks(x,y,z,phi) 
     1229      !!------------------------------------------------------------------- 
     1230      !! Function : s11ks 
     1231      !!------------------------------------------------------------------- 
     1232      REAL(wp), INTENT(in   ) ::   x,y,z,phi 
     1233 
     1234      REAL(wp) ::   s11ks, p, pih 
     1235 
     1236      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
     1237      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
     1238      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
     1239      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
     1240      REAL(wp) ::   d11, d12, d22 
     1241      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
     1242      REAL(wp) ::   Hen1t2, Hen2t1 
     1243      !!------------------------------------------------------------------- 
     1244      pih = 0.5_wp*rpi 
     1245      p = phi 
     1246 
     1247      n1t2i11 = cos(z+pih-p) * cos(z+p) 
     1248      n1t2i12 = cos(z+pih-p) * sin(z+p) 
     1249      n1t2i21 = sin(z+pih-p) * cos(z+p) 
     1250      n1t2i22 = sin(z+pih-p) * sin(z+p) 
     1251      n2t1i11 = cos(z-pih+p) * cos(z-p) 
     1252      n2t1i12 = cos(z-pih+p) * sin(z-p) 
     1253      n2t1i21 = sin(z-pih+p) * cos(z-p) 
     1254      n2t1i22 = sin(z-pih+p) * sin(z-p) 
     1255      t1t2i11 = cos(z-p) * cos(z+p) 
     1256      t1t2i12 = cos(z-p) * sin(z+p) 
     1257      t1t2i21 = sin(z-p) * cos(z+p) 
     1258      t1t2i22 = sin(z-p) * sin(z+p) 
     1259      t2t1i11 = cos(z+p) * cos(z-p) 
     1260      t2t1i12 = cos(z+p) * sin(z-p) 
     1261      t2t1i21 = sin(z+p) * cos(z-p) 
     1262      t2t1i22 = sin(z+p) * sin(z-p) 
     1263      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
     1264      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
     1265      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
     1266      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
     1267      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
     1268      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1269 
     1270      IF (-IIn1t2>=rsmall) then 
     1271      Hen1t2 = 1._wp 
     1272      ELSE 
     1273      Hen1t2 = 0._wp 
     1274      ENDIF 
     1275 
     1276      IF (-IIn2t1>=rsmall) then 
     1277      Hen2t1 = 1._wp 
     1278      ELSE 
     1279      Hen2t1 = 0._wp 
     1280      ENDIF 
     1281 
     1282      s11ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i11 + Hen2t1 * t2t1i11) 
     1283 
     1284   END FUNCTION s11ks 
     1285 
     1286   FUNCTION s12ks(x,y,z,phi) 
     1287      !!------------------------------------------------------------------- 
     1288      !! Function : s12ks 
     1289      !!------------------------------------------------------------------- 
     1290      REAL(wp), INTENT(in   ) ::   x,y,z,phi 
     1291 
     1292      REAL(wp) ::   s12ks, s12s0, s21s0, p, pih 
     1293 
     1294      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
     1295      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
     1296      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
     1297      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
     1298      REAL(wp) ::   d11, d12, d22 
     1299      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
     1300      REAL(wp) ::   Hen1t2, Hen2t1 
     1301      !!------------------------------------------------------------------- 
     1302      pih = 0.5_wp*rpi 
     1303      p =phi 
     1304 
     1305      n1t2i11 = cos(z+pih-p) * cos(z+p) 
     1306      n1t2i12 = cos(z+pih-p) * sin(z+p) 
     1307      n1t2i21 = sin(z+pih-p) * cos(z+p) 
     1308      n1t2i22 = sin(z+pih-p) * sin(z+p) 
     1309      n2t1i11 = cos(z-pih+p) * cos(z-p) 
     1310      n2t1i12 = cos(z-pih+p) * sin(z-p) 
     1311      n2t1i21 = sin(z-pih+p) * cos(z-p) 
     1312      n2t1i22 = sin(z-pih+p) * sin(z-p) 
     1313      t1t2i11 = cos(z-p) * cos(z+p) 
     1314      t1t2i12 = cos(z-p) * sin(z+p) 
     1315      t1t2i21 = sin(z-p) * cos(z+p) 
     1316      t1t2i22 = sin(z-p) * sin(z+p) 
     1317      t2t1i11 = cos(z+p) * cos(z-p) 
     1318      t2t1i12 = cos(z+p) * sin(z-p) 
     1319      t2t1i21 = sin(z+p) * cos(z-p) 
     1320      t2t1i22 = sin(z+p) * sin(z-p) 
     1321      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
     1322      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
     1323      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
     1324      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
     1325      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
     1326      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1327 
     1328      IF (-IIn1t2>=rsmall) then 
     1329      Hen1t2 = 1._wp 
     1330      ELSE 
     1331      Hen1t2 = 0._wp 
     1332      ENDIF 
     1333 
     1334      IF (-IIn2t1>=rsmall) then 
     1335      Hen2t1 = 1._wp 
     1336      ELSE 
     1337      Hen2t1 = 0._wp 
     1338      ENDIF 
     1339 
     1340      s12s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i12 + Hen2t1 * t2t1i12) 
     1341      s21s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i21 + Hen2t1 * t2t1i21) 
     1342      s12ks=0.5_wp*(s12s0+s21s0) 
     1343 
     1344   END FUNCTION s12ks 
     1345 
     1346   FUNCTION s22ks(x,y,z,phi)  
     1347      !!------------------------------------------------------------------- 
     1348      !! Function : s22ks 
     1349      !!------------------------------------------------------------------- 
     1350      REAL(wp), INTENT(in   ) ::   x,y,z,phi 
     1351 
     1352      REAL(wp) ::   s22ks, p, pih 
     1353 
     1354      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22 
     1355      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22 
     1356      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22 
     1357      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22 
     1358      REAL(wp) ::   d11, d12, d22 
     1359      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2 
     1360      REAL(wp) ::   Hen1t2, Hen2t1 
     1361      !!------------------------------------------------------------------- 
     1362      pih = 0.5_wp*rpi 
     1363      p =phi 
     1364 
     1365      n1t2i11 = cos(z+pih-p) * cos(z+p) 
     1366      n1t2i12 = cos(z+pih-p) * sin(z+p) 
     1367      n1t2i21 = sin(z+pih-p) * cos(z+p) 
     1368      n1t2i22 = sin(z+pih-p) * sin(z+p) 
     1369      n2t1i11 = cos(z-pih+p) * cos(z-p) 
     1370      n2t1i12 = cos(z-pih+p) * sin(z-p) 
     1371      n2t1i21 = sin(z-pih+p) * cos(z-p) 
     1372      n2t1i22 = sin(z-pih+p) * sin(z-p) 
     1373      t1t2i11 = cos(z-p) * cos(z+p) 
     1374      t1t2i12 = cos(z-p) * sin(z+p) 
     1375      t1t2i21 = sin(z-p) * cos(z+p) 
     1376      t1t2i22 = sin(z-p) * sin(z+p) 
     1377      t2t1i11 = cos(z+p) * cos(z-p) 
     1378      t2t1i12 = cos(z+p) * sin(z-p) 
     1379      t2t1i21 = sin(z+p) * cos(z-p) 
     1380      t2t1i22 = sin(z+p) * sin(z-p) 
     1381      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 
     1382      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 
     1383      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 
     1384      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 
     1385      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 
     1386      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
     1387 
     1388      IF (-IIn1t2>=rsmall) THEN 
     1389      Hen1t2 = 1._wp 
     1390      ELSE 
     1391      Hen1t2 = 0._wp 
     1392      ENDIF 
     1393 
     1394      IF (-IIn2t1>=rsmall) THEN 
     1395      Hen2t1 = 1._wp 
     1396      ELSE 
     1397      Hen2t1 = 0._wp 
     1398      ENDIF 
     1399 
     1400      s22ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22) 
     1401 
     1402   END FUNCTION s22ks 
     1403 
    9281404#else 
    9291405   !!---------------------------------------------------------------------- 
    9301406   !!   Default option         Empty module           NO SI3 sea-ice model 
    9311407   !!---------------------------------------------------------------------- 
     1408 
    9321409#endif 
    933  
    9341410   !!============================================================================== 
    9351411END MODULE icedyn_rhg_eap 
Note: See TracChangeset for help on using the changeset viewer.