Changeset 11920 for NEMO/branches/2019
- Timestamp:
- 2019-11-15T18:38:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
r11878 r11920 46 46 PUBLIC ice_dyn_rhg_eap ! called by icedyn_rhg.F90 47 47 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 48 55 49 56 !! * Substitutions … … 872 879 INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step 873 880 ! 874 INTEGER :: iter ! local integer881 INTEGER :: iter ! local integer 875 882 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) 876 892 !!---------------------------------------------------------------------- 877 893 ! … … 911 927 ENDIF 912 928 ! 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 913 993 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 914 994 ! ! ------------------- … … 926 1006 END SUBROUTINE rhg_eap_rst 927 1007 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 928 1404 #else 929 1405 !!---------------------------------------------------------------------- 930 1406 !! Default option Empty module NO SI3 sea-ice model 931 1407 !!---------------------------------------------------------------------- 1408 932 1409 #endif 933 934 1410 !!============================================================================== 935 1411 END MODULE icedyn_rhg_eap
Note: See TracChangeset
for help on using the changeset viewer.