diff --git a/DATA/SEMUCB_A3d/hknots2.da b/DATA/SEMUCB_A3d/hknots2.da deleted file mode 100644 index 5b9782ebf..000000000 --- a/DATA/SEMUCB_A3d/hknots2.da +++ /dev/null @@ -1,643 +0,0 @@ -642 - 0.000000E+00 0.900000E+02 4 - 0.540000E+02 0.820706E+02 4 - 0.126000E+03 0.820706E+02 4 - 0.540000E+02 0.741413E+02 4 - 0.900000E+02 0.770568E+02 4 - 0.126000E+03 0.741413E+02 4 - 0.540000E+02 0.662119E+02 4 - 0.766103E+02 0.698679E+02 4 - 0.103390E+03 0.698679E+02 4 - 0.126000E+03 0.662119E+02 4 - 0.540000E+02 0.582825E+02 4 - 0.704980E+02 0.620571E+02 4 - 0.900000E+02 0.634350E+02 4 - 0.109502E+03 0.620571E+02 4 - 0.126000E+03 0.582825E+02 4 - 0.540000E+02 0.503532E+02 4 - 0.671328E+02 0.539614E+02 4 - 0.821231E+02 0.559107E+02 4 - 0.978769E+02 0.559107E+02 4 - 0.112867E+03 0.539614E+02 4 - 0.126000E+03 0.503532E+02 4 - 0.540000E+02 0.424238E+02 4 - 0.650443E+02 0.456838E+02 4 - 0.771890E+02 0.477653E+02 4 - 0.900000E+02 0.484832E+02 4 - 0.102811E+03 0.477653E+02 4 - 0.114956E+03 0.456838E+02 4 - 0.126000E+03 0.424238E+02 4 - 0.540000E+02 0.344944E+02 4 - 0.636327E+02 0.372707E+02 4 - 0.738858E+02 0.392137E+02 4 - 0.845758E+02 0.402160E+02 4 - 0.954242E+02 0.402160E+02 4 - 0.106114E+03 0.392137E+02 4 - 0.116367E+03 0.372707E+02 4 - 0.126000E+03 0.344944E+02 4 - 0.540000E+02 0.265651E+02 4 - 0.626069E+02 0.287551E+02 4 - 0.715330E+02 0.303792E+02 4 - 0.807010E+02 0.313795E+02 4 - 0.900000E+02 0.317175E+02 4 - 0.992990E+02 0.313795E+02 4 - 0.108467E+03 0.303792E+02 4 - 0.117393E+03 0.287551E+02 4 - 0.126000E+03 0.265651E+02 4 - 0.198000E+03 0.820706E+02 4 - 0.162000E+03 0.770568E+02 4 - 0.198000E+03 0.741413E+02 4 - 0.148610E+03 0.698679E+02 4 - 0.175390E+03 0.698679E+02 4 - 0.198000E+03 0.662119E+02 4 - 0.142498E+03 0.620571E+02 4 - 0.162000E+03 0.634350E+02 4 - 0.181502E+03 0.620571E+02 4 - 0.198000E+03 0.582825E+02 4 - 0.139133E+03 0.539614E+02 4 - 0.154123E+03 0.559107E+02 4 - 0.169877E+03 0.559107E+02 4 - 0.184867E+03 0.539614E+02 4 - 0.198000E+03 0.503532E+02 4 - 0.137044E+03 0.456838E+02 4 - 0.149189E+03 0.477653E+02 4 - 0.162000E+03 0.484832E+02 4 - 0.174811E+03 0.477653E+02 4 - 0.186956E+03 0.456838E+02 4 - 0.198000E+03 0.424238E+02 4 - 0.135633E+03 0.372707E+02 4 - 0.145886E+03 0.392137E+02 4 - 0.156576E+03 0.402160E+02 4 - 0.167424E+03 0.402160E+02 4 - 0.178114E+03 0.392137E+02 4 - 0.188367E+03 0.372707E+02 4 - 0.198000E+03 0.344944E+02 4 - 0.134607E+03 0.287551E+02 4 - 0.143533E+03 0.303792E+02 4 - 0.152701E+03 0.313795E+02 4 - 0.162000E+03 0.317175E+02 4 - 0.171299E+03 0.313795E+02 4 - 0.180467E+03 0.303792E+02 4 - 0.189393E+03 0.287551E+02 4 - 0.198000E+03 0.265651E+02 4 - 0.270000E+03 0.820706E+02 4 - 0.234000E+03 0.770568E+02 4 - 0.270000E+03 0.741413E+02 4 - 0.220610E+03 0.698679E+02 4 - 0.247390E+03 0.698679E+02 4 - 0.270000E+03 0.662119E+02 4 - 0.214498E+03 0.620571E+02 4 - 0.234000E+03 0.634350E+02 4 - 0.253502E+03 0.620571E+02 4 - 0.270000E+03 0.582825E+02 4 - 0.211133E+03 0.539614E+02 4 - 0.226123E+03 0.559107E+02 4 - 0.241877E+03 0.559107E+02 4 - 0.256867E+03 0.539614E+02 4 - 0.270000E+03 0.503532E+02 4 - 0.209044E+03 0.456838E+02 4 - 0.221189E+03 0.477653E+02 4 - 0.234000E+03 0.484832E+02 4 - 0.246811E+03 0.477653E+02 4 - 0.258956E+03 0.456838E+02 4 - 0.270000E+03 0.424238E+02 4 - 0.207633E+03 0.372707E+02 4 - 0.217886E+03 0.392137E+02 4 - 0.228576E+03 0.402160E+02 4 - 0.239424E+03 0.402160E+02 4 - 0.250114E+03 0.392137E+02 4 - 0.260367E+03 0.372707E+02 4 - 0.270000E+03 0.344944E+02 4 - 0.206607E+03 0.287551E+02 4 - 0.215533E+03 0.303792E+02 4 - 0.224701E+03 0.313795E+02 4 - 0.234000E+03 0.317175E+02 4 - 0.243299E+03 0.313795E+02 4 - 0.252467E+03 0.303792E+02 4 - 0.261393E+03 0.287551E+02 4 - 0.270000E+03 0.265651E+02 4 - 0.342000E+03 0.820706E+02 4 - 0.306000E+03 0.770568E+02 4 - 0.342000E+03 0.741413E+02 4 - 0.292610E+03 0.698679E+02 4 - 0.319390E+03 0.698679E+02 4 - 0.342000E+03 0.662119E+02 4 - 0.286498E+03 0.620571E+02 4 - 0.306000E+03 0.634350E+02 4 - 0.325502E+03 0.620571E+02 4 - 0.342000E+03 0.582825E+02 4 - 0.283133E+03 0.539614E+02 4 - 0.298123E+03 0.559107E+02 4 - 0.313877E+03 0.559107E+02 4 - 0.328867E+03 0.539614E+02 4 - 0.342000E+03 0.503532E+02 4 - 0.281044E+03 0.456838E+02 4 - 0.293189E+03 0.477653E+02 4 - 0.306000E+03 0.484832E+02 4 - 0.318811E+03 0.477653E+02 4 - 0.330956E+03 0.456838E+02 4 - 0.342000E+03 0.424238E+02 4 - 0.279633E+03 0.372707E+02 4 - 0.289886E+03 0.392137E+02 4 - 0.300576E+03 0.402160E+02 4 - 0.311424E+03 0.402160E+02 4 - 0.322114E+03 0.392137E+02 4 - 0.332367E+03 0.372707E+02 4 - 0.342000E+03 0.344944E+02 4 - 0.278607E+03 0.287551E+02 4 - 0.287533E+03 0.303792E+02 4 - 0.296701E+03 0.313795E+02 4 - 0.306000E+03 0.317175E+02 4 - 0.315299E+03 0.313795E+02 4 - 0.324467E+03 0.303792E+02 4 - 0.333393E+03 0.287551E+02 4 - 0.342000E+03 0.265651E+02 4 - 0.180000E+02 0.770568E+02 4 - 0.461033E+01 0.698679E+02 4 - 0.313897E+02 0.698679E+02 4 - 0.358498E+03 0.620571E+02 4 - 0.180000E+02 0.634350E+02 4 - 0.375020E+02 0.620571E+02 4 - 0.355133E+03 0.539614E+02 4 - 0.101231E+02 0.559107E+02 4 - 0.258769E+02 0.559107E+02 4 - 0.408672E+02 0.539614E+02 4 - 0.353044E+03 0.456838E+02 4 - 0.518898E+01 0.477653E+02 4 - 0.180000E+02 0.484832E+02 4 - 0.308110E+02 0.477653E+02 4 - 0.429557E+02 0.456838E+02 4 - 0.351633E+03 0.372707E+02 4 - 0.188578E+01 0.392137E+02 4 - 0.125758E+02 0.402160E+02 4 - 0.234242E+02 0.402160E+02 4 - 0.341142E+02 0.392137E+02 4 - 0.443673E+02 0.372707E+02 4 - 0.350607E+03 0.287551E+02 4 - 0.359533E+03 0.303792E+02 4 - 0.870101E+01 0.313795E+02 4 - 0.180000E+02 0.317175E+02 4 - 0.272990E+02 0.313795E+02 4 - 0.364670E+02 0.303792E+02 4 - 0.453931E+02 0.287551E+02 4 - 0.589523E+02 0.200667E+02 4 - 0.635057E+02 0.134416E+02 4 - 0.672786E+02 0.219590E+02 4 - 0.678119E+02 0.673914E+01 4 - 0.717123E+02 0.150431E+02 4 - 0.759295E+02 0.232750E+02 4 - 0.720000E+02 0.174962E-06 4 - 0.760516E+02 0.804317E+01 4 - 0.802677E+02 0.160451E+02 4 - 0.848366E+02 0.239577E+02 4 - 0.761881E+02 -0.673914E+01 4 - 0.804220E+02 0.992240E+00 4 - 0.846759E+02 0.871818E+01 4 - 0.891097E+02 0.163948E+02 4 - 0.939068E+02 0.239710E+02 4 - 0.804943E+02 -0.134416E+02 4 - 0.849435E+02 -0.607455E+01 4 - 0.892717E+02 0.132787E+01 4 - 0.936261E+02 0.872267E+01 4 - 0.981565E+02 0.160659E+02 4 - 0.103033E+03 0.233066E+02 4 - 0.850477E+02 -0.200667E+02 4 - 0.897406E+02 -0.131138E+02 4 - 0.941730E+02 -0.608065E+01 4 - 0.984902E+02 0.987660E+00 4 - 0.102826E+03 0.805034E+01 4 - 0.107316E+03 0.150659E+02 4 - 0.112110E+03 0.219872E+02 4 - 0.900000E+02 -0.265651E+02 4 - 0.949523E+02 -0.200667E+02 4 - 0.995057E+02 -0.134416E+02 4 - 0.103812E+03 -0.673914E+01 4 - 0.108000E+03 -0.349924E-06 4 - 0.112188E+03 0.673914E+01 4 - 0.116494E+03 0.134416E+02 4 - 0.121048E+03 0.200667E+02 4 - 0.130952E+03 0.200667E+02 4 - 0.126000E+03 0.136218E+02 4 - 0.135506E+03 0.134416E+02 4 - 0.121395E+03 0.691570E+01 4 - 0.130605E+03 0.691570E+01 4 - 0.139812E+03 0.673914E+01 4 - 0.117000E+03 0.154742E-14 4 - 0.126000E+03 0.103581E-14 4 - 0.135000E+03 0.498698E-15 4 - 0.144000E+03 -0.900393E-08 4 - 0.112681E+03 -0.707892E+01 4 - 0.121559E+03 -0.725093E+01 4 - 0.130441E+03 -0.725093E+01 4 - 0.139319E+03 -0.707892E+01 4 - 0.148188E+03 -0.673914E+01 4 - 0.108293E+03 -0.142730E+02 4 - 0.117133E+03 -0.147809E+02 4 - 0.126000E+03 -0.149518E+02 4 - 0.134867E+03 -0.147809E+02 4 - 0.143707E+03 -0.142730E+02 4 - 0.152494E+03 -0.134416E+02 4 - 0.103671E+03 -0.215244E+02 4 - 0.112538E+03 -0.225217E+02 4 - 0.121501E+03 -0.230284E+02 4 - 0.130499E+03 -0.230284E+02 4 - 0.139462E+03 -0.225217E+02 4 - 0.148329E+03 -0.215244E+02 4 - 0.157048E+03 -0.200667E+02 4 - 0.986069E+02 -0.287551E+02 4 - 0.107533E+03 -0.303792E+02 4 - 0.116701E+03 -0.313795E+02 4 - 0.126000E+03 -0.317175E+02 4 - 0.135299E+03 -0.313795E+02 4 - 0.144467E+03 -0.303792E+02 4 - 0.153393E+03 -0.287551E+02 4 - 0.162000E+03 -0.265651E+02 4 - 0.139279E+03 0.219590E+02 4 - 0.143712E+03 0.150431E+02 4 - 0.147929E+03 0.232750E+02 4 - 0.148052E+03 0.804317E+01 4 - 0.152268E+03 0.160451E+02 4 - 0.156837E+03 0.239577E+02 4 - 0.152422E+03 0.992239E+00 4 - 0.156676E+03 0.871818E+01 4 - 0.161110E+03 0.163948E+02 4 - 0.165907E+03 0.239710E+02 4 - 0.156943E+03 -0.607455E+01 4 - 0.161272E+03 0.132787E+01 4 - 0.165626E+03 0.872267E+01 4 - 0.170157E+03 0.160659E+02 4 - 0.175033E+03 0.233066E+02 4 - 0.161741E+03 -0.131138E+02 4 - 0.166173E+03 -0.608065E+01 4 - 0.170490E+03 0.987660E+00 4 - 0.174826E+03 0.805034E+01 4 - 0.179316E+03 0.150659E+02 4 - 0.184110E+03 0.219872E+02 4 - 0.166952E+03 -0.200667E+02 4 - 0.171506E+03 -0.134416E+02 4 - 0.175812E+03 -0.673914E+01 4 - 0.180000E+03 -0.349924E-06 4 - 0.184188E+03 0.673914E+01 4 - 0.188494E+03 0.134416E+02 4 - 0.193048E+03 0.200667E+02 4 - 0.202952E+03 0.200667E+02 4 - 0.198000E+03 0.136218E+02 4 - 0.207506E+03 0.134416E+02 4 - 0.193395E+03 0.691570E+01 4 - 0.202605E+03 0.691570E+01 4 - 0.211812E+03 0.673914E+01 4 - 0.189000E+03 0.154742E-14 4 - 0.198000E+03 0.103581E-14 4 - 0.207000E+03 0.498698E-15 4 - 0.216000E+03 -0.900391E-08 4 - 0.184681E+03 -0.707892E+01 4 - 0.193559E+03 -0.725093E+01 4 - 0.202441E+03 -0.725093E+01 4 - 0.211319E+03 -0.707892E+01 4 - 0.220188E+03 -0.673914E+01 4 - 0.180293E+03 -0.142730E+02 4 - 0.189133E+03 -0.147809E+02 4 - 0.198000E+03 -0.149518E+02 4 - 0.206867E+03 -0.147809E+02 4 - 0.215707E+03 -0.142730E+02 4 - 0.224494E+03 -0.134416E+02 4 - 0.175671E+03 -0.215244E+02 4 - 0.184538E+03 -0.225217E+02 4 - 0.193501E+03 -0.230284E+02 4 - 0.202499E+03 -0.230284E+02 4 - 0.211462E+03 -0.225217E+02 4 - 0.220329E+03 -0.215244E+02 4 - 0.229048E+03 -0.200667E+02 4 - 0.170607E+03 -0.287551E+02 4 - 0.179533E+03 -0.303792E+02 4 - 0.188701E+03 -0.313795E+02 4 - 0.198000E+03 -0.317175E+02 4 - 0.207299E+03 -0.313795E+02 4 - 0.216467E+03 -0.303792E+02 4 - 0.225393E+03 -0.287551E+02 4 - 0.234000E+03 -0.265651E+02 4 - 0.211279E+03 0.219590E+02 4 - 0.215712E+03 0.150431E+02 4 - 0.219929E+03 0.232750E+02 4 - 0.220052E+03 0.804317E+01 4 - 0.224268E+03 0.160451E+02 4 - 0.228837E+03 0.239577E+02 4 - 0.224422E+03 0.992239E+00 4 - 0.228676E+03 0.871818E+01 4 - 0.233110E+03 0.163948E+02 4 - 0.237907E+03 0.239710E+02 4 - 0.228943E+03 -0.607455E+01 4 - 0.233272E+03 0.132787E+01 4 - 0.237626E+03 0.872267E+01 4 - 0.242157E+03 0.160659E+02 4 - 0.247033E+03 0.233066E+02 4 - 0.233741E+03 -0.131138E+02 4 - 0.238173E+03 -0.608065E+01 4 - 0.242490E+03 0.987659E+00 4 - 0.246826E+03 0.805034E+01 4 - 0.251316E+03 0.150659E+02 4 - 0.256110E+03 0.219872E+02 4 - 0.238952E+03 -0.200667E+02 4 - 0.243506E+03 -0.134416E+02 4 - 0.247812E+03 -0.673914E+01 4 - 0.252000E+03 -0.349924E-06 4 - 0.256188E+03 0.673914E+01 4 - 0.260494E+03 0.134416E+02 4 - 0.265048E+03 0.200667E+02 4 - 0.274952E+03 0.200667E+02 4 - 0.270000E+03 0.136218E+02 4 - 0.279506E+03 0.134416E+02 4 - 0.265395E+03 0.691570E+01 4 - 0.274605E+03 0.691570E+01 4 - 0.283812E+03 0.673914E+01 4 - 0.261000E+03 0.277658E-14 4 - 0.270000E+03 0.240691E-14 4 - 0.279000E+03 0.197797E-14 4 - 0.288000E+03 0.174962E-06 4 - 0.256681E+03 -0.707892E+01 4 - 0.265559E+03 -0.725093E+01 4 - 0.274441E+03 -0.725093E+01 4 - 0.283319E+03 -0.707892E+01 4 - 0.292188E+03 -0.673914E+01 4 - 0.252293E+03 -0.142730E+02 4 - 0.261133E+03 -0.147809E+02 4 - 0.270000E+03 -0.149518E+02 4 - 0.278867E+03 -0.147809E+02 4 - 0.287707E+03 -0.142730E+02 4 - 0.296494E+03 -0.134416E+02 4 - 0.247671E+03 -0.215244E+02 4 - 0.256538E+03 -0.225217E+02 4 - 0.265501E+03 -0.230284E+02 4 - 0.274499E+03 -0.230284E+02 4 - 0.283462E+03 -0.225217E+02 4 - 0.292329E+03 -0.215244E+02 4 - 0.301048E+03 -0.200667E+02 4 - 0.242607E+03 -0.287551E+02 4 - 0.251533E+03 -0.303792E+02 4 - 0.260701E+03 -0.313795E+02 4 - 0.270000E+03 -0.317175E+02 4 - 0.279299E+03 -0.313795E+02 4 - 0.288467E+03 -0.303792E+02 4 - 0.297393E+03 -0.287551E+02 4 - 0.306000E+03 -0.265651E+02 4 - 0.283279E+03 0.219590E+02 4 - 0.287712E+03 0.150431E+02 4 - 0.291929E+03 0.232750E+02 4 - 0.292052E+03 0.804317E+01 4 - 0.296268E+03 0.160451E+02 4 - 0.300837E+03 0.239577E+02 4 - 0.296422E+03 0.992240E+00 4 - 0.300676E+03 0.871818E+01 4 - 0.305110E+03 0.163948E+02 4 - 0.309907E+03 0.239710E+02 4 - 0.300943E+03 -0.607455E+01 4 - 0.305272E+03 0.132787E+01 4 - 0.309626E+03 0.872267E+01 4 - 0.314157E+03 0.160659E+02 4 - 0.319033E+03 0.233066E+02 4 - 0.305741E+03 -0.131138E+02 4 - 0.310173E+03 -0.608065E+01 4 - 0.314490E+03 0.987660E+00 4 - 0.318826E+03 0.805034E+01 4 - 0.323316E+03 0.150659E+02 4 - 0.328110E+03 0.219872E+02 4 - 0.310952E+03 -0.200667E+02 4 - 0.315506E+03 -0.134416E+02 4 - 0.319812E+03 -0.673914E+01 4 - 0.324000E+03 -0.165958E-06 4 - 0.328188E+03 0.673914E+01 4 - 0.332494E+03 0.134416E+02 4 - 0.337048E+03 0.200667E+02 4 - 0.346952E+03 0.200667E+02 4 - 0.342000E+03 0.136218E+02 4 - 0.351506E+03 0.134416E+02 4 - 0.337395E+03 0.691570E+01 4 - 0.346605E+03 0.691570E+01 4 - 0.355812E+03 0.673914E+01 4 - 0.333000E+03 0.346516E-14 4 - 0.342000E+03 0.350835E-14 4 - 0.351000E+03 0.346516E-14 4 - 0.360000E+03 0.174962E-06 4 - 0.328681E+03 -0.707892E+01 4 - 0.337559E+03 -0.725093E+01 4 - 0.346441E+03 -0.725093E+01 4 - 0.355319E+03 -0.707892E+01 4 - 0.418806E+01 -0.673914E+01 4 - 0.324293E+03 -0.142730E+02 4 - 0.333133E+03 -0.147809E+02 4 - 0.342000E+03 -0.149518E+02 4 - 0.350867E+03 -0.147809E+02 4 - 0.359707E+03 -0.142730E+02 4 - 0.849429E+01 -0.134416E+02 4 - 0.319671E+03 -0.215244E+02 4 - 0.328538E+03 -0.225217E+02 4 - 0.337501E+03 -0.230284E+02 4 - 0.346499E+03 -0.230284E+02 4 - 0.355462E+03 -0.225217E+02 4 - 0.432887E+01 -0.215244E+02 4 - 0.130477E+02 -0.200667E+02 4 - 0.314607E+03 -0.287551E+02 4 - 0.323533E+03 -0.303792E+02 4 - 0.332701E+03 -0.313795E+02 4 - 0.342000E+03 -0.317175E+02 4 - 0.351299E+03 -0.313795E+02 4 - 0.466996E+00 -0.303792E+02 4 - 0.939305E+01 -0.287551E+02 4 - 0.180000E+02 -0.265651E+02 4 - 0.355279E+03 0.219590E+02 4 - 0.359712E+03 0.150431E+02 4 - 0.392946E+01 0.232750E+02 4 - 0.405159E+01 0.804317E+01 4 - 0.826770E+01 0.160451E+02 4 - 0.128366E+02 0.239577E+02 4 - 0.842200E+01 0.992240E+00 4 - 0.126759E+02 0.871818E+01 4 - 0.171097E+02 0.163948E+02 4 - 0.219068E+02 0.239710E+02 4 - 0.129435E+02 -0.607455E+01 4 - 0.172717E+02 0.132787E+01 4 - 0.216261E+02 0.872267E+01 4 - 0.261565E+02 0.160659E+02 4 - 0.310333E+02 0.233066E+02 4 - 0.177406E+02 -0.131138E+02 4 - 0.221730E+02 -0.608065E+01 4 - 0.264902E+02 0.987660E+00 4 - 0.308259E+02 0.805034E+01 4 - 0.353156E+02 0.150659E+02 4 - 0.401104E+02 0.219872E+02 4 - 0.229523E+02 -0.200667E+02 4 - 0.275057E+02 -0.134416E+02 4 - 0.318119E+02 -0.673914E+01 4 - 0.360000E+02 -0.165958E-06 4 - 0.401881E+02 0.673914E+01 4 - 0.444943E+02 0.134416E+02 4 - 0.490477E+02 0.200667E+02 4 - 0.540000E+02 0.136218E+02 4 - 0.493949E+02 0.691570E+01 4 - 0.586051E+02 0.691570E+01 4 - 0.450000E+02 0.346516E-14 4 - 0.540000E+02 0.350835E-14 4 - 0.630000E+02 0.346516E-14 4 - 0.406806E+02 -0.707892E+01 4 - 0.495591E+02 -0.725093E+01 4 - 0.584409E+02 -0.725093E+01 4 - 0.673194E+02 -0.707892E+01 4 - 0.362931E+02 -0.142730E+02 4 - 0.451328E+02 -0.147809E+02 4 - 0.540000E+02 -0.149518E+02 4 - 0.628672E+02 -0.147809E+02 4 - 0.717069E+02 -0.142730E+02 4 - 0.316711E+02 -0.215244E+02 4 - 0.405380E+02 -0.225217E+02 4 - 0.495014E+02 -0.230284E+02 4 - 0.584986E+02 -0.230284E+02 4 - 0.674620E+02 -0.225217E+02 4 - 0.763289E+02 -0.215244E+02 4 - 0.266069E+02 -0.287551E+02 4 - 0.355330E+02 -0.303792E+02 4 - 0.447010E+02 -0.313795E+02 4 - 0.540000E+02 -0.317175E+02 4 - 0.632990E+02 -0.313795E+02 4 - 0.724670E+02 -0.303792E+02 4 - 0.813931E+02 -0.287551E+02 4 - 0.180000E+02 -0.344944E+02 4 - 0.180000E+02 -0.424238E+02 4 - 0.274537E+02 -0.367217E+02 4 - 0.180000E+02 -0.503532E+02 4 - 0.287394E+02 -0.447551E+02 4 - 0.374578E+02 -0.383360E+02 4 - 0.180000E+02 -0.582825E+02 4 - 0.307264E+02 -0.528389E+02 4 - 0.403862E+02 -0.463531E+02 4 - 0.479052E+02 -0.392315E+02 4 - 0.180000E+02 -0.662119E+02 4 - 0.339725E+02 -0.609365E+02 4 - 0.449482E+02 -0.543402E+02 4 - 0.527469E+02 -0.470314E+02 4 - 0.586161E+02 -0.393232E+02 4 - 0.180000E+02 -0.741413E+02 4 - 0.398766E+02 -0.689502E+02 4 - 0.524440E+02 -0.620982E+02 4 - 0.601856E+02 -0.545371E+02 4 - 0.654533E+02 -0.466386E+02 4 - 0.693592E+02 -0.385606E+02 4 - 0.180000E+02 -0.820706E+02 4 - 0.529139E+02 -0.765465E+02 4 - 0.657206E+02 -0.691330E+02 4 - 0.718770E+02 -0.612570E+02 4 - 0.755345E+02 -0.532149E+02 4 - 0.780252E+02 -0.450965E+02 4 - 0.798909E+02 -0.369375E+02 4 - 0.000000E+00 -0.900000E+02 4 - 0.900000E+02 -0.820706E+02 4 - 0.900000E+02 -0.741413E+02 4 - 0.900000E+02 -0.662119E+02 4 - 0.900000E+02 -0.582825E+02 4 - 0.900000E+02 -0.503532E+02 4 - 0.900000E+02 -0.424238E+02 4 - 0.900000E+02 -0.344944E+02 4 - 0.994537E+02 -0.367217E+02 4 - 0.100739E+03 -0.447551E+02 4 - 0.109458E+03 -0.383360E+02 4 - 0.102726E+03 -0.528389E+02 4 - 0.112386E+03 -0.463531E+02 4 - 0.119905E+03 -0.392315E+02 4 - 0.105972E+03 -0.609365E+02 4 - 0.116948E+03 -0.543402E+02 4 - 0.124747E+03 -0.470314E+02 4 - 0.130616E+03 -0.393232E+02 4 - 0.111877E+03 -0.689502E+02 4 - 0.124444E+03 -0.620982E+02 4 - 0.132186E+03 -0.545371E+02 4 - 0.137453E+03 -0.466386E+02 4 - 0.141359E+03 -0.385606E+02 4 - 0.124914E+03 -0.765465E+02 4 - 0.137721E+03 -0.691330E+02 4 - 0.143877E+03 -0.612570E+02 4 - 0.147535E+03 -0.532149E+02 4 - 0.150025E+03 -0.450965E+02 4 - 0.151891E+03 -0.369375E+02 4 - 0.162000E+03 -0.820706E+02 4 - 0.162000E+03 -0.741413E+02 4 - 0.162000E+03 -0.662119E+02 4 - 0.162000E+03 -0.582825E+02 4 - 0.162000E+03 -0.503532E+02 4 - 0.162000E+03 -0.424238E+02 4 - 0.162000E+03 -0.344944E+02 4 - 0.171454E+03 -0.367217E+02 4 - 0.172739E+03 -0.447551E+02 4 - 0.181458E+03 -0.383360E+02 4 - 0.174726E+03 -0.528389E+02 4 - 0.184386E+03 -0.463531E+02 4 - 0.191905E+03 -0.392315E+02 4 - 0.177972E+03 -0.609365E+02 4 - 0.188948E+03 -0.543402E+02 4 - 0.196747E+03 -0.470314E+02 4 - 0.202616E+03 -0.393232E+02 4 - 0.183877E+03 -0.689502E+02 4 - 0.196444E+03 -0.620982E+02 4 - 0.204186E+03 -0.545371E+02 4 - 0.209453E+03 -0.466386E+02 4 - 0.213359E+03 -0.385606E+02 4 - 0.196914E+03 -0.765465E+02 4 - 0.209721E+03 -0.691330E+02 4 - 0.215877E+03 -0.612570E+02 4 - 0.219535E+03 -0.532149E+02 4 - 0.222025E+03 -0.450965E+02 4 - 0.223891E+03 -0.369375E+02 4 - 0.234000E+03 -0.820706E+02 4 - 0.234000E+03 -0.741413E+02 4 - 0.234000E+03 -0.662119E+02 4 - 0.234000E+03 -0.582825E+02 4 - 0.234000E+03 -0.503532E+02 4 - 0.234000E+03 -0.424238E+02 4 - 0.234000E+03 -0.344944E+02 4 - 0.243454E+03 -0.367217E+02 4 - 0.244739E+03 -0.447551E+02 4 - 0.253458E+03 -0.383360E+02 4 - 0.246726E+03 -0.528389E+02 4 - 0.256386E+03 -0.463531E+02 4 - 0.263905E+03 -0.392315E+02 4 - 0.249972E+03 -0.609365E+02 4 - 0.260948E+03 -0.543402E+02 4 - 0.268747E+03 -0.470314E+02 4 - 0.274616E+03 -0.393232E+02 4 - 0.255877E+03 -0.689502E+02 4 - 0.268444E+03 -0.620982E+02 4 - 0.276186E+03 -0.545371E+02 4 - 0.281453E+03 -0.466386E+02 4 - 0.285359E+03 -0.385606E+02 4 - 0.268914E+03 -0.765465E+02 4 - 0.281721E+03 -0.691330E+02 4 - 0.287877E+03 -0.612570E+02 4 - 0.291534E+03 -0.532149E+02 4 - 0.294025E+03 -0.450965E+02 4 - 0.295891E+03 -0.369375E+02 4 - 0.306000E+03 -0.820706E+02 4 - 0.306000E+03 -0.741413E+02 4 - 0.306000E+03 -0.662119E+02 4 - 0.306000E+03 -0.582825E+02 4 - 0.306000E+03 -0.503532E+02 4 - 0.306000E+03 -0.424238E+02 4 - 0.306000E+03 -0.344944E+02 4 - 0.315454E+03 -0.367217E+02 4 - 0.316739E+03 -0.447551E+02 4 - 0.325458E+03 -0.383360E+02 4 - 0.318726E+03 -0.528389E+02 4 - 0.328386E+03 -0.463531E+02 4 - 0.335905E+03 -0.392315E+02 4 - 0.321972E+03 -0.609365E+02 4 - 0.332948E+03 -0.543402E+02 4 - 0.340747E+03 -0.470314E+02 4 - 0.346616E+03 -0.393232E+02 4 - 0.327877E+03 -0.689502E+02 4 - 0.340444E+03 -0.620982E+02 4 - 0.348186E+03 -0.545371E+02 4 - 0.353453E+03 -0.466386E+02 4 - 0.357359E+03 -0.385606E+02 4 - 0.340914E+03 -0.765465E+02 4 - 0.353721E+03 -0.691330E+02 4 - 0.359877E+03 -0.612570E+02 4 - 0.353451E+01 -0.532149E+02 4 - 0.602516E+01 -0.450965E+02 4 - 0.789086E+01 -0.369375E+02 4 diff --git a/EXAMPLES/regional_Berkeley/DATA/Par_file b/EXAMPLES/regional_Berkeley/DATA/Par_file index 488279249..0e8ee669a 100644 --- a/EXAMPLES/regional_Berkeley/DATA/Par_file +++ b/EXAMPLES/regional_Berkeley/DATA/Par_file @@ -81,7 +81,7 @@ FULL_GRAVITY = .false. POISSON_SOLVER = 0 # record length in minutes -RECORD_LENGTH_IN_MINUTES = 2.5d0 +RECORD_LENGTH_IN_MINUTES = 5.0d0 #----------------------------------------------------------- # @@ -269,6 +269,17 @@ USE_MONOCHROMATIC_CMT_SOURCE = .false. # print source time function PRINT_SOURCE_TIME_FUNCTION = .true. +## Berkeley source time function +STF_IS_UCB_HEAVISIDE = .true. +# UCB source frequency content (i.e., heaviside function) +SOURCE_T1 = 500.d0 +SOURCE_T2 = 400.d0 +SOURCE_T3 = 50.d0 +SOURCE_T4 = 40.d0 +# UCB source time shift +TAU = 500.d0 + + #----------------------------------------------------------- # # Seismograms @@ -347,7 +358,7 @@ APPROXIMATE_HESS_KL = .false. # forces transverse isotropy for all mantle elements # (default is to use transverse isotropy only between crust and 220) # means we allow radial anisotropy throughout the whole crust/mantle region -USE_FULL_TISO_MANTLE = .false. +USE_FULL_TISO_MANTLE = .true. # output kernel mask to zero out source region # to remove large values near the sources in the sensitivity kernels diff --git a/setup/constants.h.in b/setup/constants.h.in index 79247af01..46aadc93a 100644 --- a/setup/constants.h.in +++ b/setup/constants.h.in @@ -245,12 +245,6 @@ !! !!----------------------------------------------------------- -! source-time function is a UCB-style filtered heaviside - logical, parameter :: STF_IS_UCB_HEAVISIDE = .false. - -! Save mirror files - logical, parameter :: SAVE_MIRRORS = .false. - ! Path to A3d model ! Make sure this line ends on "/" character (len=*), parameter :: A3d_folder = "DATA/SEMUCB_A3d/" @@ -266,13 +260,19 @@ ! pathname of the topography file (un-smoothed) character (len=*), parameter :: PATHNAME_TOPO_FILE_BERKELEY = trim(A3d_folder)//"ETOPO5_1x1_filtre.dat" -!! uncomment this to use Berkeley topography as default - and comment out corresponding parameters above !! +!! uncomment this only to use Berkeley topography as default for all other Earth models as well +!! note: Berkeley topography setting is used by default for SEMUCB model (see in get_model_parameters.F90 file). +!! uncommenting the lines below here (- and commenting out corresponding parameters in the topo/bathy section above) +!! will force the Berkeley topo to be used as the default topography for all Earth models. !! Topography defaults to Berkeley ! integer, parameter :: EARTH_NX_BATHY = NX_BATHY_BERKELEY ! integer, parameter :: EARTH_NY_BATHY = NY_BATHY_BERKELEY ! double precision, parameter :: EARTH_RESOLUTION_TOPO_FILE = RESOLUTION_TOPO_FILE_BERKELEY ! character (len=*), parameter :: EARTH_PATHNAME_TOPO_FILE = PATHNAME_TOPO_FILE_BERKELEY +! Save mirror files - not implemented yet +! logical, parameter :: SAVE_MIRRORS = .false. + !!----------------------------------------------------------- !! diff --git a/src/gpu/compute_stacey_acoustic_gpu.c b/src/gpu/compute_stacey_acoustic_gpu.c index b700250b2..2f04bd763 100644 --- a/src/gpu/compute_stacey_acoustic_gpu.c +++ b/src/gpu/compute_stacey_acoustic_gpu.c @@ -76,7 +76,7 @@ void FC_FUNC_ (compute_stacey_acoustic_gpu, size_t local_work_size[2]; cl_uint idx = 0; - //daniel debug + // debug //clCheck (clFinish (mocl.command_queue)); //printf ("rank %d: stacey a %i, %i save %i num blocks x/y= %i %i nglob %i nspec2D %i nspec %i\n", // mp->myrank,interface_type,num_abs_boundary_faces,mp->save_forward,num_blocks_x,num_blocks_y, @@ -104,7 +104,7 @@ void FC_FUNC_ (compute_stacey_acoustic_gpu, clCheck (clEnqueueNDRangeKernel (mocl.command_queue, mocl.kernels.compute_stacey_acoustic_kernel, 2, NULL, global_work_size, local_work_size, 0, NULL, NULL)); - //daniel debug + // debug //clCheck (clFinish (mocl.command_queue)); //printf ("rank %d: stacey b %i, %i \n", mp->myrank,interface_type,num_abs_boundary_faces); //fflush (stdout); diff --git a/src/meshfem3D/add_topography.f90 b/src/meshfem3D/add_topography.f90 index 0a7a526ef..edc4f90e7 100644 --- a/src/meshfem3D/add_topography.f90 +++ b/src/meshfem3D/add_topography.f90 @@ -31,6 +31,10 @@ subroutine add_topography(xelm,yelm,zelm,ibathy_topo) use shared_parameters, only: REGIONAL_MESH_CUTOFF,REGIONAL_MESH_CUTOFF_DEPTH,USE_LOCAL_MESH,ELLIPTICITY use meshfem_par, only: R220,NX_BATHY,NY_BATHY,R_PLANET + ! for old version Berkeley compatibility + use constants, only: USE_OLD_VERSION_FORMAT,ICRUST_BERKELEY,THREE_D_MODEL_BERKELEY + use meshfem_models_par, only: THREE_D_MODEL,REFERENCE_CRUSTAL_MODEL + implicit none double precision,dimension(NGNOD), intent(inout) :: xelm,yelm,zelm @@ -45,6 +49,13 @@ subroutine add_topography(xelm,yelm,zelm,ibathy_topo) integer :: ia + ! for compatibility + double precision :: theta,phi + double precision :: vpvc,vphc,vsvc,vshc,etac,rhoc + double precision :: moho + double precision :: rmoho + logical :: found_crust,elem_in_crust,moho_only + ! we loop on all the points of the element do ia = 1,NGNOD @@ -76,6 +87,29 @@ subroutine add_topography(xelm,yelm,zelm,ibathy_topo) endif gamma = (r - rbottom) / (R_UNIT_SPHERE - rbottom) + ! old version compatility + if (USE_OLD_VERSION_FORMAT) then + ! Berkeley model + if (THREE_D_MODEL == THREE_D_MODEL_BERKELEY .and. & + REFERENCE_CRUSTAL_MODEL == ICRUST_BERKELEY) then + ! convert lat/lon to theta/phi + call latlon_2_thetaphi_dble(lat,lon,theta,phi) + + ! gets smoothed moho depth + elem_in_crust = .true. + moho_only = .true. + call model_berkeley_crust_aniso(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,found_crust,elem_in_crust,moho_only) + + rmoho = R_UNIT_SPHERE - moho + + ! if point above moho then move points, otherwise skip + if (r <= rmoho) cycle + + ! adjust gamma stretching to moho boundary distance + gamma = (r - rmoho) / (R_UNIT_SPHERE - rmoho) + endif + endif + ! add elevation to all the points of that element ! also make sure gamma makes sense if (gamma < -0.02 .or. gamma > 1.02) then @@ -107,7 +141,6 @@ subroutine add_topography_gll(xstore,ystore,zstore,ispec,nspec,ibathy_topo) use shared_parameters, only: REGIONAL_MESH_CUTOFF,REGIONAL_MESH_CUTOFF_DEPTH,USE_LOCAL_MESH,ELLIPTICITY use meshfem_par, only: R220,NX_BATHY,NY_BATHY - implicit none ! input parameters diff --git a/src/meshfem3D/compute_element_properties.f90 b/src/meshfem3D/compute_element_properties.f90 index 96131ceab..69af54ce3 100644 --- a/src/meshfem3D/compute_element_properties.f90 +++ b/src/meshfem3D/compute_element_properties.f90 @@ -188,7 +188,7 @@ subroutine compute_element_properties(ispec,iregion_code,idoubling,ipass, & endif ! IREGION_CRUST_MANTLE ! sets element tiso flag - call compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,ispec,nspec,idoubling) + call compute_element_tiso_flag(elem_is_tiso,elem_in_crust,elem_in_mantle,iregion_code,ispec,nspec,idoubling) ! stores as element flags ispec_is_tiso(ispec) = elem_is_tiso @@ -473,14 +473,14 @@ end subroutine compute_element_GLL_locations ! - subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,ispec,nspec,idoubling) + subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_crust,elem_in_mantle,iregion_code,ispec,nspec,idoubling) ! sets transverse isotropic flag for elements in crust/mantle use constants, only: IMAIN,myrank, & IFLAG_CRUST,IFLAG_220_80,IFLAG_80_MOHO,IFLAG_670_220,IFLAG_MANTLE_NORMAL,IREGION_CRUST_MANTLE, & REFERENCE_MODEL_1DREF,REFERENCE_MODEL_1DREF,REFERENCE_MODEL_SEMUCB, & - THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE, & + THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE,THREE_D_MODEL_BERKELEY, & USE_OLD_VERSION_FORMAT use meshfem_models_par, only: & @@ -489,7 +489,7 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is implicit none logical,intent(out) :: elem_is_tiso - logical,intent(in) :: elem_in_mantle + logical,intent(in) :: elem_in_crust,elem_in_mantle integer,intent(in) :: iregion_code,ispec integer,intent(in) :: nspec integer,dimension(nspec),intent(in) :: idoubling @@ -501,12 +501,11 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is elem_is_tiso = .false. ! checks if anything to do - if (.not. TRANSVERSE_ISOTROPY) return if (iregion_code /= IREGION_CRUST_MANTLE) return + if (.not. TRANSVERSE_ISOTROPY) return ! user output if (is_first_call) then - is_first_call = .false. if (myrank == 0) then ! only output once write(IMAIN,*) ' setting tiso flags in mantle model' @@ -516,6 +515,8 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is write(IMAIN,*) ' using formatting from old versions (7.0 to 8.0)' call flush_IMAIN() endif + ! turns off flag to show output only once + is_first_call = .false. endif ! transverse isotropic models @@ -526,6 +527,10 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is ! note: this will increase the computation time by ~ 45 % if (USE_OLD_VERSION_FORMAT) then if (elem_in_mantle) elem_is_tiso = .true. + ! adds special case for Berkeley SEMUCB model + if (THREE_D_MODEL == THREE_D_MODEL_BERKELEY) then + if (elem_in_crust) elem_is_tiso = .true. + endif else if (idoubling(ispec) == IFLAG_MANTLE_NORMAL & .or. idoubling(ispec) == IFLAG_670_220 & @@ -557,8 +562,7 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is if (idoubling(ispec) == IFLAG_670_220 & .or. idoubling(ispec) == IFLAG_220_80 & .or. idoubling(ispec) == IFLAG_80_MOHO & - .or. (idoubling(ispec) == IFLAG_CRUST & - .and. elem_in_mantle) ) then + .or. (idoubling(ispec) == IFLAG_CRUST .and. elem_in_mantle)) then elem_is_tiso = .true. endif else @@ -572,14 +576,30 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is endif case (REFERENCE_MODEL_SEMUCB) - ! SEMUCB - allows tiso for full mantle & crust - ! (same as USE_FULL_TISO_MANTLE option) - if (idoubling(ispec) == IFLAG_MANTLE_NORMAL & - .or. idoubling(ispec) == IFLAG_670_220 & - .or. idoubling(ispec) == IFLAG_220_80 & - .or. idoubling(ispec) == IFLAG_80_MOHO & - .or. idoubling(ispec) == IFLAG_CRUST) then - elem_is_tiso = .true. + ! SEMUCB + if (USE_OLD_VERSION_FORMAT) then + ! default from version 6.0 + if (idoubling(ispec) == IFLAG_220_80 & + .or. idoubling(ispec) == IFLAG_80_MOHO) then + ! models use only transverse isotropy between moho and 220 km depth + elem_is_tiso = .true. + endif + ! or additional setting for SEMUCB was to use USE_FULL_TISO_MANTLE + ! + ! note: this sets tiso for elements fully in mantle or fully in crust + ! however, this will still have tiso == .false. for elements with the moho inside going through the element, + ! having corner nodes both in mantle and crust (above and below actual moho) + !if (elem_in_crust) elem_is_tiso = .true. + !if (elem_in_mantle) elem_is_tiso = .true. + else + ! allows tiso for full mantle & crust + if (idoubling(ispec) == IFLAG_MANTLE_NORMAL & + .or. idoubling(ispec) == IFLAG_670_220 & + .or. idoubling(ispec) == IFLAG_220_80 & + .or. idoubling(ispec) == IFLAG_80_MOHO & + .or. idoubling(ispec) == IFLAG_CRUST) then + elem_is_tiso = .true. + endif endif case default @@ -639,6 +659,11 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is ! note: THREE_D_MODEL_SGLOBE_ISO ! sgloberani_iso model based on PREM, it will have tiso already set from crust down to 220 + case (THREE_D_MODEL_BERKELEY) + ! SEMUCB + ! double-check to use tiso for elements fully in crust + if (elem_in_crust) elem_is_tiso = .true. + case default ! nothing special to add continue diff --git a/src/meshfem3D/create_mass_matrices.f90 b/src/meshfem3D/create_mass_matrices.f90 index e59a342ec..4604f7b99 100644 --- a/src/meshfem3D/create_mass_matrices.f90 +++ b/src/meshfem3D/create_mass_matrices.f90 @@ -543,12 +543,6 @@ subroutine create_mass_matrices_ocean_load(ibool,xstore,ystore,zstore,NSPEC2D_TO ! checks if anything to do if (.not. OCEANS) return - ! user output - if (myrank == 0) then - write(IMAIN,*) ' updates mass matrix with ocean load' - call flush_IMAIN() - endif - ! initializes do_ocean_load = .false. @@ -556,6 +550,17 @@ subroutine create_mass_matrices_ocean_load(ibool,xstore,ystore,zstore,NSPEC2D_TO ! for 3D Earth with topography, compute local height of oceans if (TOPOGRAPHY) do_ocean_load = .true. + ! user output + if (myrank == 0) then + write(IMAIN,*) ' updates mass matrix with ocean load' + if (do_ocean_load) then + write(IMAIN,*) ' using minimum ocean thickness: ',MINIMUM_THICKNESS_3D_OCEANS,'(m)' + else + write(IMAIN,*) ' using 1D ocean PREM thickness: ',THICKNESS_OCEANS_PREM * EARTH_R,'(m)' + endif + call flush_IMAIN() + endif + ! create ocean load mass matrix for degrees of freedom at ocean bottom rmass_ocean_load(:) = 0._CUSTOM_REAL diff --git a/src/meshfem3D/meshfem3D_models.F90 b/src/meshfem3D/meshfem3D_models.F90 index 2f5363661..d6196224e 100644 --- a/src/meshfem3D/meshfem3D_models.F90 +++ b/src/meshfem3D/meshfem3D_models.F90 @@ -924,7 +924,8 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, & case (THREE_D_MODEL_BERKELEY) ! 3D Berkeley Model SEMUCB ! note: passes r/theta/phi (geocentric coordinates) - call model_berkeley_shsv(r_used,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,iregion_code,CRUSTAL) + call model_berkeley_shsv(r_used,theta,phi,vpv,vph,vsv,vsh,rho, & + dvsh,dvsv,dvph,dvpv,drho,eta_aniso,iregion_code,CRUSTAL) ! updates velocities from reference model if (TRANSVERSE_ISOTROPY) then @@ -1625,9 +1626,9 @@ subroutine meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc, & ! Berkeley crustal model ! note: passes r/theta/phi (geocentric coordinates) ! Berkeley crustal model is referencing geocentric positions in a spherical Earth frame - call model_berkeley_crust_aniso(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,found_crust) + call model_berkeley_crust_aniso(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,found_crust,elem_in_crust,moho_only) ! old version - isotropic crustal velocities - !call model_berkeley_crust(r,theta,phi,vpc,vsc,rhoc,moho,found_crust) + !call model_berkeley_crust(r,theta,phi,vpc,vsc,rhoc,moho,found_crust,elem_in_crust,moho_only) !vpvc = vpc !vphc = vpc !vsvc = vsc diff --git a/src/meshfem3D/model_1dberkeley.f90 b/src/meshfem3D/model_1dberkeley.f90 index 9d580fa17..1d6b132c5 100644 --- a/src/meshfem3D/model_1dberkeley.f90 +++ b/src/meshfem3D/model_1dberkeley.f90 @@ -243,17 +243,26 @@ subroutine model_1dberkeley(x,rho,vpv,vph,vsv,vsh,eta,Qkappa,Qmu,iregion_code,CR enddo ! make sure we stay in the right region - if (mimic_native_specfem .and. iregion_code == IREGION_INNER_CORE .and. i > NR_inner_core_berk) i = NR_inner_core_berk - - if (mimic_native_specfem .and. iregion_code == IREGION_OUTER_CORE .and. i < NR_inner_core_berk+2) i = NR_inner_core_berk+2 - if (mimic_native_specfem .and. iregion_code == IREGION_OUTER_CORE .and. i > NR_outer_core_berk) i = NR_outer_core_berk - - if (mimic_native_specfem .and. iregion_code == IREGION_CRUST_MANTLE .and. i < NR_outer_core_berk+2) i = NR_outer_core_berk+2 + if (mimic_native_specfem) then + ! inner core bounds + if (iregion_code == IREGION_INNER_CORE) then + if (i > NR_inner_core_berk) i = NR_inner_core_berk + endif + ! outer core bounds + if (iregion_code == IREGION_OUTER_CORE) then + if (i < NR_inner_core_berk+2) i = NR_inner_core_berk+2 + if (i > NR_outer_core_berk) i = NR_outer_core_berk + endif + ! crust/mantle bounds + if (iregion_code == IREGION_CRUST_MANTLE) then + if (i < NR_outer_core_berk+2) i = NR_outer_core_berk+2 + endif - ! if crustal model is used, mantle gets expanded up to surface - ! for any depth less than 24.4 km, values from mantle below moho are taken - if (mimic_native_specfem .and. CRUSTAL .and. i > modemohoberk) then - i = modemohoberk ! Warning : may need to be changed if file is modified ! + ! if crustal model is used, mantle gets expanded up to surface + ! for any depth less than 24.4 km, values from mantle below moho are taken + if (CRUSTAL) then + if (i > modemohoberk) i = modemohoberk ! Warning : may need to be changed if file is modified ! + endif endif if (i == 1) then @@ -270,6 +279,7 @@ subroutine model_1dberkeley(x,rho,vpv,vph,vsv,vsh,eta,Qkappa,Qmu,iregion_code,CR ! interpolates between one layer below to actual radius layer, ! that is from radius_ref(i-1) to r using the values at i-1 and i frac = (r-Mref_V_radius_berkeley(i-1))/(Mref_V_radius_berkeley(i)-Mref_V_radius_berkeley(i-1)) + ! interpolated model parameters rho = Mref_V_density_berkeley(i-1) + frac * (Mref_V_density_berkeley(i)- Mref_V_density_berkeley(i-1)) vpv = Mref_V_vpv_berkeley(i-1) + frac * (Mref_V_vpv_berkeley(i) - Mref_V_vpv_berkeley(i-1) ) @@ -283,11 +293,13 @@ subroutine model_1dberkeley(x,rho,vpv,vph,vsv,vsh,eta,Qkappa,Qmu,iregion_code,CR ! make sure Vs is zero in the outer core even if roundoff errors on depth ! also set fictitious attenuation to a very high value (attenuation is not used in the fluid) - if (mimic_native_specfem .and. iregion_code == IREGION_OUTER_CORE) then - vsv = 0.d0 - vsh = 0.d0 - Qkappa = 3000.d0 - Qmu = 3000.d0 + if (mimic_native_specfem) then + if (iregion_code == IREGION_OUTER_CORE) then + vsv = 0.d0 + vsh = 0.d0 + Qkappa = 3000.d0 + Qmu = 3000.d0 + endif endif ! non-dimensionalize @@ -309,21 +321,23 @@ end subroutine model_1dberkeley !! Utpal Kumar, Feb, 2022 !! subroutine to decide whether the moho node has already been computed - subroutine est_moho_node(estmohonode) - - use model_1dberkeley_par, only: modemohoberk - - implicit none - integer :: estmohonode - - if (modemohoberk < 0) then - call determine_1dberkeley_moho_node(estmohonode) - ! print *,"Determining Moho node ",estmohonode - else - estmohonode = modemohoberk - endif - - end subroutine est_moho_node +! not used yet... +! +! subroutine est_moho_node(estmohonode) +! +! use model_1dberkeley_par, only: modemohoberk +! +! implicit none +! integer :: estmohonode +! +! if (modemohoberk < 0) then +! call determine_1dberkeley_moho_node(estmohonode) +! ! print *,"Determining Moho node ",estmohonode +! else +! estmohonode = modemohoberk +! endif +! +! end subroutine est_moho_node ! !-------------------------------------------------------------------------------------------------- @@ -388,7 +402,7 @@ subroutine determine_1dberkeley_moho_node(mohonodebrk) read(FID,*) !! skip the header enddo do i = 1, num_lines - read(FID,*) radius(i), density(i), vpv(i), vsv(i), qkappa(i),qmu(i), vph(i),vsh(i),eta(i) !Read the data + read(FID,*) radius(i), density(i), vpv(i), vsv(i), qkappa(i), qmu(i), vph(i), vsh(i), eta(i) !Read the data enddo close(FID) @@ -501,7 +515,7 @@ subroutine determine_1dberkeley_moho_radius(moho_radius) read(FID,*) !! skip the header enddo do i = 1, num_lines - read(FID,*) radius(i), density(i), vpv(i), vsv(i), qkappa(i),qmu(i), vph(i),vsh(i),eta(i) !Read the data + read(FID,*) radius(i), density(i), vpv(i), vsv(i), qkappa(i), qmu(i), vph(i), vsh(i), eta(i) !Read the data enddo close(FID) diff --git a/src/meshfem3D/model_attenuation_gll.f90 b/src/meshfem3D/model_attenuation_gll.f90 index d69a51854..82cabb9e0 100644 --- a/src/meshfem3D/model_attenuation_gll.f90 +++ b/src/meshfem3D/model_attenuation_gll.f90 @@ -39,15 +39,25 @@ module model_gll_qmu_par ! GLL model_variables type model_gll_qmu_variables - sequence + !TODO: check if `sequence` is needed + ! + ! in principle, `sequence` is only needed to explicity map the memory structure, for example when + ! the variable like `MGLL_V` gets passed as a routine argument, making sure it gets accessed correctly. + ! however, here we don't pass the objects MGLL_V, MGLL_V_IC to routines, but only use the objects to store arrays + ! within this module. + ! + !sequence + ! tomographic iteration model on GLL points real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable :: qmu_new ! number of elements (crust/mantle elements) integer :: nspec - integer :: dummy_pad ! padding 4 bytes to align the structure + ! in case we use `sequence` to align memory addresses + !integer :: dummy_pad ! padding 4 bytes to align the structure end type model_gll_qmu_variables + type (model_gll_qmu_variables) :: MGLL_QMU_V end module model_gll_qmu_par diff --git a/src/meshfem3D/model_berkeley.f90 b/src/meshfem3D/model_berkeley.f90 index 8f5071e76..bc414d4a3 100644 --- a/src/meshfem3D/model_berkeley.f90 +++ b/src/meshfem3D/model_berkeley.f90 @@ -48,11 +48,19 @@ module model_berkeley_par implicit none - double precision, dimension(:), allocatable :: mdl,kntrad,aknot,oknot,aknot2,oknot2 + ! spline arrays + !double precision, dimension(:), allocatable :: aknot,oknot,aknot2,oknot2 ! old: only for getdel() + double precision, dimension(:,:), allocatable :: knot_coeff,knot_coeff2 + + double precision, dimension(:), allocatable :: mdl integer, dimension(:), allocatable :: level,level2 + ! spline node radii + real, dimension(:), allocatable :: kntrad,kntrad_hh ! real arrays - to avoid conversion in fspl() to float + + ! parameterization integer, parameter :: MAXPAR = 4 - integer, dimension(MAXPAR) :: nknotA1,nknotA2 + integer, dimension(MAXPAR) :: nknotA1 = 0,nknotA2 = 0 character(1), dimension(:), allocatable :: parblock integer :: npar,ndisc,surface,NBPARAM @@ -62,6 +70,10 @@ module model_berkeley_par ! radius of moho discontinuity from reference 1D model (in km) double precision :: moho1D_radius = -1.0 + ! work arrays + double precision, dimension(:), allocatable :: work_dh + integer, dimension(:), allocatable :: work_kindex + end module model_berkeley_par ! @@ -73,17 +85,29 @@ subroutine model_berkeley_broadcast() ! standard routine to setup model use model_berkeley_par - use constants, only: A3d_folder,myrank,EARTH_R_KM + use constants, only: A3d_folder,myrank,IMAIN,EARTH_R_KM,PI implicit none integer,parameter :: unit1 = 51,unit2 = 52 integer :: dum2,i,j,k,n,mdim,ier + integer :: size_work,size_knots + double precision, dimension(:), allocatable :: aknot,oknot,aknot2,oknot2 ! temporary for reading + double precision :: theta,phi character :: trash + character(len=100), parameter :: A3d_dat = trim(A3d_folder) // 'A3d.dat' character(len=100), parameter :: hknots_dat = trim(A3d_folder) // 'hknots.dat' character(len=100), parameter :: hknots2_dat = trim(A3d_folder) // 'hknots2.dat' + double precision, parameter :: deg2rad = PI / 180.d0 + + ! user info + if (myrank == 0) then + write(IMAIN,*) 'broadcast model: SEMUCB Berkeley' + call flush_IMAIN() + endif + ! determine moho radius from 1D reference model if (myrank == 0) then ! gets exact 1D moho radius (in km) @@ -113,12 +137,18 @@ subroutine model_berkeley_broadcast() endif read(unit1,*) nknotA2(1) - if (allocated(oknot).or.allocated(aknot).or.allocated(level)) then - print *,'A3d already initiated' + + if (allocated(oknot) .or. allocated(aknot) .or. allocated(level)) then + print *,'[model_berkeley_broadcast] A3d already initiated' return endif - allocate(oknot(nknotA2(1)),aknot(nknotA2(1)),level(nknotA2(1))) + allocate(oknot(nknotA2(1)),aknot(nknotA2(1)),level(nknotA2(1)),stat=ier) + if (ier /= 0) stop 'Error allocating oknot,.. arrays' + oknot(:) = 0.d0 + aknot(:) = 0.d0 + level(:) = 0 + do i = 1,nknotA2(1) read(unit1,*) oknot(i),aknot(i),level(i) enddo @@ -134,12 +164,18 @@ subroutine model_berkeley_broadcast() endif read(unit1,*) nknotA2(2) - if (allocated(oknot2).or.allocated(aknot2).or.allocated(level2)) then - print *,'A3d already initiated' + + if (allocated(oknot2) .or. allocated(aknot2) .or. allocated(level2)) then + print *,'[model_berkeley_broadcast] A3d hnots2 already initiated' return endif - allocate(oknot2(nknotA2(2)),aknot2(nknotA2(2)),level2(nknotA2(2))) + allocate(oknot2(nknotA2(2)),aknot2(nknotA2(2)),level2(nknotA2(2)),stat=ier) + if (ier /= 0) stop 'Error allocating oknot2,.. arrays' + oknot2(:) = 0.d0 + aknot2(:) = 0.d0 + level2(:) = 0 + do i = 1,nknotA2(2) read(unit1,*) oknot2(i),aknot2(i),level2(i) enddo @@ -159,54 +195,77 @@ subroutine model_berkeley_broadcast() NBPARAM = npar if (npar > MAXPAR) stop 'npar greater than MAXPAR' + if (npar <= 0) stop 'npar must be > 0' + + allocate(parblock(npar),stat=ier) + if (ier /= 0) stop 'Error allocating parblock array' + parblock(:) = '' - allocate(parblock(npar)) do i = 1,npar read(unit2,*) dum2,nknotA1(i),parblock(i) - if (i > 1.and.nknotA1(i) /= nknotA1(1)) then + + ! same number of splines for all parameters + if (i > 1 .and. nknotA1(i) /= nknotA1(1)) then stop 'Inconsistent A1 splines between parameters' endif + if (i > 2) then nknotA2(i) = dum2 - if (nknotA2(i) /= nknotA2(2)) stop 'Param 3 and 4 need the same A2 splines than param 2' + if (nknotA2(i) /= nknotA2(2)) stop 'Param 3 and 4 need the same A2 splines than param 2' else if (dum2 /= nknotA2(i)) then stop 'Inconsistent hknots.dat and A3d.dat' endif - if (i == 2.and.nknotA2(i) /= nknotA2(1)) then + if (i == 2 .and. nknotA2(i) /= nknotA2(1)) then unconformal = .true. - if (.not.hknots2_exist) stop 'unconformal grid requires hknots2.dat' + if (.not. hknots2_exist) stop 'unconformal grid requires hknots2.dat' endif enddo read(unit2,*) ndisc + surface = 0 if (ndisc > 0) then surface = ndisc - if (unconformal) print *,'discontinuities assumed same grid as first par' + if (unconformal) print *,'discontinuities assumed same grid as first par' do i = 1,surface read(unit2,*) dum2, trash enddo endif - allocate(kntrad(nknotA1(1))) + allocate(kntrad(nknotA1(1)), & + kntrad_hh(nknotA1(1)-1),stat=ier) + if (ier /= 0) stop 'Error allocating kntrad,.. arrays' + kntrad(:) = 0.e0 + kntrad_hh(:) = 0.e0 + read(unit2,*) (kntrad(i),i=1,nknotA1(1)) + + ! takes spacings between spline radii + ! (used for spline evaluations fspl(..) in spl_A3d.c, see routine fill_hh_A3d()) + do i = 1,nknotA1(1)-1 + kntrad_hh(i) = kntrad(i+1) - kntrad(i) + enddo + mdim = 0 do i = 1,npar - mdim = mdim+nknotA1(i)*nknotA2(i) + mdim = mdim + nknotA1(i) * nknotA2(i) enddo - mdim = mdim+ndisc*nknotA2(1) + mdim = mdim + ndisc * nknotA2(1) + + allocate(mdl(mdim),stat=ier) + if (ier /= 0) stop 'Error allocating mdl array' + mdl(:) = 0.d0 - allocate(mdl(mdim)) n = 0 do i = 1,npar do j = 1,nknotA1(i) read(unit2,*) (mdl(k+n),k=1,nknotA2(i)) - n = n+nknotA2(i) + n = n + nknotA2(i) enddo enddo do i = 1,ndisc read(unit2,*) (mdl(k+n),k=1,nknotA2(1)) - n = n+nknotA2(1) + n = n + nknotA2(1) enddo if (n /= mdim) stop 'init_A3d dimension error' @@ -214,56 +273,107 @@ subroutine model_berkeley_broadcast() endif + ! user info + if (myrank == 0) then + write(IMAIN,*) ' number of parameters: NBPARAM = ',NBPARAM + call flush_IMAIN() + endif + ! ! broadcast ! ! - call BCAST_ALL_SINGLEI(nknotA2(1)) + call bcast_all_singlei(nknotA2(1)) - if (.not.allocated(oknot)) allocate(oknot(nknotA2(1))) - if (.not.allocated(aknot)) allocate(aknot(nknotA2(1))) - if (.not.allocated(level)) allocate(level(nknotA2(1))) + if (.not. allocated(oknot)) allocate(oknot(nknotA2(1))) + if (.not. allocated(aknot)) allocate(aknot(nknotA2(1))) + if (.not. allocated(level)) allocate(level(nknotA2(1))) - call BCAST_ALL_I(level,nknotA2(1)) - call BCAST_ALL_DP(oknot,nknotA2(1)) - call BCAST_ALL_DP(aknot,nknotA2(1)) + call bcast_all_i(level,nknotA2(1)) + call bcast_all_dp(oknot,nknotA2(1)) + call bcast_all_dp(aknot,nknotA2(1)) - call BCAST_ALL_SINGLEL(hknots2_exist) + call bcast_all_singlel(hknots2_exist) if (hknots2_exist) then - call BCAST_ALL_SINGLEI(nknotA2(2)) + call bcast_all_singlei(nknotA2(2)) - if (.not.allocated(oknot2)) allocate(oknot2(nknotA2(2))) - if (.not.allocated(aknot2)) allocate(aknot2(nknotA2(2))) - if (.not.allocated(level2)) allocate(level2(nknotA2(2))) + if (.not. allocated(oknot2)) allocate(oknot2(nknotA2(2))) + if (.not. allocated(aknot2)) allocate(aknot2(nknotA2(2))) + if (.not. allocated(level2)) allocate(level2(nknotA2(2))) - call BCAST_ALL_I(level2,nknotA2(2)) - call BCAST_ALL_DP(oknot2,nknotA2(2)) - call BCAST_ALL_DP(aknot2,nknotA2(2)) + call bcast_all_i(level2,nknotA2(2)) + call bcast_all_dp(oknot2,nknotA2(2)) + call bcast_all_dp(aknot2,nknotA2(2)) else nknotA2(2) = nknotA2(1) endif - call BCAST_ALL_SINGLEI(npar) + call bcast_all_singlei(npar) NBPARAM = npar - if (.not.allocated(parblock)) allocate(parblock(npar)) - - call BCAST_ALL_CH_ARRAY(parblock,npar,1) - call BCAST_ALL_I(nknotA1,MAXPAR) - call BCAST_ALL_I(nknotA2,MAXPAR) - call BCAST_ALL_SINGLEL(unconformal) - call BCAST_ALL_SINGLEI(ndisc) + if (.not. allocated(parblock)) allocate(parblock(npar)) + + call bcast_all_ch_array(parblock,npar,1) + call bcast_all_i(nknotA1,MAXPAR) + call bcast_all_i(nknotA2,MAXPAR) + call bcast_all_singlel(unconformal) + call bcast_all_singlei(ndisc) + + ! spline radii + if (.not. allocated(kntrad)) allocate(kntrad(nknotA1(1))) + if (.not. allocated(kntrad_hh)) allocate(kntrad_hh(nknotA1(1)-1)) + + call bcast_all_r(kntrad,nknotA1(1)) + call bcast_all_r(kntrad_hh,nknotA1(1)-1) + + call bcast_all_singlei(mdim) + + if (.not. allocated(mdl)) allocate(mdl(mdim)) + + call bcast_all_dp(mdl,mdim) + + ! allocates temporary work arrays for model_berkeley_shsv() routine + ! to avoid re-allocations of the same arrays for each call + size_work = maxval(nknotA2(:)) + allocate(work_dh(size_work),work_kindex(size_work),stat=ier) + if (ier /= 0) stop 'Error allocating work arrays for Berkeley model' + work_dh(:) = 0.d0 + work_kindex(:) = 0 + + ! great-circle distance evaluations between GLL point position theta/phi and knot position + ! pre-compute knot position coefficients for getdel_pre() + size_knots = nknotA2(1) + allocate(knot_coeff(3,size_knots),stat=ier) + if (ier /= 0) stop 'Error allocating knot_coeff array for Berkeley model' + knot_coeff(:,:) = 0.d0 + do i = 1,size_knots + ! knot coefficients colat/lon in rad + theta = (90.d0 - aknot(i)) * deg2rad + phi = oknot(i) * deg2rad + knot_coeff(1,i) = sin(theta) + knot_coeff(2,i) = cos(theta) + knot_coeff(3,i) = phi + enddo + ! free unneeded arrays + deallocate(aknot,oknot) - if (.not.allocated(kntrad)) allocate(kntrad(nknotA1(1))) - - call BCAST_ALL_DP(kntrad,nknotA1(1)) - - call BCAST_ALL_SINGLEI(mdim) - - if (.not.allocated(mdl)) allocate(mdl(mdim)) - - call BCAST_ALL_DP(mdl,mdim) + if (hknots2_exist) then + size_knots = nknotA2(2) + allocate(knot_coeff2(3,size_knots),stat=ier) + if (ier /= 0) stop 'Error allocating knot_coeff2 array for Berkeley model' + knot_coeff2(:,:) = 0.d0 + do i = 1,size_knots + ! knot position colat/lon in rad + theta = (90.d0 - aknot2(i)) * deg2rad + phi = oknot2(i) * deg2rad + knot_coeff2(1,i) = sin(theta) + knot_coeff2(2,i) = cos(theta) + knot_coeff2(3,i) = phi + enddo + ! free unneeded arrays + deallocate(aknot2,oknot2) + endif end subroutine model_berkeley_broadcast @@ -272,7 +382,8 @@ end subroutine model_berkeley_broadcast !-------------------------------------------------------------------------------------------------- ! - subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,iregion_code,CRUSTAL) + subroutine model_berkeley_shsv(r,theta,phi,vpv_ref,vph_ref,vsv_ref,vsh_ref,rho_ref, & + dvsh,dvsv,dvph,dvpv,drho,eta_aniso,iregion_code,CRUSTAL) ! returns isotropic vs, vp, and rho assuming scaling dlnVs/dlnVp=2 dlnVs/dlnrho=3 ! also returns anisotropic parameters xi,fi,eta,Gc,Gs,Hc,Hs,Bc,Bs if ifanis=1 @@ -283,7 +394,9 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir implicit none double precision, intent(in) :: r,theta,phi - double precision, intent(out) :: dvsv,dvsh,dvpv,dvph,drho,eta_aniso + double precision, intent(in) :: vpv_ref,vph_ref,vsv_ref,vsh_ref,rho_ref + double precision, intent(out) :: dvsv,dvsh,dvpv,dvph,drho + double precision, intent(inout) :: eta_aniso integer, intent(in) :: iregion_code logical, intent(in) :: CRUSTAL @@ -291,30 +404,37 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir ! local parameters double precision :: x,rho1d,vpv1d,vph1d,vsv1d,vsh1d,eta1d,Qmu1d,Qkappa1d - double precision :: vs,vp,rho,Qs + double precision :: vs,vp,rho ! Qs double precision :: xi,fi,eta,Gc,Gs + double precision :: fi_inv integer :: jump,effnknot,i,j,k - integer, dimension(:), allocatable :: kindex + integer :: nknots_horiz,nknots_radial + !integer, dimension(:), allocatable :: kindex - double precision :: lat,lon double precision :: del,dr,dv double precision :: AA,CC,FF,LL,NN double precision :: eta1,adel1,r_ + ! knot great-circle distances for different spline levels + !double precision, dimension(8), parameter :: adel = (/63.4, 31.7, 15.8, 7.9, 3.95, 1.98, 0.99/) ! Include level 8 (FM, April 2021 - Mod. courtesy of Dan Frost) double precision, dimension(8), parameter :: adel = (/63.4, 31.7, 15.8, 7.9, 3.95, 1.98, 0.99, 0.495/) - !adel=(/63.4, 31.7, 15.8, 7.9, 3.95, 1.98, 0.99/) - double precision, dimension(:), allocatable :: dh - - double precision,external :: getdel - double precision,external :: spbsp + !double precision, dimension(:), allocatable :: dh + !double precision :: lat,lon ! unused + double precision :: sin_theta0,cos_theta0 double precision :: vsv,vsh,vpv,vph,scaleval double precision :: aa1, bb1 - double precision, parameter :: rad2deg = 180.d0/PI + double precision,external :: getdel_pre + double precision,external :: spbsp + + !double precision, parameter :: rad2deg = 180.d0/PI + + ! tolerance for water density check + double precision, parameter :: TOL_RHO_WATER = 1200.d0 / EARTH_RHOAV ! non-dimensionalized ! initializes model perturbations dvsv = 0.d0 @@ -322,7 +442,6 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir dvpv = 0.d0 dvph = 0.d0 drho = 0.d0 - eta_aniso = 1.d0 xi = 1.d0 fi = 1.d0 @@ -339,25 +458,54 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir ! note: r is non-dimensionalized/normalized input radius between [0,1] if (r > moho1D_radius) then r_ = moho1D_radius ! * EARTH_R_KM + ! re-evaluate reference 1D model values for this updated radius + call model_1dberkeley(r_,rho1d,vpv1d,vph1d,vsv1d,vsh1d,eta1d,Qkappa1d,Qmu1d,iregion_code,CRUSTAL) else r_ = r ! * EARTH_R_KM + ! reference model values already setup prior to this routine call + vpv1d = vpv_ref + vph1d = vph_ref + vsv1d = vsv_ref + vsh1d = vsh_ref + rho1d = rho_ref + eta1d = eta_aniso + Qkappa1d = 9999.d0 ! unused + Qmu1d = 9999.d0 ! unused endif ! non-dimensionalized radius x = r_ ! / EARTH_R_KM - call model_1dberkeley(x,rho1d,vpv1d,vph1d,vsv1d,vsh1d,eta1d,Qkappa1d,Qmu1d,iregion_code,CRUSTAL) - - if (rho1d < 1200.d0) then - ! No water in RegSEM please - x = r_ * EARTH_R_KM / 6367.999d0 - call model_1dberkeley(x,rho1d,vpv1d,vph1d,vsv1d,vsh1d,eta1d,Qkappa1d,Qmu1d,iregion_code,CRUSTAL) + if (USE_OLD_VERSION_FORMAT) then + ! old version by mistake compares non-dimensionalized rho1d with density value 1200.d0 + ! thus, this will always be evaluated since rho1d is much smaller than that. + if (rho1d < 1200.d0) then + ! No water in RegSEM please + ! + !debug + !print *,'debug: rho1d ',rho1d, rho1d * EARTH_RHOAV, rho_ref * EARTH_RHOAV,'radius ',x * EARTH_R_KM + + ! originally, this likely wanted to re-scale the radius to be below ocean at 6368.0 km? + ! however, this formula won't work - x can still be above ocean radius. + ! this scaling was moving up the radius value, thus shifting all mantle velocities from the 1D definition slightly up. + x = r_ * EARTH_R_KM / 6367.999d0 + call model_1dberkeley(x,rho1d,vpv1d,vph1d,vsv1d,vsh1d,eta1d,Qkappa1d,Qmu1d,iregion_code,CRUSTAL) + endif + else + ! fix + ! check non-dimensionalized density + if (rho1d < TOL_RHO_WATER) then + ! No water in RegSEM please + ! takes radius slightly below ocean at 6368.0 km to get mantle/crust values from the layer below ocean + x = 6367.999d0 / EARTH_R_KM + call model_1dberkeley(x,rho1d,vpv1d,vph1d,vsv1d,vsh1d,eta1d,Qkappa1d,Qmu1d,iregion_code,CRUSTAL) + endif endif ! below use r_ as radius in km r_ = r_ * EARTH_R_KM - ! non-dimensionalizes 1d model values + ! re-dimensionalizes 1d model values scaleval = EARTH_R * dsqrt(PI*GRAV*EARTH_RHOAV) rho1d = rho1d * EARTH_RHOAV @@ -372,8 +520,9 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir LL = vsv1d*vsv1d * rho NN = vsh1d*vsh1d * rho FF = eta1d * (AA - 2.d0 * LL) - qs = qmu1d + ! unused + !Qs = Qmu1d !call get_1Dmodel_TI(AA,CC,FF,LL,NN,rho,Qs,r_*1000.d0) !if (rho < 1200.d0) call get_1Dmodel_TI(AA,CC,FF,LL,NN,rho,Qs,6367999.d0) ! No water in RegSEM please @@ -384,36 +533,53 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir vs = sqrt((CC + (1.d0 - 2.d0 * eta1) * AA + (6.d0 + 4.d0 * eta1) * LL + 5.d0 * NN) / (15.d0 * rho)) eta = eta1 - xi = NN/LL - fi = CC/AA + xi = NN / LL + fi = CC / AA + + ! point lat/lon in degrees for getdel() + !lat = 90.d0 - rad2deg * theta + !lon = rad2deg * phi + + ! coefficients for getdel_pre() + sin_theta0 = sin(theta) + cos_theta0 = cos(theta) ! Vs perturbation if (r_ > kntrad(nknotA1(1)) .or. r_ < kntrad(1)) then dv = 0.d0 else jump = 0 - allocate(dh(nknotA2(1)),kindex(nknotA2(1))) - lat = 90.d0 - rad2deg * theta - lon = rad2deg * phi + + !allocate(dh(nknotA2(1)),kindex(nknotA2(1))) + work_dh(:) = 0.d0 + work_kindex(:) = 0 effnknot = 0 - do i = 1,nknotA2(1) - del = getdel(lat,lon,aknot(i),oknot(i)) + + nknots_horiz = nknotA2(1) + nknots_radial = nknotA1(1) + + do i = 1,nknots_horiz + ! gets great-circle distance between position theta/phi and knot + !old: del = getdel(lat,lon,aknot(i),oknot(i)) + del = getdel_pre(sin_theta0,cos_theta0,phi,knot_coeff(1,i),knot_coeff(2,i),knot_coeff(3,i)) if (del <= adel(level(i)) * 2.d0) then effnknot = effnknot+1 - kindex(effnknot) = i - dh(effnknot) = spbsp(del,adel(level(i))) + work_kindex(effnknot) = i + work_dh(effnknot) = spbsp(del,adel(level(i))) endif enddo dv = 0.d0 - do i = 1,nknotA1(1) - call fspl(i,nknotA1(1),kntrad,r_,dr) + do i = 1,nknots_radial + ! spline value + call fspl(i,nknots_radial,kntrad,kntrad_hh,r_,dr) + do j = 1,effnknot - dv = dv + dr * dh(j) * mdl(jump+kindex(j) + nknotA2(1) * (i-1)) + dv = dv + dr * work_dh(j) * mdl(jump + work_kindex(j) + nknots_horiz * (i-1)) enddo enddo - deallocate(dh,kindex) + !deallocate(dh,kindex) endif ! Scaling @@ -422,41 +588,109 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir rho = rho + 0.33333d0 * dv * rho ! Perturbation of other parameters - if (npar < 2) then ! no other parameters - dv = 0.d0 - else + ! note: default in A3d.dat is npar == 2 + if (npar == 2) then + ! xi perturbation (dv) + if (r_ > kntrad(nknotA1(2)) .or. r_ < kntrad(1)) then + dv = 0.d0 + else + jump = jump + nknotA1(1) * nknotA2(1) + + work_dh(:) = 0.d0 + work_kindex(:) = 0 + effnknot = 0 + + nknots_horiz = nknotA2(2) + nknots_radial = nknotA1(2) + + if (unconformal) then + ! unconformal grid + do i = 1,nknots_horiz + ! gets great-circle distance between position theta/phi and knot + !old: del = getdel(lat,lon,aknot2(i),oknot2(i)) + del = getdel_pre(sin_theta0,cos_theta0,phi,knot_coeff2(1,i),knot_coeff2(2,i),knot_coeff2(3,i)) + + adel1 = adel(level2(i)) + if (del <= adel1 * 2.d0) then + effnknot = effnknot + 1 + work_kindex(effnknot) = i + work_dh(effnknot) = spbsp(del,adel1) + endif + enddo + else + ! conformal grid + do i = 1,nknots_horiz + ! gets great-circle distance between position theta/phi and knot + !old: del = getdel(lat,lon,aknot(i),oknot(i)) + del = getdel_pre(sin_theta0,cos_theta0,phi,knot_coeff(1,i),knot_coeff(2,i),knot_coeff(3,i)) + + adel1 = adel(level(i)) + if (del <= adel1 * 2.d0) then + effnknot = effnknot + 1 + work_kindex(effnknot) = i + work_dh(effnknot) = spbsp(del,adel1) + endif + enddo + endif + + dv = 0.d0 + do i = 1,nknots_radial + ! spline value + call fspl(i,nknots_radial,kntrad,kntrad_hh,r_,dr) + + do j = 1,effnknot + dv = dv + dr * work_dh(j) * mdl(jump + work_kindex(j) + nknots_horiz * (i-1)) + enddo + enddo + endif + + ! npar == k == 2 + xi = xi + dv * xi + fi = fi - 1.5d0 * dv * fi + eta = eta - 2.5d0 * dv * eta + else if (npar > 2) then + ! npar > 2 do k = 2,npar if (r_ > kntrad(nknotA1(k)) .or. r_ < kntrad(1)) then dv = 0.d0 else jump = jump + nknotA1(k-1)*nknotA2(k-1) - allocate(dh(nknotA2(k)),kindex(nknotA2(k))) - lat = 90.d0 - rad2deg * theta - lon = rad2deg * phi + + !allocate(dh(nknotA2(k)),kindex(nknotA2(k))) + work_dh(:) = 0.d0 + work_kindex(:) = 0 + effnknot = 0 do i = 1,nknotA2(k) if (unconformal) then - del = getdel(lat,lon,aknot2(i),oknot2(i)) + ! gets great-circle distance between position theta/phi and knot + !old: del = getdel(lat,lon,aknot2(i),oknot2(i)) + del = getdel_pre(sin_theta0,cos_theta0,phi,knot_coeff2(1,i),knot_coeff2(2,i),knot_coeff2(3,i)) adel1 = adel(level2(i)) else - del = getdel(lat,lon,aknot(i),oknot(i)) + ! gets great-circle distance between position theta/phi and knot + !old: del = getdel(lat,lon,aknot(i),oknot(i)) + del = getdel_pre(sin_theta0,cos_theta0,phi,knot_coeff(1,i),knot_coeff(2,i),knot_coeff(3,i)) adel1 = adel(level(i)) endif + if (del <= adel1 * 2.d0) then effnknot = effnknot+1 - kindex(effnknot) = i - dh(effnknot) = spbsp(del,adel1) + work_kindex(effnknot) = i + work_dh(effnknot) = spbsp(del,adel1) endif enddo dv = 0.d0 do i = 1,nknotA1(k) - call fspl(i,nknotA1(k),kntrad,r_,dr) + ! spline value + call fspl(i,nknotA1(k),kntrad,kntrad_hh,r_,dr) + do j = 1,effnknot - dv = dv + dr * dh(j) * mdl(jump + kindex(j) + nknotA2(k) * (i-1)) + dv = dv + dr * work_dh(j) * mdl(jump + work_kindex(j) + nknotA2(k) * (i-1)) enddo enddo - deallocate(dh,kindex) + !deallocate(dh,kindex) endif if (k == 2) then @@ -470,6 +704,9 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir endif ! Here we can add a scaling to get Hc, Hs, Bc and Bs enddo + else + ! no other parameters + dv = 0.d0 endif ! ============= commented by < FM> on Feb 3, 2020 ===== @@ -482,22 +719,30 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir ! ==================================================== ! New conversion relationships < FM> - Feb 3, 2020 ! Auxiliar values - aa1 = 3.d0 + ( 8.d0 + 4.d0 * eta ) / fi - bb1 = 1.d0 + ( 1.d0 - 2.d0 * eta ) / fi + fi_inv = 1.d0 / fi + + aa1 = 3.d0 + ( 8.d0 + 4.d0 * eta ) * fi_inv + bb1 = 1.d0 + ( 1.d0 - 2.d0 * eta ) * fi_inv vsv = sqrt( 15.d0 * (vp*vp * bb1 - vs*vs * aa1) / & (8.d0 * (1.d0-eta) * bb1 - (6.d0 + 4.d0 * eta + 5.d0 * xi) * aa1 ) ) vsh = sqrt( xi ) * vsv vpv = sqrt( (15.d0 * vp*vp - 8.d0 * (1.d0-eta) * vsv*vsv ) / & - (3.d0 + ( 8.d0 + 4.d0 * eta ) / fi ) ) - vph = vpv * sqrt(1.d0/fi) + (3.d0 + ( 8.d0 + 4.d0 * eta ) * fi_inv ) ) + vph = vpv * sqrt( fi_inv ) rho = rho ! =================================================== ! perturbations - dvsv = vsv / vsv1d - 1.d0 - if (vsv1d == 0) dvsv = 0.d0 - dvsh = vsh / vsh1d - 1.d0 - if (vsh1d == 0) dvsh = 0.d0 + if (vsv1d == 0.d0) then + dvsv = 0.d0 + else + dvsv = vsv / vsv1d - 1.d0 + endif + if (vsh1d == 0.d0) then + dvsh = 0.d0 + else + dvsh = vsh / vsh1d - 1.d0 + endif dvpv = vpv / vpv1d - 1.d0 dvph = vph / vph1d - 1.d0 @@ -507,33 +752,64 @@ subroutine model_berkeley_shsv(r,theta,phi,dvsh,dvsv,dvph,dvpv,drho,eta_aniso,ir end subroutine model_berkeley_shsv +! +!-------------------------------------------------------------------------------------------------- +! + + double precision function getdel_pre(sin_theta0,cos_theta0,phi0,sin_theta,cos_theta,phi) + + use constants, only: PI + + implicit none + double precision,intent(in) :: sin_theta0,cos_theta0,phi0,sin_theta,cos_theta,phi + ! local parameters + double precision :: a + double precision, parameter :: rad2deg = 180.d0 / PI ! 180.d0 / (4.d0 * datan(1.d0)) + + ! great-circle distance (by law of cosines formula) in degrees + a = cos_theta * cos_theta0 + sin_theta * sin_theta0 * cos(phi - phi0) + + ! limit to +/- 1 numerical accuracy for acos(..) + if (abs(a) > 1.d0) a = sign(1.d0,a) * 1.d0 + + ! great-circle distance in degrees + getdel_pre = rad2deg * acos(a) + + end function getdel_pre + ! !-------------------------------------------------------------------------------------------------- ! double precision function getdel(a0,o0,a,o) + use constants, only: PI + implicit none double precision,intent(in) :: a0,o0,a,o ! local parameters - double precision :: q0,sq0,cq0,q,sq,cq,ff,sff,cff,arg + double precision :: q0,sq0,cq0 + double precision :: q,sq,cq + double precision :: ff,cff,arg ! sff - double precision, parameter :: deg2rad = (4.d0 * datan(1.d0)) / 180.d0 - double precision, parameter :: rad2deg = 180.d0 / (4.d0 * datan(1.d0)) + double precision, parameter :: deg2rad = PI / 180.d0 ! (4.d0 * datan(1.d0)) / 180.d0 + double precision, parameter :: rad2deg = 180.d0 / PI ! 180.d0 / (4.d0 * datan(1.d0)) - q0 = (90.d0-a0) * deg2rad + ! great-circle distance (by law of cosines formula) in degrees + q0 = (90.d0 - a0) * deg2rad sq0 = sin(q0) cq0 = cos(q0) - q = (90.d0-a) * deg2rad + q = (90.d0 - a) * deg2rad sq = sin(q) cq = cos(q) - ff = (o-o0) * deg2rad - sff = sin(ff) + ff = (o - o0) * deg2rad + !sff = sin(ff) ! not used cff = cos(ff) arg = cq * cq0 + sq * sq0 * cff + if (arg > 1.d0) arg = 1.d0 if (arg < -1.d0) arg = -1.d0 @@ -549,12 +825,22 @@ double precision function spbsp(hdel,ahdel) implicit none double precision,intent(in) :: hdel,ahdel + ! local parameters + double precision :: ratio,two_minus_ratio + + ! factor + ratio = hdel / ahdel if (hdel < ahdel) then - spbsp = 0.75d0 * (hdel/ahdel) * (hdel/ahdel) * (hdel/ahdel) - 1.5d0 * (hdel/ahdel) * (hdel/ahdel) + 1.d0 + spbsp = (0.75d0 * ratio - 1.5d0) * ratio * ratio + 1.d0 else if (hdel <= ahdel * 2.d0) then - spbsp = 0.25d0 * (2.d0 - hdel/ahdel) * (2.d0 - hdel/ahdel) * (2.d0 - hdel/ahdel) - if (spbsp < 0.d0) spbsp = 0.d0 + two_minus_ratio = 2.d0 - ratio + if (two_minus_ratio < 0.d0) then + spbsp = 0.d0 + else + spbsp = 0.25d0 * two_minus_ratio * two_minus_ratio * two_minus_ratio + endif + !if (spbsp < 0.d0) spbsp = 0.d0 ! spbsp is negative if two_minus_ratio is negative else spbsp = 0.d0 endif diff --git a/src/meshfem3D/model_crust_berkeley.f90 b/src/meshfem3D/model_crust_berkeley.f90 index 095bf13fd..8f545c50e 100644 --- a/src/meshfem3D/model_crust_berkeley.f90 +++ b/src/meshfem3D/model_crust_berkeley.f90 @@ -35,13 +35,20 @@ module model_crust_berkeley_par implicit none + ! moho model + real, dimension(:,:), allocatable :: moho_start integer :: NBT,NBP real :: drin + + ! moho filtering real, parameter :: drfiltre = 2.e0 - !double precision, parameter :: dr_ = 2.d0 - real, dimension(:,:), allocatable :: crust_array - real, dimension(:,:), allocatable :: moho_start + ! crustal model + double precision, dimension(:,:), allocatable :: crust_array + integer, dimension(:), allocatable :: crust_array_ind + ! number of crustal entries + integer :: num_crust_array = 0 + ! moho depth from 1D reference model (in km) double precision :: moho1D_depth = -1.0d0 @@ -57,7 +64,7 @@ subroutine model_berkeley_crust_broadcast() ! standard routine to setup model use model_crust_berkeley_par - use constants, only: A3d_folder,myrank,IMAIN + use constants, only: A3d_folder,myrank,IMAIN,IIN implicit none @@ -75,24 +82,24 @@ subroutine model_berkeley_crust_broadcast() endif ! moho data - open(52,file=trim(file_moho),status='old',action='read',iostat=ier) + open(IIN,file=trim(file_moho),status='old',action='read',iostat=ier) if (ier /= 0) then print *,'Error opening file: ',trim(file_moho) stop 'Error opening file crust2moho_2x2.dat' endif - call read_crustmoho_filtre(52) - close(52) + call read_crustmoho_filtre(IIN) + close(IIN) !drfiltre = dr_ ! crustal data - open(52,file=trim(file_crust),status='old',action='read',iostat=ier) + open(IIN,file=trim(file_crust),status='old',action='read',iostat=ier) if (ier /= 0) then print *,'Error opening file: ',trim(file_crust) stop 'Error opening file crust2cru2av_2x2.dat' endif - call read_crust_smooth_variable(52) - close(52) + call read_crust_smooth_variable(IIN) + close(IIN) ! gets moho depth from 1D model if (myrank == 0) then @@ -109,27 +116,32 @@ end subroutine model_berkeley_crust_broadcast !-------------------------------------------------------------------------------------------------- ! - subroutine model_berkeley_crust(x,theta,phi,vp,vs,rho,moho,found_crust) + subroutine model_berkeley_crust(x,theta,phi,vp,vs,rho,moho,found_crust,elem_in_crust,moho_only) ! returns isotropic crustal velocities use model_crust_berkeley_par - use constants, only: PI,GRAV,EARTH_RHOAV,EARTH_R,EARTH_R_KM + use constants, only: PI,GRAV,EARTH_RHOAV,EARTH_R,EARTH_R_KM, & + USE_OLD_VERSION_FORMAT implicit none double precision,intent(in) :: x,theta,phi double precision,intent(out) :: vp,vs,rho,moho logical,intent(out) :: found_crust + logical,intent(in) :: elem_in_crust,moho_only ! local parameters - double precision :: vsv,vsh - double precision :: depth - double precision :: moho_depth, scaleval + double precision :: vsv,vsh,depth,moho_depth + double precision :: scaleval,scaleval_inv !double precision, parameter :: deg2rad = PI/180.d0, rad2deg = 180.d0/PI ! initializes - found_crust = .true. + vp = 0.d0 + vs = 0.d0 + rho = 0.d0 + moho = 0.d0 + found_crust = .false. !theta = (90.d0 - lat) * deg2rad ! assumed lat range: [-90,90] !phi = lon * deg2rad ! assumed lon range: [-180,180] @@ -138,9 +150,22 @@ subroutine model_berkeley_crust(x,theta,phi,vp,vs,rho,moho,found_crust) call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth) - !TODO: daniel - check flag usage - ! sets flag if point within crust - !if (depth > moho_depth) found_crust = .false. + ! using crustal values + if (USE_OLD_VERSION_FORMAT) then + ! always uses the crustal values + found_crust = .true. + else + ! only uses the crustal values if point above moho + ! sets flag if point within crust (depth shallower than moho_depth) + if (depth <= moho_depth .or. elem_in_crust) found_crust = .true. + endif + + ! only return moho + if (moho_only) then + moho = moho_depth / EARTH_R_KM + ! all done + return + endif ! ! get equivalent isotropic vs @@ -151,9 +176,11 @@ subroutine model_berkeley_crust(x,theta,phi,vp,vs,rho,moho,found_crust) ! scale values for specfem ! scaleval = dsqrt(PI*GRAV*EARTH_RHOAV) + scaleval_inv = 1.d0 / (EARTH_R * scaleval) + + vp = vp * scaleval_inv + vs = vs * scaleval_inv - vp = vp / (EARTH_R * scaleval) - vs = vs / (EARTH_R * scaleval) rho = rho / EARTH_RHOAV moho = moho_depth / EARTH_R_KM @@ -164,28 +191,35 @@ end subroutine model_berkeley_crust !-------------------------------------------------------------------------------------------------- ! - subroutine model_berkeley_crust_aniso(x,theta,phi,vpv,vph,vsv,vsh,eta_aniso,rho,moho,found_crust) + subroutine model_berkeley_crust_aniso(x,theta,phi,vpv,vph,vsv,vsh,eta_aniso,rho,moho,found_crust,elem_in_crust,moho_only) ! returns anisotropic crustal velocities use model_crust_berkeley_par - use constants, only: PI,GRAV,EARTH_RHOAV,EARTH_R,EARTH_R_KM + use constants, only: PI,GRAV,EARTH_RHOAV,EARTH_R,EARTH_R_KM, & + USE_OLD_VERSION_FORMAT implicit none double precision,intent(in) :: x,theta,phi double precision,intent(out) :: vpv,vph,vsv,vsh,eta_aniso,rho,moho logical,intent(out) :: found_crust + logical,intent(in) :: elem_in_crust,moho_only ! local parameters - double precision :: vp - double precision :: depth - double precision :: moho_depth, scaleval + double precision :: vp,depth,moho_depth + double precision :: scaleval,scaleval_inv !double precision, parameter :: deg2rad = PI/180.d0, rad2deg = 180.d0/PI ! initializes + vpv = 0.d0 + vph = 0.d0 + vsv = 0.d0 + vsh = 0.d0 eta_aniso = 1.d0 - found_crust = .true. + rho = 0.d0 + moho = 0.d0 + found_crust = .false. !theta = (90.d0 - lat) * deg2rad ! assumed lat range: [-90,90] !phi = lon * deg2rad ! assumed lon range: [-180,180] @@ -194,19 +228,34 @@ subroutine model_berkeley_crust_aniso(x,theta,phi,vpv,vph,vsv,vsh,eta_aniso,rho, call get_crust_val_csem(theta,phi,depth,rho,vp,vsv,vsh,moho_depth) - !TODO: daniel - check flag usage - ! sets flag if point within crust - !if (depth > moho_depth) found_crust = .false. + ! using crustal values + if (USE_OLD_VERSION_FORMAT) then + ! always uses the crustal values + found_crust = .true. + else + ! only uses the crustal values if point above moho + ! sets flag if point within crust (depth shallower than moho_depth) + if (depth <= moho_depth .or. elem_in_crust) found_crust = .true. + endif + + ! only return moho + if (moho_only) then + moho = moho_depth / EARTH_R_KM + ! all done + return + endif ! ! scale values for specfem ! scaleval = dsqrt(PI*GRAV*EARTH_RHOAV) + scaleval_inv = 1.d0 / (EARTH_R * scaleval) + + vph = vp * scaleval_inv + vpv = vp * scaleval_inv + vsh = vsh * scaleval_inv + vsv = vsv * scaleval_inv - vph = vp / (EARTH_R * scaleval) - vpv = vp / (EARTH_R * scaleval) - vsh = vsh / (EARTH_R * scaleval) - vsv = vsv / (EARTH_R * scaleval) rho = rho / EARTH_RHOAV moho = moho_depth / EARTH_R_KM @@ -227,16 +276,22 @@ subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth) double precision,intent(out) :: rho,vp,vsv,vsh,moho_depth ! local parameters - double precision :: xi(5),wi(5) + ! 4-th order GLL positions + ! pre-computed by calling LGL_NODES(4,1.d-6,xi,wi)) + double precision, parameter :: xi_ref(5) = (/ -1.d0, -0.654653670707977d0, 0.d0, 0.654653670707977d0, 1.d0 /) + ! pre-computed position range in [0,1] + double precision, parameter :: xi_norm(5) = ( xi_ref(:) + 1.d0) * 0.5d0 + + double precision :: xi(5) !,wi(5) double precision :: rho_cr(5),vp_cr(5),vsv_cr(5),vsh_cr(5) double precision :: x integer :: i - real,external :: moho_filtre + double precision,external :: moho_filtre double precision,external :: lagrange ! get moho depth - moho_depth = moho1D_depth - dble(moho_filtre(theta,phi)) + moho_depth = moho1D_depth - moho_filtre(theta,phi) !debug !print *,"debug: [get_crust_val_csem] Moho depth:",moho_depth, & @@ -254,8 +309,13 @@ subroutine get_crust_val_csem(theta,phi,z,rho,vp,vsv,vsh,moho_depth) ! ! get GLL nodes position ! - call LGL_NODES(4,1.d-6,xi,wi) - xi(:) = moho_depth - (xi(:) + 1.d0)/2.d0 * moho_depth + ! 4th-order GLL node positions + !call LGL_NODES(4,1.d-6,xi,wi) + !xi(:) = moho_depth - (xi(:) + 1.d0) * 0.5d0 * moho_depth + + ! above not needed, will use pre-computed xi_norm array + xi(:) = moho_depth - xi_norm(:) * moho_depth + x = z if (x > maxval(xi)) x = maxval(xi) if (x < minval(xi)) x = minval(xi) @@ -296,13 +356,35 @@ double precision function lagrange(ind,xi,nxi,x) integer,intent(in) :: ind,nxi double precision,intent(in) :: xi(nxi),x + ! local parameters + double precision :: x0 integer :: i + ! initializes lagrange = 1.d0 - do i = 1,nxi - if (i /= ind) lagrange = lagrange * ( (x-xi(i))/(xi(ind)-xi(i)) ) - enddo + x0 = xi(ind) + + if (nxi == 5) then + ! NGLL == 5 case + select case(ind) + case (1) + lagrange = (x-xi(2))/(x0-xi(2)) * (x-xi(3))/(x0-xi(3)) * (x-xi(4))/(x0-xi(4)) * (x-xi(5))/(x0-xi(5)) + case (2) + lagrange = (x-xi(1))/(x0-xi(1)) * (x-xi(3))/(x0-xi(3)) * (x-xi(4))/(x0-xi(4)) * (x-xi(5))/(x0-xi(5)) + case (3) + lagrange = (x-xi(1))/(x0-xi(1)) * (x-xi(2))/(x0-xi(2)) * (x-xi(4))/(x0-xi(4)) * (x-xi(5))/(x0-xi(5)) + case (4) + lagrange = (x-xi(1))/(x0-xi(1)) * (x-xi(2))/(x0-xi(2)) * (x-xi(3))/(x0-xi(3)) * (x-xi(5))/(x0-xi(5)) + case (5) + lagrange = (x-xi(1))/(x0-xi(1)) * (x-xi(2))/(x0-xi(2)) * (x-xi(3))/(x0-xi(3)) * (x-xi(4))/(x0-xi(4)) + end select + else + ! general NGLL case + do i = 1,nxi + if (i /= ind) lagrange = lagrange * ( (x-xi(i))/(x0-xi(i)) ) + enddo + endif end function lagrange @@ -362,7 +444,7 @@ end subroutine LGL_NODES !-------------------------------------------------------------------------------------------------- ! - real function moho_filtre(theta,phi) + double precision function moho_filtre(theta,phi) !theta phi en radians !dr en degre @@ -370,16 +452,19 @@ real function moho_filtre(theta,phi) use model_crust_berkeley_par + use constants, only: RADIANS_TO_DEGREES implicit none double precision,intent(in) :: theta,phi ! local parameter real :: t,p - real, parameter :: pi = 3.141592653589793, deg2rad = pi/180.0, rad2deg = 180.0/pi + real, parameter :: pi = 3.141592653589793e0 ! pi = 4.0 * atan(1.0) + real, parameter :: deg2rad = pi/180.e0 - t = sngl(theta) / deg2rad - p = sngl(phi) / deg2rad - moho_filtre = gauss_filtre1(moho_start,t,p,drfiltre) + t = sngl(theta * RADIANS_TO_DEGREES) + p = sngl(phi * RADIANS_TO_DEGREES) + + moho_filtre = dble(gauss_filtre1(moho_start,t,p,drfiltre)) contains @@ -387,24 +472,34 @@ real function moho_filtre(theta,phi) real function gauss_filtre1(tin,theta,phi,dr) !----------------------------------------------------------------- implicit none - real,intent(in) :: theta,phi,dr - real, dimension(:,:),intent(in) :: tin + real, intent(in) :: theta,phi,dr + real, dimension(:,:), intent(in) :: tin ! local parameters - real :: thetar,phir,tmpnorm,inte - real :: tmp - integer :: i,ii,j,jj,LARG + real :: thetar,phir + real :: inte,rfac + real :: tmp,tmpnorm + integer :: i,ii,j,jj + + integer, parameter :: LARG = 10 + + tmp = 0.e0 + tmpnorm = 0.e0 - tmp = 0.0 - tmpnorm = 0.0 - LARG = 10 do i = 1,LARG+1 do j = 1,LARG+1 call get_indexloc(phi,theta,i,j,dr,LARG,ii,jj,phir,thetar) - inte = cos_cylindre(theta,phi,dr,thetar,phir) * (dr/real(LARG/2) * deg2rad)**2 * sin(thetar * deg2rad) + + ! integration factor + rfac = ( dr/real(LARG/2) * deg2rad )**2 + inte = cos_cylindre(theta,phi,dr,thetar,phir) * rfac * sin(thetar * deg2rad) + + ! contribution tmp = tmp + tin(ii,jj) * inte tmpnorm = tmpnorm + inte enddo enddo + + ! normalizes gauss_filtre1 = tmp / tmpnorm end function gauss_filtre1 @@ -413,30 +508,38 @@ end function gauss_filtre1 real function cos_cylindre(t0_,p0_,d0_,theta_,phi_) !---------------------------------------------------------------------- implicit none - real,intent(in) :: t0_,p0_ - real,intent(in) :: d0_,theta_,phi_ + real, intent(in) :: t0_,p0_ + real, intent(in) :: d0_,theta_,phi_ ! local parameters real :: t0,p0,d0,theta,phi,d_ang real :: cosa + ! target location t0 = t0_ * deg2rad p0 = p0_ * deg2rad + + ! raster location theta = theta_ * deg2rad phi = phi_ * deg2rad + d0 = d0_ * deg2rad - !distance angulaire au centre du cylindre: + + ! distance angulaire au centre du cylindre: cosa = cos(theta) * cos(t0) + sin(theta) * sin(t0) * cos(phi - p0) + if (cosa >= 1.0) then - d_ang = 0.0 + d_ang = 0.e0 else if (cosa <= -1.0) then - d_ang = 4.0 * atan(1.0) + d_ang = pi else - d_ang = acos(cos(theta) * cos(t0) + sin(theta) * sin(t0) * cos(phi - p0)) + !d_ang = acos(cos(theta) * cos(t0) + sin(theta) * sin(t0) * cos(phi - p0)) + d_ang = acos( cosa ) endif + if (d_ang > d0) then - cos_cylindre = 0.0 + cos_cylindre = 0.e0 else - cos_cylindre = 0.5 * (1.0 + cos(PI * d_ang/d0)) + cos_cylindre = 0.5e0 * (1.e0 + cos(pi * d_ang/d0)) endif end function cos_cylindre @@ -458,25 +561,25 @@ subroutine get_indexloc(phi,theta,i,j,dr,LARG,ii,jj,phir,thetar) integer, intent(out) :: ii,jj ! local parameters - double precision :: t,p - double precision, parameter :: eps = 1.d-8 + real :: t,p + real, parameter :: eps = 1.e-8 - p = dble(phi) + (i-1-LARG/2) * dble(dr)/dble(LARG/2) - t = dble(theta) + (j-1-LARG/2) * dble(dr)/dble(LARG/2) + p = phi + (i-1-LARG/2) * dr/real(LARG/2) + t = theta + (j-1-LARG/2) * dr/real(LARG/2) - if (p < 0.d0-eps) p = p + 360.d0 - if (p >= 360.d0-eps) p = p - 360.d0 + if (p < 0.e0 - eps) p = p + 360.e0 + if (p >= 360.e0 - eps) p = p - 360.e0 - if (t > 180.d0 - eps) then - t = t - 180.d0 - p = 360.d0 - p - else if (t < 0.d0 - eps) then - t = 180.d0 + t - p = 360.d0 - p + if (t > 180.e0 - eps) then + t = t - 180.e0 + p = 360.e0 - p + else if (t < 0.e0 - eps) then + t = 180.e0 + t + p = 360.e0 - p endif - if (p < 0.d0-eps) p = p + 360.d0 - if (p >= 360.d0-eps) p = p - 360.d0 + if (p < 0.e0 - eps) p = p + 360.e0 + if (p >= 360.e0 - eps) p = p - 360.e0 ii = nint(p/drin) + 1 if (ii > NBP) ii = NBP @@ -484,8 +587,8 @@ subroutine get_indexloc(phi,theta,i,j,dr,LARG,ii,jj,phir,thetar) jj = nint(t/drin) + 1 if (jj > NBT) jj = NBT - thetar = sngl(t) - phir = sngl(p) + thetar = t + phir = p end subroutine get_indexloc @@ -493,7 +596,7 @@ end subroutine get_indexloc !-------------------------------------------------------------------------------------------------- ! - subroutine crust_bilinear_variable(theta,phi,gll_ind,rho_cr,vp_cr,vsv_cr,vsh_cr) + subroutine crust_bilinear_variable(theta,phi,gll_ind,rho_crust,vp_crust,vsv_crust,vsh_crust) !theta phi en radians !dr en degre @@ -503,18 +606,21 @@ subroutine crust_bilinear_variable(theta,phi,gll_ind,rho_cr,vp_cr,vsv_cr,vsh_cr) use constants, only: PI implicit none -! real, dimension(2) :: crust_array - integer :: gll_ind - double precision :: theta,phi,dx,dy,t,p,D_DEG - double precision, intent(out) :: rho_cr,vp_cr,vsv_cr,vsh_cr - integer :: broj_lokacija,j,m - double precision, parameter :: deg2rad = PI/180.d0, rad2deg = 180.d0/PI + integer, intent(in) :: gll_ind + double precision, intent(in) :: theta,phi + double precision, intent(out) :: rho_crust,vp_crust,vsv_crust,vsh_crust + + ! local parameters + double precision :: dx,dy,t,p,factor + integer :: j,m + double precision, parameter :: rad2deg = 180.d0/PI ! deg2rad = PI/180.d0 ! sfrench 20110103 model grid spacing (degrees) - D_DEG = 2.d0 + double precision, parameter :: D_DEG = 2.d0 + double precision, parameter :: D_DEG_INV = 1.d0 / D_DEG ! latitude in degrees - t = 90.d0 - theta / deg2rad + t = 90.d0 - theta * rad2deg ! sfrench 20110103 new bounds on latitude in accord with new model grid spacing !if (t>89.d0) t = 89.d0 @@ -523,44 +629,47 @@ subroutine crust_bilinear_variable(theta,phi,gll_ind,rho_cr,vp_cr,vsv_cr,vsh_cr) if (t < D_DEG - 90.d0) t = D_DEG - 90.d0 ! longitude in degrees - p = phi / deg2rad + p = phi * rad2deg if (p > 180.d0) p = p - 360.d0 - broj_lokacija = size(crust_array, DIM = 1) + m = 0 - rho_cr = 0.d0 - vp_cr = 0.d0 - vsv_cr = 0.d0 - vsh_cr = 0.d0 + rho_crust = 0.d0 + vp_crust = 0.d0 + vsv_crust = 0.d0 + vsh_crust = 0.d0 - do j = 1,broj_lokacija - if (m >= 4) then - exit - else - if ( gll_ind == nint(crust_array(j,3) + 1) ) then - dy = dabs(crust_array(j,1) - t) - dy = dy / D_DEG ! sfrench 20110103 : normalized to model grid spacing - - if (dy < 1.d0) then - dx = dabs(crust_array(j,2) - p) - if (dx > 180.d0) dx = 360.d0 - dx - - dx = dx / D_DEG ! sfrench 20110103 : normalized to model grid spacing - - if (dabs(dx) < 1.d0) then - ! Increment number of found locations. Must be <= 4 - m = m + 1 - rho_cr = rho_cr + (1.d0-dx) * (1.d0-dy) * crust_array(j,4) - vp_cr = vp_cr + (1.d0-dx) * (1.d0-dy) * crust_array(j,5) - vsv_cr = vsv_cr + (1.d0-dx) * (1.d0-dy) * crust_array(j,6) - vsh_cr = vsh_cr + (1.d0-dx) * (1.d0-dy) * crust_array(j,7) - - ! print *,'m= ',m,'dx= ',dx,' dy= ',dy, 'Vsv= ',vsv_cr - endif + do j = 1,num_crust_array + ! m must be < 4 + if (m >= 4) exit + + if (gll_ind == crust_array_ind(j)) then + dy = dabs(crust_array(1,j) - t) + dy = dy * D_DEG_INV ! sfrench 20110103 : normalized to model grid spacing + + if (dy < 1.d0) then + dx = dabs(crust_array(2,j) - p) + if (dx > 180.d0) dx = 360.d0 - dx + + dx = dx * D_DEG_INV ! sfrench 20110103 : normalized to model grid spacing + + if (dabs(dx) < 1.d0) then + ! Increment number of found locations. Must be <= 4 + m = m + 1 + + factor = (1.d0-dx) * (1.d0-dy) + + rho_crust = rho_crust + factor * crust_array(3,j) + vp_crust = vp_crust + factor * crust_array(4,j) + vsv_crust = vsv_crust + factor * crust_array(5,j) + vsh_crust = vsh_crust + factor * crust_array(6,j) + + ! print *,'m= ',m,'dx= ',dx,' dy= ',dy, 'Vsv= ',vsv_crust endif endif endif + enddo end subroutine crust_bilinear_variable @@ -575,20 +684,40 @@ subroutine read_crust_smooth_variable(unit) implicit none integer, intent(in) :: unit - integer :: j,l,broj_lokacija + ! local parameters + integer :: j,ier + double precision :: t,p,rho,vp,vsv,vsh + integer :: gll_ind - read(unit,*) broj_lokacija + read(unit,*) num_crust_array - ! sfrench 20110103 note: crust_array is < t> < p> < gll_ind> < rho> < vp> < vsv> < vsh> - allocate(crust_array(broj_lokacija,7)) + ! allocates crustal array + ! note: file format is: #t #p #gll_ind #rho #vp #vsv #vsh + allocate(crust_array(6,num_crust_array), & + crust_array_ind(num_crust_array),stat=ier) + if (ier /= 0) stop 'Error allocating crust_array,..' crust_array(:,:) = -1000.0 - - do j = 1,broj_lokacija - read(unit,*) (crust_array(j,l),l=1,7) + crust_array_ind(:) = 0 + + do j = 1,num_crust_array + ! format: #t #p #gll_ind #rho #vp #vsv #vsh + read(unit,*) t,p,gll_ind,rho,vp,vsv,vsh + + ! store values + ! GLL point index + crust_array_ind(j) = gll_ind + 1 ! GLL point index in range [1,5] + + ! point values + crust_array(1,j) = t ! latitude in degree [-90,90] + crust_array(2,j) = p ! longitude in degree [-180,180] + crust_array(3,j) = rho + crust_array(4,j) = vp + crust_array(5,j) = vsv + crust_array(6,j) = vsh enddo !debug - !print *,'read_crust_smooth_variable: I have read ',broj_lokacija,' crustal inputs!' + !print *,'read_crust_smooth_variable: I have read ',num_crust_array,' crustal inputs!' end subroutine read_crust_smooth_variable @@ -605,18 +734,24 @@ subroutine read_crustmoho_filtre(unit) integer :: i,j read(unit,*) NBP,NBT,drin + + ! checks if (drin /= 2.) STOP 'read_crust_filtre: dr muste be == 2' + ! allocates moho array allocate(moho_start(NBP,NBT)) moho_start(:,:) = -1000.0 + ! reads in values do j = 1,NBP do i = 1,NBT read(unit,*) moho_start(j,i) enddo enddo + ! convert to km moho_start(:,:) = moho_start(:,:) / 1000.0 + NBT = NBT - 1 end subroutine read_crustmoho_filtre diff --git a/src/meshfem3D/model_gll.f90 b/src/meshfem3D/model_gll.f90 index 649b214b2..3814907bc 100644 --- a/src/meshfem3D/model_gll.f90 +++ b/src/meshfem3D/model_gll.f90 @@ -39,7 +39,15 @@ module model_gll_par ! GLL model_variables type model_gll_variables - sequence + !TODO: check if `sequence` is needed + ! + ! in principle, `sequence` is only needed to explicity map the memory structure, for example when + ! the variable like `MGLL_V` gets passed as a routine argument, making sure it gets accessed correctly. + ! however, here we don't pass the objects MGLL_V, MGLL_V_IC to routines, but only use the objects to store arrays + ! within this module. + ! + !sequence + ! tomographic iteration model on GLL points real(kind=CUSTOM_REAL) :: scale_velocity,scale_density,scale_GPa @@ -54,7 +62,9 @@ module model_gll_par ! number of elements (crust/mantle elements) integer :: nspec - integer :: dummy_pad ! padding 4 bytes to align the structure + + ! in case we use `sequence` to align memory addresses + !integer :: dummy_pad ! padding 4 bytes to align the structure end type model_gll_variables ! crust/mantle model @@ -967,7 +977,7 @@ subroutine model_gll_build_cij(r,theta,phi, & Gc = Gc_prime * mu0 Gs = Gs_prime * mu0 - !daniel debug: + ! debug ! test: ignore Gc,Gs contributions !Gc = 0.d0 !Gs = 0.d0 diff --git a/src/meshfem3D/moho_stretching.f90 b/src/meshfem3D/moho_stretching.f90 index 2c32aec10..a10e2db14 100644 --- a/src/meshfem3D/moho_stretching.f90 +++ b/src/meshfem3D/moho_stretching.f90 @@ -130,8 +130,8 @@ subroutine moho_stretching_honor_crust(xelm,yelm,zelm, & endif ! radii for stretching criteria - R_moho = RMOHO_FICTITIOUS_IN_MESHER/R_PLANET - R_middlecrust = RMIDDLE_CRUST/R_PLANET + R_moho = RMOHO_FICTITIOUS_IN_MESHER / R_PLANET + R_middlecrust = RMIDDLE_CRUST / R_PLANET ! loops over element's anchor points count_crust = 0 @@ -171,7 +171,7 @@ subroutine moho_stretching_honor_crust(xelm,yelm,zelm, & ! ! note: flag found_crust returns .false. for points below moho, ! nevertheless its moho depth should be set and will be used in linear stretching - if (moho < TINYVAL ) call exit_mpi(myrank,'Error moho depth to honor') + if (moho < TINYVAL) call exit_mpi(myrank,'Error moho depth to honor') ! limits moho depth to a threshold value to avoid stretching problems if (moho < MOHO_MINIMUM) then diff --git a/src/meshfem3D/save_arrays_solver.f90 b/src/meshfem3D/save_arrays_solver.f90 index 6298d88d7..c790dcc45 100644 --- a/src/meshfem3D/save_arrays_solver.f90 +++ b/src/meshfem3D/save_arrays_solver.f90 @@ -313,6 +313,8 @@ subroutine save_arrays_solver(idoubling,ibool,xstore,ystore,zstore, & write(filename,'(a,i6.6)') trim(OUTPUT_FILES)//'/ispec_is_tiso',myrank call write_VTK_data_elem_l(nspec,nglob,xstore_glob,ystore_glob,zstore_glob,ibool, & ispec_is_tiso,filename) + ! debug + print *,'debug: rank',myrank,' tiso count = ',count(ispec_is_tiso) endif endif diff --git a/src/meshfem3D/save_arrays_solver_adios.F90 b/src/meshfem3D/save_arrays_solver_adios.F90 index d59b8e2ba..8f0a37506 100644 --- a/src/meshfem3D/save_arrays_solver_adios.F90 +++ b/src/meshfem3D/save_arrays_solver_adios.F90 @@ -578,7 +578,7 @@ subroutine save_arrays_solver_adios(idoubling,ibool,xstore,ystore,zstore, & NSPEC2DMAX_XMIN_XMAX_wmax = ints_to_reduce(3) NSPEC2DMAX_YMIN_YMAX_wmax = ints_to_reduce(4) - ! debug daniel + ! debug !call synchronize_all() !print *,'debug: ',myrank,'nspec2d top :',NSPEC2D_TOP,NSPEC2D_TOP_wmax !print *,'debug: ',myrank,'nspec2d bottom :',NSPEC2D_BOTTOM,NSPEC2D_BOTTOM_wmax diff --git a/src/meshfem3D/save_model_meshfiles_adios.F90 b/src/meshfem3D/save_model_meshfiles_adios.F90 index 258f2438c..c73227380 100644 --- a/src/meshfem3D/save_model_meshfiles_adios.F90 +++ b/src/meshfem3D/save_model_meshfiles_adios.F90 @@ -126,7 +126,7 @@ subroutine save_model_meshfiles_adios() write(region_name_scalar,"('reg',i1)") iregion_code write(group_name,"('SPECFEM3D_GLOBE_MODEL_reg',i1)") iregion_code -!daniel + ! note: on Mac OsX with OpenMPI 3.1, an error can occur at the end when finalizing MPI: ! !-------------------------------------------------------------------------- diff --git a/src/meshfem3D/setup_model.f90 b/src/meshfem3D/setup_model.f90 index 6bc18d235..7c4f0eff5 100644 --- a/src/meshfem3D/setup_model.f90 +++ b/src/meshfem3D/setup_model.f90 @@ -130,7 +130,8 @@ end subroutine setup_model subroutine sm_output_info() - use constants, only: IMAIN,NGLLX,NGLLY,NGLLZ,NGNOD,NGNOD2D,N_SLS,ASSUME_PERFECT_SPHERE + use constants, only: IMAIN,NGLLX,NGLLY,NGLLZ,NGNOD,NGNOD2D,N_SLS,ASSUME_PERFECT_SPHERE, & + USE_OLD_VERSION_FORMAT use shared_parameters, only: R_PLANET_KM use meshfem_models_par @@ -164,6 +165,12 @@ subroutine sm_output_info() write(IMAIN,*) 'Surface shape functions defined by NGNOD2D = ',NGNOD2D,' control nodes' write(IMAIN,*) + ! backward compatibility + if (USE_OLD_VERSION_FORMAT) then + write(IMAIN,*) 'using old version backward compatibility (versions 7.0 to 8.0)' + write(IMAIN,*) + endif + ! model user parameters write(IMAIN,*) 'model: ',trim(MODEL) if (OCEANS) then diff --git a/src/shared/auto_ner.f90 b/src/shared/auto_ner.f90 index 732cb5867..d63e47944 100644 --- a/src/shared/auto_ner.f90 +++ b/src/shared/auto_ner.f90 @@ -272,7 +272,7 @@ end subroutine auto_time_stepping ! subroutine auto_attenuation_periods(WIDTH, NEX_MAX, MIN_ATTENUATION_PERIOD, MAX_ATTENUATION_PERIOD) - use constants, only: N_SLS + use constants, only: N_SLS,USE_OLD_VERSION_FORMAT use shared_parameters, only: T_min_period @@ -318,6 +318,13 @@ subroutine auto_attenuation_periods(WIDTH, NEX_MAX, MIN_ATTENUATION_PERIOD, MAX_ ! maximum period MAX_ATTENUATION_PERIOD = tmp + ! older version uses integer values for min/max attenuation period + if (USE_OLD_VERSION_FORMAT) then + ! convert to integer (and back to double) + MIN_ATTENUATION_PERIOD = dble(int(MIN_ATTENUATION_PERIOD)) + MAX_ATTENUATION_PERIOD = dble(int(MAX_ATTENUATION_PERIOD)) + endif + ! check if (MIN_ATTENUATION_PERIOD <= 0.d0) then print *,'Error: invalid attenuation minimum: ',MIN_ATTENUATION_PERIOD diff --git a/src/shared/broadcast_computed_parameters.f90 b/src/shared/broadcast_computed_parameters.f90 index 0ba221f5e..1ce5aeb98 100644 --- a/src/shared/broadcast_computed_parameters.f90 +++ b/src/shared/broadcast_computed_parameters.f90 @@ -27,7 +27,7 @@ subroutine broadcast_computed_parameters() - use constants, only: myrank,STF_IS_UCB_HEAVISIDE + use constants, only: myrank use shared_parameters implicit none @@ -206,6 +206,7 @@ subroutine broadcast_computed_parameters() call bcast_all_singlei(NZ_DOUBLING_5) ! (optional) Berkeley UCB stf + call bcast_all_singlel(STF_IS_UCB_HEAVISIDE) if (STF_IS_UCB_HEAVISIDE) then call bcast_all_singledp(UCB_SOURCE_T1) call bcast_all_singledp(UCB_SOURCE_T2) diff --git a/src/shared/define_all_layers.f90 b/src/shared/define_all_layers.f90 index e8aea5803..e3538dc84 100644 --- a/src/shared/define_all_layers.f90 +++ b/src/shared/define_all_layers.f90 @@ -536,7 +536,7 @@ subroutine define_all_layers(NUMBER_OF_MESH_LAYERS,layer_offset,last_doubling_la ! simulations (larger time step), 1D models can be run with just one average crustal ! layer instead of two. - !daniel debug + ! debug !print *,'one_crust case in define_all_layers' ! check with define_all_layers_number_and_offset() @@ -769,7 +769,7 @@ subroutine define_all_layers(NUMBER_OF_MESH_LAYERS,layer_offset,last_doubling_la ! contains the crustal layers ! doubling at the base of the crust - !daniel debug + ! debug !print *,'default case in define_all_layers' ! check with define_all_layers_number_and_offset() diff --git a/src/shared/get_model_parameters.F90 b/src/shared/get_model_parameters.F90 index 097e65f0d..fc87ab5c3 100644 --- a/src/shared/get_model_parameters.F90 +++ b/src/shared/get_model_parameters.F90 @@ -350,6 +350,11 @@ subroutine get_model_parameters_flags() HONOR_1D_SPHERICAL_MOHO = .true. REFERENCE_1D_MODEL = REFERENCE_MODEL_CCREM + case ('1d_berkeley') + REFERENCE_1D_MODEL = REFERENCE_MODEL_SEMUCB + TRANSVERSE_ISOTROPY = .true. + HONOR_1D_SPHERICAL_MOHO = .true. + ! Mars 1D models case ('1d_sohl') ! Mars model A @@ -1153,6 +1158,17 @@ subroutine get_model_planet_constants() HONOR_DEEP_MOHO = EARTH_HONOR_DEEP_MOHO end select + ! special cases - additional overwrites + ! Berkeley model + if (REFERENCE_1D_MODEL == REFERENCE_MODEL_SEMUCB) then + ! SEMUCB uses a homogenized crustal model with its own smooth topo + ! uses default Berkeley topo + PATHNAME_TOPO_FILE = PATHNAME_TOPO_FILE_BERKELEY + RESOLUTION_TOPO_FILE = RESOLUTION_TOPO_FILE_BERKELEY + NX_BATHY = NX_BATHY_BERKELEY + NY_BATHY = NY_BATHY_BERKELEY + endif + end subroutine get_model_planet_constants ! @@ -1229,7 +1245,11 @@ subroutine get_model_parameters_radii() R771 = PREM_R771 RTOPDDOUBLEPRIME = PREM_RTOPDDOUBLEPRIME RCMB = PREM_RCMB - RICB = PREM_RICB + if (USE_OLD_VERSION_FORMAT) then + RICB = PREM_RICB_OLD + else + RICB = PREM_RICB + endif ! non-dimensionalizes densities ! density ocean diff --git a/src/shared/get_timestep_and_layers.f90 b/src/shared/get_timestep_and_layers.f90 index 133b29e29..8b003d48b 100644 --- a/src/shared/get_timestep_and_layers.f90 +++ b/src/shared/get_timestep_and_layers.f90 @@ -949,7 +949,7 @@ end subroutine get_timestep_limit_significant_digit subroutine get_minimum_period_estimate() - use constants, only: NGLLX,PI,NPTS_PER_WAVELENGTH,REFERENCE_MODEL_CASE65TAY + use constants, only: NGLLX,PI,NPTS_PER_WAVELENGTH,REFERENCE_MODEL_CASE65TAY,USE_OLD_VERSION_FORMAT use shared_parameters, only: T_min_period,estimated_min_wavelength, & ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, & @@ -1062,6 +1062,10 @@ subroutine get_minimum_period_estimate() ! example Earth: radius 6371 km -> 2 * PI * R / 360 ~ 111.19 km deg2km = R_PLANET / 1000.d0 * 2.d0 * PI / 360.d0 + ! for compatibility w/ older versions (7.0, 8.0,..) + ! (fixed value taken from old version auto_attenuation_periods()) + if (USE_OLD_VERSION_FORMAT) deg2km = 111.00d0 + ! computes Min Period ! ! width of element in km = (Angular width in degrees / NEX_MAX) * degrees to km @@ -1074,7 +1078,7 @@ subroutine get_minimum_period_estimate() tmp = tmp * dble(NPTS_PER_WAVELENGTH-1) ! minimum period resolved = (minimum wavelength) / V_min - tmp = tmp/S_VELOCITY_MIN + tmp = tmp / S_VELOCITY_MIN ! estimated minimum period resolved T_min_period = tmp diff --git a/src/shared/make_ellipticity.f90 b/src/shared/make_ellipticity.f90 index 2d29507f5..e1c529dc2 100644 --- a/src/shared/make_ellipticity.f90 +++ b/src/shared/make_ellipticity.f90 @@ -85,7 +85,7 @@ subroutine make_ellipticity_epsilon_eta(nspl,rspl,ellipicity_spline,ellipicity_s ! for mars, creates a spline for the ellipticity profile in Mars (Sohl $ Spohn) - use constants, only: NR_DENSITY,TWO_PI,PI,GRAV,R_UNIT_SPHERE,myrank + use constants, only: NR_DENSITY,TWO_PI,PI,GRAV,R_UNIT_SPHERE,myrank,USE_OLD_VERSION_FORMAT use shared_parameters, only: PLANET_TYPE,IPLANET_EARTH,IPLANET_MARS,IPLANET_MOON, & R_PLANET,R_PLANET_KM,RHOAV, & @@ -137,9 +137,9 @@ subroutine make_ellipticity_epsilon_eta(nspl,rspl,ellipicity_spline,ellipicity_s ! please check the effect. ! Earth - ! PREM radius of the Earth for gravity calculation + ! PREM radius of the Earth for calculation double precision, parameter :: R_EARTH_ELLIPTICITY = 6371000.d0 - ! PREM radius of the ocean floor for gravity calculation + ! PREM radius of the ocean floor for calculation double precision, parameter :: ROCEAN_ELLIPTICITY = PREM_ROCEAN ! initializes @@ -162,7 +162,7 @@ subroutine make_ellipticity_epsilon_eta(nspl,rspl,ellipicity_spline,ellipicity_s case (IPLANET_EARTH) ! Earth ! default PREM - ! radius of the planet for gravity calculation + ! radius of the planet for calculation RSURFACE = R_EARTH_ELLIPTICITY ! physical surface (Earth: 6371000, ..) ROCEAN = ROCEAN_ELLIPTICITY RMIDDLE_CRUST = PREM_RMIDDLE_CRUST @@ -175,7 +175,11 @@ subroutine make_ellipticity_epsilon_eta(nspl,rspl,ellipicity_spline,ellipicity_s R771 = PREM_R771 RTOPDDOUBLEPRIME = PREM_RTOPDDOUBLEPRIME RCMB = PREM_RCMB - RICB = PREM_RICB + if (USE_OLD_VERSION_FORMAT) then + RICB = PREM_RICB_OLD + else + RICB = PREM_RICB + endif case (IPLANET_MARS) ! Mars diff --git a/src/shared/model_prem.f90 b/src/shared/model_prem.f90 index be8e94f95..d78882131 100644 --- a/src/shared/model_prem.f90 +++ b/src/shared/model_prem.f90 @@ -74,7 +74,7 @@ module model_prem_par double precision, parameter :: PREM_RCMB = 3480000.d0 ! 2891 km depth ! note: SPECFEM versions up to 8.0 (Aug, 2021) used an inner core radius of 1221 km; ! based on Table 1 in Dziewonski&Anderson's PREM paper, the inner core radius is at 1221.5 km - !double precision, parameter :: PREM_RICB = 1221000.d0 ! old versions + double precision, parameter :: PREM_RICB_OLD = 1221000.d0 ! old versions double precision, parameter :: PREM_RICB = 1221500.d0 ! PREM2 additional radii for modifications @@ -114,6 +114,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! local parameters double precision :: r,scaleval + double precision :: ROCEAN,RMIDDLE_CRUST,RMOHO,R80,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB ! compute real physical radius in meters r = x * R_PLANET @@ -123,11 +124,49 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! note: using stop statements, not exit_mpi() calls to avoid the need for MPI libraries when linking xcreate_header_file - if (check_doubling_flag) then + ! default setting + ROCEAN = PREM_ROCEAN + RMIDDLE_CRUST = PREM_RMIDDLE_CRUST + RMOHO = PREM_RMOHO + R80 = PREM_R80 + R220 = PREM_R220 + R400 = PREM_R400 + R600 = PREM_R600 + R670 = PREM_R670 + R771 = PREM_R771 + RTOPDDOUBLEPRIME = PREM_RTOPDDOUBLEPRIME + RCMB = PREM_RCMB + RICB = PREM_RICB + + ! note: different versions use different radii when calling this routine + if (USE_OLD_VERSION_FORMAT) then + ! version compatibility with older versions 7.0 - 8.0 + RICB = PREM_RICB_OLD + + ! version compatibility with special case Berkeley model + ! note: old Berkeley implementation used radii different to PREM when calling this routine for gravity evaluation. + ! this was inconsistent with the the spline evaluation in the old make_gravity.f90. + ! the change here below is to re-introduce this inconsistency to match the old Berkeley model output. + ! + ! new(er) versions use the "original" PREM model for calculating gravity (and ellipticity). + ! that is, gravity (acceleration g) and ellipticity (based on density model) is based on PREM, since we know that + ! PREM was constructed with corresponding constraints. + ! other reference models might better adapt to velocities, but might violate for example total mass constraints. + ! the different radii of the reference models are used for meshing purposes only. + ! + if (REFERENCE_1D_MODEL == REFERENCE_MODEL_SEMUCB) then + RICB = PREM_RICB + RMOHO = 6341000.d0 + R400 = 5961000.d0 + R670 = 5721000.d0 + endif + endif + + if (check_doubling_flag) then ! !--- inner core ! - if (r >= 0.d0 .and. r < PREM_RICB) then + if (r >= 0.d0 .and. r < RICB) then if (idoubling /= IFLAG_INNER_CORE_NORMAL .and. & idoubling /= IFLAG_MIDDLE_CENTRAL_CUBE .and. & idoubling /= IFLAG_BOTTOM_CENTRAL_CUBE .and. & @@ -137,31 +176,31 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! !--- outer core ! - else if (r > PREM_RICB .and. r < PREM_RCMB) then + else if (r > RICB .and. r < RCMB) then if (idoubling /= IFLAG_OUTER_CORE_NORMAL) & stop 'wrong doubling flag for outer core point in model_prem_iso()' ! !--- D" at the base of the mantle ! - else if (r > PREM_RCMB .and. r < PREM_RTOPDDOUBLEPRIME) then + else if (r > RCMB .and. r < RTOPDDOUBLEPRIME) then if (idoubling /= IFLAG_MANTLE_NORMAL) & stop 'wrong doubling flag for D" point in model_prem_iso()' ! !--- mantle: from top of D" to d670 ! - else if (r > PREM_RTOPDDOUBLEPRIME .and. r < PREM_R670) then + else if (r > RTOPDDOUBLEPRIME .and. r < R670) then if (idoubling /= IFLAG_MANTLE_NORMAL) & stop 'wrong doubling flag for top D" to d670 point in model_prem_iso()' ! !--- mantle: from d670 to d220 ! - else if (r > PREM_R670 .and. r < PREM_R220) then + else if (r > R670 .and. r < R220) then if (idoubling /= IFLAG_670_220) & stop 'wrong doubling flag for d670 to d220 point in model_prem_iso()' ! !--- mantle and crust: from d220 to MOHO and then to surface ! - else if (r > PREM_R220) then + else if (r > R220) then if (idoubling /= IFLAG_220_80 .and. idoubling /= IFLAG_80_MOHO .and. idoubling /= IFLAG_CRUST) & stop 'wrong doubling flag for d220 to Moho to surface point in model_prem_iso()' endif @@ -171,7 +210,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! !--- inner core ! - if (r >= 0.d0 .and. r <= PREM_RICB) then + if (r >= 0.d0 .and. r <= RICB) then drhodr = -2.0d0*8.8381d0*x rho = 13.0885d0 - 8.8381d0*x*x if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then @@ -193,7 +232,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! !--- outer core ! - else if (r > PREM_RICB .and. r <= PREM_RCMB) then + else if (r > RICB .and. r <= RCMB) then drhodr = -1.2638d0 - 2.0d0*3.6426d0*x - 3.0d0*5.5281d0*x*x rho = 12.5815d0 - 1.2638d0*x - 3.6426d0*x*x - 5.5281d0*x*x*x if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then @@ -215,7 +254,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! !--- D" at the base of the mantle ! - else if (r > PREM_RCMB .and. r <= PREM_RTOPDDOUBLEPRIME) then + else if (r > RCMB .and. r <= RTOPDDOUBLEPRIME) then drhodr = -6.4761d0 + 2.0d0*5.5283d0*x - 3.0d0*3.0807d0*x*x rho = 7.9565d0 - 6.4761d0*x + 5.5283d0*x*x - 3.0807d0*x*x*x if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then @@ -231,7 +270,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! !--- mantle: from top of D" to d670 ! - else if (r > PREM_RTOPDDOUBLEPRIME .and. r <= PREM_R771) then + else if (r > RTOPDDOUBLEPRIME .and. r <= R771) then drhodr = -6.4761d0 + 2.0d0*5.5283d0*x - 3.0d0*3.0807d0*x*x rho = 7.9565d0 - 6.4761d0*x + 5.5283d0*x*x - 3.0807d0*x*x*x if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then @@ -250,7 +289,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ vs = 11.1671d0 - 13.7818d0*x + 17.4575d0*x*x - 9.2777d0*x*x*x Qmu = 312.0d0 Qkappa = 57827.0d0 - else if (r > PREM_R771 .and. r <= PREM_R670) then + else if (r > R771 .and. r <= R670) then drhodr = -6.4761d0 + 2.0d0*5.5283d0*x - 3.0d0*3.0807d0*x*x rho = 7.9565d0 - 6.4761d0*x + 5.5283d0*x*x - 3.0807d0*x*x*x vp = 29.2766d0 - 23.6027d0*x + 5.5242d0*x*x - 2.5514d0*x*x*x @@ -260,28 +299,28 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ ! !--- mantle: above d670 ! - else if (r > PREM_R670 .and. r <= PREM_R600) then + else if (r > R670 .and. r <= R600) then drhodr = -1.4836d0 rho = 5.3197d0 - 1.4836d0*x vp = 19.0957d0 - 9.8672d0*x vs = 9.9839d0 - 4.9324d0*x Qmu = 143.0d0 Qkappa = 57827.0d0 - else if (r > PREM_R600 .and. r <= PREM_R400) then + else if (r > R600 .and. r <= R400) then drhodr = -8.0298d0 rho = 11.2494d0 - 8.0298d0*x vp = 39.7027d0 - 32.6166d0*x vs = 22.3512d0 - 18.5856d0*x Qmu = 143.0d0 Qkappa = 57827.0d0 - else if (r > PREM_R400 .and. r <= PREM_R220) then + else if (r > R400 .and. r <= R220) then drhodr = -3.8045d0 rho = 7.1089d0 - 3.8045d0*x vp = 20.3926d0 - 12.2569d0*x vs = 8.9496d0 - 4.4597d0*x Qmu = 143.0d0 Qkappa = 57827.0d0 - else if (r > PREM_R220 .and. r <= PREM_R80) then + else if (r > R220 .and. r <= R80) then drhodr = 0.6924d0 rho = 2.6910d0 + 0.6924d0*x vp = 4.1875d0 + 3.9382d0*x @@ -291,7 +330,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ else if (CRUSTAL .and. .not. SUPPRESS_CRUSTAL_MESH) then ! fill with PREM mantle and later add CRUST2.0 - if (r > PREM_R80) then + if (r > R80) then ! density/velocity from mantle just below moho drhodr = 0.6924d0 rho = 2.6910d0 + 0.6924d0*x @@ -303,7 +342,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ endif else ! use PREM crust - if (r > PREM_R80 .and. r <= PREM_RMOHO) then + if (r > R80 .and. r <= RMOHO) then drhodr = 0.6924d0 rho = 2.6910d0 + 0.6924d0*x vp = 4.1875d0 + 3.9382d0*x @@ -314,13 +353,13 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ else if (SUPPRESS_CRUSTAL_MESH) then ! extend the Moho up to the surface instead of the crust drhodr = 0.6924d0 - rho = 2.6910d0 + 0.6924d0*(PREM_RMOHO / R_PLANET) - vp = 4.1875d0 + 3.9382d0*(PREM_RMOHO / R_PLANET) - vs = 2.1519d0 + 2.3481d0*(PREM_RMOHO / R_PLANET) + rho = 2.6910d0 + 0.6924d0*(RMOHO / R_PLANET) + vp = 4.1875d0 + 3.9382d0*(RMOHO / R_PLANET) + vs = 2.1519d0 + 2.3481d0*(RMOHO / R_PLANET) Qmu = 600.0d0 Qkappa = 57827.0d0 - else if (r > PREM_RMOHO .and. r <= PREM_RMIDDLE_CRUST) then + else if (r > RMOHO .and. r <= RMIDDLE_CRUST) then drhodr = 0.0d0 rho = 2.9d0 vp = 6.8d0 @@ -338,7 +377,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ Qkappa = 57827.0d0 endif - else if (r > PREM_RMIDDLE_CRUST .and. r <= PREM_ROCEAN) then + else if (r > RMIDDLE_CRUST .and. r <= ROCEAN) then drhodr = 0.0d0 rho = 2.6d0 vp = 5.8d0 @@ -346,7 +385,7 @@ subroutine model_prem_iso(x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL,check_ Qmu = 600.0d0 Qkappa = 57827.0d0 ! for density profile for gravity, we do not check that r <= R_PLANET - else if (r > PREM_ROCEAN) then + else if (r > ROCEAN) then drhodr = 0.0d0 rho = 2.6d0 vp = 5.8d0 @@ -394,6 +433,7 @@ subroutine model_prem_aniso(x,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling ! local parameters double precision :: r double precision :: scaleval + double precision :: RICB ! compute real physical radius in meters r = x * R_PLANET @@ -401,11 +441,18 @@ subroutine model_prem_aniso(x,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling ! check flags to make sure we correctly honor the discontinuities ! we use strict inequalities since r has been slightly changed in mesher + ! version compatibility + if (USE_OLD_VERSION_FORMAT) then + RICB = PREM_RICB_OLD + else + RICB = PREM_RICB + endif + if (check_doubling_flag) then ! !--- inner core ! - if (r >= 0.d0 .and. r < PREM_RICB) then + if (r >= 0.d0 .and. r < RICB) then if (idoubling /= IFLAG_INNER_CORE_NORMAL .and. & idoubling /= IFLAG_MIDDLE_CENTRAL_CUBE .and. & idoubling /= IFLAG_BOTTOM_CENTRAL_CUBE .and. & @@ -415,7 +462,7 @@ subroutine model_prem_aniso(x,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling ! !--- outer core ! - else if (r > PREM_RICB .and. r < PREM_RCMB) then + else if (r > RICB .and. r < PREM_RCMB) then if (idoubling /= IFLAG_OUTER_CORE_NORMAL) & stop 'wrong doubling flag for outer core point in model_prem_aniso()' ! @@ -453,7 +500,7 @@ subroutine model_prem_aniso(x,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling ! !--- inner core ! - if (r >= 0.d0 .and. r <= PREM_RICB) then + if (r >= 0.d0 .and. r <= RICB) then rho = 13.0885d0 - 8.8381d0*x*x if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then ! PREM2 @@ -476,7 +523,7 @@ subroutine model_prem_aniso(x,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling ! !--- outer core ! - else if (r > PREM_RICB .and. r <= PREM_RCMB) then + else if (r > RICB .and. r <= PREM_RCMB) then rho = 12.5815d0 - 1.2638d0*x - 3.6426d0*x*x - 5.5281d0*x*x*x if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM2) then ! PREM2 @@ -853,14 +900,22 @@ subroutine prem_density(x,rho) ! local parameters double precision :: r + double precision :: RICB ! compute real physical radius in meters r = x * R_PLANET + ! version compatibility + if (USE_OLD_VERSION_FORMAT) then + RICB = PREM_RICB_OLD + else + RICB = PREM_RICB + endif + ! calculates density according to radius - if (r <= PREM_RICB) then + if (r <= RICB) then rho = 13.0885d0 - 8.8381d0*x*x - else if (r > PREM_RICB .and. r <= PREM_RCMB) then + else if (r > RICB .and. r <= PREM_RCMB) then rho = 12.5815d0 - 1.2638d0*x - 3.6426d0*x*x - 5.5281d0*x*x*x else if (r > PREM_RCMB .and. r <= PREM_RTOPDDOUBLEPRIME) then rho = 7.9565d0 - 6.4761d0*x + 5.5283d0*x*x - 3.0807d0*x*x*x diff --git a/src/shared/model_topo_bathy.f90 b/src/shared/model_topo_bathy.f90 index 9946e7522..847f82c8f 100644 --- a/src/shared/model_topo_bathy.f90 +++ b/src/shared/model_topo_bathy.f90 @@ -47,13 +47,20 @@ module smooth_etopo5_par private + ! unused !integer, parameter :: NBT_TOPO5 = 2160, NBP_TOPO5 = 4320, INITR_TOPO5 = 12 !integer, dimension(:,:), allocatable :: topo_start - !integer :: NBT,NBP,INITR - + ! values taken from constants.h integer, parameter :: NBT = NY_BATHY_BERKELEY - 1 ! 181 - 1 == 180 to match gauss_filtre integer, parameter :: NBP = NX_BATHY_BERKELEY ! 360 + + ! resolution + ! given in minutes -> 1/60 to get resolution in degrees + ! default: + ! RESOLUTION_TOPO_FILE_BERKELEY == 60 -> INITR = int( 1.0 / ( 60 / 60.0) ) == 1 + ! or finer example: (would need NX/NY_BATHY adapted also) + ! RESOLUTION_TOPO_FILE_BERKELEY == 30 -> INITR = int( 1.0 / ( 30 / 60.0) ) == 2 integer, parameter :: INITR = int(1.0 / (RESOLUTION_TOPO_FILE_BERKELEY / 60.0)) ! resolution contains @@ -125,64 +132,70 @@ double precision function get_etopo5_filtre(xlat,xlon,ibathy_topo) integer, dimension(NX_BATHY_BERKELEY,NY_BATHY_BERKELEY),intent(in) :: ibathy_topo ! local parameters - real :: value - real :: theta,phi - !real, parameter :: pi = 3.141592653589793, deg2rad = pi / 180. - - real, parameter :: drfiltre = 2.e0 + double precision :: value + double precision :: colat,lon ! colatitude/longitude - theta = sngl(90.d0 - xlat) - phi = sngl(xlon) - - if (phi < -180.0) phi = phi + 360.0 - if (phi > 180.0) phi = phi - 360.0 + colat = 90.d0 - xlat ! range [0,180] + lon = xlon - ! in rad - !theta = theta * deg2rad - !phi = phi * deg2rad - ! elevation (in meters) - !value = dble(etopo_filtre(theta,phi)) + ! puts lon in range [-180,180] + if (lon < -180.d0) lon = lon + 360.d0 + if (lon > 180.d0) lon = lon - 360.d0 - ! or - ! in degrees - value = gauss_filtre(ibathy_topo,theta,phi,drfiltre) + ! filtered topo value + value = gauss_filtre(ibathy_topo,colat,lon) ! return value - get_etopo5_filtre = dble(value) + get_etopo5_filtre = value end function get_etopo5_filtre !----------------------------------------------------------------- - real function gauss_filtre(tin,theta,phi,dr) + double precision function gauss_filtre(tin,colat,lon) + + use constants, only: USE_OLD_VERSION_FORMAT implicit none - integer, dimension(:,:),intent(in) :: tin - real,intent(in) :: theta,phi,dr + integer, dimension(NX_BATHY_BERKELEY,NY_BATHY_BERKELEY),intent(in) :: tin + double precision, intent(in) :: colat,lon ! local parameters - real :: tmp,thetar,phir,tmpnorm,int_val + real :: t,p,thetar,phir + real :: tmp,tmpnorm,int_val integer :: i,ii,j,jj - real, parameter :: pi = 3.141592653589793, deg2rad = pi / 180. + real, parameter :: pi = 3.141592653589793, deg2rad = pi / 180.0 integer, parameter :: LARG = 2 - tmp = 0. - tmpnorm = 0. + real, parameter :: drfiltre = 2.e0 + + ! converts to real values + t = sngl(colat) + p = sngl(lon) + + tmp = 0.0 + tmpnorm = 0.0 - do i = 1,int(LARG*dr*INITR + 1) - do j = 1,int(LARG*dr*INITR + 1) - call get_indexloc(phi,theta,i,j,dr,LARG,ii,jj,phir,thetar) + do i = 1,int(LARG * drfiltre * INITR + 1) + do j = 1,int(LARG * drfiltre * INITR + 1) + call get_indexloc(p,t,i,j,drfiltre,LARG,ii,jj,phir,thetar) - int_val = cos_cylindre(theta,phi,dr,thetar,phir) * (1./real(INITR))**2 * sin(thetar * deg2rad) + int_val = cos_cylindre(t,p,drfiltre,thetar,phir) * (1./real(INITR))**2 * sin(thetar * deg2rad) tmp = tmp + tin(ii,jj) * int_val tmpnorm = tmpnorm + int_val enddo enddo - gauss_filtre = tmp / tmpnorm + gauss_filtre = dble(tmp / tmpnorm) + + ! old version uses integer for elevation values from gauss_filtre() routine + if (USE_OLD_VERSION_FORMAT) then + ! convert elevation to integer value before returning as a real/double number + gauss_filtre = dble(int(gauss_filtre)) + endif end function gauss_filtre @@ -191,10 +204,11 @@ end function gauss_filtre real function cos_cylindre(t0_,p0_,d0_,theta_,phi_) implicit none - real :: t0,p0,d0,theta,phi, d_ang - real :: t0_,p0_,d0_,theta_,phi_ + real, intent(in) :: t0_,p0_,d0_,theta_,phi_ + ! local parameters + real :: t0,p0,d0,theta,phi,d_ang real :: d_arg_clean - real, parameter :: pi = 3.141592653589793, deg2rad = pi / 180. + real, parameter :: pi = 3.141592653589793e0, deg2rad = pi / 180.e0 t0 = t0_ * deg2rad p0 = p0_ * deg2rad @@ -204,13 +218,13 @@ real function cos_cylindre(t0_,p0_,d0_,theta_,phi_) !distance angulaire au centre du cylindre: !d_ang=acos(cos(theta)*cos(t0)+sin(theta)*sin(t0)*cos(phi-p0)) - d_arg_clean = min( 1.0, max( 0.0, real( cos(theta)*cos(t0)+sin(theta)*sin(t0)*cos(phi-p0) ) ) ) + d_arg_clean = min( 1.e0, max( 0.e0, real( cos(theta) * cos(t0) + sin(theta) * sin(t0) * cos(phi-p0) ) ) ) d_ang = acos( d_arg_clean ) if (d_ang > d0) then - cos_cylindre = 0.0 + cos_cylindre = 0.e0 else - cos_cylindre = 0.5 * (1.0 + cos(PI*d_ang/d0)) + cos_cylindre = 0.5e0 * (1.e0 + cos(pi * d_ang/d0)) endif end function cos_cylindre @@ -233,20 +247,20 @@ subroutine get_indexloc(phi,theta,i,j,dr,LARG,ii,jj,phir,thetar) t = theta + (j - 1 - (LARG * dr * INITR)/2) / real(INITR) ! longitude in range [0,360] - if (p < 0.0 - eps) p = p + 360.0 - if (p >= 360.0 - eps) p = p - 360.0 + if (p < 0.e0 - eps) p = p + 360.e0 + if (p >= 360.e0 - eps) p = p - 360.e0 ! colatitude in range [0,180] - if (t > 180.0 - eps) then - t = t - 180.0 - p = 360.0 - p - else if (t < 0.0 - eps) then - t = 180.0 + t - p = 360.0 - p + if (t > 180.e0 - eps) then + t = t - 180.e0 + p = 360.e0 - p + else if (t < 0.e0 - eps) then + t = 180.e0 + t + p = 360.e0 - p endif - if (p < 0.0 - eps) p = p + 360.0 - if (p >= 360.0 - eps) p = p - 360.0 + if (p < 0.e0 - eps) p = p + 360.e0 + if (p >= 360.e0 - eps) p = p - 360.e0 ii = nint(p * INITR) + 1 @@ -419,6 +433,7 @@ subroutine model_topo_bathy_broadcast(ibathy_topo,LOCAL_PATH) double precision :: time1,tCPU double precision, external :: wtime + ! gets topo file ending ('.bin' or '.dat') ending = '' if (len_trim(PATHNAME_TOPO_FILE) > 4) ending = PATHNAME_TOPO_FILE(len_trim(PATHNAME_TOPO_FILE)-3:len_trim(PATHNAME_TOPO_FILE)) diff --git a/src/shared/read_parameter_file.F90 b/src/shared/read_parameter_file.F90 index 37c479fe9..bbbb9477e 100644 --- a/src/shared/read_parameter_file.F90 +++ b/src/shared/read_parameter_file.F90 @@ -221,6 +221,7 @@ subroutine read_parameter_file() if (ier /= 0) stop 'an error occurred while reading the parameter file: USE_MONOCHROMATIC_CMT_SOURCE' ! (optional) Berkeley UCB STF parameters + call read_value_logical(STF_IS_UCB_HEAVISIDE, 'STF_IS_UCB_HEAVISIDE', ier); ier = 0 if (STF_IS_UCB_HEAVISIDE) then call read_value_double_precision(UCB_SOURCE_T1,'SOURCE_T1', ier); ier = 0 call read_value_double_precision(UCB_SOURCE_T2,'SOURCE_T2', ier); ier = 0 diff --git a/src/shared/shared_par.f90 b/src/shared/shared_par.f90 index 30e943739..ab85596e7 100644 --- a/src/shared/shared_par.f90 +++ b/src/shared/shared_par.f90 @@ -209,6 +209,7 @@ module shared_input_parameters double precision :: FILESYSTEM_IO_BANDWIDTH = 0.d0 ! UCB Source time function parameters + logical :: STF_IS_UCB_HEAVISIDE = .false. ! source-time function is a UCB-style filtered heaviside double precision :: UCB_SOURCE_T1 = 400.d0, UCB_SOURCE_T2 = 250.d0, UCB_SOURCE_T3 = 53.d0, UCB_SOURCE_T4 = 40.d0 double precision :: UCB_TAU = 400.d0 diff --git a/src/shared/spl_A3d.c b/src/shared/spl_A3d.c index 251d92592..501f0a95a 100644 --- a/src/shared/spl_A3d.c +++ b/src/shared/spl_A3d.c @@ -41,14 +41,14 @@ #include // for malloc/free #define U (unsigned) -#define TOL 1.5 +#define TOL 1.5f float *farray1(int, int); void free_farray1(float *, int); void stop(char *); void fill_hh_A3d(float *, float *, int); -float spl_A3d(int, int, float *, float); +float spl_A3d(int, int, float *, float *, float); /* ----------------------------------------------------------------------------- */ @@ -61,27 +61,35 @@ float spl_A3d(int, int, float *, float); //#endif void -FC_FUNC(fspl,FSPL)(int *ord, int *nknots, double *knot, double *xi, double *rho) +FC_FUNC(fspl,FSPL)(int *ord, int *nknots, float *knot_s, float *knot_hh, double *xi, double *rho) { - float *knot_s,xi_s,val; - int ord_s,i; + float xi_s,val; + int ord_s; - ord_s = *ord-1; /* change from fortran to c array convention */ + ord_s = *ord - 1; /* change from fortran to c array convention */ xi_s = (float) (*xi); - knot_s = farray1(0,*nknots-1); - for(i=0; i < *nknots; i++) - knot_s[i] = (float) knot[i]; + // allocate temporary float array + //float *knot_s; + //knot_s = farray1(0,*nknots - 1); + //for(int i=0; i < *nknots; i++) + // knot_s[i] = (float) knot[i]; - val = spl_A3d(ord_s,*nknots,knot_s,xi_s); + // spline value + val = spl_A3d(ord_s, *nknots, knot_s, knot_hh, xi_s); + + // return value *rho = (double) val; - free_farray1(knot_s,0); + + // free temporary array + //free_farray1(knot_s,0); + return; } /* ----------------------------------------------------------------------------- */ -float spl_A3d(int ord, int nknots, float *knot, float xi) +float spl_A3d(int ord, int nknots, float *knot, float *hh, float xi) { /* ord: number of rho(x) nknots : # of knkots = Nx+1 (Nx=index of highest spline) @@ -89,248 +97,250 @@ float spl_A3d(int ord, int nknots, float *knot, float xi) splines defined as : f_i(x) = a_i(x-x_i)^3 + b_i(x-x_i)^2 + c_i(x-x_i) + d_i */ - int ii,Nx; - float *hh,rho_x; + int Nx; + float rho_x; float coefa,coefb,coefc,coefd; + float dxi; - Nx = nknots-1; - /* Compute vector hh of spacings */ - hh= farray1(0,Nx-1); - fill_hh_A3d(hh,knot,Nx); + Nx = nknots - 1; + + /* Compute vector hh of spacings + note: this is done already in the fortran routine. + the spline node radii won't change and so do their spacings hh + between different spline evaluation calls. + thus, there is no need to re-allocate and compute this array for each spline call. + */ + //float *hh; + //hh = farray1(0,Nx - 1); + //fill_hh_A3d(hh,knot,Nx); /* Consistency checks */ - if ((xi-(float)TOL)>knot[Nx]) { + if ((xi - TOL) > knot[Nx]) { printf("xi=%g / knot[%d]=%g",xi,Nx,knot[Nx]); stop("spl: xi>knot[Nx]"); } - else if ((xi+(float)TOL)Nx) + else if (ord > Nx) stop("spl: order > Nx"); - if (ord==0) { /* LHS */ - float denom; - denom = 3.*hh[ord]*hh[ord]+3.*hh[ord]*hh[ord+1]+hh[ord+1]*hh[ord+1]; + if (ord == 0) { /* LHS */ + float hh0 = hh[ord]; + float hhp1 = hh[ord+1]; + + float denom = 3.0f * hh0 * hh0 + 3.0f * hh0 * hhp1 + hhp1 * hhp1; + if (xi >= knot[ord] && xi <= knot[ord+1]) { /* x0<=x<=x1 */ - coefa = 4./(hh[ord]*(hh[ord]+hh[ord+1])*denom); - coefb = 0.0; - coefc = -12/denom; - coefd = 4*(2*hh[ord]+hh[ord+1])/denom; - - rho_x = coefa*(xi-knot[ord])*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefb*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefc*(xi-knot[ord]); - rho_x += coefd; + coefa = 4.0f / (hh0 * (hh0 + hhp1) * denom); + coefb = 0.0f; + coefc = -12.0f / denom; + coefd = 4.0f * (2.0f * hh0 + hhp1) / denom; + + dxi = xi - knot[ord]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi > knot[ord+1] && xi <= knot[ord+2]){ /* x1<=x<=x2 */ - coefa = -4./(hh[ord+1]*(hh[ord]+hh[ord+1])*denom); - coefb = 12/((hh[ord]+hh[ord+1])*denom); - coefc = -12.*hh[ord+1]/((hh[ord]+hh[ord+1])*denom); - coefd = 4.*hh[ord+1]*hh[ord+1]/((hh[ord]+hh[ord+1])*denom); - - rho_x = coefa*(xi-knot[ord+1])*(xi-knot[ord+1])*(xi-knot[ord+1]); - rho_x += coefb*(xi-knot[ord+1])*(xi-knot[ord+1]); - rho_x += coefc*(xi-knot[ord+1]); - rho_x += coefd; + coefa = -4.0f / (hhp1 * (hh0 + hhp1) * denom); + coefb = 12.0f / ((hh0 + hhp1) * denom); + coefc = -12.0f * hhp1 / ((hh0 + hhp1) * denom); + coefd = 4.0f * hhp1 * hhp1 / ((hh0 + hhp1) * denom); + + dxi = xi - knot[ord+1]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else /* x>x2 */ - rho_x = 0.0; + rho_x = 0.0f; } - else if (ord==1) { /* LHS+1 */ - float denom,denomsum,dd; - denom = (3.*hh[ord-1]*hh[ord-1]+4.*hh[ord-1]*hh[ord]+hh[ord]*hh[ord]+ - 2.*hh[ord-1]*hh[ord+1]+hh[ord]*hh[ord+1]); - denomsum = hh[ord-1]+hh[ord]+hh[ord+1]; - dd = denomsum*denom; + else if (ord == 1) { /* LHS+1 */ + float hh0 = hh[ord]; + float hhp1 = hh[ord+1]; + float hhm1 = hh[ord-1]; + + float denom = 3.0f * hhm1 * hhm1 + 4.0f * hhm1 * hh0 + hh0 * hh0 + 2.0f * hhm1 * hhp1 + hh0 * hhp1; + float denomsum = hhm1 + hh0 + hhp1; + float dd = denomsum * denom; + if (xi >= knot[ord-1] && xi <= knot[ord]) { /* x0<=x<=x1 */ - coefa = -4.*(3.*hh[ord-1]+2.*hh[ord]+hh[ord+1])/ - (hh[ord-1]*(hh[ord-1]+hh[ord])*dd); - coefb = 0.; - coefc = 12./denom; - coefd = 0.; - - rho_x = coefa*(xi-knot[ord-1])*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefb*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefc*(xi-knot[ord-1]); - rho_x += coefd; + coefa = -4.0f * (3.0f * hhm1 + 2.0f * hh0 + hhp1) / + (hhm1 * (hhm1 + hh0) * dd); + coefb = 0.0f; + coefc = 12.0f / denom; + coefd = 0.0f; + + dxi = xi - knot[ord-1]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord] && xi <= knot[ord+1]) { /* x1<=x<=x2 */ - coefa = 4.*(2.*hh[ord-1]*hh[ord-1]+6.*hh[ord-1]*hh[ord]+3.*hh[ord]*hh[ord]+3.*hh[ord-1]*hh[ord+1]+ - 3.*hh[ord]*hh[ord+1]+hh[ord+1]*hh[ord+1])/ - (hh[ord]*(hh[ord-1]+hh[ord])*(hh[ord]+hh[ord+1])*dd); - coefb = -12.*(3.*hh[ord-1]+2.*hh[ord]+hh[ord+1])/ - ((hh[ord-1]+hh[ord])*dd); - coefc = 12.*(-2.*hh[ord-1]*hh[ord-1]+hh[ord]*hh[ord]+hh[ord]*hh[ord+1])/ - ((hh[ord-1]+hh[ord])*dd); - coefd = 4.*hh[ord-1]*(4.*hh[ord-1]*hh[ord]+3.*hh[ord]*hh[ord]+2.*hh[ord-1]*hh[ord+1]+3.*hh[ord]*hh[ord+1])/ - ((hh[ord-1]+hh[ord])*dd); - - rho_x = coefa*(xi-knot[ord])*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefb*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefc*(xi-knot[ord]); - rho_x += coefd; + coefa = 4.0f * (2.0f * hhm1 * hhm1 + 6.0f * hhm1 * hh0 + 3.0f * hh0 * hh0 + 3.0f * hhm1 * hhp1 + + 3.0f * hh0 * hhp1 + hhp1 * hhp1) / + (hh0 * (hhm1 + hh0) * (hh0 + hhp1) * dd); + coefb = -12.0f * (3.0f * hhm1 + 2.0f * hh0 + hhp1) / + ((hhm1 + hh0) * dd); + coefc = 12.0f * (-2.0f * hhm1 * hhm1 + hh0 * hh0 + hh0 * hhp1) / + ((hhm1 + hh0) * dd); + coefd = 4.0f * hhm1 * (4.0f * hhm1 * hh0 + 3.0f * hh0 * hh0 + 2.0f * hhm1 * hhp1 + 3.0f * hh0 * hhp1) / + ((hhm1 + hh0) * dd); + + dxi = xi - knot[ord]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord+1] && xi <= knot[ord+2]) { /* x2<=x<=x3 */ - dd *= (hh[ord]+hh[ord+1]); - coefa = -4.*(2.*hh[ord-1]+hh[ord])/(hh[ord+1]*dd); - coefb = 12.*(2.*hh[ord-1]+hh[ord])/dd; - coefc = -12.*(2.*hh[ord-1]+hh[ord])*hh[ord+1]/dd; - coefd = 4.*(2.*hh[ord-1]+hh[ord])*hh[ord+1]*hh[ord+1]/dd; - - rho_x = coefa*(xi-knot[ord+1])*(xi-knot[ord+1])*(xi-knot[ord+1]); - rho_x += coefb*(xi-knot[ord+1])*(xi-knot[ord+1]); - rho_x += coefc*(xi-knot[ord+1]); - rho_x += coefd; + dd *= (hh0 + hhp1); + coefa = -4.0f * (2.0f * hhm1 + hh0) / (hhp1 * dd); + coefb = 12.0f * (2.0f * hhm1 + hh0) / dd; + coefc = -12.0f * (2.0f * hhm1 + hh0) * hhp1 / dd; + coefd = 4.0f * (2.0f * hhm1 + hh0) * hhp1 * hhp1 / dd; + + dxi = xi - knot[ord+1]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else /* x>x3 */ - rho_x = 0.0; + rho_x = 0.0f; } else if (ord==Nx-1) { /* RHS-1 */ - float denom,denomsum,dd; - denom = hh[ord-2]*hh[ord-1]+hh[ord-1]*hh[ord-1]+2.*hh[ord-2]*hh[ord]+4.*hh[ord-1]*hh[ord]+3.*hh[ord]*hh[ord]; - denomsum = hh[ord-2]+hh[ord-1]+hh[ord]; - dd = denomsum*denom; + float hh0 = hh[ord]; + float hhm1 = hh[ord-1]; + float hhm2 = hh[ord-2]; + + float denom = hhm2 * hhm1 + hhm1 * hhm1 + 2.0f * hhm2 * hh0 + 4.0f * hhm1 * hh0 + 3.0f * hh0 * hh0; + float denomsum = hhm2 + hhm1 + hh0; + float dd = denomsum * denom; + if (xi >= knot[ord-2] && xi <= knot[ord-1]) { /* x0<=x<=x1 */ - coefa = 4.*(hh[ord-1]+2.*hh[ord])/(hh[ord-2]*(hh[ord-2]+hh[ord-1])*dd); - coefb = coefc = coefd = 0.0; + coefa = 4.0f * (hhm1 + 2.0f * hh0) / (hhm2 * (hhm2 + hhm1) * dd); + coefb = coefc = coefd = 0.0f; - rho_x = coefa*(xi-knot[ord-2])*(xi-knot[ord-2])*(xi-knot[ord-2]); - rho_x += coefb*(xi-knot[ord-2])*(xi-knot[ord-2]); - rho_x += coefc*(xi-knot[ord-2]); - rho_x += coefd; + dxi = xi - knot[ord-2]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord-1] && xi <= knot[ord]) { /* x1<=x<=x2 */ - coefa = -4.*(hh[ord-2]*hh[ord-2]+3.*hh[ord-2]*hh[ord-1]+3.*hh[ord-1]*hh[ord-1]+3.*hh[ord-2]*hh[ord]+ - 6.*hh[ord-1]*hh[ord]+2.*hh[ord]*hh[ord])/ - (hh[ord-1]*(hh[ord-2]+hh[ord-1])*(hh[ord-1]+hh[ord])*dd); - coefb = 12.*(hh[ord-1]+2.*hh[ord])/((hh[ord-2]+hh[ord-1])*dd); - coefc = 12.*hh[ord-2]*(hh[ord-1]+2.*hh[ord])/((hh[ord-2]+hh[ord-1])*dd); - coefd = 4.*hh[ord-2]*hh[ord-2]*(hh[ord-1]+2.*hh[ord])/((hh[ord-2]+hh[ord-1])*dd); - - rho_x = coefa*(xi-knot[ord-1])*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefb*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefc*(xi-knot[ord-1]); - rho_x += coefd; + coefa = -4.0f * (hhm2 * hhm2 + 3.0f * hhm2 * hhm1 + 3.0f * hhm1 * hhm1 + 3.0f * hhm2 * hh0 + + 6.0f * hhm1 * hh0 + 2.0f * hh0 * hh0) / + (hhm1 * (hhm2 + hhm1) * (hhm1 + hh0) * dd); + coefb = 12.0f * (hhm1 + 2.0f * hh0) / ((hhm2 + hhm1) * dd); + coefc = 12.0f * hhm2 * (hhm1 + 2.0f * hh0) / ((hhm2 + hhm1) * dd); + coefd = 4.0f * hhm2 * hhm2 * (hhm1 + 2.0f * hh0) / ((hhm2 + hhm1) * dd); + + dxi = xi - knot[ord-1]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord] && xi <= knot[ord+1]) { /* x2<=x<=x3 */ - dd *= (hh[ord-1]+hh[ord]); - coefa = 4.*(hh[ord-2]+2.*hh[ord-1]+3.*hh[ord])/(hh[ord]*dd); - coefb = -12.*(hh[ord-2]+2.*hh[ord-1]+3.*hh[ord])/dd; - coefc = 12.*(-hh[ord-2]*hh[ord-1]-hh[ord-1]*hh[ord-1]+2.*hh[ord]*hh[ord])/dd; - coefd = 4.*hh[ord]*(3.*hh[ord-2]*hh[ord-1]+3.*hh[ord-1]*hh[ord-1]+2.*hh[ord-2]*hh[ord]+4.*hh[ord-1]*hh[ord])/dd; - - rho_x = coefa*(xi-knot[ord])*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefb*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefc*(xi-knot[ord]); - rho_x += coefd; + dd *= (hhm1 + hh0); + coefa = 4.0f * (hhm2 + 2.0f * hhm1 + 3.0f * hh0) / (hh0 * dd); + coefb = -12.0f * (hhm2 + 2.0f * hhm1 + 3.0f * hh0) / dd; + coefc = 12.0f * (-hhm2 * hhm1 - hhm1 * hhm1 + 2.0f * hh0 * hh0) / dd; + coefd = 4.0f * hh0 * (3.0f * hhm2 * hhm1 + 3.0f * hhm1 * hhm1 + 2.0f * hhm2 * hh0 + 4.0f * hhm1 * hh0) / dd; + + dxi = xi - knot[ord]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else /* x>x4 */ - rho_x = 0.0; + rho_x = 0.0f; } else if (ord==Nx) { /* RHS */ - float denom; - denom = (hh[ord-2]+hh[ord-1])*(hh[ord-2]*hh[ord-2]+3.*hh[ord-2]*hh[ord-1]+3.*hh[ord-1]*hh[ord-1]); + float hhm1 = hh[ord-1]; + float hhm2 = hh[ord-2]; + + float denom = (hhm2 + hhm1) * (hhm2 * hhm2 + 3.0f * hhm2 * hhm1 + 3.0f * hhm1 * hhm1); + if (xi >= knot[ord-2] && xi <= knot[ord-1]) { /* x0<=x<=x1 */ - coefa = 4./(hh[ord-2]*denom); - coefb = coefc = coefd = 0.0; + coefa = 4.0f / (hhm2 * denom); + coefb = coefc = coefd = 0.0f; - rho_x = coefa*(xi-knot[ord-2])*(xi-knot[ord-2])*(xi-knot[ord-2]); - rho_x += coefb*(xi-knot[ord-2])*(xi-knot[ord-2]); - rho_x += coefc*(xi-knot[ord-2]); - rho_x += coefd; + dxi = xi - knot[ord-2]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord-1] && xi <= knot[ord]) { /* x1<=x<=x2 */ - coefa = -4./(hh[ord-1]*denom); - coefb = 12/denom; - coefc = 12*hh[ord-2]/denom; - coefd = 4.*hh[ord-2]*hh[ord-2]/denom; - - rho_x = coefa*(xi-knot[ord-1])*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefb*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefc*(xi-knot[ord-1]); - rho_x += coefd; + coefa = -4.0f / (hhm1 * denom); + coefb = 12.0f / denom; + coefc = 12.0f * hhm2 / denom; + coefd = 4.0f * hhm2 * hhm2 / denom; + + dxi = xi - knot[ord-1]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else /* x>x2 */ - rho_x = 0.0; + rho_x = 0.0f; } - else { /* Away from borders */ - float denom1,denom2,denom; - denom1 = hh[ord-2]+hh[ord-1]+hh[ord]+hh[ord+1]; + else { + /* Away from borders */ + float hh0 = hh[ord]; + float hhp1 = hh[ord+1]; + float hhm1 = hh[ord-1]; + float hhm2 = hh[ord-2]; + + float denom1 = hhm2 + hhm1 + hh0 + hhp1; + if (xi >= knot[ord-2] && xi <= knot[ord-1]) { /* x0<=x<=x1 */ - coefa = 4./(hh[ord-2]*(hh[ord-2]+hh[ord-1])*(hh[ord-2]+hh[ord-1]+hh[ord])*denom1); - coefb = coefc = coefd = 0.; + coefa = 4.0f / (hhm2 * (hhm2 + hhm1) * (hhm2 + hhm1 + hh0) * denom1); + coefb = coefc = coefd = 0.0f; - rho_x = coefa*(xi-knot[ord-2])*(xi-knot[ord-2])*(xi-knot[ord-2]); - rho_x += coefb*(xi-knot[ord-2])*(xi-knot[ord-2]); - rho_x += coefc*(xi-knot[ord-2]); - rho_x += coefd; + dxi = xi - knot[ord-2]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord-1] && xi <= knot[ord]) { /* x1<=x<=x2 */ - denom2 = (hh[ord-2]+hh[ord-1])*(hh[ord-2]+hh[ord-1]+hh[ord]); - denom = denom1*denom2; - - coefa = -4.*(hh[ord-2]*hh[ord-2]+3.*hh[ord-2]*hh[ord-1]+3.*hh[ord-1]*hh[ord-1]+2.*hh[ord-2]*hh[ord]+ - 4.*hh[ord-1]*hh[ord]+hh[ord]*hh[ord]+hh[ord-2]*hh[ord+1]+2.*hh[ord-1]*hh[ord+1]+hh[ord]*hh[ord+1])/ - (hh[ord-1]*(hh[ord-1]+hh[ord])*(hh[ord-1]+hh[ord]+hh[ord+1])*denom); - coefb = 12./denom; - coefc = 12.*hh[ord-2]/denom; - coefd = 4.*hh[ord-2]*hh[ord-2]/denom; - - rho_x = coefa*(xi-knot[ord-1])*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefb*(xi-knot[ord-1])*(xi-knot[ord-1]); - rho_x += coefc*(xi-knot[ord-1]); - rho_x += coefd; + float denom2 = (hhm2 + hhm1) * (hhm2 + hhm1 + hh0); + float denom = denom1 * denom2; + + coefa = -4.0f * (hhm2 * hhm2 + 3.0f * hhm2 * hhm1 + 3.0f * hhm1 * hhm1 + 2.0f * hhm2 * hh0 + + 4.0f * hhm1 * hh0 + hh0 * hh0 + hhm2 * hhp1 + 2.0f * hhm1 * hhp1 + hh0 * hhp1) / + (hhm1 * (hhm1 + hh0) * (hhm1 + hh0 + hhp1) * denom); + coefb = 12.0f / denom; + coefc = 12.0f * hhm2 / denom; + coefd = 4.0f * hhm2 * hhm2 / denom; + + dxi = xi - knot[ord-1]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord] && xi <= knot[ord+1]) { /* x2<=x<=x3 */ - denom2 = (hh[ord-1]+hh[ord])*(hh[ord-2]+hh[ord-1]+hh[ord])*(hh[ord-1]+hh[ord]+hh[ord+1]); - denom = denom1*denom2; - - coefa = 4.*(hh[ord-2]*hh[ord-1]+hh[ord-1]*hh[ord-1]+2.*hh[ord-2]*hh[ord]+4.*hh[ord-1]*hh[ord]+3.*hh[ord]*hh[ord]+ - hh[ord-2]*hh[ord+1]+2.*hh[ord-1]*hh[ord+1]+3.*hh[ord]*hh[ord+1]+hh[ord+1]*hh[ord+1])/ - (hh[ord]*(hh[ord]+hh[ord+1])*denom); - coefb = -12.*(hh[ord-2]+2.*hh[ord-1]+2.*hh[ord]+hh[ord+1])/denom; - coefc = 12.*(-hh[ord-2]*hh[ord-1]-hh[ord-1]*hh[ord-1]+hh[ord]*hh[ord]+hh[ord]*hh[ord+1])/denom; - coefd = 4.*(2.*hh[ord-2]*hh[ord-1]*hh[ord]+2.*hh[ord-1]*hh[ord-1]*hh[ord]+hh[ord-2]*hh[ord]*hh[ord]+ - 2.*hh[ord-1]*hh[ord]*hh[ord]+hh[ord-2]*hh[ord-1]*hh[ord+1]+hh[ord-1]*hh[ord-1]*hh[ord+1]+ - hh[ord-2]*hh[ord]*hh[ord+1]+2.*hh[ord-1]*hh[ord]*hh[ord+1])/denom; - - rho_x = coefa*(xi-knot[ord])*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefb*(xi-knot[ord])*(xi-knot[ord]); - rho_x += coefc*(xi-knot[ord]); - rho_x += coefd; + float denom2 = (hhm1 + hh0) * (hhm2 + hhm1 + hh0) * (hhm1 + hh0 + hhp1); + float denom = denom1 * denom2; + + coefa = 4.0f * (hhm2 * hhm1 + hhm1 * hhm1 + 2.0f * hhm2 * hh0 + 4.0f * hhm1 * hh0 + 3.0f * hh0 * hh0 + + hhm2 * hhp1 + 2.0f * hhm1 * hhp1 + 3.0f * hh0 * hhp1 + hhp1 * hhp1) / + (hh0 * (hh0 + hhp1) * denom); + coefb = -12.0f * (hhm2 + 2.0f * hhm1 + 2.0f * hh0 + hhp1) / denom; + coefc = 12.0f * (-hhm2 * hhm1 - hhm1 * hhm1 + hh0 * hh0 + hh0 * hhp1) / denom; + coefd = 4.0f * (2.0f * hhm2 * hhm1 * hh0 + 2.0f * hhm1 * hhm1 * hh0 + hhm2 * hh0 * hh0 + + 2.0f * hhm1 * hh0 * hh0 + hhm2 * hhm1 * hhp1 + hhm1 * hhm1 * hhp1 + + hhm2 * hh0 * hhp1 + 2.0f * hhm1 * hh0 * hhp1) / denom; + + dxi = xi - knot[ord]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else if (xi >= knot[ord+1] && xi <= knot[ord+2]) { /* x3<=x<=x4 */ - denom2 = (hh[ord]+hh[ord+1])*(hh[ord-1]+hh[ord]+hh[ord+1]); - denom = denom1*denom2; - - coefa = -4./(hh[ord+1]*denom); - coefb = 12/denom; - coefc = -12*hh[ord+1]/denom; - coefd = 4.*hh[ord+1]*hh[ord+1]/denom; - - rho_x = coefa*(xi-knot[ord+1])*(xi-knot[ord+1])*(xi-knot[ord+1]); - rho_x += coefb*(xi-knot[ord+1])*(xi-knot[ord+1]); - rho_x += coefc*(xi-knot[ord+1]); - rho_x += coefd; + float denom2 = (hh0 + hhp1) * (hhm1 + hh0 + hhp1); + float denom = denom1 * denom2; + + coefa = -4.0f / (hhp1 * denom); + coefb = 12.0f / denom; + coefc = -12.0f * hhp1 / denom; + coefd = 4.0f * hhp1 * hhp1 / denom; + + dxi = xi - knot[ord+1]; + rho_x = ((coefa * dxi + coefb) * dxi + coefc) * dxi + coefd; } else /* x>x4 */ - rho_x=0.0; + rho_x = 0.0f; } - free_farray1(hh,0); - return(rho_x); + + // frees temporary array + //free_farray1(hh,0); + + return (rho_x); } /* ----------------------------------------------------------------------------- */ @@ -348,11 +358,11 @@ void fill_hh_A3d(float *hh, float *knot, int Nx) float *farray1(int n11, int n12) { float *m; - m = (float *) malloc( U (n12-n11+1)*sizeof(float) ); + m = (float *) malloc( U (n12-n11+1) * sizeof(float) ); if(!m) { stop("allocation error in farray1"); } - return(m-n11); + return (m-n11); } /* ----------------------------------------------------------------------------- */ diff --git a/src/specfem3D/SIEM_prepare_solver.F90 b/src/specfem3D/SIEM_prepare_solver.F90 index d27871e20..25344f4f9 100644 --- a/src/specfem3D/SIEM_prepare_solver.F90 +++ b/src/specfem3D/SIEM_prepare_solver.F90 @@ -31,13 +31,27 @@ subroutine SIEM_prepare_solver() use specfem_par_full_gravity implicit none - + integer :: ier ! timing double precision :: tstart,tstart0,tCPU double precision, external :: wtime ! check if anything to do - if (.not. FULL_GRAVITY) return + if (.not. FULL_GRAVITY) then + ! allocate dummy arrays + ! (needed for function arguments even in case full gravity is turned off) + allocate(pgrav_cm(1),pgrav_ic(1),pgrav_oc(1),stat=ier) + if (ier /= 0) stop 'Error allocating dummy pgrav_cm,.. arrays' + pgrav_cm(:) = 0._CUSTOM_REAL; pgrav_ic(:) = 0._CUSTOM_REAL; pgrav_oc(:) = 0._CUSTOM_REAL + if (SIMULATION_TYPE == 3 ) then + allocate(b_pgrav_cm(1),b_pgrav_ic(1),b_pgrav_oc(1),stat=ier) + if (ier /= 0) stop 'Error allocating dummy b_pgrav_cm,.. arrays' + b_pgrav_cm(:) = 0._CUSTOM_REAL; b_pgrav_ic(:) = 0._CUSTOM_REAL; b_pgrav_oc(:) = 0._CUSTOM_REAL + endif + + ! nothing left to do + return + endif ! user output if (myrank == 0) then diff --git a/src/specfem3D/comp_source_time_function.f90 b/src/specfem3D/comp_source_time_function.f90 index f12180054..6f5c11e44 100644 --- a/src/specfem3D/comp_source_time_function.f90 +++ b/src/specfem3D/comp_source_time_function.f90 @@ -27,10 +27,10 @@ double precision function comp_source_time_function(t,hdur,it_index) - use constants, only: EXTERNAL_SOURCE_TIME_FUNCTION, & - STF_IS_UCB_HEAVISIDE + use constants, only: EXTERNAL_SOURCE_TIME_FUNCTION ! for berkeley source time function + use shared_parameters, only: STF_IS_UCB_HEAVISIDE use ucb_heaviside, only: comp_source_time_function_ucb_stf implicit none diff --git a/src/specfem3D/compute_element_att_memory.F90 b/src/specfem3D/compute_element_att_memory.F90 index e54a0a2d0..5ba2eacc4 100644 --- a/src/specfem3D/compute_element_att_memory.F90 +++ b/src/specfem3D/compute_element_att_memory.F90 @@ -515,7 +515,7 @@ end subroutine compute_element_att_memory_ic_lddrk ! helper functions ! ! -!daniel debug: att - debug update +! debug: att - debug update ! ! subroutine compute_element_att_mem_up_cm(ispec,i,j,k, & ! R_xx_loc,R_yy_loc,R_xy_loc,R_xz_loc,R_yz_loc, & @@ -524,7 +524,7 @@ end subroutine compute_element_att_memory_ic_lddrk !! update memory variables based upon the Runge-Kutta scheme ! ! -!!daniel: att - debug update +!! debug: att - debug update ! use specfem_par, only: tau_sigma_dble,deltat,b_deltat ! ! use specfem_par_crustmantle, only: factor_common => factor_common_crust_mantle diff --git a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 index bcdfd309a..85a13f822 100644 --- a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 +++ b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 @@ -386,7 +386,7 @@ subroutine compute_forces_viscoelastic_backward() ! checks if (SIMULATION_TYPE /= 3 ) return -!daniel debug: att - debug +! debug: att - debug ! integer :: iglob ! logical,parameter :: DEBUG = .false. ! if (DEBUG) then @@ -719,7 +719,7 @@ subroutine compute_forces_viscoelastic_backward() call update_veloc_elastic_newmark_backward() endif -!daniel debug: att - debug +! debug: att - debug ! if (DEBUG) then ! if (SIMULATION_TYPE == 1) then ! if (it > NSTEP - 1000 .and. myrank == 0) then diff --git a/src/specfem3D/get_attenuation.f90 b/src/specfem3D/get_attenuation.f90 index fa2c5a15e..5810b1eec 100644 --- a/src/specfem3D/get_attenuation.f90 +++ b/src/specfem3D/get_attenuation.f90 @@ -201,7 +201,10 @@ subroutine get_attenuation_property_values(tau_s, tau_e, factor_common, one_minu ! local parameters double precision, dimension(N_SLS) :: tauinv,beta - double precision :: omsb1,omsb2 + + ! debugging + !double precision :: omsb1,omsb2 + integer :: i tauinv(:) = 0.d0 @@ -223,20 +226,20 @@ subroutine get_attenuation_property_values(tau_s, tau_e, factor_common, one_minu ! M_U = M_R / [ 1 - sum( (tau_eps - tau_sigma)/tau_eps ) ] = M_R / [ 1 - sum(1 - tau_sigma/tau_eps) ] ! with a scaling factor = 1 / [ 1 - sum(1 - tau_sigma/tau_eps) ] ! - ! to compare these two factors: - if (.false.) then - ! factor to scale from relaxed to unrelaxed moduli: see Komatitsch - omsb1 = 1.0d0 - do i = 1,N_SLS - omsb1 = omsb1 - (1.d0 - tau_e(i)/tau_s(i)) - enddo - ! factor to scale from relaxed to unrelaxed moduli: see Liu - omsb2 = 1.0d0 - do i = 1,N_SLS - omsb2 = omsb2 - (1.d0 - tau_s(i)/tau_e(i)) - enddo - print *,'debug: one_minus_sum_beta1,2 = ',omsb1,1.d0/omsb2 - endif + ! debugging - to compare these two factors: + !if (.false.) then + ! ! factor to scale from relaxed to unrelaxed moduli: see Komatitsch + ! omsb1 = 1.0d0 + ! do i = 1,N_SLS + ! omsb1 = omsb1 - (1.d0 - tau_e(i)/tau_s(i)) + ! enddo + ! ! factor to scale from relaxed to unrelaxed moduli: see Liu + ! omsb2 = 1.0d0 + ! do i = 1,N_SLS + ! omsb2 = omsb2 - (1.d0 - tau_s(i)/tau_e(i)) + ! enddo + ! print *,'debug: one_minus_sum_beta1,2 = ',omsb1,1.d0/omsb2 + !endif ! results: ! Q ~ 65: one_minus_sum_beta1,2 = 1.0648735275165400 1.0677815840121743 -> ratio = 0.997276 ! Q ~ 184: one_minus_sum_beta1,2 = 1.0233224020801992 1.0236877302481358 -> ratio = 0.999643 diff --git a/src/specfem3D/initialize_simulation.F90 b/src/specfem3D/initialize_simulation.F90 index 90cd641e5..7f45c3a72 100644 --- a/src/specfem3D/initialize_simulation.F90 +++ b/src/specfem3D/initialize_simulation.F90 @@ -168,6 +168,12 @@ subroutine initialize_simulation() tiny(1._CUSTOM_REAL),huge(1._CUSTOM_REAL) write(IMAIN,*) + ! backward compatibility + if (USE_OLD_VERSION_FORMAT) then + write(IMAIN,*) 'using old version backward compatibility (versions 7.0 to 8.0)' + write(IMAIN,*) + endif + ! model user parameters write(IMAIN,*) 'model: ',trim(MODEL) if (OCEANS) then diff --git a/src/specfem3D/iterate_time_undoatt.F90 b/src/specfem3D/iterate_time_undoatt.F90 index 5c23708ba..40609a2d6 100644 --- a/src/specfem3D/iterate_time_undoatt.F90 +++ b/src/specfem3D/iterate_time_undoatt.F90 @@ -473,7 +473,7 @@ subroutine iterate_time_undoatt() endif call register_host_array(NDIM*NGLOB_CRUST_MANTLE_ADJOINT, b_displ_cm_store_buffer(:,:, it_of_buffer)) #endif - ! daniel debug: check if these transfers could be made async to overlap + ! daniel TODO: check if these transfers could be made async to overlap call transfer_ofs_b_displ_cm_from_device(NDIM*NGLOB_CRUST_MANTLE_ADJOINT,it_of_buffer, & b_displ_cm_store_buffer,Mesh_pointer) call transfer_ofs_b_displ_ic_from_device(NDIM*NGLOB_INNER_CORE_ADJOINT,it_of_buffer, & @@ -521,7 +521,7 @@ subroutine iterate_time_undoatt() ! crust/mantle ! transfers wavefields from CPU to GPU if (GPU_MODE) then - ! daniel debug: check if these transfers could be made async to overlap + ! daniel TODO: check if these transfers could be made async to overlap call transfer_ofs_b_displ_cm_to_device(NDIM*NGLOB_CRUST_MANTLE_ADJOINT,it_of_buffer, & b_displ_cm_store_buffer,Mesh_pointer) call transfer_ofs_b_displ_ic_to_device(NDIM*NGLOB_INNER_CORE_ADJOINT,it_of_buffer, & diff --git a/src/specfem3D/locate_sources.f90 b/src/specfem3D/locate_sources.f90 index 8d7ce18f3..bfd1151b8 100644 --- a/src/specfem3D/locate_sources.f90 +++ b/src/specfem3D/locate_sources.f90 @@ -60,7 +60,7 @@ subroutine locate_sources() use specfem_par_movie, only: vtkdata_source_x,vtkdata_source_y,vtkdata_source_z ! for Berkeley stf - use shared_parameters, only: UCB_SOURCE_T1,UCB_SOURCE_T2,UCB_SOURCE_T3,UCB_SOURCE_T4,UCB_TAU + use shared_parameters, only: STF_IS_UCB_HEAVISIDE,UCB_SOURCE_T1,UCB_SOURCE_T2,UCB_SOURCE_T3,UCB_SOURCE_T4,UCB_TAU implicit none diff --git a/src/specfem3D/make_gravity.f90 b/src/specfem3D/make_gravity.f90 index 814758aa2..11f9dd592 100644 --- a/src/specfem3D/make_gravity.f90 +++ b/src/specfem3D/make_gravity.f90 @@ -30,7 +30,7 @@ subroutine make_gravity(nspl,rspl,gravity_spline,gravity_spline2) ! creates a spline for the gravity profile in PREM ! radius and density are non-dimensional - use constants, only: NR_DENSITY,myrank + use constants, only: NR_DENSITY,myrank,USE_OLD_VERSION_FORMAT use shared_parameters, only: PLANET_TYPE,IPLANET_EARTH,IPLANET_MARS,IPLANET_MOON,R_PLANET,R_PLANET_KM ! reference models @@ -99,7 +99,12 @@ subroutine make_gravity(nspl,rspl,gravity_spline,gravity_spline2) R771 = PREM_R771 RTOPDDOUBLEPRIME = PREM_RTOPDDOUBLEPRIME RCMB = PREM_RCMB - RICB = PREM_RICB + ! version compatibility + if (USE_OLD_VERSION_FORMAT) then + RICB = PREM_RICB_OLD + else + RICB = PREM_RICB + endif case (IPLANET_MARS) ! Mars diff --git a/src/specfem3D/prepare_gravity.f90 b/src/specfem3D/prepare_gravity.f90 index 1c73201c4..848c48260 100644 --- a/src/specfem3D/prepare_gravity.f90 +++ b/src/specfem3D/prepare_gravity.f90 @@ -130,11 +130,15 @@ subroutine prepare_gravity() ! non-dimensionalizes sampling point locations do int_radius = 1,NRAD_GRAVITY - ! old: assuming R_PLANET_KM = R_EARTH_KM = 6371.d0, ranges from 1/10 * 1/R_PLANET_KM ~ 1.e-5 to 70000/10 * 1/R_PLANET ~ 1.09 - ! r(int_radius) = dble(int_radius) / (R_PLANET_KM * 10.d0) - ! - ! setting sampling points up to range_max - r(int_radius) = dble(int_radius) / dble(NRAD_GRAVITY) * range_max + ! sets radius value + if (USE_OLD_VERSION_FORMAT) then + ! old version increments in 1/10 equals steps of 100 m (see constants.h: NRAD_GRAVITY = 70000) + ! old: assuming R_PLANET_KM = R_EARTH_KM = 6371.d0, ranges from 1/10 * 1/R_PLANET_KM ~ 1.e-5 to 70000/10 * 1/R_PLANET ~ 1.09 + r(int_radius) = dble(int_radius) / (R_PLANET_KM * 10.d0) + else + ! setting sampling points up to range_max + r(int_radius) = dble(int_radius) / dble(NRAD_GRAVITY) * range_max + endif ! r in range [1.4e-5,1.001] ! @@ -226,33 +230,56 @@ subroutine prepare_gravity() endif enddo + ! debug - file output + if (.false.) then + open(21,file='tmp_grav.dat',status='unknown') + write(21,*) "# int_radius # radius # rho # minus_g # minus_dg" + do int_radius = 1,NRAD_GRAVITY + radius = r(int_radius) + rho = density_table(int_radius) + minus_g = minus_gravity_table(int_radius) + minus_dg = minus_deriv_gravity_table(int_radius) + write(21,*) int_radius,radius,rho,minus_g,minus_dg + enddo + close(21) + endif + ! make sure fluid array is only assigned in outer core between 1222 and 3478 km ! lookup table is defined every 100 m do int_radius = 1,NRAD_GRAVITY - radius_km = r(int_radius) * R_PLANET_KM + ! radius in km + if (USE_OLD_VERSION_FORMAT) then + radius_km = dble(int_radius) / 10.d0 + else + radius_km = r(int_radius) * R_PLANET_KM + endif + ! upper limit for fluid core if (radius_km > RCMB/1000.d0 - 3.d0) then ! gets index for fluid - ! - ! old: based on int_radius = r * R_PLANET/1000 * 10.d0 = (RCMB/1000 - 3)/R_PLANET_KM * R_PLANET/1000 * 10.d0 - ! = (RCMB/1000 - 3) * 10.d0 - !index_fluid = nint((RCMB/1000.d0 - 3.d0)*10.d0) - ! - ! new: based on int_radius = r * NRAD_GRAVITY / range_max (with r being non-dimensionalized between [0,1.001]) - index_fluid = nint( (RCMB/1000.d0 - 3.d0)/R_PLANET_KM * NRAD_GRAVITY / range_max) + if (USE_OLD_VERSION_FORMAT) then + ! old: based on int_radius = r * R_PLANET/1000 * 10.d0 = (RCMB/1000 - 3)/R_PLANET_KM * R_PLANET/1000 * 10.d0 + ! = (RCMB/1000 - 3) * 10.d0 + index_fluid = nint((RCMB/1000.d0 - 3.d0)*10.d0) + else + ! new: based on int_radius = r * NRAD_GRAVITY / range_max (with r being non-dimensionalized between [0,1.001]) + index_fluid = nint( (RCMB/1000.d0 - 3.d0)/R_PLANET_KM * NRAD_GRAVITY / range_max) + endif ! stays with fluid properties minus_rho_g_over_kappa_fluid(int_radius) = minus_rho_g_over_kappa_fluid(index_fluid) endif + ! lower limit for fluid core if (radius_km < RICB/1000.d0 + 3.d0) then ! gets index for fluid - ! - ! old: based on int_radius = r * R_PLANET/1000 * 10.d0 = (RICB/1000 + 3)/R_PLANET_KM * R_PLANET/1000 * 10.d0 - ! = (RICB/1000 + 3) * 10.d0 - !index_fluid = nint((RICB/1000.d0 + 3.d0)*10.d0) - ! - ! new: based on int_radius = r * NRAD_GRAVITY / range_max (with r being non-dimensionalized between [0,1.001]) - index_fluid = nint( (RICB/1000.d0 + 3.d0)/R_PLANET_KM * NRAD_GRAVITY / range_max) + if (USE_OLD_VERSION_FORMAT) then + ! old: based on int_radius = r * R_PLANET/1000 * 10.d0 = (RICB/1000 + 3)/R_PLANET_KM * R_PLANET/1000 * 10.d0 + ! = (RICB/1000 + 3) * 10.d0 + index_fluid = nint((RICB/1000.d0 + 3.d0)*10.d0) + else + ! new: based on int_radius = r * NRAD_GRAVITY / range_max (with r being non-dimensionalized between [0,1.001]) + index_fluid = nint( (RICB/1000.d0 + 3.d0)/R_PLANET_KM * NRAD_GRAVITY / range_max) + endif ! stays with fluid properties minus_rho_g_over_kappa_fluid(int_radius) = minus_rho_g_over_kappa_fluid(index_fluid) endif @@ -286,13 +313,19 @@ subroutine prepare_gravity() call revert_ellipticity_rtheta(r_table,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2) ! integrated and multiply by rho / Kappa - ! old: int_radius = nint(10.d0 * radius * R_PLANET_KM) - ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) - ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max - int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + if (USE_OLD_VERSION_FORMAT) then + ! old: int_radius = nint(10.d0 * radius * R_PLANET_KM) + ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) + int_radius = nint(radius * R_PLANET_KM * 10.d0) + else + ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max + int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + endif + ! limits range if (int_radius < 1) int_radius = 1 if (int_radius > NRAD_GRAVITY) int_radius = NRAD_GRAVITY + ! debug if (DEBUG) then if (myrank == 0) print *,'debug: prepare gravity radius = ',radius,r_table,'r = ',r(int_radius),int_radius @@ -301,9 +334,9 @@ subroutine prepare_gravity() fac = minus_rho_g_over_kappa_fluid(int_radius) ! gravitational acceleration (integrated and multiply by rho / Kappa) in Cartesian coordinates - gxl = dsin(theta) * dcos(phi) * fac - gyl = dsin(theta) * dsin(phi) * fac - gzl = dcos(theta) * fac + gxl = (dsin(theta) * dcos(phi)) * fac + gyl = (dsin(theta) * dsin(phi)) * fac + gzl = (dcos(theta)) * fac gravity_pre_store_outer_core(1,iglob) = real(gxl,kind=CUSTOM_REAL) gravity_pre_store_outer_core(2,iglob) = real(gyl,kind=CUSTOM_REAL) gravity_pre_store_outer_core(3,iglob) = real(gzl,kind=CUSTOM_REAL) @@ -396,13 +429,19 @@ subroutine prepare_gravity() call revert_ellipticity_rtheta(r_table,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2) ! radius index - ! old: int_radius = nint(10.d0 * radius * R_PLANET_KM) - ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) - ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max - int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + if (USE_OLD_VERSION_FORMAT) then + ! old: int_radius = nint(10.d0 * radius * R_PLANET_KM) + ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) + int_radius = nint(radius * R_PLANET_KM * 10.d0) + else + ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max + int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + endif + ! limits range if (int_radius < 1) int_radius = 1 if (int_radius > NRAD_GRAVITY) int_radius = NRAD_GRAVITY + ! debug if (DEBUG) then if (myrank == 0) print *,'debug: prepare no gravity radius = ',radius,r_table,'r = ',r(int_radius),int_radius @@ -412,9 +451,9 @@ subroutine prepare_gravity() fac = d_ln_density_dr_table(int_radius) ! gradient of d ln(rho)/dr in Cartesian coordinates - gxl = dsin(theta) * dcos(phi) * fac - gyl = dsin(theta) * dsin(phi) * fac - gzl = dcos(theta) * fac + gxl = (dsin(theta) * dcos(phi)) * fac + gyl = (dsin(theta) * dsin(phi)) * fac + gzl = (dcos(theta)) * fac gravity_pre_store_outer_core(1,iglob) = real(gxl,kind=CUSTOM_REAL) gravity_pre_store_outer_core(2,iglob) = real(gyl,kind=CUSTOM_REAL) gravity_pre_store_outer_core(3,iglob) = real(gzl,kind=CUSTOM_REAL) @@ -502,10 +541,15 @@ subroutine prepare_gravity() call revert_ellipticity_rtheta(r_table,theta,nspl_ellip,rspl_ellip,ellipicity_spline,ellipicity_spline2) ! for efficiency replace with lookup table every 100 m in radial direction - ! old: int_radius = nint(10.d0 * radius * R_PLANET_KM) - ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) - ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max - int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + if (USE_OLD_VERSION_FORMAT) then + ! old: int_radius = nint(10.d0 * radius * R_PLANET_KM) + ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) + int_radius = nint(radius * R_PLANET_KM * 10.d0) + else + ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max + int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + endif + ! limits range if (int_radius < 1) int_radius = 1 if (int_radius > NRAD_GRAVITY) int_radius = NRAD_GRAVITY @@ -516,9 +560,9 @@ subroutine prepare_gravity() ! Cartesian components of the gravitational acceleration ! multiplied by common factor rho - gxl = minus_g * sin_theta * cos_phi * rho - gyl = minus_g * sin_theta * sin_phi * rho - gzl = minus_g * cos_theta * rho + gxl = (minus_g * sin_theta * cos_phi) * rho + gyl = (minus_g * sin_theta * sin_phi) * rho + gzl = (minus_g * cos_theta) * rho gravity_pre_store_crust_mantle(1,iglob) = real(gxl,kind=CUSTOM_REAL) gravity_pre_store_crust_mantle(2,iglob) = real(gyl,kind=CUSTOM_REAL) gravity_pre_store_crust_mantle(3,iglob) = real(gzl,kind=CUSTOM_REAL) @@ -644,10 +688,15 @@ subroutine prepare_gravity() ! for efficiency replace with lookup table every 100 m in radial direction ! make sure we never use zero for point exactly at the center of the Earth - ! old: int_radius = max(1,nint(10.d0 * radius * R_EARTH_KM)) - ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) - ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max - int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + if (USE_OLD_VERSION_FORMAT) then + ! old: int_radius = max(1,nint(10.d0 * radius * R_EARTH_KM)) + ! based on radius = dble(int_radius) / (R_PLANET_KM * 10.d0) + int_radius = nint(radius * R_EARTH_KM * 10.d0) + else + ! new: radius = dble(int_radius) / dble(NRAD_GRAVITY) * range_max + int_radius = nint( r_table / range_max * dble(NRAD_GRAVITY) ) + endif + ! limits range if (int_radius < 1) int_radius = 1 if (int_radius > NRAD_GRAVITY) int_radius = NRAD_GRAVITY @@ -658,9 +707,9 @@ subroutine prepare_gravity() ! Cartesian components of the gravitational acceleration ! multiplied by common factor rho - gxl = minus_g * sin_theta * cos_phi * rho - gyl = minus_g * sin_theta * sin_phi * rho - gzl = minus_g * cos_theta * rho + gxl = (minus_g * sin_theta * cos_phi) * rho + gyl = (minus_g * sin_theta * sin_phi) * rho + gzl = (minus_g * cos_theta) * rho gravity_pre_store_inner_core(1,iglob) = real(gxl,kind=CUSTOM_REAL) gravity_pre_store_inner_core(2,iglob) = real(gyl,kind=CUSTOM_REAL) diff --git a/src/specfem3D/prepare_oceans.f90 b/src/specfem3D/prepare_oceans.f90 index 46cdc6ca7..4a02721e3 100644 --- a/src/specfem3D/prepare_oceans.f90 +++ b/src/specfem3D/prepare_oceans.f90 @@ -106,22 +106,30 @@ subroutine prepare_oceans() npoin_oceans = ipoin ! Assemble normals - call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, & - normx, & - num_interfaces_crust_mantle,max_nibool_interfaces_cm, & - nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle, & - my_neighbors_crust_mantle) - call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, & - normy, & - num_interfaces_crust_mantle,max_nibool_interfaces_cm, & - nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle, & - my_neighbors_crust_mantle) - call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, & - normz, & - num_interfaces_crust_mantle,max_nibool_interfaces_cm, & - nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle, & - my_neighbors_crust_mantle) - + ! note: having multiple MPI slices, the normal on a halo point that is shared between MPI slices might be slightly different + ! depending from which slice it is taken. to avoid such a jump, here we count its valence and average the normal + ! across different MPI slices. + if (USE_OLD_VERSION_FORMAT) then + ! no assembly between different slices done + continue + else + ! assembles normal values from halo points + call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, & + normx, & + num_interfaces_crust_mantle,max_nibool_interfaces_cm, & + nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle, & + my_neighbors_crust_mantle) + call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, & + normy, & + num_interfaces_crust_mantle,max_nibool_interfaces_cm, & + nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle, & + my_neighbors_crust_mantle) + call assemble_MPI_scalar(NPROCTOT_VAL,NGLOB_CRUST_MANTLE, & + normz, & + num_interfaces_crust_mantle,max_nibool_interfaces_cm, & + nibool_interfaces_crust_mantle,ibool_interfaces_crust_mantle, & + my_neighbors_crust_mantle) + endif ! user info if (myrank == 0) then @@ -150,21 +158,35 @@ subroutine prepare_oceans() do i = 1,NGLLX ! get global point number iglob = ibool_crust_mantle(i,j,k,ispec) + ! updates once if (.not. updated_dof_ocean_load(iglob)) then ipoin = ipoin + 1 + ! fills arrays ibool_ocean_load(ipoin) = iglob rmass_ocean_load_selected(ipoin) = rmass_ocean_load(iglob) - ! Make normal continuous - if (valence(iglob) > 1.) then - norm = sqrt(normx(iglob)**2 + normy(iglob)**2 + normz(iglob)**2) - normal_ocean_load(1,ipoin) = normx(iglob) / norm - normal_ocean_load(2,ipoin) = normy(iglob) / norm - normal_ocean_load(3,ipoin) = normz(iglob) / norm + + ! normals + if (USE_OLD_VERSION_FORMAT) then + ! old version compatibility + ! uses single point normal + ! note: this might lead to small jumps across MPI slices, i.e., + ! having different normals on shared halo points depending from which slice it is taken + normal_ocean_load(:,ipoin) = normal_top_crust_mantle(:,i,j,ispec2D) else - normal_ocean_load(:,ipoin) = normal_top_crust_mantle(:,i,j,ispec2D) + ! considers assembled normals on halo points + ! Make normal continuous + if (valence(iglob) > 1.) then + norm = sqrt(normx(iglob)**2 + normy(iglob)**2 + normz(iglob)**2) + normal_ocean_load(1,ipoin) = normx(iglob) / norm + normal_ocean_load(2,ipoin) = normy(iglob) / norm + normal_ocean_load(3,ipoin) = normz(iglob) / norm + else + normal_ocean_load(:,ipoin) = normal_top_crust_mantle(:,i,j,ispec2D) + endif endif + ! masks this global point updated_dof_ocean_load(iglob) = .true. endif diff --git a/src/specfem3D/prepare_timerun.F90 b/src/specfem3D/prepare_timerun.F90 index d09df9384..78adbdc34 100644 --- a/src/specfem3D/prepare_timerun.F90 +++ b/src/specfem3D/prepare_timerun.F90 @@ -64,7 +64,7 @@ subroutine prepare_timerun() call prepare_gravity() ! full gravity preparation - if (FULL_GRAVITY) call SIEM_prepare_solver() + call SIEM_prepare_solver() ! precomputes attenuation factors call prepare_attenuation() diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90 index 295bcef40..9af62a78b 100644 --- a/src/specfem3D/setup_sources_receivers.f90 +++ b/src/specfem3D/setup_sources_receivers.f90 @@ -582,6 +582,7 @@ subroutine setup_sources() use specfem_par_movie ! for berkeley ucb stf + !use shared_parameters, only: STF_IS_UCB_HEAVISIDE use ucb_heaviside, only: init_ucb_heaviside implicit none diff --git a/src/specfem3D/specfem3D_par.F90 b/src/specfem3D/specfem3D_par.F90 index 9bbadf92e..77cadf84c 100644 --- a/src/specfem3D/specfem3D_par.F90 +++ b/src/specfem3D/specfem3D_par.F90 @@ -33,7 +33,7 @@ module constants_solver implicit none -! daniel debug: todo + #ifdef USE_STATIC_COMPILATION ! static compilation diff --git a/src/specfem3D/ucb_heaviside.f90 b/src/specfem3D/ucb_heaviside.f90 index 0ff9a80ba..3757e771e 100644 --- a/src/specfem3D/ucb_heaviside.f90 +++ b/src/specfem3D/ucb_heaviside.f90 @@ -90,6 +90,7 @@ end function comp_source_time_function_ucb_stf subroutine init_ucb_heaviside(nstep_,dt_) use constants, only: myrank + use shared_parameters, only: STF_IS_UCB_HEAVISIDE implicit none @@ -107,6 +108,9 @@ subroutine init_ucb_heaviside(nstep_,dt_) double precision :: pi = 4.d0 * datan(1.d0) + ! check if anything to do + if (.not. STF_IS_UCB_HEAVISIDE) return + ! frequencies f1h = 1.d0 / UCB_SOURCE_T1 f2h = 1.d0 / UCB_SOURCE_T2 @@ -118,11 +122,12 @@ subroutine init_ucb_heaviside(nstep_,dt_) ! user info if (myrank == 0) then write(*,*) "*****************************************" - write(*,*) "Heaviside Source T1:", UCB_SOURCE_T1 - write(*,*) "Heaviside Source T2:", UCB_SOURCE_T2 - write(*,*) "Heaviside Source T3:", UCB_SOURCE_T3 - write(*,*) "Heaviside Source T4:", UCB_SOURCE_T4 - write(*,*) "Source time-shift:", t0 + write(*,*) "STF is UCB Heaviside:" + write(*,*) " Heaviside Source T1:", UCB_SOURCE_T1 + write(*,*) " Heaviside Source T2:", UCB_SOURCE_T2 + write(*,*) " Heaviside Source T3:", UCB_SOURCE_T3 + write(*,*) " Heaviside Source T4:", UCB_SOURCE_T4 + write(*,*) " Source time-shift:", t0 write(*,*) "*****************************************" endif