From 19997ba027a0fa6a3b5e97eaaa6bcdcf74081e2c Mon Sep 17 00:00:00 2001 From: ibanos90 Date: Wed, 18 Feb 2026 10:32:21 -0700 Subject: [PATCH 1/5] Convert ERA5 model level analyses to MPAS initial conditions --- bin/ExternalAnalysisToMPAS.csh | 7 + bin/GetERA5AnalysisFromGDEX.csh | 93 +++ config/mpas/initic/ecmwf_coeffs | 138 +++++ config/tools.csh | 1 + initialize/data/ExternalAnalyses.py | 34 +- scenarios/GenerateERA5Analyses.yaml | 18 + scenarios/defaults/externalanalyses.yaml | 12 + tools/WPSUtils.py | 139 +++++ tools/era5_to_int.py | 725 +++++++++++++++++++++++ tools/fortran_io.py | 137 +++++ 10 files changed, 1298 insertions(+), 6 deletions(-) create mode 100755 bin/GetERA5AnalysisFromGDEX.csh create mode 100644 config/mpas/initic/ecmwf_coeffs create mode 100644 scenarios/GenerateERA5Analyses.yaml create mode 100644 tools/WPSUtils.py create mode 100755 tools/era5_to_int.py create mode 100644 tools/fortran_io.py diff --git a/bin/ExternalAnalysisToMPAS.csh b/bin/ExternalAnalysisToMPAS.csh index 89355899..1e8c9805 100755 --- a/bin/ExternalAnalysisToMPAS.csh +++ b/bin/ExternalAnalysisToMPAS.csh @@ -127,6 +127,13 @@ sed -i 's@nCells@'${ArgNCells}'@' $NamelistFileInit sed -i 's@{{meshRatio}}@'${ArgRatio}'@' $NamelistFileInit sed -i 's@{{UngribPrefix}}@'${externalanalyses__UngribPrefix}'@' $NamelistFileInit +# ERA5 adjustments +if ( "${externalanalyses__UngribPrefix}" == "ERA5" ) then + sed -i "s@config_met_prefix = .*@config_met_prefix = 'ERA5'@g" $NamelistFileInit + sed -i "s@config_nfglevels = .*@config_nfglevels = 138@g" $NamelistFileInit + sed -i "s@config_use_spechumd = .*@config_use_spechumd = true@g" $NamelistFileInit +endif + # Run the executable # ================== rm ./${InitEXE} diff --git a/bin/GetERA5AnalysisFromGDEX.csh b/bin/GetERA5AnalysisFromGDEX.csh new file mode 100755 index 00000000..5af0753b --- /dev/null +++ b/bin/GetERA5AnalysisFromGDEX.csh @@ -0,0 +1,93 @@ +#!/bin/csh -f + +# (C) Copyright 2023 UCAR +# +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + +# Get GFS analysis (0-h forecast) for cold start initial conditions + +# Process arguments +# ================= +## args +# ArgDT: int, valid time offset beyond CYLC_TASK_CYCLE_POINT in hours +set ArgDT = "$1" + +# ArgWorkDir: my location +set ArgWorkDir = "$2" + +set test = `echo $ArgDT | grep '^[0-9]*$'` +set isNotInt = ($status) +if ( $isNotInt ) then + echo "ERROR in $0 : ArgDT must be an integer, not $ArgDT" + exit 1 +endif + +date + +# Setup environment +# ================= +source config/auto/build.csh +source config/auto/experiment.csh +source config/tools.csh +set ccyymmdd = `echo ${CYLC_TASK_CYCLE_POINT} | cut -c 1-8` +set hh = `echo ${CYLC_TASK_CYCLE_POINT} | cut -c 10-11` +set thisCycleDate = ${ccyymmdd}${hh} +set thisValidDate = `$advanceCYMDH ${thisCycleDate} ${ArgDT}` +set prevValidDate = `$advanceCYMDH ${thisCycleDate} -6` + +source ./bin/getCycleVars.csh + +set cyyyy = `echo ${thisValidDate} | cut -c1-4` +set cmm = `echo ${thisValidDate} | cut -c5-6` +set cdd = `echo ${thisValidDate} | cut -c7-8` +set chh = `echo ${thisValidDate} | cut -c9-10` +set datestr = "${cyyyy}-${cmm}-${cdd}_${chh}" + +set WorkDir = ${ExperimentDirectory}/`echo "$ArgWorkDir" \ + | sed 's@{{thisValidDate}}@'${thisValidDate}'@' \ + ` +echo "WorkDir = ${WorkDir}" +mkdir -p ${WorkDir} +cd ${WorkDir} + +# ================================================================================================ + +if ( -e GETSUCCESS ) then + echo "$0 (INFO): GETSUCCESS file already exists, exiting with success" + echo "$0 (INFO): if regenerating the output files is desired, delete GETSUCCESS" + + date + + exit 0 +endif + +# ================================================================================================ + +echo "Getting ERA5 analysis from GDEX" + + +# Link ECMWF coefficients +ln -sf ${ModelConfigDir}/initic/ecmwf_coeffs ./ecmwf_coeffs + +# Convert ERA5 NetCDF to WPS intermediate +setenv myCommand `$era5_to_int ${datestr}` +echo "$myCommand" +${myCommand} + +# Build ERA5 intermediate filename (YYYY-MM-DD_HH) +set era5file = ERA5:${datestr} + +# Check that ERA5 intermediate file exists +if ( ! -e ${era5file} ) then + echo "${era5file} not found -- exiting" + exit 1 +endif + +echo "ERA5 intermediate file created successfully: ${era5file}" + +date + +touch GETSUCCESS + +exit 0 diff --git a/config/mpas/initic/ecmwf_coeffs b/config/mpas/initic/ecmwf_coeffs new file mode 100644 index 00000000..236ec06e --- /dev/null +++ b/config/mpas/initic/ecmwf_coeffs @@ -0,0 +1,138 @@ +0 0.000000 0.000000 0.000000 - - - - - +1 2.000365 0.000000 0.0200 0.0100 79301.79 80301.65 198.05 0.000018 +2 3.102241 0.000000 0.0310 0.0255 73721.58 74584.91 209.21 0.000042 +3 4.666084 0.000000 0.0467 0.0388 71115.75 71918.79 214.42 0.000063 +4 6.827977 0.000000 0.0683 0.0575 68618.43 69365.77 221.32 0.000090 +5 9.746966 0.000000 0.0975 0.0829 66210.99 66906.53 228.06 0.000127 +6 13.605424 0.000000 0.1361 0.1168 63890.03 64537.43 234.56 0.000173 +7 18.608931 0.000000 0.1861 0.1611 61651.77 62254.39 240.83 0.000233 +8 24.985718 0.000000 0.2499 0.2180 59492.50 60053.46 246.87 0.000308 +9 32.985710 0.000000 0.3299 0.2899 57408.61 57930.78 252.71 0.000400 +10 42.879242 0.000000 0.4288 0.3793 55396.62 55882.68 258.34 0.000512 +11 54.955463 0.000000 0.5496 0.4892 53453.20 53905.62 263.78 0.000646 +12 69.520576 0.000000 0.6952 0.6224 51575.15 51996.21 269.04 0.000806 +13 86.895882 0.000000 0.8690 0.7821 49767.41 50159.36 270.65 0.001007 +14 107.415741 0.000000 1.0742 0.9716 48048.70 48413.94 270.65 0.001251 +15 131.425507 0.000000 1.3143 1.1942 46416.22 46756.98 269.02 0.001546 +16 159.279404 0.000000 1.5928 1.4535 44881.17 45199.69 264.72 0.001913 +17 191.338562 0.000000 1.9134 1.7531 43440.23 43738.55 260.68 0.002343 +18 227.968948 0.000000 2.2797 2.0965 42085.00 42364.93 256.89 0.002843 +19 269.539581 0.000000 2.6954 2.4875 40808.05 41071.20 253.31 0.003421 +20 316.420746 0.000000 3.1642 2.9298 39602.76 39850.56 249.94 0.004084 +21 368.982361 0.000000 3.6898 3.4270 38463.25 38696.94 246.75 0.004838 +22 427.592499 0.000000 4.2759 3.9829 37384.22 37604.95 243.73 0.005693 +23 492.616028 0.000000 4.9262 4.6010 36360.94 36569.72 240.86 0.006655 +24 564.413452 0.000000 5.6441 5.2851 35389.15 35586.89 238.14 0.007731 +25 643.339905 0.000000 6.4334 6.0388 34465.00 34652.52 235.55 0.008931 +26 729.744141 0.000000 7.2974 6.8654 33585.02 33763.05 233.09 0.010261 +27 823.967834 0.000000 8.2397 7.7686 32746.04 32915.27 230.74 0.011729 +28 926.344910 0.000000 9.2634 8.7516 31945.53 32106.57 228.60 0.013337 +29 1037.201172 0.000000 10.3720 9.8177 31177.59 31330.96 227.83 0.015012 +30 1156.853638 0.000000 11.5685 10.9703 30438.54 30584.71 227.09 0.016829 +31 1285.610352 0.000000 12.8561 12.2123 29726.69 29866.09 226.38 0.018793 +32 1423.770142 0.000000 14.2377 13.5469 29040.48 29173.50 225.69 0.020910 +33 1571.622925 0.000000 15.7162 14.9770 28378.46 28505.47 225.03 0.023186 +34 1729.448975 0.000000 17.2945 16.5054 27739.29 27860.64 224.39 0.025624 +35 1897.519287 0.000000 18.9752 18.1348 27121.74 27237.73 223.77 0.028232 +36 2076.095947 0.000000 20.7610 19.8681 26524.63 26635.56 223.17 0.031013 +37 2265.431641 0.000000 22.6543 21.7076 25946.90 26053.04 222.60 0.033972 +38 2465.770508 0.000000 24.6577 23.6560 25387.55 25489.15 222.04 0.037115 +39 2677.348145 0.000000 26.7735 25.7156 24845.63 24942.93 221.50 0.040445 +40 2900.391357 0.000000 29.0039 27.8887 24320.28 24413.50 220.97 0.043967 +41 3135.119385 0.000000 31.3512 30.1776 23810.67 23900.02 220.46 0.047685 +42 3381.743652 0.000000 33.8174 32.5843 23316.04 23401.71 219.97 0.051604 +43 3640.468262 0.000000 36.4047 35.1111 22835.68 22917.85 219.49 0.055727 +44 3911.490479 0.000000 39.1149 37.7598 22368.91 22447.75 219.02 0.060059 +45 4194.930664 0.000000 41.9493 40.5321 21915.16 21990.82 218.57 0.064602 +46 4490.817383 0.000000 44.9082 43.4287 21473.98 21546.62 218.12 0.069359 +47 4799.149414 0.000000 47.9915 46.4498 21045.00 21114.77 217.70 0.074330 +48 5119.895020 0.000000 51.1990 49.5952 20627.87 20694.90 217.28 0.079516 +49 5452.990723 0.000000 54.5299 52.8644 20222.24 20286.66 216.87 0.084916 +50 5798.344727 0.000000 57.9834 56.2567 19827.95 19889.88 216.65 0.090458 +51 6156.074219 0.000000 61.5607 59.7721 19443.55 19503.09 216.65 0.096110 +52 6526.946777 0.000000 65.2695 63.4151 19068.35 19125.61 216.65 0.101968 +53 6911.870605 0.000000 69.1187 67.1941 18701.27 18756.34 216.65 0.108045 +54 7311.869141 0.000000 73.1187 71.1187 18341.27 18394.25 216.65 0.114355 +55 7727.412109 0.000007 77.2810 75.1999 17987.41 18038.35 216.65 0.120917 +56 8159.354004 0.000024 81.6182 79.4496 17638.78 17687.77 216.65 0.127751 +57 8608.525391 0.000059 86.1450 83.8816 17294.53 17341.62 216.65 0.134877 +58 9076.400391 0.000112 90.8774 88.5112 16953.83 16999.08 216.65 0.142321 +59 9562.682617 0.000199 95.8280 93.3527 16616.09 16659.55 216.65 0.150106 +60 10065.978516 0.000340 101.0047 98.4164 16281.10 16322.83 216.65 0.158248 +61 10584.631836 0.000562 106.4153 103.7100 15948.85 15988.88 216.65 0.166760 +62 11116.662109 0.000890 112.0681 109.2417 15619.30 15657.70 216.65 0.175655 +63 11660.067383 0.001353 117.9714 115.0198 15292.44 15329.24 216.65 0.184946 +64 12211.547852 0.001992 124.1337 121.0526 14968.24 15003.50 216.65 0.194646 +65 12766.873047 0.002857 130.5637 127.3487 14646.68 14680.44 216.65 0.204770 +66 13324.668945 0.003971 137.2703 133.9170 14327.75 14360.05 216.65 0.215331 +67 13881.331055 0.005378 144.2624 140.7663 14011.41 14042.30 216.65 0.226345 +68 14432.139648 0.007133 151.5493 147.9058 13697.65 13727.18 216.65 0.237825 +69 14975.615234 0.009261 159.1403 155.3448 13386.45 13414.65 216.65 0.249786 +70 15508.256836 0.011806 167.0450 163.0927 13077.79 13104.70 216.65 0.262244 +71 16026.115234 0.014816 175.2731 171.1591 12771.64 12797.30 216.65 0.275215 +72 16527.322266 0.018318 183.8344 179.5537 12467.99 12492.44 216.65 0.288713 +73 17008.789063 0.022355 192.7389 188.2867 12166.81 12190.10 216.65 0.302755 +74 17467.613281 0.026964 201.9969 197.3679 11868.08 11890.24 216.65 0.317357 +75 17901.621094 0.032176 211.6186 206.8078 11571.79 11592.86 216.65 0.332536 +76 18308.433594 0.038026 221.6146 216.6166 11277.92 11297.93 216.65 0.348308 +77 18685.718750 0.044548 231.9954 226.8050 10986.70 11005.69 216.74 0.364545 +78 19031.289063 0.051773 242.7719 237.3837 10696.22 10714.22 218.62 0.378253 +79 19343.511719 0.059728 253.9549 248.3634 10405.61 10422.64 220.51 0.392358 +80 19620.042969 0.068448 265.5556 259.7553 10114.89 10130.98 222.40 0.406868 +81 19859.390625 0.077958 277.5852 271.5704 9824.08 9839.26 224.29 0.421790 +82 20059.931641 0.088286 290.0548 283.8200 9533.20 9547.49 226.18 0.437130 +83 20219.664063 0.099462 302.9762 296.5155 9242.26 9255.70 228.08 0.452897 +84 20337.863281 0.111505 316.3607 309.6684 8951.30 8963.90 229.97 0.469097 +85 20412.308594 0.124448 330.2202 323.2904 8660.32 8672.11 231.86 0.485737 +86 20442.078125 0.138313 344.5663 337.3932 8369.35 8380.36 233.75 0.502825 +87 20425.718750 0.153125 359.4111 351.9887 8078.41 8088.67 235.64 0.520367 +88 20361.816406 0.168910 374.7666 367.0889 7787.51 7797.04 237.53 0.538370 +89 20249.511719 0.185689 390.6450 382.7058 7496.68 7505.51 239.42 0.556842 +90 20087.085938 0.203491 407.0583 398.8516 7205.93 7214.09 241.31 0.575790 +91 19874.025391 0.222333 424.0190 415.5387 6915.29 6922.80 243.20 0.595219 +92 19608.572266 0.242244 441.5395 432.7792 6624.76 6631.66 245.09 0.615138 +93 19290.226563 0.263242 459.6321 450.5858 6334.38 6340.68 246.98 0.635553 +94 18917.460938 0.285354 478.3096 468.9708 6044.15 6049.89 248.86 0.656471 +95 18489.707031 0.308598 497.5845 487.9470 5754.10 5759.30 250.75 0.677899 +96 18006.925781 0.332939 517.4198 507.5021 5464.60 5469.30 252.63 0.699815 +97 17471.839844 0.358254 537.7195 527.5696 5176.77 5180.98 254.50 0.722139 +98 16888.687500 0.384363 558.3430 548.0312 4892.26 4896.02 256.35 0.744735 +99 16262.046875 0.411125 579.1926 568.7678 4612.58 4615.92 258.17 0.767472 +100 15596.695313 0.438391 600.1668 589.6797 4338.77 4341.73 259.95 0.790242 +101 14898.453125 0.466003 621.1624 610.6646 4071.80 4074.41 261.68 0.812937 +102 14173.324219 0.493800 642.0764 631.6194 3812.53 3814.82 263.37 0.835453 +103 13427.769531 0.521619 662.8084 652.4424 3561.70 3563.69 265.00 0.857686 +104 12668.257813 0.549301 683.2620 673.0352 3319.94 3321.67 266.57 0.879541 +105 11901.339844 0.576692 703.3467 693.3043 3087.75 3089.25 268.08 0.900929 +106 11133.304688 0.603648 722.9795 713.1631 2865.54 2866.83 269.52 0.921768 +107 10370.175781 0.630036 742.0855 732.5325 2653.58 2654.69 270.90 0.941988 +108 9617.515625 0.655736 760.5996 751.3426 2452.04 2452.99 272.21 0.961527 +109 8880.453125 0.680643 778.4661 769.5329 2260.99 2261.80 273.45 0.980334 +110 8163.375000 0.704669 795.6396 787.0528 2080.41 2081.09 274.63 0.998368 +111 7470.343750 0.727739 812.0847 803.8622 1910.19 1910.76 275.73 1.015598 +112 6804.421875 0.749797 827.7756 819.9302 1750.14 1750.63 276.77 1.032005 +113 6168.531250 0.770798 842.6959 835.2358 1600.04 1600.44 277.75 1.047576 +114 5564.382813 0.790717 856.8376 849.7668 1459.58 1459.91 278.66 1.062310 +115 4993.796875 0.809536 870.2004 863.5190 1328.43 1328.70 279.52 1.076209 +116 4457.375000 0.827256 882.7910 876.4957 1206.21 1206.44 280.31 1.089286 +117 3955.960938 0.843881 894.6222 888.7066 1092.54 1092.73 281.05 1.101558 +118 3489.234375 0.859432 905.7116 900.1669 987.00 987.15 281.73 1.113047 +119 3057.265625 0.873929 916.0815 910.8965 889.17 889.29 282.37 1.123777 +120 2659.140625 0.887408 925.7571 920.9193 798.62 798.72 282.96 1.133779 +121 2294.242188 0.899900 934.7666 930.2618 714.94 715.02 283.50 1.143084 +122 1961.500000 0.911448 943.1399 938.9532 637.70 637.76 284.00 1.151724 +123 1659.476563 0.922096 950.9082 947.0240 566.49 566.54 284.47 1.159733 +124 1387.546875 0.931881 958.1037 954.5059 500.91 500.95 284.89 1.167147 +125 1143.250000 0.940860 964.7584 961.4311 440.58 440.61 285.29 1.173999 +126 926.507813 0.949064 970.9046 967.8315 385.14 385.16 285.65 1.180323 +127 734.992188 0.956550 976.5737 973.7392 334.22 334.24 285.98 1.186154 +128 568.062500 0.963352 981.7968 979.1852 287.51 287.52 286.28 1.191523 +129 424.414063 0.969513 986.6036 984.2002 244.68 244.69 286.56 1.196462 +130 302.476563 0.975078 991.0230 988.8133 205.44 205.44 286.81 1.201001 +131 202.484375 0.980072 995.0824 993.0527 169.50 169.51 287.05 1.205168 +132 122.101563 0.984542 998.8081 996.9452 136.62 136.62 287.26 1.208992 +133 62.781250 0.988500 1002.2250 1000.5165 106.54 106.54 287.46 1.212498 +134 22.835938 0.991984 1005.3562 1003.7906 79.04 79.04 287.64 1.215710 +135 3.757813 0.995003 1008.2239 1006.7900 53.92 53.92 287.80 1.218650 +136 0.000000 0.997630 1010.8487 1009.5363 30.96 30.96 287.95 1.221341 +137 0.000000 1.000000 1013.2500 1012.0494 10.00 10.00 288.09 1.223803 diff --git a/config/tools.csh b/config/tools.csh index 83c0c816..3a6dc312 100644 --- a/config/tools.csh +++ b/config/tools.csh @@ -22,6 +22,7 @@ set pyTools = ( \ update_sensorScanPosition \ updateXTIME \ concatenate \ + era5_to_int \ ) foreach tool ($pyTools) setenv ${tool} "python ${pyDir}/${tool}.py" diff --git a/initialize/data/ExternalAnalyses.py b/initialize/data/ExternalAnalyses.py index 402e5e67..a85cc051 100644 --- a/initialize/data/ExternalAnalyses.py +++ b/initialize/data/ExternalAnalyses.py @@ -116,8 +116,11 @@ def __init__(self, config:Config, hpc:HPC, meshes:dict): 'queue': {'def': hpc['SharedQueue']}, 'account': {'def': hpc['CriticalAccount']}, } - ungribjob = Resource(self._conf, attr, ('job', 'ungrib')) - self.__ungribtask = TaskLookup[hpc.system](ungribjob) + if self['externalanalyses__UngribPrefix'] == 'GFS': + ungribjob = Resource(self._conf, attr, ('job', 'ungrib')) + self.__ungribtask = TaskLookup[hpc.system](ungribjob) + else: + self.__ungribtask = None ######### # outputs @@ -169,6 +172,25 @@ def export(self, dtOffsets:list=[0]): [['''+base+''']] inherit = '''+base+zeroHR] + # ERA5 GDEX + base = 'GetERA5AnalysisFromGDEX' + queue = 'GetExternalAnalyses' + if base in self['PrepareExternalAnalysisOuter']: + subqueues.append(queue) + taskNames[base] = base+dtLen + self._tasks += [''' + [['''+taskNames[base]+''']] + inherit = '''+queue+''', SingleBatch + script = $origin/bin/'''+base+'''.csh '''+dt_work_Args+''' + execution time limit = PT20M + execution retry delays = '''+self.__getRetry] + + # generic 0hr task name for external classes/tasks to grab + if dt == 0: + self._tasks += [''' + [['''+base+''']] + inherit = '''+base+zeroHR] + # GFS RDA base = 'GetGFSAnalysisFromRDA' queue = 'GetExternalAnalyses' @@ -210,13 +232,13 @@ def export(self, dtOffsets:list=[0]): # ungrib base = 'UngribExternalAnalysis' queue = 'UngribExternalAnalyses' - if base in self['PrepareExternalAnalysisOuter']: + if base in self['PrepareExternalAnalysisOuter'] and self['externalanalyses__UngribPrefixOuter'] != 'ERA5': subqueues.append(queue) taskNames[base] = base+dtLen self._tasks += [''' - [['''+taskNames[base]+''']] - inherit = '''+queue+''', BATCH - script = $origin/bin/'''+base+'''.csh '''+dt_work_Args+''' +[['''+taskNames[base]+''']] + inherit = '''+queue+''', BATCH + script = $origin/bin/'''+base+'''.csh '''+dt_work_Args+''' '''+self.__ungribtask.job()+self.__ungribtask.directives()] # generic 0hr task name for external classes/tasks to grab diff --git a/scenarios/GenerateERA5Analyses.yaml b/scenarios/GenerateERA5Analyses.yaml new file mode 100644 index 00000000..20fea5b8 --- /dev/null +++ b/scenarios/GenerateERA5Analyses.yaml @@ -0,0 +1,18 @@ +suite: GenerateExternalAnalyses + +externalanalyses: + resource: "ERA5.GDEX" +experiment: + name: 'GenerateERA5Analyses' + prefix: '' +model: + outerMesh: 120km +hpc: + CriticalPriority: economy + NonCriticalPriority: economy + #CriticalQueue: main + #NonCriticalQueue: main +workflow: + first cycle point: 20220606T00 + final cycle point: 20220606T06 + max active cycle points: 60 diff --git a/scenarios/defaults/externalanalyses.yaml b/scenarios/defaults/externalanalyses.yaml index ed68a063..9b764850 100644 --- a/scenarios/defaults/externalanalyses.yaml +++ b/scenarios/defaults/externalanalyses.yaml @@ -80,3 +80,15 @@ externalanalyses: - "GetGDASAnalysisFromFTP" #Vtable: /glade/campaign/mmm/parc/liuz/pandac_common/vtables/Vtable.GFS.O3MR retry: '40*PT10M' + + ERA5: + GDEX: + common: + UngribPrefix: ERA5 + PrepareExternalAnalysisTasks: + - "GetERA5AnalysisFromGDEX" + - "ExternalAnalysisToMPAS-{{mesh}}" + - "ExternalAnalysisReady__" + job: + GetAnalysisFrom: + retry: '40*PT10M' diff --git a/tools/WPSUtils.py b/tools/WPSUtils.py new file mode 100644 index 00000000..a0d6de69 --- /dev/null +++ b/tools/WPSUtils.py @@ -0,0 +1,139 @@ +from enum import Enum +import fortran_io as f_io + +class Projections( Enum ) : + LATLON = 0 # Cylindrical equidistant + LC = 1 # Lambert conformal + PS = 2 # Polar stereographic + PS_WGS84 = 102 + MERC = 3 # Mercator + GAUSS = 4 # Gaussian + CYL = 5 + CASSINI = 6 + ALBERS_NAD83 = 105 + ROTLL = 203 + +class IntermediateFile( object ): + def __init__( self, prefix, datestr ) : + self.prefix_ = prefix + self.datestr_ = datestr + self.filename_ = self.prefix_.strip() + ":" + self.datestr_.strip() + self.file_ = open( self.filename_, "wb" ) + + def close( self ) : + self.file_.close() + + def write_next_met_field( + self, + version, + nx, ny, + iproj, + xfcst, + xlvl, + startlat, startlon, starti, startj, + deltalat, deltalon, dx, dy, + xlonc, truelat1, truelat2, + earth_radius, + is_wind_grid_rel, + field, + hdate, + units, + map_source, + desc, + slab + ) : + + f_io.unfmt_ftn_rec_write( + [ version ], + fmt=f_io.StructFormats.INT32, + file=self.file_ + ) + + if field == "GHT" : + field = "HGT" + + if version == 5 : + # FP check okay here? + startloc = None + if starti == 1.0 and startj == 1.0 : + startloc = "SWCORNER" + else : + startloc = "CENTER " + + f_io.unfmt_ftn_rec_write( + [ + hdate.ljust(24)[0:24], + xfcst, + map_source.ljust(32)[0:32], + field.ljust(9)[0:9], + units.ljust(25)[0:25], + desc.ljust(46)[0:46], + xlvl, nx, ny, iproj.value ], + fmt=( [ f_io.StructFormats.ARRCHAR ] + + [ f_io.StructFormats.FP32 ] + + [ f_io.StructFormats.ARRCHAR ] * 4 + + [ f_io.StructFormats.FP32 ] + + [ f_io.StructFormats.INT32 ] * 3 ), + file=self.file_ + ) + if iproj == Projections.LATLON or iproj == Projections.GAUSS: + f_io.unfmt_ftn_rec_write( + [ + startloc.ljust(8)[0:8], + startlat, + startlon, + deltalat, + deltalon, + earth_radius ], + fmt=[ f_io.StructFormats.ARRCHAR ] + [ f_io.StructFormats.FP32 ] * 5, + file=self.file_ + ) + elif iproj == Projections.MERC : + f_io.unfmt_ftn_rec_write( + [ + startloc.ljust(8)[0:8], + startlat, + startlon, + dx, + dy, + truelat1, + earth_radius ], + fmt=[ f_io.StructFormats.ARRCHAR ] + [ f_io.StructFormats.FP32 ] * 6, + file=self.file_ + ) + elif iproj == Projections.LC : + f_io.unfmt_ftn_rec_write( + [ + startloc.ljust(8)[0:8], + startlat, + startlon, + dx, + dy, + xlonc, + truelat1, + truelat2, + earth_radius ], + fmt=[ f_io.StructFormats.ARRCHAR ] + [ f_io.StructFormats.FP32 ] * 8, + file=self.file_ + ) + elif iproj == Projections.PS : + f_io.unfmt_ftn_rec_write( + [ + startloc.ljust(8)[0:8], + startlat, + startlon, + dx, + dy, + xlonc, + truelat1, + earth_radius ], + fmt=[ f_io.StructFormats.ARRCHAR ] + [ f_io.StructFormats.FP32 ] * 7, + file=self.file_ + ) + + f_io.unfmt_ftn_rec_write( [ is_wind_grid_rel ], fmt=f_io.StructFormats.INT32, file=self.file_ ) + f_io.unfmt_ftn_rec_write( slab, fmt=f_io.StructFormats.FP32, file=self.file_ ) + return 0 + else : + print( "Didn't recognize format number " + str( version ) ) + return 1 diff --git a/tools/era5_to_int.py b/tools/era5_to_int.py new file mode 100755 index 00000000..dac7dd9f --- /dev/null +++ b/tools/era5_to_int.py @@ -0,0 +1,725 @@ +#!/usr/bin/env python3 + +class SnowDiags: + """ Implements the computation of SNOW (water equivalent snow depth) and + SNOWH (physical snow depth) from ERA5 RSN (water equivalent snow depth) + and SD (snow density) fields. + """ + + def __init__(self): + self.snow_den = None + self.snow_ec = None + + def consider(self, field, xlvl, proj, hdate, slab, intfile): + """ Considers whether a given field may be used in the computation + of SNOW or SNOWH. When all available information has been acquired, this + method computes these two fields and writes them to the output + intermediate file. + """ + + if field == 'SNOW_DEN' and xlvl == 200100.0: + self.snow_den = slab + elif field == 'SNOW_EC' and xlvl == 200100.0: + self.snow_ec = slab + else: + return + + if self.snow_den is not None and self.snow_ec is not None: + print('Computing SNOWH and SNOW') + snow = self.snow_ec * 1000.0 + snowh = snow / self.snow_den + + write_slab(intfile, snow, 200100.0, proj, 'SNOW', hdate, 'kg m**-2', + 'ERA5 reanalysis grid', 'Water equivalent snow depth') + write_slab(intfile, snowh, 200100.0, proj, 'SNOWH', hdate, 'm', + 'ERA5 reanalysis grid', 'Physical snow depth') + + self.snow_den = None + self.snow_ec = None + + +class RH2mDiags: + """ Implements the computation of RH (relative humidity in % at the surface) + from ERA5 VAR_2T (2m temperature) and VAR_2D (2m dewpoint) fields. + """ + + def __init__(self): + self.t = None + self.td = None + + def consider(self, field, xlvl, proj, hdate, slab, intfile): + """ Considers whether a given field may be used in the computation + of RH at the surface. When all available information has been acquired, + this method computes the RH field and writes it to the output + intermediate file. + """ + + if field == 'TT' and xlvl == 200100.0: + self.t = slab + elif field == 'DEWPT' and xlvl == 200100.0: + self.td = slab + else: + return + + if self.t is not None and self.td is not None: + print('Computing RH at 200100.0') + + Xlv = 2.5e6 + Rv = 461.5 + rh2 = np.exp(Xlv/Rv*(1.0/self.t - 1.0/self.td)) * 1.0e2 + + write_slab(intfile, rh2, 200100.0, proj, 'RH', hdate, '%', + 'ERA5 reanalysis grid', 'Relative humidity') + + self.t = None + self.td = None + + +class RHDiags: + """ Implements the computation of RH (relative humidity in % at pressure level) + from ERA5 T (temperature) and Q (specific humidity) fields at a corresponding + pressure level. + """ + + def __init__(self): + self.savefields = {} + + def liquidSaturationVaporMixRatio(self, t, p): + """ Based on the RSLF function from the Thompson microphysics scheme in WRF + Compute the saturation vapor mixing ratio with respect to liquid water. + """ + + t0 = 273.16 + C0= .611583699E03 + C1= .444606896E02 + C2= .143177157E01 + C3= .264224321E-1 + C4= .299291081E-3 + C5= .203154182E-5 + C6= .702620698E-8 + C7= .379534310E-11 + C8=-.321582393E-13 + + t_bounded =np.maximum(-80.0, t - t0) + # Simplified calc of saturation vapor pressure + # esL = 612.2 * np.exp(17.67 * t_bounded / (self.t-29.65)) + # esL with coefficients used instead + esL = (C0 + t_bounded + * (C1 + t_bounded + * (C2 + t_bounded + * (C3 + t_bounded + * (C4 + t_bounded + * (C5 + t_bounded + * (C6 + t_bounded + * (C7 + t_bounded * C8) + ) + ) + ) + ) + ) + ) + ) + + # Even with P=1050mb and T=55C, the sat. vap. pres only contributes to + # ~15% of total pres. + esL = np.minimum(esL, p * 0.15) + mixRatioESL = .622 * esL / (p - esL) + return mixRatioESL + + def consider(self, field, xlvl, proj, hdate, slab, intfile): + """ Considers whether a given field may be used in the computation + of RH at the given pressure level. When all available information has + been acquired, this method computes the RH field and writes it to the + output intermediate file. + """ + + if field == 'SPECHUMD': + self.savefields[('q', xlvl )] = slab + # Differentiate from surface pressure TT denoted by 200100.0 + elif field == 'TT' and xlvl != 200100.0: + self.savefields[('tt', xlvl)] = slab + else: + return + + if ('tt', xlvl) in self.savefields and ('q', xlvl) in self.savefields: + print('Computing RH at ' + str(xlvl)) + mixRatioESL = self.liquidSaturationVaporMixRatio( self.savefields[('tt', xlvl)], xlvl) + mixRatioWaterVapor = self.savefields[('q', xlvl)] / (1.0 - self.savefields[('q', xlvl)]) + rh = 100.0 * mixRatioWaterVapor / mixRatioESL + + write_slab(intfile, rh, xlvl, proj, 'RH', hdate, '%', + 'ERA5 reanalysis grid', 'Relative humidity') + + del self.savefields[('tt', xlvl)] + del self.savefields[('q', xlvl)] + + +class GeopotentialHeightDiags: + """ Implements the computation of SOILHGT and GHT (geopotential height at + the surface and at upper-air levels, respectively) from ERA5 Z + (geopotential), whose WPS name is either GEOPT at upper-air levels or + SOILGEO at the surface. + """ + + def consider(self, field, xlvl, proj, hdate, slab, intfile): + """ Considers whether the given field represents geopotential, and if + so, computes geopotential height, writing the result to the output + intermediate file. + """ + + if field == 'GEOPT': + print('Computing GHT at ', xlvl) + + write_slab(intfile, slab / 9.81, xlvl, proj, 'GHT', hdate, 'm', + 'ERA5 reanalysis grid', 'Geopotential height') + + elif field == 'SOILGEO': + print('Computing SOILHGT') + + write_slab(intfile, slab / 9.81, 200100.0, proj, 'SOILHGT', hdate, 'm', + 'ERA5 reanalysis grid', 'Geopotential height') + + else: + return + + +class ModelLevelPresGhtDiags: + """ Implements the computation of the 3D pressure field and geopotential height + from ECMWF model-level data utilizing the A and B coefficients provided from + a coefficients array. + This takes the place of the utility program calc_ecmwf_p.exe and thus the + utility program need not be run on the resulting intermediate file. + """ + + def __init__(self, coeff_a, coeff_b): + from collections import defaultdict + + self.coeff_a = coeff_a + self.coeff_b = coeff_b + self.savefields = defaultdict(dict) + + def consider(self, field, xlvl, proj, hdate, slab, intfile): + import numpy as np + + if field == 'SOILGEO': + self.savefields['hgtsfc'] = slab / 9.81 + elif field == 'PSFC': + self.savefields['psfc'] = slab + elif field == 'TT' and xlvl != 200100.0: + self.savefields['tt'][xlvl] = slab + elif field == 'SPECHUMD': + self.savefields['q'][xlvl] = slab + + if ('hgtsfc' in self.savefields and + 'psfc' in self.savefields and + len(self.savefields['tt']) == self.coeff_a.shape[0] - 1 and + len(self.savefields['q']) == self.coeff_a.shape[0] - 1): + + print('Computing pressure and geopotential height from model levels') + levels = sorted(list(self.savefields['tt'].keys())) + + q = self.savefields['q'] + tt = self.savefields['tt'] + + # psfc and hgtsfc are the first at bottom + psfc = self.savefields['psfc'] + hgtprev = self.savefields['hgtsfc'] + print('Computing pressure and geopotential height for surface') + write_slab(intfile, psfc, 200100.0, proj, 'PRESSURE', hdate, 'Pa', 'ERA5 reanalysis grid', 'Pressure') + write_slab(intfile, hgtprev, 200100.0, proj, 'GHT', hdate, 'm', 'ERA5 reanalysis grid', 'Geopotential height') + + # Loop from the surface to top, where level 1 is the top and max level is surface + # Note that there is an extra blank at the very top (index 0) + for idx, level in reversed(list(enumerate(levels))): + print(f'Computing pressure and geopotential height for {level}') + + # Offset by one to "match" levels starting at zero and the way calc_ecmwf_p + # indexes into the coeffs + i = idx + 1 + + # Compute using the half level pressure at the current half-level + a = 0.5 * (self.coeff_a[i-1] + self.coeff_a[i]) + b = 0.5 * (self.coeff_b[i-1] + self.coeff_b[i]) + + # Initial values at the bottom of this level + pres = a + b * psfc + pres_start = psfc + q_half = q[level] + tt_half = tt[level] + + if i != len(levels): + # if not the surface level use half level pressure, temperature, and specific humidity + a = 0.5 * (self.coeff_a[i+1] + self.coeff_a[i]) + b = 0.5 * (self.coeff_b[i+1] + self.coeff_b[i]) + pres_start = a + b * psfc + next_level = levels[levels.index(level)+1] + q_half = 0.5 * (q_half + q[next_level]) + tt_half = 0.5 * (tt_half + tt[next_level]) + + # Compute virtual temperature + tv = 287.05 * tt_half * (1.0 + 0.61 * q_half) + hgtprev = hgtprev + tv * np.log(pres_start / pres) / 9.81 + + write_slab(intfile, pres, level, proj, 'PRESSURE', hdate, 'Pa', 'ERA5 reanalysis grid', '3-d pressure') + write_slab(intfile, hgtprev, level, proj, 'GHT', hdate, 'm', 'ERA5 reanalysis grid', '3-d geopotential height') + + for key in list(self.savefields.keys()): + del self.savefields[key] + + +def days_in_month(year, month): + """ Returns the number of days in a month, depending on the year. + A Gregorian calendar is assumed for the purposes of determining leap + years. + """ + non_leap_year = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ] + leap_year = [ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ] + + if (year % 4 == 0) and ((year % 100 != 0) or (year % 400 == 0)): + return leap_year[month - 1] + else: + return non_leap_year[month - 1] + + +def intdate_to_string(utc_int): + """ Converts an integer date representation into a string. + The digits of utc_int represent yyyymmddhh, and the returned string + is of the form 'yyyy-mm-dd_hh:00:00'. + """ + year = utc_int // 1000000 + month = (utc_int // 10000) % 100 + day = (utc_int // 100) % 100 + hour = utc_int % 100 + return '{:04d}-{:02d}-{:02d}_{:02d}:00:00'.format(year, month, day, hour) + + +def datetime_to_string(dt): + """ Converts a Python datetime instance into a string of the form + 'yyyy-mm-dd_hh'. + """ + return '{:04d}-{:02d}-{:02d}_{:02d}'.format(dt.year, dt.month, dt.day, dt.hour) + + +def string_to_yyyymmddhh(str): + """ Given a string of the form yyyy-mm-dd_hh, returns the component year, + month, day, and hour as a tuple of integers. + """ + tmp = str.split('-') + yyyy = int(tmp[0]) + mm = int(tmp[1]) + ddhh = tmp[2] + dd = int(ddhh.split('_')[0]) + hh = int(ddhh.split('_')[1]) + + return yyyy, mm, dd, hh + + +def begin_6hourly(yyyy, mm, dd, hh): + """ Returns a date-time string of the form yyyymmddhh for the year, month, + day, and hour with the hours rounded down the the beginning of a six-hour + interval, i.e., 0, 6, 12, and 18. + """ + return f'{yyyy:04d}{mm:02d}{dd:02d}{((hh // 6) * 6):02d}' + + +def end_6hourly(yyyy, mm, dd, hh): + """ Returns a date-time string of the form yyyymmddhh for the year, month, + day, and hour with the hours rounded up the the last hour of a six-hour + interval, i.e., 5, 11, 17, and 23. + """ + return f'{yyyy:04d}{mm:02d}{dd:02d}{((hh // 6) * 6 + 5):02d}' + + +def begin_daily(yyyy, mm, dd, hh): + """ Returns a date-time string of the form yyyymmddhh for the year, month, + day, and hour with the hours rounded down the beginning of a day (i.e., to 00). + """ + return f'{yyyy:04d}{mm:02d}{dd:02d}{0:02d}' + + +def end_daily(yyyy, mm, dd, hh): + """ Returns a date-time string of the form yyyymmddhh for the year, month, + day, and hour with the hours rounded up the end of a day (i.e., to 23). + """ + return f'{yyyy:04d}{mm:02d}{dd:02d}{23:02d}' + + +def begin_monthly(yyyy, mm, dd, hh): + """ Returns a date-time string of the form yyyymmddhh for the year, month, + day, and hour with the day and hour rounded down the the beginning of a + monthly interval. + """ + return f'{yyyy:04d}{mm:02d}{1:02d}{0:02d}' + + +def end_monthly(yyyy, mm, dd, hh): + """ Returns a date-time string of the form yyyymmddhh for the year, month, + day, and hour with the day and hour rounded up the the last day and hour + of the month. + """ + return f'{yyyy:04d}{mm:02d}{days_in_month(yyyy,mm):02d}{23:02d}' + + +class MapProjection: + """ Stores parameters of map projections as used in the WPS intermediate + file format. + """ + + def __init__(self, projType, startLat, startLon, startI, startJ, + deltaLat, deltaLon, + dx=0.0, dy=0.0, truelat1=0.0, truelat2=0.0, xlonc=0.0): + self.projType = projType + self.startLat = startLat + self.startLon = startLon + self.startI = startI + self.startJ = startJ + self.deltaLat = deltaLat + self.deltaLon = deltaLon + self.dx = dx + self.dy = dy + self.truelat1 = truelat1 + self.truelat2 = truelat2 + self.xlonc = xlonc + + +class MetVar: + """ Describes a variable to be converted from netCDF to intermediate file + format. + """ + + def __init__(self, WPSname, ERA5name, ERA5file, beginDateFn, endDateFn, + mapProj, isInvariant=False): + self.WPSname = WPSname + self.ERA5name = ERA5name + self.ERA5file = ERA5file + self.beginDateFn = beginDateFn + self.endDateFn = endDateFn + self.mapProj = mapProj + self.isInvariant = isInvariant + + +def find_era5_file(var, validtime, localpaths=None): + """ Returns a filename with path information of an ERA5 netCDF file given + the ERA5 netCDF variable name and the valid time of the variable. + The localpaths argument may be used to specify an array of local + directories to search for ERA5 netCDF files. If localpaths is not provided + or is None, known paths on the NWSC Glade filesystem are searched. + If no file can be found containing the specified variable at the specified + time, a RuntimeError exception is raised. + """ + import os.path + + glade_paths = [ + # New Glade paths as of March 2025 + '/glade/campaign/collections/rda/data/d633006/e5.oper.an.ml/', + '/glade/campaign/collections/rda/data/d633000/e5.oper.an.pl/', + '/glade/campaign/collections/rda/data/d633000/e5.oper.an.sfc/', + '/glade/campaign/collections/rda/data/d633000/e5.oper.invariant/197901/', + '/glade/campaign/collections/rda/data/d633006/e5.oper.invariant/', + + # Old Glade paths prior to March 2025 + '/glade/campaign/collections/rda/data/ds633.6/e5.oper.an.ml/', + '/glade/campaign/collections/rda/data/ds633.0/e5.oper.an.pl/', + '/glade/campaign/collections/rda/data/ds633.0/e5.oper.an.sfc/', + '/glade/campaign/collections/rda/data/ds633.0/e5.oper.invariant/197901/', + '/gpfs/csfs1/collections/rda/decsdata/ds630.0/P/e5.oper.invariant/201601/' + ] + + if localpaths != None: + file_paths = localpaths + else: + file_paths = glade_paths + + yyyy, mm, dd, hh = string_to_yyyymmddhh(validtime) + + begin_date = var.beginDateFn(yyyy, mm, dd, hh) + end_date = var.endDateFn(yyyy, mm, dd, hh) + + for p in file_paths: + if var.isInvariant or localpaths != None: + # For time-invariant fields, or if we are searching local paths, + # assume that there is no need to append yyyymm to the end of + # directories. + base_path = p + else: + base_path = p + f'{yyyy:04d}{mm:02d}/' + + filename = base_path + var.ERA5file.format(begin_date, end_date) + if os.path.isfile(filename): + return filename + + raise RuntimeError('Error: Could not find file ' + + var.ERA5file.format(begin_date, end_date) + + ' needed for ERA5 variable ' + var.ERA5name) + + +def find_time_index(ncfilename, validtime): + """ Returns a 0-based offset into the unlimited dimension of the specified + valid time within the specified etCDF file. + If the specified time is not found in the file, an offset of -1 is + returned. + The 'utc_date' variable is assumed to exist in the netCDF file, and it is + through this variable that the search for the valid time takes place. + """ + from netCDF4 import Dataset + import numpy as np + + yyyy, mm, dd, hh = string_to_yyyymmddhh(validtime) + + needed_date = yyyy * 1000000 + mm * 10000 + dd * 100 + hh + + with Dataset(ncfilename) as f: + utc_date = f.variables['utc_date'][:] + + idx = np.where(needed_date == utc_date) + if idx[0].size == 0: + return -1 + else: + return idx[0][0] + + +def write_slab( intfile, slab, xlvl, proj, WPSname, hdate, units, map_source, desc): + """ Writes a 2-d array (a 'slab' of data) to an opened intermediate file + using the provided level, projection, and other metadata. + This routine assumes that the intermediate file has already been created + through a previous call to the WPSUtils.intermediate.write_met_init + routine. + """ + missing_value = -1.0e30 + stat = intfile.write_next_met_field( + 5, slab.shape[1], slab.shape[0], proj.projType, 0.0, xlvl, + proj.startLat, proj.startLon, proj.startI, proj.startJ, + proj.deltaLat, proj.deltaLon, proj.dx, proj.dy, proj.xlonc, + proj.truelat1, proj.truelat2, 6371229.0, 0, WPSname, + hdate, units, map_source, desc, slab.filled(missing_value)) + + +def add_trailing_slash(str): + """ Returns str with a forward slash appended to it if the str argument does + not end with a forward slash character. Otherwise, return str unmodified if + it already ends with a slash. + """ + if str[-1] == '/': + return str + else: + return str + '/' + + +def handle_datetime_args(args): + """ Returns starting and ending datetime objects along with an interval + timedelta object given an argparse namespace with members 'datetime', + 'until_datetime', and 'interval_hours'. If 'until_datetime' is None, + the returned ending datetime is the same as the starting datetime. + If 'interval_hours' is 0, it defaults to a interval of six hours. + """ + + yyyy, mm, dd, hh = string_to_yyyymmddhh(args.datetime) + startDate = datetime.datetime(yyyy, mm, dd, hh) + + if args.until_datetime != None: + yyyy, mm, dd, hh = string_to_yyyymmddhh(args.until_datetime) + endDate = datetime.datetime(yyyy, mm, dd, hh) + else: + endDate = startDate + + if endDate < startDate: + raise ValueError('until_datetime precedes datetime') + + if args.interval_hours > 0: + intvH = datetime.timedelta(hours=args.interval_hours) + else: + raise ValueError('interval_hours is not a positive integer') + + return startDate, endDate, intvH + + +if __name__ == '__main__': + from netCDF4 import Dataset + import numpy as np + import WPSUtils + import argparse + import datetime + import sys + + parser = argparse.ArgumentParser() + parser.add_argument('datetime', help='the date-time to convert in YYYY-MM-DD_HH format') + parser.add_argument('until_datetime', nargs='?', default=None, help='the date-time in YYYY-MM-DD_HH format until which records are converted (Default: datetime)') + parser.add_argument('interval_hours', type=int, nargs='?', default=6, help='the interval in hours between records to be converted (Default: %(default)s)') + parser.add_argument('-p', '--path', help='the local path to search for ERA5 netCDF files') + parser.add_argument('-i', '--isobaric', action='store_true', help='use ERA5 pressure-level data rather than model-level data') + parser.add_argument('-v', '--variables', help='a comma-separated list, without spaces, of WPS variable names to process') + args = parser.parse_args() + + try: + startDate, endDate, intvH = handle_datetime_args(args) + except ValueError as e: + print('Error in argument list: ' + e.args[0]) + sys.exit(1) + + print('datetime = ', startDate) + print('until_datetime = ', endDate) + print('interval_hours = ', intvH) + + # Set up the two map projections used in the ERA5 fields to be converted + Gaussian = MapProjection(WPSUtils.Projections.GAUSS, + 89.7848769072, 0.0, 1.0, 1.0, 640.0 / 2.0, 360.0 / 1280.0) + LatLon = MapProjection(WPSUtils.Projections.LATLON, + 90.0, 0.0, 1.0, 1.0, -0.25, 0.25) + + diagnostics = [] + diagnostics.append(SnowDiags()) + diagnostics.append(RH2mDiags()) + diagnostics.append(GeopotentialHeightDiags()) + + if args.isobaric: + diagnostics.append(RHDiags()) + else: + coeffs = np.loadtxt('./ecmwf_coeffs', usecols=(1,2)) + diagnostics.append(ModelLevelPresGhtDiags(coeffs[:,0], coeffs[:,1])) + + int_vars = [] + if args.isobaric: + int_vars.append(MetVar('GEOPT', 'Z', 'e5.oper.an.pl.128_129_z.ll025sc.{}_{}.nc', begin_daily, end_daily, LatLon)) + int_vars.append(MetVar('SPECHUMD', 'Q', 'e5.oper.an.pl.128_133_q.ll025sc.{}_{}.nc', begin_daily, end_daily, LatLon)) + int_vars.append(MetVar('TT', 'T', 'e5.oper.an.pl.128_130_t.ll025sc.{}_{}.nc', begin_daily, end_daily, LatLon)) + int_vars.append(MetVar('UU', 'U', 'e5.oper.an.pl.128_131_u.ll025uv.{}_{}.nc', begin_daily, end_daily, LatLon)) + int_vars.append(MetVar('VV', 'V', 'e5.oper.an.pl.128_132_v.ll025uv.{}_{}.nc', begin_daily, end_daily, LatLon)) + else: + int_vars.append(MetVar('SPECHUMD', 'Q', 'e5.oper.an.ml.0_5_0_1_0_q.regn320sc.{}_{}.nc', begin_6hourly, end_6hourly, Gaussian)) + int_vars.append(MetVar('TT', 'T', 'e5.oper.an.ml.0_5_0_0_0_t.regn320sc.{}_{}.nc', begin_6hourly, end_6hourly, Gaussian)) + int_vars.append(MetVar('UU', 'U', 'e5.oper.an.ml.0_5_0_2_2_u.regn320uv.{}_{}.nc', begin_6hourly, end_6hourly, Gaussian)) + int_vars.append(MetVar('VV', 'V', 'e5.oper.an.ml.0_5_0_2_3_v.regn320uv.{}_{}.nc', begin_6hourly, end_6hourly, Gaussian)) + int_vars.append(MetVar('LANDSEA', 'LSM', 'e5.oper.invariant.128_172_lsm.ll025sc.1979010100_1979010100.nc', begin_monthly, end_monthly, LatLon, isInvariant=True)) + int_vars.append(MetVar('SST', 'SSTK', 'e5.oper.an.sfc.128_034_sstk.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SKINTEMP', 'SKT', 'e5.oper.an.sfc.128_235_skt.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SM000007', 'SWVL1', 'e5.oper.an.sfc.128_039_swvl1.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SM007028', 'SWVL2', 'e5.oper.an.sfc.128_040_swvl2.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SM028100', 'SWVL3', 'e5.oper.an.sfc.128_041_swvl3.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SM100289', 'SWVL4', 'e5.oper.an.sfc.128_042_swvl4.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('ST000007', 'STL1', 'e5.oper.an.sfc.128_139_stl1.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('ST007028', 'STL2', 'e5.oper.an.sfc.128_170_stl2.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('ST028100', 'STL3', 'e5.oper.an.sfc.128_183_stl3.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('ST100289', 'STL4', 'e5.oper.an.sfc.128_236_stl4.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SEAICE', 'CI', 'e5.oper.an.sfc.128_031_ci.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('TT', 'VAR_2T', 'e5.oper.an.sfc.128_167_2t.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('DEWPT', 'VAR_2D', 'e5.oper.an.sfc.128_168_2d.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('UU', 'VAR_10U', 'e5.oper.an.sfc.128_165_10u.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('VV', 'VAR_10V', 'e5.oper.an.sfc.128_166_10v.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SNOW_DEN', 'RSN', 'e5.oper.an.sfc.128_033_rsn.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('SNOW_EC', 'SD', 'e5.oper.an.sfc.128_141_sd.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('PMSL', 'MSL', 'e5.oper.an.sfc.128_151_msl.ll025sc.{}_{}.nc', begin_monthly, end_monthly, LatLon)) + int_vars.append(MetVar('PSFC', 'SP', 'e5.oper.an.ml.128_134_sp.regn320sc.{}_{}.nc', begin_6hourly, end_6hourly, Gaussian)) + int_vars.append(MetVar('SOILGEO', 'Z', 'e5.oper.invariant.128_129_z.regn320sc.2016010100_2016010100.nc', begin_monthly, end_monthly, Gaussian, isInvariant=True)) + + # WPS variable names in dont_output cause the corresponding field to not be + # written to the output intermediate file. This is useful when fields are + # read from ERA5 netCDF files only to be used in diagnosing other fields. + dont_output = [] + dont_output.append('GEOPT') # Only used to diagnose GHT + dont_output.append('SOILGEO') # Only used to diagnose SOILHGT + dont_output.append('DEWPT') # Only used to diagnose RH + dont_output.append('SNOW_EC') # Only used in diagnosing SNOW and SNOWH + dont_output.append('SNOW_DEN') # Only used in diagnosing SNOW and SNOWH + + if args.path != None: + paths = [ add_trailing_slash(p) for p in args.path.split(',') ] + else: + paths = None + + # Prepare a list of variables to process. By default, this list contains + # all WPS variables from the int_vars list, but the user may have supplied + # a subset of these variables. + var_set = [ int_var.WPSname for int_var in int_vars ] + if args.variables != None: + user_var_set = args.variables.split(',') + + errcount = 0 + for v in user_var_set: + if v not in var_set: + errcount = errcount + 1 + print('Error: ' + v + ' is not a known WPS variable') + + if errcount > 0: + print('The following are known WPS variables:', var_set) + sys.exit(1) + + var_set = user_var_set + + currDate = startDate + while currDate <= endDate: + initdate = datetime_to_string(currDate) + print('Processing time record ' + initdate) + + intfile = WPSUtils.IntermediateFile('ERA5', initdate) + + for v in int_vars: + + if v.WPSname not in var_set: + continue + + try: + e5filename = find_era5_file(v, initdate, localpaths=paths) + except RuntimeError as e: + print(e.args[0]) + intfile.close() + sys.exit(1) + + idx = find_time_index(e5filename, initdate) + if idx == -1: + idx = 0 + proj = v.mapProj + + print(e5filename) + with Dataset(e5filename) as f: + hdate = intdate_to_string(f.variables['utc_date'][idx]) + print('Converting ' + v.WPSname + ' at ' + hdate) + map_source = 'ERA5 reanalysis grid ' + units = f.variables[v.ERA5name].units + desc = f.variables[v.ERA5name].long_name + field_arr = f.variables[v.ERA5name][idx,:] + + if field_arr.ndim == 2: + slab = field_arr + xlvl = 200100.0 + + # There are some special cases in which 2-d fields are not + # surface (xlvl = 200100.0) variables + if v.WPSname == 'SOILGEO': + xlvl = 1.0 + elif v.WPSname == 'PMSL': + xlvl = 201300.0 + + if v.WPSname in dont_output: + print(v.WPSname + ' is NOT being written to the ' + + 'intermediate file at level', xlvl) + else: + write_slab(intfile, slab, xlvl, proj, v.WPSname, hdate, + units, map_source, desc) + + for diag in diagnostics: + diag.consider(v.WPSname, xlvl, proj, hdate, slab, intfile) + else: + for k in range(f.dimensions['level'].size): + slab = field_arr[k,:,:] + + # For isobaric levels, the WPS intermediate file xlvl + # value should be in Pa, while the ERA5 netCDF files + # provide units of hPa + if args.isobaric: + xlvl = f.variables['level'][k] * 100.0 # Convert hPa to Pa + else: + xlvl = float(f.variables['level'][k]) # Level index + + if v.WPSname in dont_output: + print(v.WPSname + ' is NOT being written to the ' + + 'intermediate file at level', xlvl) + else: + write_slab(intfile, slab, xlvl, proj, + v.WPSname, hdate, units, map_source, desc) + + for diag in diagnostics: + diag.consider(v.WPSname, xlvl, proj, hdate, slab, intfile) + + intfile.close() + + currDate += intvH diff --git a/tools/fortran_io.py b/tools/fortran_io.py new file mode 100644 index 00000000..df45fdae --- /dev/null +++ b/tools/fortran_io.py @@ -0,0 +1,137 @@ +from collections.abc import Iterable +import struct +from enum import Enum +from typing import Union +from io import TextIOWrapper +import numpy as np + +################################################################################ +# +# Wrapper classes around struct formatting options just to constrain inputs +# +class StructFormats( Enum ) : + # # C Type Python Type Standard size in bytes (note) + NULL = "x" # pad byte no value (7) + CHAR = "c" # char bytes of length 1 1 + SCHAR = "b" # signed char integer 1 (1), (2) + UCHAR = "B" # unsigned char integer 1 (2) + BOOL = "?" # _Bool bool 1 (1) + INT16 = "h" # short integer 2 (2) + UINT16 = "H" # unsigned short integer 2 (2) + INT32 = "i" # int integer 4 (2) + UINT32 = "I" # unsigned int integer 4 (2) + # = "l" # long integer 4 (2) + # = "L" # unsigned long integer 4 (2) + INT64 = "q" # long long integer 8 (2) + UINT64 = "Q" # unsigned long long integer 8 (2) + # = "n" # ssize_t integer (3) + # = "N" # size_t integer (3) + FP16 = "e" # *IEEE754 half float(6) float 2 (4) + FP32 = "f" # float float 4 (4) + FP64 = "d" # double float 8 (4) + ARRCHAR = "s" # char[] bytes (9) + PCHAR = "p" # char[] bytes (8) + # = "P" # void* integer (5) + # Notes : refer to https://docs.python.org/3/library/struct.html#format-characters + + def to_dtype( self ) : + if "INT" in self.name or "FP" in self.name : + return self.name[0].lower() + str(int(int(self.name[len(self.name.rstrip( "0123456789" )):]) / 8)) + +class StructEndian( Enum ) : + BIG = ">" + LITTLE = "<" + + +def calcsize( data : Iterable, fmt : Union[StructFormats, Iterable[StructFormats]] ) : + """Calcuclates size of iterable data using format + data - iterable of primitive data (including str) + fmt - a single StructFormat to use on all data elements or iterable index matched to each data element + """ + sizeBytes = 0 + isIter = isinstance( fmt, Iterable ) + ithFmt = None + + if isIter : + if len( fmt ) != len( data ) : + raise Exception( "Mismatch size of struct format and data" ) + else : + ithFmt = fmt.value + + if ithFmt == StructFormats.ARRCHAR.value or ithFmt == StructFormats.PCHAR.value or isIter : + for i, elem in enumerate( data ) : + if isIter : + ithFmt = fmt[i].value + + count = "" + if ithFmt == StructFormats.ARRCHAR.value or ithFmt == StructFormats.PCHAR.value : + count = str( len( elem ) ) + + sizeBytes += struct.calcsize( count + ithFmt ) + else : + # we know they are all the same + totalElems = len( data ) + if isinstance( data, np.ndarray ) : + totalElems = data.size + sizeBytes += struct.calcsize( ithFmt ) * totalElems + + return sizeBytes + + + + +def unfmt_ftn_rec_write( + data : Iterable, + filename : str = None, + file : TextIOWrapper = None, + endian : StructEndian = StructEndian.BIG, + fmt : Union[StructFormats, Iterable[StructFormats]] = StructFormats.INT32, + append : bool = True + ) : + """Writes out a Fortran unformatted record using byte amount record labels + filename - file to write to if provided + file - file handle to write to if provided + data - iterable of primitive data (including str) + endian - endianness of data and record label + fmt - a single StructFormats to use on all data elements or iterable index matched to each data element + append - append to file or overwrite existing data + """ + if filename is not None : + file = open( filename, ( "a" if append else "w" ) + "b" ) + elif file is None : + print( "Error : No filename or file handle provided!" ) + return 1 + + sizeBytes = calcsize( data, fmt ) + file.write( struct.pack( endian.value + StructFormats.UINT32.value, sizeBytes ) ) + isIter = isinstance( fmt, Iterable ) + ithFmt = None + + if not isIter : + ithFmt = fmt.value + + if ithFmt == StructFormats.ARRCHAR.value or ithFmt == StructFormats.PCHAR.value or isIter : + for i, elem in enumerate( data ) : + if isIter : + ithFmt = fmt[i].value + + if isinstance( elem, str ) : + file.write( struct.pack( endian.value + str( len( elem ) ) + ithFmt, elem.encode( 'utf-8' ) ) ) + else : + file.write( struct.pack( endian.value + ithFmt, elem ) ) + else : + # we know they are all the same + if isinstance( data, np.ndarray ) : + file.write( np.asfortranarray( data, dtype=endian.value + fmt.to_dtype() ).tobytes() ) + else : + baseStruct = struct.Struct( endian.value + str( len( data ) ) + ithFmt ) + file.write( baseStruct.pack( *data ) ) + + + file.write( struct.pack( endian.value + StructFormats.UINT32.value, sizeBytes ) ) + + if filename is not None : + file.close() + return 0 + + From bab0a2b32fe309c656e4c34e4de9e2ab1eecd44c Mon Sep 17 00:00:00 2001 From: ibanos90 Date: Wed, 18 Feb 2026 11:07:54 -0700 Subject: [PATCH 2/5] Increase walltime slightly --- scenarios/defaults/initic.yaml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scenarios/defaults/initic.yaml b/scenarios/defaults/initic.yaml index 91eac6d6..fc8be233 100644 --- a/scenarios/defaults/initic.yaml +++ b/scenarios/defaults/initic.yaml @@ -15,18 +15,18 @@ initic: retry: '0*PT30S' 60-3km: - seconds: 320 + seconds: 350 nodes: 2 PEPerNode: 128 30km: - seconds: 320 + seconds: 350 nodes: 2 PEPerNode: 128 60km: - seconds: 160 + seconds: 190 nodes: 1 PEPerNode: 128 120km: - seconds: 120 + seconds: 150 nodes: 1 PEPerNode: 128 From fa7e210a71db58b25a91937e9eabf18dd106b2bf Mon Sep 17 00:00:00 2001 From: ibanos90 Date: Wed, 18 Feb 2026 11:30:11 -0700 Subject: [PATCH 3/5] Remove empty line --- bin/GetERA5AnalysisFromGDEX.csh | 1 - 1 file changed, 1 deletion(-) diff --git a/bin/GetERA5AnalysisFromGDEX.csh b/bin/GetERA5AnalysisFromGDEX.csh index 5af0753b..0d6450b6 100755 --- a/bin/GetERA5AnalysisFromGDEX.csh +++ b/bin/GetERA5AnalysisFromGDEX.csh @@ -66,7 +66,6 @@ endif echo "Getting ERA5 analysis from GDEX" - # Link ECMWF coefficients ln -sf ${ModelConfigDir}/initic/ecmwf_coeffs ./ecmwf_coeffs From 6ca65fcabd790d9220fd5e2de8df7c2c004e2d60 Mon Sep 17 00:00:00 2001 From: ibanos90 Date: Fri, 20 Feb 2026 16:33:22 -0700 Subject: [PATCH 4/5] Attend Jim's comments --- scenarios/GenerateERA5Analyses.yaml | 1 - scenarios/defaults/externalanalyses.yaml | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/scenarios/GenerateERA5Analyses.yaml b/scenarios/GenerateERA5Analyses.yaml index 20fea5b8..44b9a461 100644 --- a/scenarios/GenerateERA5Analyses.yaml +++ b/scenarios/GenerateERA5Analyses.yaml @@ -15,4 +15,3 @@ hpc: workflow: first cycle point: 20220606T00 final cycle point: 20220606T06 - max active cycle points: 60 diff --git a/scenarios/defaults/externalanalyses.yaml b/scenarios/defaults/externalanalyses.yaml index 9b764850..be147bef 100644 --- a/scenarios/defaults/externalanalyses.yaml +++ b/scenarios/defaults/externalanalyses.yaml @@ -91,4 +91,4 @@ externalanalyses: - "ExternalAnalysisReady__" job: GetAnalysisFrom: - retry: '40*PT10M' + retry: '3*PT5S' From 686901c46f509e753727253f616b0a6cd62db71b Mon Sep 17 00:00:00 2001 From: ibanos90 Date: Fri, 20 Feb 2026 16:37:26 -0700 Subject: [PATCH 5/5] Remove old paths to era5 data --- tools/era5_to_int.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/tools/era5_to_int.py b/tools/era5_to_int.py index dac7dd9f..bf0a1efb 100755 --- a/tools/era5_to_int.py +++ b/tools/era5_to_int.py @@ -410,19 +410,12 @@ def find_era5_file(var, validtime, localpaths=None): import os.path glade_paths = [ - # New Glade paths as of March 2025 + # Glade paths as of March 2025 '/glade/campaign/collections/rda/data/d633006/e5.oper.an.ml/', '/glade/campaign/collections/rda/data/d633000/e5.oper.an.pl/', '/glade/campaign/collections/rda/data/d633000/e5.oper.an.sfc/', '/glade/campaign/collections/rda/data/d633000/e5.oper.invariant/197901/', '/glade/campaign/collections/rda/data/d633006/e5.oper.invariant/', - - # Old Glade paths prior to March 2025 - '/glade/campaign/collections/rda/data/ds633.6/e5.oper.an.ml/', - '/glade/campaign/collections/rda/data/ds633.0/e5.oper.an.pl/', - '/glade/campaign/collections/rda/data/ds633.0/e5.oper.an.sfc/', - '/glade/campaign/collections/rda/data/ds633.0/e5.oper.invariant/197901/', - '/gpfs/csfs1/collections/rda/decsdata/ds630.0/P/e5.oper.invariant/201601/' ] if localpaths != None: