MODULE CSM_Report ! Code imported from FiniteMap & Map_Tools (and new code) ! to give SHELLS_for_CSM the capability to report ! model stresses in Community Stress Model* format, ! at lon/lat/elevation points specified by template ! file CSM_grid.txt. ! *A project of the Southern California Earthquake Center, ! led by Dr. Jeanne Hardebeck of USGS. ! This module by Peter Bird, UCLA, 2012.09.04+. USE Sphere ! Spherical-geometry, by Peter Bird; provided as Sphere.f90 IMPLICIT NONE REAL, DIMENSION(3,2,2,2):: dG REAL, DIMENSION(3,2,2) :: G REAL, DIMENSION(3) :: eps_dot INTEGER, DIMENSION(:,:), ALLOCATABLE :: neighbor ! 3 neighboring elements(?) for each element REAL, DIMENSION(:), ALLOCATABLE :: a_ ! area, per element REAL, DIMENSION(:,:), ALLOCATABLE :: center ! uvec, per element REAL, DIMENSION(:,:), ALLOCATABLE :: node_uvec ! uvec, per node CONTAINS SUBROUTINE CSM_Reporter(INPUT,ACREEP,ALPHAT, & & BCREEP,BIOT, & & CCREEP,CFRIC,CONDUC,cooling_curvature,& & DCREEP,density_anomaly,DQDTDA, & & ECREEP,ELEV, & & GMEAN,GRADIE, & & MXDOF,MXEL,MXFEL,MXNODE, & & NDOF,NFL,NODEF,NODES,NUMEL,NUMNOD, & & ONEKM, & & PULLED, & & RADIO,RADIUS,RHOAST,RHOBAR,RHOH2O, & & SIGHB,SIGZZN, & & TADIAB, TEMLIM,TLNODE,TRHMAX,TSURF, & & V,VISMAX, & & XNODE, & & YNODE, & & ZMNODE) IMPLICIT NONE INTEGER, INTENT(IN) :: INPUT, MXDOF, MXEL, MXFEL, MXNODE, & & NDOF, NFL, NUMEL, NUMNOD INTEGER, DIMENSION(3,MXEL), INTENT(IN) :: NODES INTEGER, DIMENSION(4,MXFEL),INTENT(IN) :: NODEF REAL, DIMENSION(2), INTENT(IN) :: ACREEP, ALPHAT, BCREEP, CCREEP, CONDUC,& & DCREEP, RADIO, RHOBAR, TEMLIM REAL, INTENT(IN) :: BIOT, CFRIC, ECREEP, GMEAN, GRADIE, ONEKM, & & RADIUS, RHOAST, RHOH2O, TADIAB, TRHMAX,& & TSURF, VISMAX REAL, DIMENSION(MXNODE), INTENT(IN) :: cooling_curvature, density_anomaly, & & DQDTDA, ELEV, SIGZZN, TLNODE, XNODE, YNODE, ZMNODE REAL, DIMENSION(2,7,MXEL), INTENT(IN) :: SIGHB LOGICAL, DIMENSION(7,MXEL), INTENT(IN) :: PULLED DOUBLE PRECISION, DIMENSION(2,MXNODE),INTENT(IN) :: V !N.B. These variables contain everything needed to compute ! strain-rates and stresses at any point in the continuum ! elements, which make up "100% of the volume" of any ! SHELLS domain. ! Most values are scalars (BIOT,CFRIC,GMEAN,ONEKM,RADIUS,...) ! two-vectors (crust/mantle: ACREEP, BCREEP, ...) ! or lists of values at nodes (ELEV,DQDTDA,V,...); ! the only quantities provided at integration points are ! PULLED (logical, T/F) and SIGHB (a horizontal 2-vector) ! for each integration point of each continuum element). ! The new code will be responsible for ! reading in the desired grid points, interpolating, ! and writing the new output file SHELLS_for_CSM.txt. INTEGER :: i, iLayer, ios, j, l_, M_of_closest_IP, n1, n2, n3, points_done LOGICAL :: cold_start, isotropic_stress, new_surface_point, success REAL :: air_pressure_Pa, ANGAT2,ANGAT3, ARGUME, asthenosphere_below_MSL_in_km, & & crust, curviness, & & delta_quadratic, depth_below_MSL_in_km, DTDZC, DTDZM, & & E1AT1,E2AT1,E1AT2,E2AT2,E1AT3,E2AT3,E1AT4,E2AT4, & & elevation, ELon, equat, err, e1h, e2h, & & FRAC, & & GAMMA, GEOTH1,GEOTH2,GEOTH3,GEOTH4,GEOTH5,GEOTH6,GEOTH7,GEOTH8, & & hard_surface_below_MSL_in_km, hard_surface_pore_pressure_Pa, & & local_density_anomaly, & & m_from_top, mantle_lithosphere, Moho_below_MSL_in_km, Moho_T, & & NLat, & & old_lat, old_lon, & & Q, & & rho_mean_for_layer, & & s1, S1EFF, s2, S2EFF, s3, SIGHBI, SIGMA1,SIGMA2, & & sigma_1h_MPa, sigma_1h_Pa, sigma_2h_MPa, sigma_2h_Pa, & & sigma_1h_creep_Pa, sigma_2h_creep_Pa, sigma_1h_friction_Pa, sigma_2h_friction_Pa, & & sigma_rr_MPa, sigma_rr_Pa, STFRIC, SZEFF, & & sigma_EE_MPa, sigma_EN_MPa, sigma_Er_MPa, sigma_NN_MPa, sigma_Nr_MPa, & & sigma_mean_horizontal_MPa, sigma_radius_horizontal_MPa, & & T_K, T_mean_for_layer_K, TASTHK, & & target_azimuth_degrees, target_pore_pressure_Pa, target_T, TERR0R, TEST, theta_, & & u1theta,u1phi, u2theta,u2phi, & & VIS, VISDCR, VISINF, VISSHB REAL, DIMENSION(3) :: uvec, uvec1, uvec2, uvec3 DOUBLE PRECISION :: SECINV ! to prevent squared strain-rates from underflowing, NOT for great precision! !Learn the topology of the SHELLS .FEG: ALLOCATE ( node_uvec(3, NUMNOD) )! N.B. This array is global in MODULE CSM_Report DO i = 1, NUMNOD ELon = degrees_per_radian * YNODE(i) NLat = 90.0 - degrees_per_radian * XNODE(i) NLat = MAX(MIN(NLat, 90.0), -90.0) CALL LonLat_2_Uvec(ELon, NLat, uvec) node_uvec(1:3, i) = uvec(1:3) END DO ALLOCATE ( a_(numel) ) ! N.B. This array is global in MODULE CSM_Report ALLOCATE ( center(3, numel) ) ! N.B. This array is global in MODULE CSM_Report ALLOCATE ( neighbor(3, numel) ) ! N.B. This array is global in MODULE CSM_Report CALL Learn_Spherical_Triangles (NUMEL, NODES, node_uvec, .TRUE., & & a_, center, neighbor) 10 WRITE (*, "(' ')") OPEN (UNIT = 11, FILE = "CSM_grid.txt", STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: CSM_grid.txt not found (in this folder). Please supply it.')") CALL Pause() GO TO 10 END IF OPEN (UNIT = 12, FILE = "SHELLS_for_CSM.txt") ! unconditional; overwrites any existing WRITE (12, "('# Stress model SHELLS_for_CSM, version 1.0, from SHELLS_for_CSM.f90.')") WRITE (12, "('# General characteristics of SHELLS models:')") WRITE (12, "('# -Domain includes crust, and often mantle lithosphere (if any, depending on heat-flow).')") WRITE (12, "('# -Rheology is frictional plasticity at low T and dislocation creep at high T (non-elastic).')") WRITE (12, "('# -Vertical integrals of stress through plate satisfy 2-D quasi-static equilibrium on a sphere.')") WRITE (12, "('# -No variation in horizontal velocity with depth, unless there is basal shear traction.')") WRITE (12, "('# -Models NOT constrained by internal GPS velocities or geologic slip rates or stress directions;')") WRITE (12, "('# so prediction errors in any of these quantities are tests of model quality.')") WRITE (12, "('# -Codes, citations, examples, and tutorial available at my web site: http://peterbird.name')") WRITE (12, "('# Specific characteristics of this model:')") WRITE (12, "('# -Elevations with 5-minute resolution from ETOPO5.')") WRITE (12, "('# -Heat-flow map (smooth, hand-contoured) of Blackwell & Steele [1992].')") WRITE (12, "('# -Lithosphere assumed isostatic, with all geotherms in steady-state.')") WRITE (12, "('# -No basal shear tractions applied to the lithosphere.')") WRITE (12, "('# -Inland boundary velocities from NeoKinema model UCERF3073 (constrained by GPS).')") WRITE (12, "('# -Fault traces and dips from UCERF3 Fault Model 3.1; slip-rates and rakes are free.')") WRITE (12, "('# -Quartz-like dislocation creep in crust, and olivine-like creep in mantle lithosphere.')") WRITE (12, "('# -Effective friction coefficient 0.15 on faults, but 0.85 between the faults.')") WRITE (12, "('# By Peter Bird, UCLA, 2012.09.05, for the SCEC Community Stress Model.')") WRITE (12, "('# Note that fault friction has not yet been formally optimized by scoring prediction errors.')") points_done = 0 WRITE (*, *) old_lat = 999. old_lon = 999. WRITE (*, "(' ')") processing: DO ! indefinite loop, defined by length of template grid file on unit 1 READ (11, *, IOSTAT = ios) ELon, NLat, depth_below_MSL_in_km IF (ios == -1) EXIT processing ! EOF condition new_surface_point = (ELon /= old_lon).OR.(NLat /= old_lat) IF (new_surface_point) THEN CALL LonLat_2_Uvec(ELon, NLat, uvec) cold_start = (ABS(ELon - old_lon) > 0.1).OR.(ABS(NLat - old_lat) > 0.1) CALL Which_Spherical_Triangle (uvec, cold_start, & & NUMEL, NODES, node_uvec, center, a_, neighbor, & & success, l_, s1, s2, s3) END IF IF (success) THEN !Decide which layer the test point is in, according to its depth (in km, below MSL): n1 = NODES(1, l_) n2 = NODES(2, l_) n3 = NODES(3, l_) elevation = s1 * ELEV(n1) + s2 * ELEV(n2) + s3 * ELEV(n3) ! of hard surface, in m, + above sea level hard_surface_below_MSL_in_km = -elevation/1000. crust = s1 * ZMNODE(n1) + s2 * ZMNODE(n2) + s3 * ZMNODE(n3) ! crustal thickness, in m Moho_below_MSL_in_km = hard_surface_below_MSL_in_km + crust/1000. mantle_lithosphere = s1 * TLNODE(n1) + s2 * TLNODE(n2) + s3 * TLNODE(n3) asthenosphere_below_MSL_in_km = Moho_below_MSL_in_km + mantle_lithosphere/1000. IF (depth_below_MSL_in_km < hard_surface_below_MSL_in_km) THEN ! in atmosphere or ocean IF (depth_below_MSL_in_km < 0.0) THEN iLayer = -1 ! atmosphere ELSE iLayer = 0 ! ocean END IF ELSE ! depth is inside the solid-Earth IF (depth_below_MSL_in_km < Moho_below_MSL_in_km) THEN iLayer = 1 ! crust ELSE IF (depth_below_MSL_in_km <= asthenosphere_below_MSL_in_km) THEN iLayer = 2 ! mantle lithosphere ELSE iLayer = 3 ! asthenosphere END IF END IF ! in fluid layers, or in solid layers? IF ((iLayer == 1).OR.(iLayer == 2).OR.(iLayer == 3)) THEN ! only report from hard-surface down !(not for points in atmosphere or ocean) !Compute the strain-rate tensor in this element: uvec1(1:3) = node_uvec(1:3, NODES(1, l_)) uvec2(1:3) = node_uvec(1:3, NODES(2, l_)) uvec3(1:3) = node_uvec(1:3, NODES(3, l_)) equat = SQRT(uvec(1)**2 + uvec(2)**2) theta_ = ATAN2(equat, uvec(3)) CALL Gjxy (l_, uvec1, & & uvec2, & & uvec3, & & uvec, G) CALL Del_Gjxy_del_thetaphi (l_, uvec1, & & uvec2, & & uvec3, & & uvec, dG) CALL E_rate(RADIUS, l_, NODES, G, dG, theta_, V, eps_dot) ! eps_dot(1) = epsilon_dot_sub_theta_theta = epsilon_dot_sub_S_S = epsilon_dot_sub_N_N ! eps_dot(2) = epsilon_dot_sub_theta_phi = epsilon_dot_sub_S_E = -epsilon_dot_sub_N_E ! eps_dot(3) = epsilon_dot_sub_phi_phi = epsilon_dot_sub_E_E !find the principal axes and values: CALL Principal_Axes_22 (eps_dot(1),eps_dot(2),eps_dot(3), & & e1h, e2h, u1theta,u1phi, u2theta,u2phi) !where principal axis #1 is the more-compressive of the two. target_azimuth_degrees = 180.0 - degrees_per_radian * ATAN2F(u1phi, u1theta) IF (target_azimuth_degrees < 0.0) target_azimuth_degrees = target_azimuth_degrees + 180.0 IF (target_azimuth_degrees < 0.0) target_azimuth_degrees = target_azimuth_degrees + 180.0 IF (target_azimuth_degrees > 180.0) target_azimuth_degrees = target_azimuth_degrees - 180.0 IF (target_azimuth_degrees > 180.0) target_azimuth_degrees = target_azimuth_degrees - 180.0 !compute (and adjust) geotherm, including target_T and Moho_T (both in Kelvin). ! (Note that this logic comes from subprogram FILLIN in SHELLS.) TASTHK = TADIAB + GRADIE * 100.E3 GEOTH1=TSURF GEOTH3= -0.5*RADIO(1)/CONDUC(1) GEOTH4=0. GEOTH7= -0.5*RADIO(2)/CONDUC(2) GEOTH8=0. Q = s1 * DQDTDA(n1) + s2 * DQDTDA(n2) + s3 * DQDTDA(n3) ! heat-flow, in W/m^2 GEOTH2 = Q / CONDUC(1) Moho_T = GEOTH1 + GEOTH2 * crust + GEOTH3 * crust**2 + GEOTH4 * crust**3 GEOTH5 = Moho_T DTDZC = GEOTH2 + 2. * GEOTH3 * crust + 3. * GEOTH4 * crust**2 DTDZM = DTDZC * CONDUC(1) / CONDUC(2) GEOTH6 = DTDZM !Now, correct geotherm to hit TASTHK: IF (mantle_lithosphere.GT.0.) THEN TEST=GEOTH5+ & & GEOTH6*mantle_lithosphere+ & & GEOTH7*mantle_lithosphere**2+ & & GEOTH8*mantle_lithosphere**3 ELSE TEST=GEOTH1+ & & GEOTH2*crust+ & & GEOTH3*crust**2+ & & GEOTH4*crust**3 END IF TERR0R=TEST-TASTHK delta_quadratic= -TERR0R/(crust+mantle_lithosphere)**2 curviness = -2. * delta_quadratic GEOTH3=GEOTH3 + delta_quadratic GEOTH7=GEOTH7 + delta_quadratic GEOTH5=GEOTH1+ & & GEOTH2*crust+ & & GEOTH3*crust**2+ & & GEOTH4*crust**3 DTDZC= GEOTH2+ & & 2.*GEOTH3*crust+ & & 3.*GEOTH4*crust**2 DTDZM=DTDZC*CONDUC(1)/CONDUC(2) GEOTH6=DTDZM Moho_T = GEOTH1 + GEOTH2 * crust + GEOTH3 * crust**2 + GEOTH4 * crust**3 IF (iLayer == 1) THEN m_from_top = 1000.0 * (depth_below_MSL_in_km - hard_surface_below_MSL_in_km) target_T = GEOTH1 + GEOTH2 * m_from_top + GEOTH3 * m_from_top**2 + GEOTH4 * m_from_top**3 ELSE IF (iLayer == 2) THEN m_from_top = 1000.0 * (depth_below_MSL_in_km - Moho_below_MSL_in_km) target_T = GEOTH5 + GEOTH6 * m_from_top + GEOTH7 * m_from_top**2 + GEOTH8 * m_from_top**3 ELSE IF (iLayer == 3) THEN target_T = TASTHK END IF !air pressure at the bottom of the atmosphere: IF (elevation <= 0) THEN ! at sea surface: air_pressure_Pa = 1.01325E5 ELSE ! positive elevation; land surface: air_pressure_Pa = 1.01325E5 * EXP(-0.693 * elevation/5300.) ! 5300 m gives half of sea-level P END IF !pore pressure and vertical compressive stress at hard surface (seafloor or land): IF (elevation <= 0) THEN ! at sea-floor: hard_surface_pore_pressure_Pa = air_pressure_Pa + RHOH2O * GMEAN * (-elevation) ELSE ! positive elevation; land surface: hard_surface_pore_pressure_Pa = air_pressure_Pa END IF sigma_rr_Pa = -hard_surface_pore_pressure_Pa !pore pressure and vertical compressive stress at target point within (or below) lithosphere: target_pore_pressure_Pa = hard_surface_pore_pressure_Pa + RHOH2O * GMEAN * & & (depth_below_MSL_in_km - hard_surface_below_MSL_in_km) * 1000.0 local_density_anomaly = s1 * density_anomaly(n1) + & ! of chemical origin, in kg/m^3 & s2 * density_anomaly(n2) + & & s3 * density_anomaly(n3) IF (iLayer == 1) THEN ! target in crust T_mean_for_layer_K = (TSURF + target_T) / 2.0 rho_mean_for_layer = (RHOBAR(1) + local_density_anomaly) * (1.0 - ALPHAT(1) * T_mean_for_layer_K) sigma_rr_Pa = sigma_rr_Pa - rho_mean_for_layer * GMEAN * & & (depth_below_MSL_in_km - hard_surface_below_MSL_in_km) * 1000.0 ELSE IF (iLayer == 2) THEN ! target in mantle lithosphere !first, add the contribution of the whole crust: T_mean_for_layer_K = (TSURF + Moho_T) / 2.0 rho_mean_for_layer = (RHOBAR(1) + local_density_anomaly) * (1.0 - ALPHAT(1) * T_mean_for_layer_K) sigma_rr_Pa = sigma_rr_Pa - rho_mean_for_layer * GMEAN * & & (Moho_below_MSL_in_km - hard_surface_below_MSL_in_km) * 1000.0 ! (this line = crust) !second, add the contribution of part of the mantle lithosphere T_mean_for_layer_K = (Moho_T + target_T) / 2.0 rho_mean_for_layer = (RHOBAR(2) + local_density_anomaly) * (1.0 - ALPHAT(2) * T_mean_for_layer_K) sigma_rr_Pa = sigma_rr_Pa - rho_mean_for_layer * GMEAN * & & (depth_below_MSL_in_km - Moho_below_MSL_in_km) * 1000.0 ELSE IF (iLayer == 3) THEN ! target in asthenosphere !first, add the contribution of the whole crust: T_mean_for_layer_K = (TSURF + Moho_T) / 2.0 rho_mean_for_layer = (RHOBAR(1) + local_density_anomaly) * (1.0 - ALPHAT(1) * T_mean_for_layer_K) sigma_rr_Pa = sigma_rr_Pa - rho_mean_for_layer * GMEAN * & & (Moho_below_MSL_in_km - hard_surface_below_MSL_in_km) * 1000.0 ! (this line = crust) !second, add the contribution of part of the whole mantle lithosphere T_mean_for_layer_K = (Moho_T + TASTHK) / 2.0 rho_mean_for_layer = (RHOBAR(2) + local_density_anomaly) * (1.0 - ALPHAT(2) * T_mean_for_layer_K) sigma_rr_Pa = sigma_rr_Pa - rho_mean_for_layer * GMEAN * & & (asthenosphere_below_MSL_in_km - Moho_below_MSL_in_km) * 1000.0 !third, add contribution of asthenosphere above target: sigma_rr_Pa = sigma_rr_Pa - RHOAST * GMEAN * & & (depth_below_MSL_in_km - asthenosphere_below_MSL_in_km) * 1000.0 END IF isotropic_stress = (iLayer == 3) ! This needs to be improved to use PULLED and SIGHB, in future! IF (isotropic_stress) THEN ! in asthenosphere (IF no basal traction). sigma_rr_MPa = sigma_rr_Pa / 1.0E6 WRITE (12, "(F9.4, F9.4, F7.2, ' 1 NaN',4F8.1,' 0 0',F8.1,' 0',F8.1,' 0 0 1')") & & ELon, NLat, depth_below_MSL_in_km, sigma_rr_MPa, sigma_rr_MPa, sigma_rr_MPa, & & sigma_rr_MPa, sigma_rr_MPa, sigma_rr_MPa WRITE (12, "(F9.4, F9.4, F7.2, ' 2 NaN NaN NaN NaN 0 NaN 0 NaN NaN NaN NaN NaN NaN')") & & ELon, NLat, depth_below_MSL_in_km ELSE ! normal case (within lithosphere, or with basal traction in asthenosphere) !Compute & report deviatoric stress at target, based on lesser of frictional or dislocation-creep strength. !This code is adapted from that of complex subprogram DIAMND of SHELLS (revised 1998). !Further characterize the strain-rate tensor: err = -(e1h + e2h) ! because non-elastic deformation is incompressible SECINV= -((1.D0*e1h)*e2h + (1.D0*e1h)*err + (1.D0*e2h)*err) !(One possible form for the second invariant of the matrix. ! Note that the double-precision is just to prevent underflows ! from squaring small strain rates, not for precision.) VISINF=0.5*ACREEP(MIN(iLayer,2))*(2.*SQRT(SECINV))**(ECREEP-1.) !VISINF is the viscosity for dislocation creep, lacking only !the exponential term; therefore, as a mathematical abstraction, !we can say that it is the viscosity at infinite temperature. !Characterize the continuum friction: STFRIC=SIN(ATAN(CFRIC)) GAMMA=(1+STFRIC)/(1-STFRIC) !Note: For thrusting, effective-sigma1h is effective-sigma1z ! times GAMMA. For normal faulting, effective-sigma2h ! is effective-sigmaz/GAMMA. For very small CFRIC, GAMMA ! would be approximately equal to (1.+2.*CFRIC). !Compute the effective vertical stress at the target depth, ! ASSUMING that it is within the frictional layer: SZEFF = sigma_rr_Pa + BIOT * target_pore_pressure_Pa ! (less negative than sigma_rr_Pa) !Compute effective horizontal principal stresses, ! assuming that the target is in the frictional layer, according ! to the methods in Bird (1989), pages 3973-3977 ! (except, correcting the typos in the caption for Figure 4): !Define the corner points of the diamond in the ! ordered principal strain-rate plane: E1AT1=((1./GAMMA)-1.)*SZEFF/(6.*VISMAX) E2AT1=E1AT1 E1AT2=(1.-(1./GAMMA))*SZEFF/(6.*VISMAX) E2AT2=((2./GAMMA)-2.)*SZEFF/(6.*VISMAX) E1AT3=(2.*GAMMA-2.)*SZEFF/(6.*VISMAX) E2AT3=(1.-GAMMA)*SZEFF/(6.*VISMAX) E1AT4=(GAMMA-1.)*SZEFF/(6.*VISMAX) E2AT4=E1AT4 ANGAT2=ATAN2F((e2h-E2AT2),(e1h-E1AT2)) ANGAT3=ATAN2F((e2h-E2AT3),(e1h-E1AT3)) !Select proper segment of diagram and assign effective principal stresses. IF (e1h.GT.E1AT1) THEN !Region N/N: two conjugate sets of normal faults S1EFF=SZEFF/GAMMA S2EFF=S1EFF ELSE IF ((e1h.GE.E1AT2).AND. & & (ANGAT2.GT.ATAN2F((E2AT1-E2AT2),(E1AT1-E1AT2)))) THEN !Region N: single conjugate set of normal faults S2EFF=SZEFF/GAMMA FRAC=(e1h-E1AT1)/(E1AT2-E1AT1) !fraction increases in -E1 direction, from point 1 -> 2 S1EFF=SZEFF*((1/GAMMA)+FRAC*(1.-(1./GAMMA))) ELSE IF ((ANGAT2.LE.1.9635).AND.(ANGAT2.GE.1.5707)) THEN !Region N/S: transtension, dominantly normal. S1EFF=SZEFF S2EFF=SZEFF/GAMMA ELSE IF ((ANGAT2.LE.2.3562).AND.(ANGAT2.GE.1.9635)) THEN !Region S/N: transtension, dominantly strike-slip. S1EFF=SZEFF S2EFF=SZEFF/GAMMA ELSE IF ((ANGAT3.LE.2.3562).AND. & & (ANGAT3.GE.ATAN2F((E2AT2-E2AT3),(E1AT2-E1AT3)))) THEN !Region S: single set of conjugate strike-slip faults FRAC=((e1h+e2h)-(E1AT2+E2AT2))/((E1AT3+E2AT3)-(E1AT2+E2AT2)) !FRAC increases across band from the S/N (point 2) side !toward the S/T (point 3) side; contours of FRAC are !parallel to the band sides, not normal to the diamond. S1EFF=SZEFF*(1.+FRAC*(GAMMA-1.)) S2EFF=SZEFF*((1./GAMMA)+FRAC*(1.-(1./GAMMA))) !Notes: The equation of this line is S2EFF=S1EFF/GAMMA. !I used algebra to check (1998.04.21) that the !pure strike-slip stress (S1EFF,S2EFF)= !SZZEFF*(1.+STFRIC,1.-STFRIC) correctly falls on !this line, at the correct point (E1= -E2). ELSE IF ((ANGAT3.LE.2.7489).AND.(ANGAT3.GE.2.3562)) THEN !Region S/T: transpression; strike-slip dominant. S1EFF=SZEFF*GAMMA S2EFF=SZEFF ELSE IF ((e2h.GE.E2AT3).AND.(ANGAT3.GE.2.7489)) THEN !Region T/S: transpression; thrusting dominant. S1EFF=SZEFF*GAMMA S2EFF=SZEFF ELSE IF ((e2h.GE.E2AT4).AND. & & (ANGAT3.LE.ATAN2F((E2AT4-E2AT3),(E1AT4-E1AT3)))) THEN !Region T: single conjugate thrust fault set. S1EFF=SZEFF*GAMMA FRAC=(e2h-E2AT3)/(E2AT4-E2AT3) !FRAC increases in the -E2 direction across the band. S2EFF=SZEFF*(1.+FRAC*(GAMMA-1.)) ELSE IF (e2h.LE.E2AT4) THEN !Region T/T: Two set of conjugate thrust faults. S1EFF=SZEFF*GAMMA S2EFF=S1EFF ELSE !Region V: linear viscosity !Note that equations are now for SIGMA1,2 and no !longer for S1EFF and S2EFF. However, we can easily compute both: SIGMA1=sigma_rr_Pa+VISMAX*(4.*e1h+2.*e2h) SIGMA2=sigma_rr_Pa+VISMAX*(2.*e1h+4.*e2h) S1EFF=SIGMA1+BIOT*target_pore_pressure_Pa S2EFF=SIGMA2+BIOT*target_pore_pressure_Pa END IF ! different regions of the polygon of Bird [1998, JGR]. sigma_1h_friction_Pa = sigma_rr_Pa + (S1EFF - SZEFF) sigma_2h_friction_Pa = sigma_rr_Pa + (S2EFF - SZEFF) !because differences between principal stresses equal differences between !effective principal stresses (both being corrected by the same BIOT*target_pore_pressure_Pa). !Evaluate alternative estimates of principal stresses by !ASSUMING that the target depth is in the dislocation-creep regime: !Precompute the maximum viscosity limit imposed by the !requirement that creep shear stress never exceeds !DCREEP(MIN(iLayer,2)) on any plane: VISDCR=DCREEP(MIN(iLayer,2))/(MAX(e1h,e2h,err)-MIN(e1h,e2h,err)) !Effective dislocation-creep viscosity is now: T_K=MIN(target_T,TEMLIM(MIN(iLayer,2))) ARGUME=(BCREEP(MIN(iLayer,2))+CCREEP(MIN(iLayer,2))*(depth_below_MSL_in_km*1000.))/T_K ARGUME=MAX(MIN(ARGUME,87.),-87.) VIS=VISINF*EXP(ARGUME) VIS=MIN(VIS,VISDCR) sigma_1h_creep_Pa = sigma_rr_Pa + 2.0 * VIS * (e1h - err) sigma_2h_creep_Pa = sigma_rr_Pa + 2.0 * VIS * (e2h - err) !Select estimate closest to vertical stress (for each axis); minimize deviatoric stress. IF (ABS(sigma_1h_friction_Pa - sigma_rr_Pa) < ABS(sigma_1h_creep_Pa - sigma_rr_Pa)) THEN sigma_1h_Pa = sigma_1h_friction_Pa ELSE sigma_1h_Pa = sigma_1h_creep_Pa END IF IF (ABS(sigma_2h_friction_Pa - sigma_rr_Pa) < ABS(sigma_2h_creep_Pa - sigma_rr_Pa)) THEN sigma_2h_Pa = sigma_2h_friction_Pa ELSE sigma_2h_Pa = sigma_2h_creep_Pa END IF !GPBhere !FINALLY, report the full stress tensor! sigma_1h_MPa = sigma_1h_Pa / 1.E6 sigma_2h_MPa = sigma_2h_Pa / 1.E6 sigma_rr_MPa = sigma_rr_Pa / 1.E6 sigma_mean_horizontal_MPa = (sigma_1h_MPa + sigma_2h_MPa)/2.0 ! always negative sigma_radius_horizontal_MPa = ABS(sigma_1h_MPa - sigma_2h_MPa)/2.0 ! always positive sigma_EE_MPa = sigma_mean_horizontal_MPa + & & sigma_radius_horizontal_MPa * COS(2.0 * target_azimuth_degrees * radians_per_degree) sigma_EN_MPa = -sigma_radius_horizontal_MPa * SIN(2.0 * target_azimuth_degrees * radians_per_degree) sigma_NN_MPa = sigma_mean_horizontal_MPa - & & sigma_radius_horizontal_MPa * COS(2.0 * target_azimuth_degrees * radians_per_degree) sigma_Er_MPa = 0.0 ! this will need to change if there is basal shear traction! sigma_Nr_MPa = 0.0 ! this will need to change if there is basal shear traction! WRITE (12, "(F9.4, F9.4, F7.2, ' 1',F6.1,9F8.1,' 1 1 1')") & & ELon, NLat, depth_below_MSL_in_km, target_azimuth_degrees, & & sigma_1h_MPa, sigma_2h_MPa, sigma_rr_MPa, & & sigma_EE_MPa, sigma_EN_MPa, sigma_Er_MPa, sigma_NN_MPa, sigma_Nr_MPa, sigma_rr_MPa WRITE (12, "(F9.4, F9.4, F7.2, ' 2 NaN NaN NaN NaN 0 NaN 0 NaN NaN NaN NaN NaN NaN')") & & ELon, NLat, depth_below_MSL_in_km END IF END IF ! iLayer = 1 or 2 or 3 (crust, mantle lithosphere, or asthenosphere); no reporting for air or sea! END IF ! "success" in finding the finite element that contains the target (ELon, NLat); no reporting otherwise! old_lat = NLat old_lon = Elon points_done = points_done + 1 IF (MOD(points_done, 10000) == 0) THEN WRITE (*, "('+points processed = ',I8)") points_done END IF END DO processing ! indefinite loop, defined by length of grid file on unit 2 WRITE (*, "('+points processed = DONE')") WRITE (*, *) CLOSE (UNIT = 11) ! CSM_grid.txt (input) CLOSE (UNIT = 12) ! SHELLS_for_CSM.txt (output) END SUBROUTINE CSM_Reporter !============================================================================================ ! SUPPORTING SUBPROGRAMS IN ALPHABETICAL ORDER (most from FiniteMap and/or Map_Tools): REAL FUNCTION ATAN2F(y, x) !Corrects for problem of abend due to ATAN2(0.,0.): IMPLICIT NONE REAL, INTENT(IN) :: y, x IF ((y /= 0.).OR.(x /= 0.)) THEN ATAN2F=ATAN2(y,x) ELSE ATAN2F=0. END IF END FUNCTION ATAN2F SUBROUTINE Del_Gjxy_del_thetaphi (l_, uvec1, uvec2, uvec3, r_, dG) INTEGER, INTENT(IN) :: l_ ! element number or code REAL, DIMENSION(3), INTENT(IN) :: uvec1, uvec2, uvec3, & & r_ ! position vector REAL, DIMENSION (3,2,2,2), INTENT(OUT) :: dG ! Computes array of 2 derivitives of each of the 2 components of ! each of the 6 nodal functions for element l_, ! whose corners are uvec1, uvec2, and uvec3 (counterclockwise!), ! evaluated at position r_ (Cartesian unit vector). ! Results are in 1./radian (dimensionless), NOT 1./degree ! ! Also note that derivitives with respect to phi are algebraic, ! not spatial; that is, they have not been corrected for the ! changing length of a unit of phi with latitude, so these ! derivitives tend to become small near the poles, where ! a unit of phi does not cover very much distance. ! It is user's responsibility that element l_ contains r_. SAVE ! allows fast re-entry when l_ is unchanged. INTEGER :: l_last = 0 ! remembers l_ from previous invocation INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2)? INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function? INTEGER :: m ! 1:2 = theta (S-ward) or phi (E-ward) derivitive? REAL, DIMENSION(3,2) :: del_r_ ! theta- and phi-derivitives of r_ (in 3-D) REAL, DIMENSION(3,2) :: local ! local Theta, Phi unit vectors at r_ (xyz, SE) REAL, DIMENSION(3,2,2) :: del_local ! theta-, phi- derivitives of local REAL, DIMENSION(3,3) :: corner ! positions vector of corner nodes (xyz, 123) REAL, DIMENSION(3,3,2) :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) REAL, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, tv3, vfa, vfb ! temporary vector factors REAL :: cos_phi, cos_theta, phi, sin_phi, sin_theta INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ corner(1:3, 1) = uvec1(1:3) corner(1:3, 2) = uvec2(1:3) corner(1:3, 3) = uvec3(1:3) DO j = 1, 3 tvi = corner(1:3, j) CALL Local_Theta(tvi, tvo) post(1:3, j, 1) = tvo CALL Local_Phi (tvi, tvo) post(1:3, j, 2) = tvo END DO END IF ! begin calculations which depend on r_ CALL Local_Theta(r_, tv) local(1:3,1) = tv CALL Local_Phi(r_, tv) local(1:3,2) = tv ! Note: these functions will catch polar points; don't test again phi = ATAN2(r_(2), r_(1)) cos_phi = COS(phi) sin_phi = SIN(phi) cos_theta = r_(3) sin_theta = SQRT(r_(1)**2 + r_(2)**2) del_r_(1:3,1) = local(1:3,1) ! d.r_/d.theta = Theta del_r_(1:3,2) = local(1:3,2) * sin_theta ! d.r_/d.phi = Phi * SIN(theta) del_local(1:3,1,1) = - r_(1:3) ! d.Theta/d.theta = - r_ del_local(1:3,1,2) = local(1:3,2) * cos_theta ! d.Theta/d.phi = Phi * COS(theta) del_local(1:3,2,1) = (/ 0., 0., 0. /) ! d.Phi/d.theta = 0 del_local(1:3,2,2) = (/ -cos_phi, -sin_phi, 0. /) ! d.Phi/d.phi = (-COS(phi),-SIN(phi,0) DO j = 1, 3 ! 3 corner nodes of element i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1 = corner(1:3, i1) tv2 = corner(1:3, i2) tv3 = corner(1:3, i3) CALL Cross(tv2, tv3, vfa) vfb = vfa / Dot (tv1, vfa) DO x = 1, 2 ! unit velocity at node is S or E DO y = 1, 2 ! S- or E- component of nodal function tv1 = post(1:3, j, x) tvi = local(1:3, y) DO m = 1, 2 ! theta- or phi-derivitive tv = del_r_(1:3, m) tvo = del_local(1:3, y, m) dG(j, x, y, m) = & & (Dot(tv,vfb)*Dot(tv1,tvi)) + & & (Dot(r_,vfb)*Dot(tv1,tvo)) END DO END DO END DO END DO END SUBROUTINE Del_Gjxy_del_thetaphi SUBROUTINE Dumb_s123 (element, vector, node, xyz_nod, center, a_, & & s1, s2, s3) ! Finds s1, s2, s3 coordinates of position vector "in" element ! (whether the point is actually in the element or NOT). IMPLICIT NONE INTEGER, INTENT(IN) :: element ! element # REAL, DIMENSION(3), INTENT(IN) :: vector ! uvec to point INTEGER, DIMENSION(:,:), INTENT(IN) :: node ! element definitions REAL, DIMENSION(:,:), INTENT(IN) :: xyz_nod ! uvecs of nodes REAL, DIMENSION(:,:), INTENT(IN) :: center ! uvec of each element (uvec) REAL, DIMENSION(:), INTENT(IN) :: a_ ! element areas (plane; R == 1.0) REAL, INTENT(OUT) :: s1, s2, s3 ! results !- - - - - - - - - - - - - - - - - - - - - - - - INTEGER :: i1, i2, i3 REAL, DIMENSION(3) :: tv, tvi, tvo, tv1, tv2, v1 REAL :: d1, dc, t IF (element == 0) THEN WRITE (*,"(' ERROR: element = 0 in Dumb_s123')") CALL Traceback END IF i1 = node(1, element) i2 = node(2, element) i3 = node(3, element) !shorten(?) vector to just touch plane element -> v1 tv1 = center(1:3, element) dc = DOT_PRODUCT(vector, tv1) IF (dc <= 0.) THEN WRITE (*,"(' ERROR: Internal vector >= 90 deg. from element in Dumb_s123')") CALL Traceback END IF tv2 = xyz_nod(1:3, i1) d1 = DOT_PRODUCT(tv2, tv1) t = d1 / dc v1 = t * vector tvi = xyz_nod(1:3,i3) - xyz_nod(1:3,i2) tvo = v1(1:3) - xyz_nod(1:3,i3) CALL Cross(tvi, tvo, tv) s1 = 0.5 * DOT_PRODUCT(tv1, tv) / a_(element) tvi = xyz_nod(1:3,i1) - xyz_nod(1:3,i3) tvo = v1(1:3) - xyz_nod(1:3,i1) CALL Cross(tvi, tvo, tv) s2 = 0.5 * DOT_PRODUCT(tv1, tv) / a_(element) s3 = 1.00 - s1 - s2 END SUBROUTINE Dumb_s123 SUBROUTINE E_rate(R, l_, nodes, G, dG, theta_, V, eps_dot) ! evaluate strain-rate in spherical continuum element REAL, INTENT(IN) :: R ! radius of planet, in m INTEGER, INTENT(IN) :: l_ ! element number INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes REAL, DIMENSION(3,2,2) :: G ! nodal functions @ selected point REAL, DIMENSION(3,2,2,2):: dG ! derivitives of nodal functions REAL, INTENT(IN) :: theta_ ! colatitude, radians DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: V REAL, DIMENSION(3), INTENT(OUT) :: eps_dot INTEGER :: iv, iw, j REAL :: cott, csct, prefix eps_dot = 0. ! (1..3) cott = 1. / TAN(theta_) csct = 1. / SIN(theta_) prefix = 1./R DO j = 1, 3 ! epsilon_dot_sub_theta_theta eps_dot(1) = eps_dot(1) + & & V(1, nodes(j, l_)) * prefix * dG(j,1,1,1) + & & V(2, nodes(j, l_)) * prefix * dG(j,2,1,1) ! epsilon_dot_sub_theta_phi eps_dot(2) = eps_dot(2) + & & V(1, nodes(j, l_)) * prefix * 0.5 * (csct * dG(j,1,1,2) + dG(j,1,2,1) - cott * G(j,1,2)) + & & V(2, nodes(j, l_)) * prefix * 0.5 * (csct * dG(j,2,1,2) + dG(j,2,2,1) - cott * G(j,2,2)) ! epsilon_dot_sub_phi_phi eps_dot(3) = eps_dot(3) + & & V(1, nodes(j, l_)) * prefix * (csct * dG(j,1,2,2) + cott * G(j,1,1)) + & & V(2, nodes(j, l_)) * prefix * (csct * dG(j,2,2,2) + cott * G(j,2,1)) END DO ! 3 local nodes END SUBROUTINE E_rate ! eps_dot(1) = epsilon_dot_sub_theta_theta = epsilon_dot_sub_S_S = epsilon_dot_sub_N_N ! eps_dot(2) = epsilon_dot_sub_theta_phi = epsilon_dot_sub_S_E = -epsilon_dot_sub_N_E ! eps_dot(3) = epsilon_dot_sub_phi_phi = epsilon_dot_sub_E_E SUBROUTINE Gjxy (l_, uvec1, uvec2, uvec3, r_, G) INTEGER, INTENT(IN) :: l_ ! element number or code REAL, DIMENSION(3), INTENT(IN) :: uvec1, uvec2, uvec3, & ! corners & r_ ! position vector REAL, DIMENSION (3,2,2), INTENT(OUT) :: G ! Computes matrix of 6 vector nodal functions for element l_, ! whose corners are at uvec1, uvec2, and uvec3 (counterclockwise!), ! evaluated at position r_ (Cartesian unit vector). ! It is user's responsibility that element l_ contains r_. SAVE ! allows fast re-entry when l_ is unchanged. INTEGER :: l_last = -999 ! remembers l_ from previous invocation INTEGER :: j ! 1:3 = local node numbering in element l_ INTEGER :: x ! 1:2 = node j has unit velocity to South(1) or East(2) INTEGER :: y ! 1:2 = South(1) or East(2) component of vector nodal function REAL, DIMENSION(3,2) :: local ! local unit vectors at r_ (xyz, SE) REAL, DIMENSION(3,3) :: corner ! positions vector of corner nodes (xyz, 123) REAL, DIMENSION(3,3,2) :: post ! unit coordinate vectors at corner nodes: ! (xyz, 123, SE) REAL, DIMENSION(3) :: tvi, tvo, tv1, tv2, tv3, vf ! temporary vector factor REAL :: f_sup_j ! as in Kong and Bird (1995) [ j == k ] INTEGER :: i1, i2, i3 ! 1, 2, or 3 in cyclic rotation (depends on j) IF (l_ /= l_last) THEN ! new finite element l_last = l_ corner(1:3, 1) = uvec1(1:3) corner(1:3, 2) = uvec2(1:3) corner(1:3, 3) = uvec3(1:3) DO j = 1, 3 tvi = corner(1:3, j) CALL Local_Theta(tvi, tvo) post(1:3, j, 1) = tvo CALL Local_Phi (tvi, tvo) post(1:3, j, 2) = tvo END DO END IF ! begin computations which depend on r_ CALL Local_Theta(r_, tvo) local(1:3,1) = tvo CALL Local_Phi(r_, tvo) local(1:3,2) = tvo DO j = 1, 3 i1 = j i2 = 1 + MOD(j, 3) i3 = 1 + MOD(i2,3) tv1 = corner(1:3, i1) tv2 = corner(1:3, i2) tv3 = corner(1:3, i3) CALL Cross(tv2, tv3, vf) f_sup_j = Dot(r_, vf) / Dot (tv1, vf) DO x = 1, 2 tv1 = post(1:3, j, x) DO y = 1, 2 tv2 = local(1:3, y) G(j, x, y) = f_sup_j * Dot(tv1, tv2) END DO END DO END DO END SUBROUTINE Gjxy SUBROUTINE Learn_Spherical_Triangles (numel, nodes, node_uvec, chatty, & & a_, center, neighbor) !Creates arrays needed by lookup subr. Which_Spherical_Triangle: ! a_ = area of plane triangle below element, when radius == 1.0 ! center = uvec pointing to center of element ! neighbor = neighboring spherical triangular elements on each ! side (or zero if none); note that algorithm ! depends on node-location match, not on node-number ! match, and therefore ignores intevening faults. !These arrays are only meaningful for finite element grids used ! with SHELLS and/or RESTORE. IMPLICIT NONE INTEGER, INTENT(IN) :: numel ! number of spherical triangle elements INTEGER, DIMENSION(:,:), INTENT(IN) :: nodes ! element definitions REAL, DIMENSION(:,:), INTENT(IN) :: node_uvec ! uvecs of nodes LOGICAL, INTENT(IN) :: chatty REAL, DIMENSION(:), INTENT(OUT) :: a_ REAL, DIMENSION(:,:), INTENT(OUT) :: center INTEGER, DIMENSION(:,:), INTENT(OUT) :: neighbor INTEGER :: furthest, i, ia, ib, i1, i2, i3, j, j1, j2, j3, k, l_, m, step_aside REAL, DIMENSION(3) :: a, b, c, t, u IF (chatty) WRITE (*,"(' Learning the spherical triangles...')") furthest = (numel + 1) / 2 neighbor = 0 ! whole array, initialized to "no neighbor on this side" homes: DO l_ = 1, numel !first, a_ i1 = nodes(1,l_) i2 = nodes(2,l_) i3 = nodes(3,l_) a = node_uvec(1:3,i2) - node_uvec(1:3,i1) b = node_uvec(1:3,i3) - node_uvec(1:3,i2) CALL Cross (a, b, c) a_(l_) = 0.5 * Magnitude(c) !second, compute center t(1:3) = (node_uvec(1:3,i1)+node_uvec(1:3,i2)+node_uvec(1:3,i3))/3.0 CALL Make_Uvec(t, u) center(1:3, l_) = u(1:3) !third, find neighbor(?) for each side of element sides: DO j = 1, 3 ! 3 sides k = 1 + MOD (j, 3) ia = nodes(k, l_) ! 1st node along side ib = nodes(1 + MOD (k, 3), l_) ! 2nd node along side strangers: DO step_aside = 1, furthest m = l_ + step_aside ! I also try -step_aside, below m = 1 + MOD(m-1, numel) ! wraps around j1 = nodes(1, m) j2 = nodes(2, m) j3 = nodes(3, m) IF (Same(j1, ib) .AND. Same(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j2, ib) .AND. Same(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j3, ib) .AND. Same(j1, ia)) THEN neighbor(j, l_) = m EXIT strangers END IF m = l_ - step_aside ! I also try +step_aside, above m = 1 + MOD(m-1+numel, numel) ! wraps around j1 = nodes(1, m) j2 = nodes(2, m) j3 = nodes(3, m) IF (Same(j1, ib) .AND. Same(j2, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j2, ib) .AND. Same(j3, ia)) THEN neighbor(j, l_) = m EXIT strangers ELSE IF (Same(j3, ib) .AND. Same(j1, ia)) THEN neighbor(j, l_) = m EXIT strangers END IF END DO strangers END DO sides IF (chatty) WRITE (*,"('+Learning the spherical triangles...',I6)") l_ END DO homes IF (chatty) WRITE (*,"('+Learning the spherical triangles...DONE ')") CONTAINS LOGICAL FUNCTION Same(i,j) ! Are node_uvec #i and #j the same vector? INTEGER, INTENT(IN) :: i, j !the logic is: !Same = (node_uvec(1,i) == node_uvec(1,j)).AND. & ! & (node_uvec(2,i) == node_uvec(2,j)).AND. & ! & (node_uvec(3,i) == node_uvec(3,j)) !But, it is written this way for speed: IF (node_uvec(1,i) == node_uvec(1,j)) THEN IF (node_uvec(2,i) == node_uvec(2,j)) THEN IF (node_uvec(3,i) == node_uvec(3,j)) THEN Same = .TRUE. ELSE Same = .FALSE. END IF ELSE Same = .FALSE. END IF ELSE Same = .FALSE. END IF END FUNCTION Same END SUBROUTINE Learn_Spherical_Triangles SUBROUTINE Principal_Axes_22 (e11, e12, e22, & & e1, e2, u1x,u1y, u2x,u2y) ! Find principal values (e1,e2) of the symmetric 2x2 tensor ! e11 e12 ! e12 e22 ! and also the associated unit eigenvectors: ! #1 = (u1x, u1y); #2 = (u2x, u2y). ! The convention is that e1 <= e2. IMPLICIT NONE REAL, INTENT(IN) :: e11, e12, e22 REAL, INTENT(OUT) :: e1, e2, u1x, u1y, u2x, u2y REAL :: c, f1, f11, f12, f2, f22, r, scale, smallest, test, theta ! Smallest number that can be squared without underflowing: smallest = 1.1 * SQRT(TINY(f1)) ! First, check for trivial isotropic case: IF ((e11 == e22).AND.(e12 == 0.)) THEN ! In this case, directions are arbitrary: e1 = e11 u1x = 1. u1y = 0. e2 = e22 u2x = 0. u2y = 1. RETURN END IF ! Else, re-scale matrix to prevent under/overflows: scale = MAX(ABS(e11), ABS(e12), ABS(e22)) f11 = e11 / scale IF (ABS(f11) < smallest) f11 = 0. f12 = e12 / scale IF (ABS(f12) < smallest) f12 = 0. f22 = e22 / scale IF (ABS(f22) < smallest) f22 = 0. ! Find eigenvalues and eigenvectors of scaled matrix: r = SQRT(((f11 - f22)*0.5)**2 + f12**2) c = (f11 + f22)*0.5 f1 = c - r f2 = c + r test = 0.01 * MAX (ABS(f1), ABS(f2)) IF ((ABS(f12) > test).OR.(ABS(f11 - f1) > test)) THEN theta = Atan2F((f11 - f1), -f12) ELSE theta = Atan2F(f12, (f1 - f22)) END IF u1x = COS(theta) u1y = SIN(theta) u2x = u1y u2y = -u1x ! Undo the scaling e1 = scale * f1 e2 = scale * f2 END SUBROUTINE Principal_Axes_22 SUBROUTINE Pull_in(s) ! If necessary, adjusts internal coordinates s(1..3) so ! that none is negative. IMPLICIT NONE REAL, DIMENSION(3), INTENT(INOUT) :: s INTEGER, DIMENSION(1) :: array ! stupid, to satisfy MINLOC REAL factor, lowest, highest, medium INTEGER :: side, sidea, sideb lowest = MINVAL(s) IF (lowest < 0.) THEN highest = MAXVAL(s) medium = 1.00 - lowest - highest IF (medium > 0.) THEN ! correct to nearest edge array = MINLOC(s) side = array(1) s(side) = 0. sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) factor = 1.00 / (1.00 - lowest) s(sidea) = factor * s(sidea) ! s(sideb) = factor * s(sideb) would be logical s(sideb) = 1.00 - s(sidea) ! is safer ELSE ! correct to nearest vertex array = MAXLOC(s) side = array(1) s(side) = 1.00 sidea = 1 + MOD(side, 3) sideb = 1 + MOD(sidea, 3) s(sidea) = 0. s(sideb) = 0. END IF END IF END SUBROUTINE Pull_in SUBROUTINE Which_Spherical_Triangle (b_, cold_start, & & num_ele, node, xyz_node, center, a_, neighbor, & & success, iele, s1, s2, s3) !Locates a point (b_, a uvec) in element iele with internal !coordinates (s1, s2, s3) in a SHELLS or RESTORE .feg. !and reports success. !If (cold_start), makes no use of input iele, s1, s2, s3 !If not, uses these values to initialize the search. ! !Note that Learn_Spherical_Triangles can be used to initialize !necessary arrays a_, center, and neighbor. ! !Beware of variable name changes (numel = num_ele, nodes = node, etc.). IMPLICIT NONE REAL, DIMENSION(3), INTENT(IN) :: b_ ! uvec of unknown point LOGICAL, INTENT(IN) :: cold_start ! mode switch INTEGER, INTENT(IN) :: num_ele ! count of elements INTEGER, DIMENSION(:,:),INTENT(IN) :: node ! element definitions REAL, DIMENSION(:,:),INTENT(IN) :: xyz_node ! uvecs of nodes REAL, DIMENSION(:,:),INTENT(IN) :: center ! center uvecs of elements REAL, DIMENSION(:), INTENT(IN) :: a_ ! (plane) areas of elements INTEGER, DIMENSION(:,:),INTENT(IN) :: neighbor ! neighbors of each element LOGICAL, INTENT(OUT) :: success ! OUTPUT INTEGER, INTENT(INOUT) :: iele ! OUTPUT REAL, INTENT(INOUT) :: s1, s2, s3 ! OUTPUT !- - - - - - - - - - - - - - - - - - - - - - - - - - - - INTEGER :: back1, back2, back3, i, iet, l_ REAL :: r2, r2min, s1t, s2t, s3t REAL, DIMENSION(3) :: s_temp, tv ! establish defaults (not found) in case of quick exit success = .FALSE. IF (cold_start) THEN iele = 0 ! default s1 = 0.0; s2 = 0.0; s3 = 0.0 ! default !find closest element center to initialize search r2min = 4.01 ! radians DO l_ = 1, num_ele r2 = (b_(1) - center(1,l_))**2 +(b_(2) - center(2,l_))**2 +(b_(3) - center(3,l_))**2 IF (r2 < r2min) THEN r2min = r2 iet = l_ END IF END DO ! If closest element center is more than 8 degrees away, give up. tv = center(1:3, iet) IF (DOT_PRODUCT(b_, tv) < 0.99) RETURN END IF ! cold_start ! initialize search memory (with impossible numbers) back1 = -1 back2 = -2 back3 = -3 is_it_here: DO ! first, check for infinite loop between 2 elements! IF (iet == back2) THEN ! in loop; force location in one or the other! CALL Dumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL Pull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ! then, check for infinite loop between 3 elements! ELSE IF (iet == back3) THEN ! in loop; force location in one or the other! CALL Dumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) s_temp(1) = s1t; s_temp(2) = s2t; s_temp(3) = s3t CALL Pull_in(s_temp) s1t = s_temp(1); s2t = s_temp(2); s3t = s_temp(3) EXIT is_it_here ELSE ! normal operation CALL Dumb_s123 (iet, b_, node, xyz_node, center, a_, & & s1t, s2t, s3t) IF ((s1t < s2t) .AND. (s1t < s3t)) THEN ! s1 is most negative; most critical IF (s1t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(1, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE IF ((s2t < s1t) .AND. (s2t < s3t)) THEN ! s2 is most negative; most critical IF (s2t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(2, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF ELSE ! s3 is most negative; most critical IF (s3t >= 0.) THEN EXIT is_it_here ! success ELSE i = neighbor(3, iet) IF (i > 0) THEN back3 = back2 back2 = back1 back1 = iet iet = i CYCLE is_it_here ELSE RETURN ! fell off edge of grid ENDIF ENDIF END IF END IF ! in/not in a loop END DO is_it_here ! successful completion iele = iet s1 = s1t s2 = s2t s3 = s3t success = .TRUE. END SUBROUTINE Which_Spherical_Triangle END MODULE CSM_Report