Changeset 14117 for NEMO/trunk/src/ICE
- Timestamp:
- 2020-12-07T11:21:42+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_rhg_eap.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_rhg_eap.F90
r14072 r14117 16 16 !! CICE code (Tsamados, Heorton) 17 17 !!---------------------------------------------------------------------- 18 #if defined key_si3 && ! defined key_agrif18 #if defined key_si3 19 19 !!---------------------------------------------------------------------- 20 20 !! 'key_si3' SI3 sea-ice model … … 66 66 INTEGER :: ncvgid ! netcdf file id 67 67 INTEGER :: nvarid ! netcdf variable id 68 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 68 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aimsk00 69 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: eap_res , aimsk15 69 70 !!---------------------------------------------------------------------- 70 71 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 202 203 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' 203 204 ! 204 ! for diagnostics and convergence tests 205 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 205 IF( kt == nit000 ) THEN 206 ! 207 ! for diagnostics 208 ALLOCATE( aimsk00(jpi,jpj) ) 209 ! for convergence tests 210 IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) ) 211 ENDIF 212 ! 206 213 DO_2D( 1, 1, 1, 1 ) 207 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 208 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 214 aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 209 215 END_2D 216 IF( nn_rhg_chkcvg > 0 ) THEN 217 DO_2D( 1, 1, 1, 1 ) 218 aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 219 END_2D 220 ENDIF 210 221 ! 211 222 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... … … 349 360 zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 350 361 ! ice-bottom stress at U points 351 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 362 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 352 363 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 353 364 ! ice-bottom stress at V points 354 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 365 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 355 366 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 356 367 ! ice_bottom stress at T points 357 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 368 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 358 369 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 359 370 END_2D … … 819 830 & ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) 820 831 ! 821 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )822 CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )823 CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )824 CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )825 CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )826 CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )832 CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 ) 833 CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 ) 834 CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 ) 835 CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 ) 836 CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 ) 837 CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 ) 827 838 ENDIF 828 839 829 840 ! --- divergence, shear and strength --- ! 830 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence831 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear832 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta833 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength841 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * aimsk00 ) ! divergence 842 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * aimsk00 ) ! shear 843 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * aimsk00 ) ! delta 844 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * aimsk00 ) ! strength 834 845 835 846 ! --- Stress tensor invariants (SIMIP diags) --- ! … … 856 867 ! 857 868 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 858 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress859 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress869 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * aimsk00(:,:) ) ! Normal stress 870 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress 860 871 861 872 DEALLOCATE ( zsig_I, zsig_II ) … … 903 914 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp ) 904 915 905 CALL iom_put( 'yield11', zyield11 * zmsk00 )906 CALL iom_put( 'yield22', zyield22 * zmsk00 )907 CALL iom_put( 'yield12', zyield12 * zmsk00 )916 CALL iom_put( 'yield11', zyield11 * aimsk00 ) 917 CALL iom_put( 'yield22', zyield22 * aimsk00 ) 918 CALL iom_put( 'yield12', zyield12 * aimsk00 ) 908 919 ENDIF 909 920 … … 911 922 IF( iom_use('aniso') ) THEN 912 923 CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp ) 913 CALL iom_put( 'aniso' , paniso_11 * zmsk00 )924 CALL iom_put( 'aniso' , paniso_11 * aimsk00 ) 914 925 ENDIF 915 926 … … 922 933 & zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp ) 923 934 924 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)925 CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)926 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x)927 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y)928 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)929 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)935 CALL iom_put( 'dssh_dx' , zspgU * aimsk00 ) ! Sea-surface tilt term in force balance (x) 936 CALL iom_put( 'dssh_dy' , zspgV * aimsk00 ) ! Sea-surface tilt term in force balance (y) 937 CALL iom_put( 'corstrx' , zCorU * aimsk00 ) ! Coriolis force term in force balance (x) 938 CALL iom_put( 'corstry' , zCorV * aimsk00 ) ! Coriolis force term in force balance (y) 939 CALL iom_put( 'intstrx' , zfU * aimsk00 ) ! Internal force term in force balance (x) 940 CALL iom_put( 'intstry' , zfV * aimsk00 ) ! Internal force term in force balance (y) 930 941 ENDIF 931 942 … … 938 949 DO_2D( 0, 0, 0, 0 ) 939 950 ! 2D ice mass, snow mass, area transport arrays (X, Y) 940 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)941 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)951 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj) 952 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj) 942 953 943 954 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component … … 973 984 IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 974 985 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 975 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )986 & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) ) 976 987 ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 977 988 CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 978 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )989 & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) ) 979 990 ENDIF 980 991 ENDIF 981 992 ENDIF 982 !983 DEALLOCATE( zmsk00, zmsk15 )984 993 ! 985 994 END SUBROUTINE ice_dyn_rhg_eap … … 1006 1015 REAL(wp) :: zresm ! local real 1007 1016 CHARACTER(len=20) :: clname 1008 REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence1009 1017 !!---------------------------------------------------------------------- 1010 1018 … … 1037 1045 ELSE 1038 1046 DO_2D( 1, 1, 1, 1 ) 1039 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &1040 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)1047 eap_res(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 1048 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * aimsk15(ji,jj) 1041 1049 END_2D 1042 zresm = MAXVAL( zres ) 1050 1051 zresm = MAXVAL( eap_res ) 1043 1052 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain 1044 1053 ENDIF … … 1080 1089 REAL(wp) :: zsig11, zsig12, zsig22 1081 1090 REAL(wp) :: zsgprm11, zsgprm12, zsgprm22 1082 REAL(wp) :: zinvstressconviso1083 1091 REAL(wp) :: zAngle_denom_gamma, zAngle_denom_alpha 1084 1092 REAL(wp) :: zTany_1, zTany_2 1085 REAL(wp) :: zx, zy, z dx, zdy, zda, zkxw, kyw, kaw1093 REAL(wp) :: zx, zy, zkxw, kyw, kaw 1086 1094 REAL(wp) :: zinvdx, zinvdy, zinvda 1087 REAL(wp) :: zdtemp1, zdtemp2, zatempprime, zinvsin 1088 1089 REAL(wp), PARAMETER :: kfriction = 0.45_wp 1090 !!--------------------------------------------------------------------- 1095 REAL(wp) :: zdtemp1, zdtemp2, zatempprime 1096 1097 REAL(wp), PARAMETER :: ppkfriction = 0.45_wp 1091 1098 ! Factor to maintain the same stress as in EVP (see Section 3) 1092 1099 ! Can be set to 1 otherwise 1093 ! zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction)1094 zinvstressconviso = 1._wp1095 1096 zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso1097 !now uses phi as set in higher code1100 ! REAL(wp), PARAMETER :: ppinvstressconviso = 1._wp/(1._wp+ppkfriction*ppkfriction) 1101 REAL(wp), PARAMETER :: ppinvstressconviso = 1._wp 1102 1103 ! next statement uses pphi set in main module (icedyn_rhg_eap) 1104 REAL(wp), PARAMETER :: ppinvsin = 1._wp/sin(2._wp*pphi) * ppinvstressconviso 1098 1105 1099 1106 ! compute eigenvalues, eigenvectors and angles for structure tensor, strain … … 1175 1182 1176 1183 ! 3) update anisotropic stress tensor 1177 zdx = rpi/real(nx_yield-1,kind=wp) 1178 zdy = rpi/real(ny_yield-1,kind=wp) 1179 zda = 0.5_wp/real(na_yield-1,kind=wp) 1180 zinvdx = 1._wp/zdx 1181 zinvdy = 1._wp/zdy 1182 zinvda = 1._wp/zda 1184 zinvdx = real(nx_yield-1,kind=wp)/rpi 1185 zinvdy = real(ny_yield-1,kind=wp)/rpi 1186 zinvda = 2._wp*real(na_yield-1,kind=wp) 1183 1187 1184 1188 ! % need 8 coords and 8 weights … … 1258 1262 ! Tsamados 2013) 1259 1263 1260 zsig11 = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin1261 zsig12 = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin1262 zsig22 = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin1264 zsig11 = pstrength*(zstemp11r + ppkfriction*zstemp11s) * ppinvsin 1265 zsig12 = pstrength*(zstemp12r + ppkfriction*zstemp12s) * ppinvsin 1266 zsig22 = pstrength*(zstemp22r + ppkfriction*zstemp22s) * ppinvsin 1263 1267 1264 1268 ! Back - rotation of the stress from principal axes into general coordinates … … 1319 1323 REAL (wp) :: zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 1320 1324 1321 !!$ REAL (wp), PARAMETER ::kfrac = 0.0001_wp ! rate of fracture formation1322 REAL (wp), PARAMETER :: kfrac = 1.e-3_wp! rate of fracture formation1323 REAL (wp), PARAMETER :: threshold = 0.3_wp ! critical confinement ratio1325 !!$ REAL (wp), PARAMETER :: ppkfrac = 0.0001_wp ! rate of fracture formation 1326 REAL (wp), PARAMETER :: ppkfrac = 1.e-3_wp ! rate of fracture formation 1327 REAL (wp), PARAMETER :: ppthreshold = 0.3_wp ! critical confinement ratio 1324 1328 !!--------------------------------------------------------------- 1325 1329 ! … … 1363 1367 ! which leads to the loss of their shape, so we again model it through diffusion 1364 1368 ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp)) THEN 1365 pmresult11 = - kfrac * (pa11 - zQ12Q12)1366 pmresult12 = - kfrac * (pa12 + zQ11Q12)1369 pmresult11 = - ppkfrac * (pa11 - zQ12Q12) 1370 pmresult12 = - ppkfrac * (pa12 + zQ11Q12) 1367 1371 1368 1372 ! Shear faulting … … 1370 1374 pmresult11 = 0.0_wp 1371 1375 pmresult12 = 0.0_wp 1372 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN1373 pmresult11 = - kfrac * (pa11 - zQ12Q12)1374 pmresult12 = - kfrac * (pa12 + zQ11Q12)1376 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= ppthreshold)) THEN 1377 pmresult11 = - ppkfrac * (pa11 - zQ12Q12) 1378 pmresult12 = - ppkfrac * (pa12 + zQ11Q12) 1375 1379 1376 1380 ! Horizontal spalling … … 1405 1409 !!clem 1406 1410 REAL(wp) :: zw1, zw2, zfac, ztemp 1407 REAL(wp) :: idx, idy, idz 1411 REAL(wp) :: zidx, zidy, zidz 1412 REAL(wp) :: zsaak(6) ! temporary array 1408 1413 1409 1414 REAL(wp), PARAMETER :: eps6 = 1.0e-6_wp … … 1522 1527 zw2 = w2(ainit+ia*da) 1523 1528 DO iz = 1, nz 1524 idz = zinit+iz*dz1529 zidz = zinit+iz*dz 1525 1530 ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz)) 1526 1531 DO iy = 1, ny_yield 1527 idy = yinit+iy*dy1532 zidy = yinit+iy*dy 1528 1533 DO ix = 1, nx_yield 1529 idx = xinit+ix*dx 1530 s11r(ix,iy,ia) = s11r(ix,iy,ia) + ztemp * s11kr(idx,idy,idz)*zfac 1531 s12r(ix,iy,ia) = s12r(ix,iy,ia) + ztemp * s12kr(idx,idy,idz)*zfac 1532 s22r(ix,iy,ia) = s22r(ix,iy,ia) + ztemp * s22kr(idx,idy,idz)*zfac 1533 s11s(ix,iy,ia) = s11s(ix,iy,ia) + ztemp * s11ks(idx,idy,idz)*zfac 1534 s12s(ix,iy,ia) = s12s(ix,iy,ia) + ztemp * s12ks(idx,idy,idz)*zfac 1535 s22s(ix,iy,ia) = s22s(ix,iy,ia) + ztemp * s22ks(idx,idy,idz)*zfac 1534 zidx = xinit+ix*dx 1535 call all_skr_sks(zidx,zidy,zidz,zsaak) 1536 zsaak = ztemp*zsaak*zfac 1537 s11r(ix,iy,ia) = s11r(ix,iy,ia) + zsaak(1) 1538 s12r(ix,iy,ia) = s12r(ix,iy,ia) + zsaak(2) 1539 s22r(ix,iy,ia) = s22r(ix,iy,ia) + zsaak(3) 1540 s11s(ix,iy,ia) = s11s(ix,iy,ia) + zsaak(4) 1541 s12s(ix,iy,ia) = s12s(ix,iy,ia) + zsaak(5) 1542 s22s(ix,iy,ia) = s22s(ix,iy,ia) + zsaak(6) 1536 1543 END DO 1537 1544 END DO 1538 1545 END DO 1539 1546 END DO 1540 1541 1547 zfac = 1._wp/sin(2._wp*pphi) 1542 1548 ia = na_yield 1543 1549 DO iy = 1, ny_yield 1544 idy = yinit+iy*dy1550 zidy = yinit+iy*dy 1545 1551 DO ix = 1, nx_yield 1546 idx = xinit+ix*dx 1547 s11r(ix,iy,ia) = 0.5_wp*s11kr(idx,idy,0._wp)*zfac 1548 s12r(ix,iy,ia) = 0.5_wp*s12kr(idx,idy,0._wp)*zfac 1549 s22r(ix,iy,ia) = 0.5_wp*s22kr(idx,idy,0._wp)*zfac 1550 s11s(ix,iy,ia) = 0.5_wp*s11ks(idx,idy,0._wp)*zfac 1551 s12s(ix,iy,ia) = 0.5_wp*s12ks(idx,idy,0._wp)*zfac 1552 s22s(ix,iy,ia) = 0.5_wp*s22ks(idx,idy,0._wp)*zfac 1552 zidx = xinit+ix*dx 1553 call all_skr_sks(zidx,zidy,0._wp,zsaak) 1554 zsaak = 0.5_wp*zsaak*zfac 1555 s11r(ix,iy,ia) = zsaak(1) 1556 s12r(ix,iy,ia) = zsaak(2) 1557 s22r(ix,iy,ia) = zsaak(3) 1558 s11s(ix,iy,ia) = zsaak(4) 1559 s12s(ix,iy,ia) = zsaak(5) 1560 s22s(ix,iy,ia) = zsaak(6) 1553 1561 ENDDO 1554 1562 ENDDO … … 1616 1624 END FUNCTION w2 1617 1625 1618 FUNCTION s11kr(px,py,pz) 1619 !!------------------------------------------------------------------- 1620 !! Function : s11kr 1621 !!------------------------------------------------------------------- 1626 SUBROUTINE all_skr_sks( px, py, pz, allsk ) 1622 1627 REAL(wp), INTENT(in ) :: px,py,pz 1623 1624 REAL(wp) :: s11kr, zpih 1625 1628 REAL(wp), INTENT(out ) :: allsk(6) 1629 1630 REAL(wp) :: zs12r0, zs21r0 1631 REAL(wp) :: zs12s0, zs21s0 1632 1633 REAL(wp) :: zpih 1626 1634 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1627 1635 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 … … 1673 1681 ENDIF 1674 1682 1675 s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 1676 1677 END FUNCTION s11kr 1678 1679 FUNCTION s12kr(px,py,pz) 1683 !!------------------------------------------------------------------- 1684 !! Function : s11kr 1685 !!------------------------------------------------------------------- 1686 allsk(1) = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 1680 1687 !!------------------------------------------------------------------- 1681 1688 !! Function : s12kr 1682 1689 !!------------------------------------------------------------------- 1683 REAL(wp), INTENT(in ) :: px,py,pz1684 1685 REAL(wp) :: s12kr, zs12r0, zs21r0, zpih1686 1687 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i221688 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i221689 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i221690 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i221691 REAL(wp) :: zd11, zd12, zd221692 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t21693 REAL(wp) :: zHen1t2, zHen2t11694 !!-------------------------------------------------------------------1695 zpih = 0.5_wp*rpi1696 1697 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)1698 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)1699 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)1700 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)1701 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)1702 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)1703 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)1704 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)1705 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)1706 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)1707 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)1708 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)1709 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)1710 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)1711 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)1712 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)1713 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))1714 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))1715 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))1716 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd221717 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd221718 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd221719 1720 IF (-zIIn1t2>=rsmall) THEN1721 zHen1t2 = 1._wp1722 ELSE1723 zHen1t2 = 0._wp1724 ENDIF1725 1726 IF (-zIIn2t1>=rsmall) THEN1727 zHen2t1 = 1._wp1728 ELSE1729 zHen2t1 = 0._wp1730 ENDIF1731 1732 1690 zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 1733 1691 zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21) 1734 s12kr=0.5_wp*(zs12r0+zs21r0) 1735 1736 END FUNCTION s12kr 1737 1738 FUNCTION s22kr(px,py,pz) 1692 allsk(2)=0.5_wp*(zs12r0+zs21r0) 1739 1693 !!------------------------------------------------------------------- 1740 1694 !! Function : s22kr 1741 1695 !!------------------------------------------------------------------- 1742 REAL(wp), INTENT(in ) :: px,py,pz 1743 1744 REAL(wp) :: s22kr, zpih 1745 1746 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1747 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1748 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1749 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1750 REAL(wp) :: zd11, zd12, zd22 1751 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1752 REAL(wp) :: zHen1t2, zHen2t1 1753 !!------------------------------------------------------------------- 1754 zpih = 0.5_wp*rpi 1755 1756 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1757 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1758 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1759 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1760 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1761 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1762 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1763 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1764 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1765 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1766 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1767 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1768 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1769 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1770 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1771 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1772 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1773 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1774 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1775 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1776 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1777 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1778 1779 IF (-zIIn1t2>=rsmall) THEN 1780 zHen1t2 = 1._wp 1781 ELSE 1782 zHen1t2 = 0._wp 1783 ENDIF 1784 1785 IF (-zIIn2t1>=rsmall) THEN 1786 zHen2t1 = 1._wp 1787 ELSE 1788 zHen2t1 = 0._wp 1789 ENDIF 1790 1791 s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 1792 1793 END FUNCTION s22kr 1794 1795 FUNCTION s11ks(px,py,pz) 1696 allsk(3) = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 1796 1697 !!------------------------------------------------------------------- 1797 1698 !! Function : s11ks 1798 1699 !!------------------------------------------------------------------- 1799 REAL(wp), INTENT(in ) :: px,py,pz 1800 1801 REAL(wp) :: s11ks, zpih 1802 1803 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1804 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1805 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1806 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1807 REAL(wp) :: zd11, zd12, zd22 1808 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1809 REAL(wp) :: zHen1t2, zHen2t1 1810 !!------------------------------------------------------------------- 1811 zpih = 0.5_wp*rpi 1812 1813 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1814 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1815 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1816 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1817 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1818 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1819 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1820 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1821 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1822 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1823 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1824 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1825 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1826 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1827 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1828 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1829 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1830 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1831 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1832 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1833 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1834 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1835 1836 IF (-zIIn1t2>=rsmall) THEN 1837 zHen1t2 = 1._wp 1838 ELSE 1839 zHen1t2 = 0._wp 1840 ENDIF 1841 1842 IF (-zIIn2t1>=rsmall) THEN 1843 zHen2t1 = 1._wp 1844 ELSE 1845 zHen2t1 = 0._wp 1846 ENDIF 1847 1848 s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 1849 1850 END FUNCTION s11ks 1851 1852 FUNCTION s12ks(px,py,pz) 1700 1701 allsk(4) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 1853 1702 !!------------------------------------------------------------------- 1854 1703 !! Function : s12ks 1855 1704 !!------------------------------------------------------------------- 1856 REAL(wp), INTENT(in ) :: px,py,pz1857 1858 REAL(wp) :: s12ks, zs12s0, zs21s0, zpih1859 1860 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i221861 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i221862 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i221863 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i221864 REAL(wp) :: zd11, zd12, zd221865 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t21866 REAL(wp) :: zHen1t2, zHen2t11867 !!-------------------------------------------------------------------1868 zpih = 0.5_wp*rpi1869 1870 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)1871 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)1872 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)1873 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)1874 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)1875 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)1876 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)1877 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)1878 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)1879 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)1880 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)1881 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)1882 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)1883 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)1884 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)1885 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)1886 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))1887 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))1888 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))1889 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd221890 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd221891 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd221892 1893 IF (-zIIn1t2>=rsmall) THEN1894 zHen1t2 = 1._wp1895 ELSE1896 zHen1t2 = 0._wp1897 ENDIF1898 1899 IF (-zIIn2t1>=rsmall) THEN1900 zHen2t1 = 1._wp1901 ELSE1902 zHen2t1 = 0._wp1903 ENDIF1904 1905 1705 zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 1906 1706 zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21) 1907 s12ks=0.5_wp*(zs12s0+zs21s0) 1908 1909 END FUNCTION s12ks 1910 1911 FUNCTION s22ks(px,py,pz) 1707 allsk(5)=0.5_wp*(zs12s0+zs21s0) 1912 1708 !!------------------------------------------------------------------- 1913 1709 !! Function : s22ks 1914 1710 !!------------------------------------------------------------------- 1915 REAL(wp), INTENT(in ) :: px,py,pz 1916 1917 REAL(wp) :: s22ks, zpih 1918 1919 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1920 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1921 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1922 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1923 REAL(wp) :: zd11, zd12, zd22 1924 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1925 REAL(wp) :: zHen1t2, zHen2t1 1926 !!------------------------------------------------------------------- 1927 zpih = 0.5_wp*rpi 1928 1929 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1930 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1931 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1932 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1933 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1934 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1935 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1936 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1937 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1938 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1939 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1940 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1941 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1942 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1943 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1944 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1945 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1946 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1947 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1948 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1949 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1950 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1951 1952 IF (-zIIn1t2>=rsmall) THEN 1953 zHen1t2 = 1._wp 1954 ELSE 1955 zHen1t2 = 0._wp 1956 ENDIF 1957 1958 IF (-zIIn2t1>=rsmall) THEN 1959 zHen2t1 = 1._wp 1960 ELSE 1961 zHen2t1 = 0._wp 1962 ENDIF 1963 1964 s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 1965 1966 END FUNCTION s22ks 1711 allsk(6) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 1712 END SUBROUTINE all_skr_sks 1967 1713 1968 1714 #else
Note: See TracChangeset
for help on using the changeset viewer.
