Changeset 11933
- Timestamp:
- 2019-11-19T20:10:30+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
r11920 r11933 13 13 !! 3.7 ! 2017 (C. Rousset) add aEVP (Kimmritz 2016-2017) 14 14 !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] 15 !! ! 2019 (S. Rynders) change into eap rheology from 16 !! CICE code (Tsamados, Heorton) 15 17 !!---------------------------------------------------------------------- 16 18 #if defined key_si3 … … 46 48 PUBLIC ice_dyn_rhg_eap ! called by icedyn_rhg.F90 47 49 PUBLIC rhg_eap_rst ! called by icedyn_rhg.F90 50 51 REAL(wp), PARAMETER :: phi = 3.141592653589793_wp/12._wp ! diamond shaped floe smaller angle (default phi = 30 deg) 48 52 49 53 ! look-up table for calculating structure tensor … … 126 130 REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling 127 131 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 128 REAL(wp) :: zalph1, z1_alph1 , zalph2, z1_alph2! alpha coef from Bouillon 2009 or Kimmritz 2017132 REAL(wp) :: zalph1, z1_alph1 ! alpha coef from Bouillon 2009 or Kimmritz 2017 129 133 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV, zvU, zvV ! ice/snow mass and volume 130 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 134 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2, zdsT ! temporary scalars 131 135 REAL(wp) :: zTauO, zTauB, zRHS, zvel ! temporary scalars 132 136 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice … … 137 141 REAL(wp) :: zfac_x, zfac_y 138 142 REAL(wp) :: zshear, zdum1, zdum2 143 REAL(wp) :: stressptmp, stressmtmp ! anisotropic stress tensor components 144 REAL(wp) :: alphar, alphas ! for mechanical redistribution 139 145 ! 146 REAL(wp), DIMENSION(jpi,jpj) :: stress12tmp ! anisotropic stress tensor component for regridding 140 147 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 141 148 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 … … 237 244 IF( .NOT. ln_aEVP ) THEN 238 245 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 239 zalph2 = zalph1 * z1_ecc2240 241 246 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 242 z1_alph2 = 1._wp / ( zalph2 + 1._wp )243 247 ENDIF 244 248 … … 376 380 DO ji = 2, jpi ! no vector loop 377 381 382 ! shear at T points 383 zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 384 & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 385 & ) * 0.25_wp * r1_e1e2t(ji,jj) 386 378 387 ! shear**2 at T points (doc eq. A16) 379 388 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & … … 393 402 zdt2 = zdt * zdt 394 403 404 ! --- anisotropic stress calculation --- ! 405 CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, & 406 zdsT, aniso_11(ji,jj), aniso_12(ji,jj), & 407 stressptmp, stressmtmp, & 408 stress12tmp(ji,jj), strength(ji,jj), & 409 alphar, alphas) 410 395 411 ! delta at T points 396 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )412 !zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 397 413 398 414 ! P/delta at T points 399 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )415 !zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 400 416 401 417 ! alpha & beta for aEVP 402 418 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 403 419 ! alpha = beta = sqrt(4*gamma) 404 IF( ln_aEVP ) THEN 405 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 406 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 407 zalph2 = zalph1 408 z1_alph2 = z1_alph1 409 ENDIF 420 !IF( ln_aEVP ) THEN 421 ! zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 422 ! z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 423 !ENDIF 410 424 411 425 ! stress at T points (zkt/=0 if landfast) 412 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 413 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 426 !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 427 !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 428 zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1 ) * z1_alph1 429 zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1 ) * z1_alph1 414 430 415 431 END DO 416 432 END DO 417 CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 433 !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 434 CALL lbc_lnk( 'icedyn_rhg_eap', stress12tmp, 'T', 1. ) 418 435 419 436 DO jj = 1, jpjm1 420 437 DO ji = 1, jpim1 438 ! stress12tmp at F points 439 zdsT = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) & 440 & + stress12tmp(ji,jj ) * e1e2t(ji,jj ) + stress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) & 441 & ) * 0.25_wp * r1_e1e2f(ji,jj) 421 442 422 443 ! alpha & beta for aEVP 423 IF( ln_aEVP ) THEN 424 zalph2 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 425 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 426 zbeta(ji,jj) = zalph2 427 ENDIF 428 429 ! P/delta at F points 430 zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 444 !IF( ln_aEVP ) THEN 445 ! zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 446 ! z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 447 !ENDIF 431 448 432 449 ! stress at F points (zkt/=0 if landfast) 433 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 450 !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 451 zs12(ji,jj) = ( zs12(ji,jj) + stress12tmp(ji,jj) * zalph1 ) * z1_alph1 434 452 435 453 END DO … … 867 885 END SUBROUTINE ice_dyn_rhg_eap 868 886 887 SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu, tension, & 888 shear, a11, a12, & 889 stressp, stressm, & 890 stress12, strength, & 891 alphar, alphas) 892 !!--------------------------------------------------------------------- 893 !! *** SUBROUTINE rhg_eap_rst *** 894 !! Updates the stress depending on values of strain rate and structure 895 !! tensor and for ksub=ndte it computes closing and sliding rate 896 !!--------------------------------------------------------------------- 897 INTEGER, INTENT(in ) :: ksub, ndte 898 REAL(wp), INTENT(in ) :: phi, strength 899 REAL(wp), INTENT(in ) :: divu, tension, shear 900 REAL(wp), INTENT(in ) :: a11, a12 901 REAL(wp), INTENT( out) :: stressp, stressm, stress12 902 REAL(wp), INTENT( out) :: alphar, alphas 903 904 INTEGER :: kx ,ky, ka 905 906 REAL(wp) :: stemp11r, stemp12r, stemp22r 907 REAL(wp) :: stemp11s, stemp12s, stemp22s 908 REAL(wp) :: a22, Qd11Qd11, Qd11Qd12, Qd12Qd12 909 REAL(wp) :: Q11Q11, Q11Q12, Q12Q12 910 REAL(wp) :: dtemp11, dtemp12, dtemp22 911 REAL(wp) :: rotstemp11r, rotstemp12r, rotstemp22r 912 REAL(wp) :: rotstemp11s, rotstemp12s, rotstemp22s 913 REAL(wp) :: sig11, sig12, sig22 914 REAL(wp) :: sgprm11, sgprm12, sgprm22 915 REAL(wp) :: invstressconviso 916 REAL(wp) :: Angle_denom_gamma, Angle_denom_alpha 917 REAL(wp) :: Tany_1, Tany_2 918 REAL(wp) :: x, y, dx, dy, da, kxw, kyw, kaw 919 REAL(wp) :: invdx, invdy, invda 920 REAL(wp) :: dtemp1, dtemp2, atempprime, a, invsin 921 922 REAL(wp), PARAMETER :: kfriction = 0.45_wp 923 !!--------------------------------------------------------------------- 924 ! 1) structure tensor 925 a22 = 1._wp-a11 926 Q11Q11 = 1._wp 927 Q12Q12 = rsmall 928 Q11Q12 = rsmall 929 930 ! gamma: angle between general coordiantes and principal axis of A 931 ! here Tan2gamma = 2 a12 / (a11 - a22) 932 if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then 933 Angle_denom_gamma = 1._wp/sqrt( ( a11 - a22 )*( a11 - a22) + & 934 4._wp*a12*a12 ) 935 936 Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma !Cos^2 937 Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma !Sin^2 938 Q11Q12 = a12*Angle_denom_gamma !CosSin 939 endif 940 941 ! rotation Q*atemp*Q^T 942 atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22 943 944 ! make first principal value the largest 945 atempprime = max(atempprime, 1.0_wp - atempprime) 946 947 ! 2) strain rate 948 dtemp11 = 0.5_wp*(divu + tension) 949 dtemp12 = shear*0.5_wp 950 dtemp22 = 0.5_wp*(divu - tension) 951 952 ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 953 954 Qd11Qd11 = 1.0_wp 955 Qd12Qd12 = rsmall 956 Qd11Qd12 = rsmall 957 958 if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then 959 960 Angle_denom_alpha = 1.0_wp/sqrt( ( dtemp11 - dtemp22 )* & 961 ( dtemp11 - dtemp22 ) + 4.0_wp*dtemp12*dtemp12) 962 963 Qd11Qd11 = 0.5_wp + ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Cos^2 964 Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^2 965 Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin 966 endif 967 968 dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp22 969 dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp22 970 ! In cos and sin values 971 x = 0._wp 972 if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then 973 x = atan2(dtemp2,dtemp1) 974 endif 975 976 ! to ensure the angle lies between pi/4 and 9 pi/4 977 if (x < rpi*0.25_wp) x = x + rpi*2.0_wp 978 979 ! y: angle between major principal axis of strain rate and structure 980 ! tensor 981 ! y = gamma - alpha 982 ! Expressesed componently with 983 ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 984 985 Tany_1 = Q11Q12 - Qd11Qd12 986 Tany_2 = Q11Q11 - Qd12Qd12 987 988 y = 0._wp 989 990 if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then 991 y = atan2(Tany_1,Tany_2) 992 endif 993 994 ! to make sure y is between 0 and pi 995 if (y > rpi) y = y - rpi 996 if (y < 0) y = y + rpi 997 998 ! 3) update anisotropic stress tensor 999 dx = rpi/real(nx_yield-1,kind=wp) 1000 dy = rpi/real(ny_yield-1,kind=wp) 1001 da = 0.5_wp/real(na_yield-1,kind=wp) 1002 invdx = 1._wp/dx 1003 invdy = 1._wp/dy 1004 invda = 1._wp/da 1005 1006 ! % need 8 coords and 8 weights 1007 ! % range in kx 1008 kx = int((x-rpi*0.25_wp-rpi)*invdx) + 1 1009 kxw = kx - (x-rpi*0.25_wp-rpi)*invdx 1010 1011 ky = int(y*invdy) + 1 1012 kyw = ky - y*invdy 1013 1014 ka = int((atempprime-0.5_wp)*invda) + 1 1015 kaw = ka - (atempprime-0.5_wp)*invda 1016 1017 ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 1018 stemp11r = kxw * kyw * kaw * s11r(kx ,ky ,ka ) & 1019 & + (1._wp-kxw) * kyw * kaw * s11r(kx+1,ky ,ka ) & 1020 & + kxw * (1._wp-kyw) * kaw * s11r(kx ,ky+1,ka ) & 1021 & + kxw * kyw * (1._wp-kaw) * s11r(kx ,ky ,ka+1) & 1022 & + (1._wp-kxw) * (1._wp-kyw) * kaw * s11r(kx+1,ky+1,ka ) & 1023 & + (1._wp-kxw) * kyw * (1._wp-kaw) * s11r(kx+1,ky ,ka+1) & 1024 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s11r(kx ,ky+1,ka+1) & 1025 & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 1026 stemp12r = kxw * kyw * kaw * s12r(kx ,ky ,ka ) & 1027 & + (1._wp-kxw) * kyw * kaw * s12r(kx+1,ky ,ka ) & 1028 & + kxw * (1._wp-kyw) * kaw * s12r(kx ,ky+1,ka ) & 1029 & + kxw * kyw * (1._wp-kaw) * s12r(kx ,ky ,ka+1) & 1030 & + (1._wp-kxw) * (1._wp-kyw) * kaw * s12r(kx+1,ky+1,ka ) & 1031 & + (1._wp-kxw) * kyw * (1._wp-kaw) * s12r(kx+1,ky ,ka+1) & 1032 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s12r(kx ,ky+1,ka+1) & 1033 & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 1034 stemp22r = kxw * kyw * kaw * s22r(kx ,ky ,ka ) & 1035 & + (1._wp-kxw) * kyw * kaw * s22r(kx+1,ky ,ka ) & 1036 & + kxw * (1._wp-kyw) * kaw * s22r(kx ,ky+1,ka ) & 1037 & + kxw * kyw * (1._wp-kaw) * s22r(kx ,ky ,ka+1) & 1038 & + (1._wp-kxw) * (1._wp-kyw) * kaw * s22r(kx+1,ky+1,ka ) & 1039 & + (1._wp-kxw) * kyw * (1._wp-kaw) * s22r(kx+1,ky ,ka+1) & 1040 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) & 1041 & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 1042 1043 stemp11s = kxw * kyw * kaw * s11s(kx ,ky ,ka ) & 1044 & + (1._wp-kxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & 1045 & + kxw * (1._wp-kyw) * kaw * s11s(kx ,ky+1,ka ) & 1046 & + kxw * kyw * (1._wp-kaw) * s11s(kx ,ky ,ka+1) & 1047 & + (1._wp-kxw) * (1._wp-kyw) * kaw * s11s(kx+1,ky+1,ka ) & 1048 & + (1._wp-kxw) * kyw * (1._wp-kaw) * s11s(kx+1,ky ,ka+1) & 1049 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s11s(kx ,ky+1,ka+1) & 1050 & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 1051 stemp12s = kxw * kyw * kaw * s12s(kx ,ky ,ka ) & 1052 & + (1._wp-kxw) * kyw * kaw * s12s(kx+1,ky ,ka ) & 1053 & + kxw * (1._wp-kyw) * kaw * s12s(kx ,ky+1,ka ) & 1054 & + kxw * kyw * (1._wp-kaw) * s12s(kx ,ky ,ka+1) & 1055 & + (1._wp-kxw) * (1._wp-kyw) * kaw * s12s(kx+1,ky+1,ka ) & 1056 & + (1._wp-kxw) * kyw * (1._wp-kaw) * s12s(kx+1,ky ,ka+1) & 1057 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s12s(kx ,ky+1,ka+1) & 1058 & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 1059 stemp22s = kxw * kyw * kaw * s22s(kx ,ky ,ka ) & 1060 & + (1._wp-kxw) * kyw * kaw * s22s(kx+1,ky ,ka ) & 1061 & + kxw * (1._wp-kyw) * kaw * s22s(kx ,ky+1,ka ) & 1062 & + kxw * kyw * (1._wp-kaw) * s22s(kx ,ky ,ka+1) & 1063 & + (1._wp-kxw) * (1._wp-kyw) * kaw * s22s(kx+1,ky+1,ka ) & 1064 & + (1._wp-kxw) * kyw * (1._wp-kaw) * s22s(kx+1,ky ,ka+1) & 1065 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s22s(kx ,ky+1,ka+1) & 1066 & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 1067 1068 ! Calculate mean ice stress over a collection of floes (Equation 3 in 1069 ! Tsamados 2013) 1070 1071 sig11 = strength*(stemp11r + kfriction*stemp11s) * invsin 1072 sig12 = strength*(stemp12r + kfriction*stemp12s) * invsin 1073 sig22 = strength*(stemp22r + kfriction*stemp22s) * invsin 1074 1075 ! Back - rotation of the stress from principal axes into general coordinates 1076 1077 ! Update stress 1078 sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 - 2._wp*Q11Q12 *sig12 1079 sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12 1080 sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 + 2._wp*Q11Q12 *sig12 1081 1082 stressp = sgprm11 + sgprm22 1083 stress12 = sgprm12 1084 stressm = sgprm11 - sgprm22 1085 1086 ! Compute ridging and sliding functions in general coordinates 1087 ! (Equation 11 in Tsamados 2013) 1088 if (ksub == ndte) then 1089 rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12* stemp12r & 1090 + Q12Q12*stemp22r 1091 rotstemp12r = Q11Q11*stemp12r + Q11Q12*(stemp11r-stemp22r) & 1092 - Q12Q12*stemp12r 1093 rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12* stemp12r & 1094 + Q11Q11*stemp22r 1095 1096 rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12* stemp12s & 1097 + Q12Q12*stemp22s 1098 rotstemp12s = Q11Q11*stemp12s + Q11Q12*(stemp11s-stemp22s) & 1099 - Q12Q12*stemp12s 1100 rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12* stemp12s & 1101 + Q11Q11*stemp22s 1102 1103 alphar = rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 & 1104 + rotstemp22r*dtemp22 1105 alphas = rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 & 1106 + rotstemp22s*dtemp22 1107 endif 1108 END SUBROUTINE update_stress_rdg 869 1109 870 1110 SUBROUTINE rhg_eap_rst( cdrw, kt ) … … 873 1113 !! 874 1114 !! ** Purpose : Read or write RHG file in restart file 875 !! 1115 !! 876 1116 !! ** Method : use of IOM library 877 1117 !!---------------------------------------------------------------------- … … 889 1129 890 1130 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)892 1131 !!---------------------------------------------------------------------- 893 1132 !
Note: See TracChangeset
for help on using the changeset viewer.
