From f7a1b91f5ff438f51deb6e0d0d5b5e759db932b6 Mon Sep 17 00:00:00 2001 From: William Putman Date: Tue, 30 Jul 2024 12:46:07 -0400 Subject: [PATCH 01/15] added more protections on QSAT tables --- GEOS_Shared/GEOS_Utilities.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOS_Shared/GEOS_Utilities.F90 b/GEOS_Shared/GEOS_Utilities.F90 index 137db95d..d9aa04fe 100644 --- a/GEOS_Shared/GEOS_Utilities.F90 +++ b/GEOS_Shared/GEOS_Utilities.F90 @@ -546,7 +546,7 @@ function QSAT0(TL,PL,RAMP,PASCALS,DQSAT) result(QSAT) end if TI = (TI - TMINTBL)*DEGSUBS+1 - IT = int(TI) + IT = min(max(1,int(TI)),TABLESIZE-1) if(URAMP==TMIX) then DQ = ESTBLX(IT+1) - ESTBLX(IT) @@ -764,7 +764,7 @@ function DQSAT0(TL,PL,RAMP,PASCALS,QSAT) result(DQSAT) end if TT = (TI - TMINTBL)*DEGSUBS+1 - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) if(URAMP==TMIX) then DQQ = ESTBLX(IT+1) - ESTBLX(IT) From 20092d4d9980ba615750da376c75a0a8aba14747 Mon Sep 17 00:00:00 2001 From: William Putman Date: Tue, 13 Aug 2024 12:49:42 -0400 Subject: [PATCH 02/15] more protections on bounds for QSAT tables --- GEOS_Shared/qsat_fast.code | 6 +++--- GEOS_Shared/qsat_new.code | 10 +++++----- GEOS_Shared/qsatice.code | 2 +- GEOS_Shared/qsatlqu.code | 2 +- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/GEOS_Shared/qsat_fast.code b/GEOS_Shared/qsat_fast.code index 8026d666..9f06696c 100644 --- a/GEOS_Shared/qsat_fast.code +++ b/GEOS_Shared/qsat_fast.code @@ -16,7 +16,7 @@ end if TT = TT - (TMIN-1) - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) Qsat = (TT-IT)*(ESTBLX(IT+1) - ESTBLX(IT)) + ESTBLX(IT) DD = (PL - ERFAC*Qsat) @@ -33,7 +33,7 @@ end if TT = TT - (TMIN-1) - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) Qsat(I) = (TT-IT)*(ESTBLX(IT+1) - ESTBLX(IT)) + ESTBLX(IT) DD = (PL(I) - ERFAC*Qsat(I)) @@ -52,7 +52,7 @@ end if TT = TT - (TMIN-1) - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) Qsat(I,J) = (TT-IT)*(ESTBLX(IT+1) - ESTBLX(IT)) + ESTBLX(IT) DD = (PL(I,J) - ERFAC*Qsat(I,J)) diff --git a/GEOS_Shared/qsat_new.code b/GEOS_Shared/qsat_new.code index 230d34a0..00b599a4 100644 --- a/GEOS_Shared/qsat_new.code +++ b/GEOS_Shared/qsat_new.code @@ -28,7 +28,7 @@ TI = MAX(MIN(TL,TMAX),TMIN) TT = TI - TMIN - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) TT = TT-IT TI = TI - ZEROC @@ -67,7 +67,7 @@ TI = MAX(MIN(TL(I),TMAX),TMIN) TT = TI - TMIN - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) TT = TT-IT TI = TI - ZEROC @@ -109,7 +109,7 @@ TI = MAX(MIN(TL(I,J),TMAX),TMIN) TT = TI - TMIN - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) TT = TT-IT TI = TI - ZEROC @@ -152,7 +152,7 @@ TI = MAX(MIN(TL(I,J,K),TMAX),TMIN) TT = TI - TMIN - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) TT = TT-IT TI = TI - ZEROC @@ -197,7 +197,7 @@ TI = MAX(MIN(TL(I,J,K,L),TMAX),TMIN) TT = TI - TMIN - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) TT = TT-IT TI = TI - ZEROC diff --git a/GEOS_Shared/qsatice.code b/GEOS_Shared/qsatice.code index 083993a3..5c699a3f 100644 --- a/GEOS_Shared/qsatice.code +++ b/GEOS_Shared/qsatice.code @@ -14,7 +14,7 @@ if(UTBL) then if(present(DQ)) DDQ = 0.0 else TT = (TX - TMINTBL)*DEGSUBS+1 - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) DDQ = ESTBLE(IT+1) - ESTBLE(IT) EX = ((TT-IT)*DDQ + ESTBLE(IT)) end if diff --git a/GEOS_Shared/qsatlqu.code b/GEOS_Shared/qsatlqu.code index 767974c3..7583309a 100644 --- a/GEOS_Shared/qsatlqu.code +++ b/GEOS_Shared/qsatlqu.code @@ -14,7 +14,7 @@ if(UTBL) then if(present(DQ)) DDQ = 0.0 else TT = (TX - TMINTBL)*DEGSUBS+1 - IT = int(TT) + IT = min(max(1,int(TT)),TABLESIZE-1) DDQ = ESTBLW(IT+1) - ESTBLW(IT) EX = ((TT-IT)*DDQ + ESTBLW(IT)) end if From e13522ecff2f9750fd609d6a7a8ac3a1dbb4d1af Mon Sep 17 00:00:00 2001 From: William Putman Date: Tue, 13 Aug 2024 12:50:24 -0400 Subject: [PATCH 03/15] added L186 extended model top --- GMAO_hermes/m_set_eta.F90 | 88 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 87 insertions(+), 1 deletion(-) diff --git a/GMAO_hermes/m_set_eta.F90 b/GMAO_hermes/m_set_eta.F90 index 149fd2d8..eb676bd0 100644 --- a/GMAO_hermes/m_set_eta.F90 +++ b/GMAO_hermes/m_set_eta.F90 @@ -91,7 +91,10 @@ subroutine set_eta_r8_(km, ks, ptop, pint, ak, bk) real(r8) a91(92), b91(92) real(r8) a137(138), b137(138) real(r8) a181(182), b181(182) - + +! Extended model top + real(r8) a186(187), b186(187) + ! NCEP real(r8) a64(65),b64(65), a64_sig(65),b64_sig(65) real(r8) a127(128),b127(128) @@ -884,6 +887,84 @@ subroutine set_eta_r8_(km, ks, ptop, pint, ak, bk) 0.985649, 0.988688, 0.991432, 0.993898, 0.996109, & 0.998114, 1.00000 / +! Extended model top + data a186 /0.00025000, 0.00600000, 0.0500000, 0.20000, 0.5000000,& + 1.0000000, 1.6461494, 2.3988098, 3.3373925, 4.5180100, & + 5.9846950, 7.7937295, 10.010172, 12.710444, 15.973860, & + 19.886346, 24.538486, 30.024769, 36.447639, 43.915803, & + 52.545737, 62.453836, 73.764708, 86.607943, 101.11680, & + 117.42552, 135.67413, 156.00783, 178.57081, 203.51057, & + 230.97658, 261.12225, 294.09498, 330.04762, 369.13346, & + 411.50830, 457.32002, 506.71984, 559.85860, 616.88715, & + 677.95229, 743.20095, 812.77976, 886.82963, 965.49068, & + 1048.9029, 1137.2056, 1230.5308, 1329.0094, 1432.7708, & + 1541.9456, 1656.6555, 1777.0232, 1903.1678, 2035.2088, & + 2173.2571, 2317.4255, 2467.8252, 2624.5658, 2787.7489, & + 2957.4768, 3133.8507, 3316.9714, 3506.9325, 3703.8235, & + 3907.7214, 4118.6851, 4336.7484, 4561.9281, 4794.2279, & + 5033.6431, 5280.1592, 5533.7681, 5794.5055, 6062.4980, & + 6337.9868, 6621.3268, 6912.9118, 7213.1229, 7522.2804, & + 7840.7098, 8168.7861, 8506.9408, 8855.5426, 9214.7142, & + 9584.2399, 9963.5983, 10352.111, 10748.987, 11153.318, & + 11564.007, 11979.737, 12399.133, 12820.830, 13243.570, & + 13665.922, 14086.428, 14503.591, 14916.111, 15322.602, & + 15721.715, 16112.059, 16492.407, 16861.516, 17218.218, & + 17561.331, 17889.768, 18202.451, 18498.379, 18776.566, & + 19036.070, 19275.971, 19495.400, 19693.530, 19869.553, & + 20022.686, 20152.202, 20257.406, 20337.620, 20392.130, & + 20420.286, 20421.458, 20395.044, 20340.425, 20257.007, & + 20144.234, 20001.538, 19828.378, 19624.186, 19388.455, & + 19120.682, 18820.454, 18487.599, 18122.389, 17725.721, & + 17299.086, 16844.410, 16363.835, 15859.595, 15334.131, & + 14790.135, 14230.490, 13658.113, 13075.906, 12486.888, & + 11894.071, 11300.423, 10708.651, 10121.401, 9541.1803, & + 8970.3675, 8411.0514, 7865.0798, 7334.1450, 6819.7266, & + 6323.0938, 5845.2210, 5386.9040, 4948.7590, 4531.2215, & + 4134.5581, 3758.8643, 3404.0995, 3070.1029, 2756.5979, & + 2463.2467, 2189.6016, 1935.1705, 1699.3937, 1481.6793, & + 1281.4331, 1098.0062, 930.81258, 779.20329, 642.64526, & + 520.55054, 412.42879, 317.77547, 236.18368, 167.34212, & + 111.00464, 67.106120, 35.426839, 15.342204, 4.8543294, & + 0.89721680, 0.0000000 / + data b186 /0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, & + 0.00000, 0.00000, 0.00000, 8.66281e-07, 4.30872e-06, & + 1.23622e-05, 2.78373e-05, 5.24551e-05, 8.85466e-05, 0.000139698, & + 0.000213463, 0.000319596, 0.000468744, 0.000671591, 0.000940946, & + 0.00129107, 0.00173889, 0.00230345, 0.00300141, 0.00384951, & + 0.00486404, 0.00606625, 0.00747397, 0.00910526, 0.0109764, & + 0.0131071, 0.0155140, 0.0182150, 0.0212261, 0.0245646, & + 0.0282452, 0.0322847, 0.0366982, 0.0415011, 0.0467072, & + 0.0523312, 0.0583875, 0.0648898, 0.0718511, 0.0792848, & + 0.0872044, 0.0956224, 0.104550, 0.114000, 0.123985, & + 0.134518, 0.145609, 0.157270, 0.169514, 0.182353, & + 0.195797, 0.209857, 0.224545, 0.239873, 0.255851, & + 0.272491, 0.289803, 0.307785, 0.326419, 0.345659, & + 0.365435, 0.385673, 0.406299, 0.427248, 0.448454, & + 0.469840, 0.491328, 0.512844, 0.534317, 0.555677, & + 0.576847, 0.597762, 0.618360, 0.638587, 0.658385, & + 0.677706, 0.696509, 0.714758, 0.732423, 0.749476, & + 0.765898, 0.781674, 0.796796, 0.811257, 0.825057, & + 0.838199, 0.850689, 0.862535, 0.873752, 0.884352, & + 0.894353, 0.903769, 0.912623, 0.920933, 0.928720, & + 0.936005, 0.942807, 0.949150, 0.955052, 0.960536, & + 0.965618, 0.970318, 0.974655, 0.978647, 0.982307, & + 0.985649, 0.988688, 0.991432, 0.993898, 0.996109, & + 0.998114, 1.00000 / + ! GFS 127-levels data a127 /1.0000000, 1.6100000, 2.5300000, 3.9200001, 5.9800000, & 8.9499998, 13.179999, 19.100000, 27.240000, 38.280003, & @@ -1165,6 +1246,11 @@ subroutine set_eta_r8_(km, ks, ptop, pint, ak, bk) bk(k) = b181(k) end do + case (186) + do k=1,km+1 + ak(k) = a186(k) + bk(k) = b186(k) + end do end select From 434cd55f672bbeced529d18ae6dce694e6f89c5b Mon Sep 17 00:00:00 2001 From: William Putman Date: Mon, 5 May 2025 12:35:03 -0400 Subject: [PATCH 04/15] added 41 ML levels for PYMLINC --- GMAO_hermes/m_set_eta.F90 | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/GMAO_hermes/m_set_eta.F90 b/GMAO_hermes/m_set_eta.F90 index eb676bd0..8230fbff 100644 --- a/GMAO_hermes/m_set_eta.F90 +++ b/GMAO_hermes/m_set_eta.F90 @@ -92,6 +92,9 @@ subroutine set_eta_r8_(km, ks, ptop, pint, ak, bk) real(r8) a137(138), b137(138) real(r8) a181(182), b181(182) +! ML Increment Levels + integer l181tol41(42) + ! Extended model top real(r8) a186(187), b186(187) @@ -887,6 +890,9 @@ subroutine set_eta_r8_(km, ks, ptop, pint, ak, bk) 0.985649, 0.988688, 0.991432, 0.993898, 0.996109, & 0.998114, 1.00000 / +! For ML levels + data l181tol41 /0,2,5,8,11,14,17,20,23,26,29,32,35,38,41,44,47,50,53,56,60,64,69,74,80,85,90,95,100,105,110,115,120,125,130,135,140,145,152,160,170,181/ + ! Extended model top data a186 /0.00025000, 0.00600000, 0.0500000, 0.20000, 0.5000000,& 1.0000000, 1.6461494, 2.3988098, 3.3373925, 4.5180100, & @@ -1246,6 +1252,13 @@ subroutine set_eta_r8_(km, ks, ptop, pint, ak, bk) bk(k) = b181(k) end do + ! ML levels as a subset of L181 + case (41) + do k=1,km+1 + ak(k) = a181(l181tol41(k)+1) + bk(k) = b181(l181tol41(k)+1) + end do + case (186) do k=1,km+1 ak(k) = a186(k) From 583d816de5a3bdbbb37896061a7d59597ae57ea6 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Tue, 13 May 2025 10:50:28 -0400 Subject: [PATCH 05/15] Update CI --- .circleci/config.yml | 6 +++--- .github/workflows/workflow.yml | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index e6526f6d..a89e348f 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,11 +1,11 @@ version: 2.1 # Anchors in case we need to override the defaults from the orb -#baselibs_version: &baselibs_version v7.33.0 -#bcs_version: &bcs_version v11.6.0 +#baselibs_version: &baselibs_version v8.14.0 +#bcs_version: &bcs_version v12.0.0 orbs: - ci: geos-esm/circleci-tools@4 + ci: geos-esm/circleci-tools@5 workflows: build-test: diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index d84b5092..c94ed43a 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -21,7 +21,7 @@ jobs: if: "!contains(github.event.pull_request.labels.*.name, '0 diff trivial')" runs-on: ubuntu-24.04 container: - image: gmao/ubuntu24-geos-env:v7.33.0-intelmpi_2021.13-ifort_2021.13 + image: gmao/ubuntu24-geos-env:v8.14.0-intelmpi_2021.13-ifort_2021.13 # Per https://github.com/actions/virtual-environments/issues/1445#issuecomment-713861495 # It seems like we might not need secrets on GitHub Actions which is good for forked # pull requests From fa0cef0c0c60f7500709b60ed37e195be127c170 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 20 Oct 2025 12:08:54 -0400 Subject: [PATCH 06/15] v12: Updates for ifx builds --- GMAO_gfio/CMakeLists.txt | 11 +++------- GMAO_gfio/{gfio.f => gfio.F} | 0 GMAO_gfio/{gfioutil.f => gfioutil.F} | 0 .../{m_GFIO_PutFld.f => m_GFIO_PutFld.F} | 0 GMAO_hermes/CMakeLists.txt | 20 +++++++++---------- .../{GFIO_Shave.f90 => GFIO_Shave.F90} | 0 .../{GFIO_remap.f90 => GFIO_remap.F90} | 0 7 files changed, 12 insertions(+), 19 deletions(-) rename GMAO_gfio/{gfio.f => gfio.F} (100%) rename GMAO_gfio/{gfioutil.f => gfioutil.F} (100%) rename GMAO_gfio/{m_GFIO_PutFld.f => m_GFIO_PutFld.F} (100%) rename GMAO_hermes/{GFIO_Shave.f90 => GFIO_Shave.F90} (100%) rename GMAO_hermes/{GFIO_remap.f90 => GFIO_remap.F90} (100%) diff --git a/GMAO_gfio/CMakeLists.txt b/GMAO_gfio/CMakeLists.txt index 5fe73bf6..93c2c641 100644 --- a/GMAO_gfio/CMakeLists.txt +++ b/GMAO_gfio/CMakeLists.txt @@ -3,21 +3,16 @@ string (REPLACE GMAO_gfio_ "" precision ${this}) set (srcs - gfio.f + gfio.F gfioCF.F diffdate.f getdate.f timesubs.f - gfioutil.f + gfioutil.F m_GFIO_GetFld.f - m_GFIO_PutFld.f + m_GFIO_PutFld.F ) -set_source_files_properties (gfio.f PROPERTIES COMPILE_FLAGS ${PP}) -set_source_files_properties (gfioutil.f PROPERTIES COMPILE_FLAGS ${PP}) -set_source_files_properties (m_GFIO_PutFld.f PROPERTIES COMPILE_FLAGS ${PP}) - - # Need to build r4 and r8 versions. See details at bottom esma_add_library (${this} SRCS ${srcs} DEPENDENCIES GMAO_mpeu esmf NetCDF::NetCDF_Fortran) set_property(TARGET ${this} PROPERTY POSITION_INDEPENDENT_CODE ON) diff --git a/GMAO_gfio/gfio.f b/GMAO_gfio/gfio.F similarity index 100% rename from GMAO_gfio/gfio.f rename to GMAO_gfio/gfio.F diff --git a/GMAO_gfio/gfioutil.f b/GMAO_gfio/gfioutil.F similarity index 100% rename from GMAO_gfio/gfioutil.f rename to GMAO_gfio/gfioutil.F diff --git a/GMAO_gfio/m_GFIO_PutFld.f b/GMAO_gfio/m_GFIO_PutFld.F similarity index 100% rename from GMAO_gfio/m_GFIO_PutFld.f rename to GMAO_gfio/m_GFIO_PutFld.F diff --git a/GMAO_hermes/CMakeLists.txt b/GMAO_hermes/CMakeLists.txt index 16a47e3c..7bad2f29 100644 --- a/GMAO_hermes/CMakeLists.txt +++ b/GMAO_hermes/CMakeLists.txt @@ -9,7 +9,7 @@ if ( NOT HERMES_LIGHT ) include_directories (${include_GEOS_Shared}) include_directories (${include_GMAO_gfio_r8}) - + set (srcs m_ana2dyn.F m_const.f90 m_dyn.f90 m_dyn2dyn.f90 m_dynp.F90 m_insitu.F m_interp.F90 m_lp.F m_maph.F90 m_maph_pert.F90 m_mapz.F m_mapz_pert.F90 @@ -17,10 +17,10 @@ if ( NOT HERMES_LIGHT ) m_interpack.F m_interpack_ad.F m_ec_set_eta.F90 m_ecdyn.f90 m_tick.f90 m_dyn_util.F90 m_spline.f90 m_nc_JEDIinc.f90 ) - + esma_add_library(${this} SRCS ${srcs} DEPENDENCIES GMAO_mpeu GMAO_gfio_r8) target_compile_options(${this} PRIVATE ${flags}) - + if (EXTENDED_SOURCE) set_source_files_properties (m_interpack.F PROPERTIES COMPILE_FLAGS ${EXTENDED_SOURCE}) set_source_files_properties (m_interpack_ad.F PROPERTIES COMPILE_FLAGS ${EXTENDED_SOURCE}) @@ -31,13 +31,11 @@ if ( NOT HERMES_LIGHT ) set_source_files_properties (gcmbkg2ana.F PROPERTIES COMPILE_FLAGS ${EXTENDED_SOURCE}) set_source_files_properties (diag2dyn.F PROPERTIES COMPILE_FLAGS ${EXTENDED_SOURCE}) endif () - set_source_files_properties (GFIO_Shave.f90 PROPERTIES COMPILE_FLAGS ${PP}) - set_source_files_properties (GFIO_remap.f90 PROPERTIES COMPILE_FLAGS ${PP}) - - ecbuild_add_executable (TARGET GFIO_remap.x SOURCES GFIO_remap.f90 LIBS ${this} MAPL_cfio_r8) - ecbuild_add_executable (TARGET GFIO_Shave.x SOURCES GFIO_Shave.f90 LIBS ${this} MAPL_cfio_r8) + + ecbuild_add_executable (TARGET GFIO_remap.x SOURCES GFIO_remap.F90 LIBS ${this} MAPL_cfio_r8) + ecbuild_add_executable (TARGET GFIO_Shave.x SOURCES GFIO_Shave.F90 LIBS ${this} MAPL_cfio_r8) ecbuild_add_executable (TARGET lcv2prs.x SOURCES lcv2prs.F90 LIBS ${this} MAPL_cfio_r8 MAPL GEOS_Shared) - + ecbuild_add_executable (TARGET ana5sfc.x SOURCES ana5sfc.F90 LIBS ${this}) ecbuild_add_executable (TARGET diag2dyn.x SOURCES diag2dyn.F LIBS ${this}) ecbuild_add_executable (TARGET drs2dyn.x SOURCES drs2dyn.f90 LIBS ${this}) @@ -76,7 +74,7 @@ if ( NOT HERMES_LIGHT ) ecbuild_add_executable (TARGET dyn_inflate.x SOURCES dyn_inflate.f90 LIBS ${this}) ecbuild_add_executable (TARGET dyn_blob.x SOURCES dyn_blob.f90 LIBS ${this}) ecbuild_add_executable (TARGET dyn_fsens_conv.x SOURCES dyn_fsens_conv.f90 LIBS ${this}) - + install(PROGRAMS cnv2prs.pl ana5sfc echorc.pl DESTINATION bin) else () @@ -84,7 +82,7 @@ else () set (srcs m_set_eta.F90 shared_topo_remap.F90 ) - + esma_add_library(${this} SRCS ${srcs} ) target_compile_options(${this} PRIVATE ${flags}) diff --git a/GMAO_hermes/GFIO_Shave.f90 b/GMAO_hermes/GFIO_Shave.F90 similarity index 100% rename from GMAO_hermes/GFIO_Shave.f90 rename to GMAO_hermes/GFIO_Shave.F90 diff --git a/GMAO_hermes/GFIO_remap.f90 b/GMAO_hermes/GFIO_remap.F90 similarity index 100% rename from GMAO_hermes/GFIO_remap.f90 rename to GMAO_hermes/GFIO_remap.F90 From 3b308af5e84b9da01e16f43164af54b6c5bc71cc Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Fri, 2 Jan 2026 13:09:41 -0500 Subject: [PATCH 07/15] v12: Merge main, update CI --- .circleci/config.yml | 5 ++++- .github/workflows/workflow.yml | 2 ++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index 08b8daaf..31a4f50c 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -21,7 +21,10 @@ workflows: #baselibs_version: *baselibs_version repo: GEOSgcm checkout_fixture: true - mepodevelop: true + # V12 code uses a special branch for now. + fixture_branch: feature/sdrabenh/gcm_v12 + # We comment out this as it will "undo" the fixture_branch + #mepodevelop: true persist_workspace: true # Needs to be true to run fv3/gcm experiment, costs extra # Run AMIP GCM (1 hour, no ExtData) diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 41fd37c7..7696a7d1 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -26,6 +26,7 @@ jobs: compiler: ${{ matrix.compiler }} cmake-build-type: ${{ matrix.build-type }} fixture-repo: GEOS-ESM/GEOSgcm + fixture-ref: feature/sdrabenh/gcm_v12 spack_build: uses: GEOS-ESM/CI-workflows/.github/workflows/spack_gcc_build.yml@project/geosgcm @@ -34,4 +35,5 @@ jobs: BUILDCACHE_TOKEN: ${{ secrets.BUILDCACHE_TOKEN }} with: fixture-repo: GEOS-ESM/GEOSgcm + fixture-ref: feature/sdrabenh/gcm_v12 From 6669be4cd0b12fdfd7d04371a8485f29612b155d Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 4 Mar 2026 15:10:24 -0500 Subject: [PATCH 08/15] Fix for Open MPI 5 --- GMAO_etc/esma_mpirun | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GMAO_etc/esma_mpirun b/GMAO_etc/esma_mpirun index 2860c402..69e83596 100755 --- a/GMAO_etc/esma_mpirun +++ b/GMAO_etc/esma_mpirun @@ -331,12 +331,12 @@ sub get_xtraflags { if ($mpi_type eq "openmpi") { if ($siteID eq "gmao") { - $xtraflags .= " -oversubscribe"; + $xtraflags .= " --oversubscribe"; } if ($perhost) { - $xtraflags .= " -map-by ppr:$perhost:node -bind-to core"; + $xtraflags .= " --map-by ppr:$perhost:node --bind-to core"; } else { - $xtraflags .= " -map-by node -bind-to core"; + $xtraflags .= " --bind-to core"; } } From 27961c6688ad58024f37b1d752f5b045b59b1bea Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Mon, 9 Mar 2026 13:41:02 -0400 Subject: [PATCH 09/15] Fix the openmpi5 fix --- GMAO_etc/esma_mpirun | 2 -- 1 file changed, 2 deletions(-) diff --git a/GMAO_etc/esma_mpirun b/GMAO_etc/esma_mpirun index 69e83596..07732359 100755 --- a/GMAO_etc/esma_mpirun +++ b/GMAO_etc/esma_mpirun @@ -335,8 +335,6 @@ sub get_xtraflags { } if ($perhost) { $xtraflags .= " --map-by ppr:$perhost:node --bind-to core"; - } else { - $xtraflags .= " --bind-to core"; } } From 6ef9e05e6892f8fc2bcb89868cd04a67238026e9 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 11 Mar 2026 10:27:58 -0400 Subject: [PATCH 10/15] v12: Fixes for Flang 22 --- GMAO_mpeu/m_random.F | 78 ++++++++++++++++++++++---------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/GMAO_mpeu/m_random.F b/GMAO_mpeu/m_random.F index 7f0a08cd..f9a8687e 100644 --- a/GMAO_mpeu/m_random.F +++ b/GMAO_mpeu/m_random.F @@ -44,93 +44,93 @@ module m_random * for: uniform (u in [0,1)), normal ( = 0, = 1), and * Poisson distributions. The basic module, the uniform generator, * uses a lagged Fibonacci series generator: -* +* * t = u(n-273) + u(n-607) * u(n) = t - float(int(t)) -* +* * where each number generated, u(k), is floating point. Since * the numbers are floating point, the left end boundary of the * range contains zero. This package is nearly portable except * for the following. (1) It is written in lower case, (2) the * test package contains a timer (second) which is not portable, -* and (3) there are cycle times (in seconds) in data statements -* for NEC SX-3, Fujitsu VP2200, and Cray Y-MP. Select your -* favorite and comment out the others. Replacement functions -* for 'second' are included - comment out the others. Otherwise -* the package is portable and returns the same set of floating -* point numbers up to word precision on any machine. There are -* compiler directives ($cdir for Cray, *vdir for SX-3, and VOCL +* and (3) there are cycle times (in seconds) in data statements +* for NEC SX-3, Fujitsu VP2200, and Cray Y-MP. Select your +* favorite and comment out the others. Replacement functions +* for 'second' are included - comment out the others. Otherwise +* the package is portable and returns the same set of floating +* point numbers up to word precision on any machine. There are +* compiler directives ($cdir for Cray, *vdir for SX-3, and VOCL * for Fujitsu VP2200) which should be otherwise ignored. -* +* * To compile this beast, note that all floating point numbers * are declared 'double precision'. On Cray X-MP, Y-MP, and C-90 * machines, use the cft77 (cf77) option -dp to run this in 64 * bit mode (not 128 bit double). -* +* * External documentation, "Lagged Fibonacci Random Number Generators * for the NEC SX-3," is to be published in the International * Journal of High Speed Computing (1994). Otherwise, ask the -* author: -* -* W. P. Petersen +* author: +* +* W. P. Petersen * IPS, RZ F-5 * ETHZ * CH 8092, Zurich * Switzerland -* +* * e-mail: wpp@ips.ethz.ch. -* +* * The package contains the following routines: -* +* * ------------------------------------------------------ * UNIFORM generator routines: -* +* * subroutine zufalli(seed) * integer seed * c initializes common block containing seeds. if seed=0, * c the default value is 1802. -* +* * subroutine zufall(n,u) * integer n * double precision u(n) * c returns set of n uniforms u(1), ..., u(n). -* +* * subroutine zufallsv(zusave) * double precision zusave(608) * c saves buffer and pointer in zusave, for later restarts -* +* * subroutine zufallrs(zusave) * double precision zusave(608) * c restores seed buffer and pointer from zusave * ------------------------------------------------------ -* +* * NORMAL generator routines: -* +* * subroutine normalen(n,g) * integer n * double precision g(n) * c returns set of n normals g(1), ..., g(n) such that * c mean = 0, and variance = 1. -* +* * subroutine normalsv(normsv) * double precision normsv(1634) * c saves zufall seed buffer and pointer in normsv * c buffer/pointer for normalen restart also in normsv -* +* * subroutine normalrs(normsv) * double precision normsv(1634) -* c restores zufall seed buffer/pointer and +* c restores zufall seed buffer/pointer and * c buffer/pointer for normalen restart from normsv * ------------------------------------------------------ -* +* * POISSON generator routine: -* +* * subroutine fische(n,mu,q) * integer n,q(n) * double precision mu * c returns set of n integers q, with poisson * c distribution, density p(q,mu) = exp(-mu) mu**q/q! -* c +* c * c USE zufallsv and zufallrs for stop/restart sequence * c #ifdef _MAIN_ @@ -545,7 +545,7 @@ double precision function mysecond() c c Generic Unix version real xx(2) -#if defined(sysAIX) || defined(__GFORTRAN__) +#if defined(sysAIX) || defined(__GFORTRAN__) || defined (__flang__) call get_zeits(xx(1)) #else call etime(xx) @@ -631,16 +631,16 @@ subroutine zufall(n,a) else c c more than 1 full segment -c +c do 6 i=1,left a(i+aptr) = buff(ptr+i) 6 continue nn = nn - left ptr = 0 aptr = aptr+left -c +c c buff -> a(aptr0) -c +c VL = 273 k273 = 334 k607 = 0 @@ -765,7 +765,7 @@ end subroutine zufalli subroutine zufallsv(svblk) implicit none c -c saves common blocks klotz0, containing seeds and +c saves common blocks klotz0, containing seeds and c pointer to position in seed block. IMPORTANT: svblk must be c dimensioned at least 608 in driver. The entire contents c of klotz0 (pointer in buff, and buff) must be saved. @@ -810,7 +810,7 @@ subroutine normalen(n,x) double precision x(*) double precision xbuff(1024) integer i,ptr,xptr,first - integer buffsz,nn,n,left + integer buffsz,nn,n,left common /klotz1/xbuff,first,xptr data buffsz/1024/ c @@ -867,8 +867,8 @@ subroutine normalsv(svbox) implicit none c c saves common block klotz0 containing buffers -c and pointers. IMPORTANT: svbox must be dimensioned at -c least 1634 in driver. The entire contents of blocks +c and pointers. IMPORTANT: svbox must be dimensioned at +c least 1634 in driver. The entire contents of blocks c klotz0 (via zufallsv) and klotz1 must be saved. c double precision buff(607) @@ -900,7 +900,7 @@ subroutine normalrs(svbox) implicit none c c restores common blocks klotz0, klotz1 containing buffers -c and pointers. IMPORTANT: svbox must be dimensioned at +c and pointers. IMPORTANT: svbox must be dimensioned at c least 1634 in driver. The entire contents c of klotz0 and klotz1 must be restored. c @@ -946,7 +946,7 @@ subroutine fische(n,mu,p) pmu = exp(-mu) p0 = 0 c - nsegs = (n-1)/1024 + nsegs = (n-1)/1024 left = n - nsegs*1024 nsegs = nsegs + 1 nl0 = left From 8362d17d0809fcf13582b5d250eda45fecda8fa5 Mon Sep 17 00:00:00 2001 From: William Putman Date: Tue, 14 Apr 2026 11:13:47 -0400 Subject: [PATCH 11/15] updated iliqskin to apply to just snow or snow & graupel as in SFE25 --- GEOS_Shared/calcdbz.F | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/GEOS_Shared/calcdbz.F b/GEOS_Shared/calcdbz.F index 88cbd3a6..8061bed3 100644 --- a/GEOS_Shared/calcdbz.F +++ b/GEOS_Shared/calcdbz.F @@ -25,7 +25,10 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, c calculated in Thompson, Rasmussen, and Manning (2004, Monthly c Weather Review, Vol. 132, No. 2, pp. 519-542.) c -c 4. If iliqskin=1, frozen particles that are at a temperature above +c 4. If iliqskin=1, snow particles that are at a temperature above +c freezing are assumed to scatter as a liquid particle. +c +c 4. If iliqskin=2, snow & graupel particles that are at a temperature above c freezing are assumed to scatter as a liquid particle. c c More information on the derivation of simulated reflectivity in @@ -51,7 +54,7 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, REAL TEMP_C,VIRTUAL_T REAL GONV,RONV,SONV REAL FACTOR_G,FACTOR_R,FACTOR_S - REAL FACTORB_G,FACTORB_R,FACTORB_S + REAL FACTORB_G,FACTORB_S REAL RHOAIR,Z_E c Constants used to calculate variable intercepts @@ -150,6 +153,8 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, IF (ILIQSKIN.EQ.1 .AND. TMK(I,J,K).GT.CELKEL) THEN FACTORB_S = FACTOR_S/ALPHA + ELSEIF (ILIQSKIN.EQ.2 .AND. TMK(I,J,K).GT.CELKEL) THEN + FACTORB_S = FACTOR_S/ALPHA FACTORB_G = FACTOR_G/ALPHA ELSE FACTORB_S = FACTOR_S From adaf9a73da8c56b14a64db964ade239fbce18670 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Fri, 1 May 2026 08:37:54 -0400 Subject: [PATCH 12/15] Update spack CI --- .github/workflows/workflow.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 7696a7d1..df000662 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -18,6 +18,7 @@ concurrency: jobs: build_gcm: strategy: + fail-fast: false matrix: compiler: [ifort, gfortran-14, gfortran-15] build-type: [Debug] @@ -36,4 +37,5 @@ jobs: with: fixture-repo: GEOS-ESM/GEOSgcm fixture-ref: feature/sdrabenh/gcm_v12 + load-fms: true From 9e2f1cfbc33b9fe72ba3bd9e0bc969cdeb6a8526 Mon Sep 17 00:00:00 2001 From: William Putman Date: Wed, 6 May 2026 16:54:31 -0400 Subject: [PATCH 13/15] minimal zero-diff for L72 merge from latest v12-SFE26 code --- GEOS_Shared/CMakeLists.txt | 2 +- GEOS_Shared/calcdbz.F | 58 +++++++++++++++++++++++++++++--------- 2 files changed, 45 insertions(+), 15 deletions(-) diff --git a/GEOS_Shared/CMakeLists.txt b/GEOS_Shared/CMakeLists.txt index e5cb6a83..9e76e46b 100644 --- a/GEOS_Shared/CMakeLists.txt +++ b/GEOS_Shared/CMakeLists.txt @@ -4,7 +4,7 @@ set (srcs GEOS_Utilities.F90 DragCoefficients.F90 GEOS.F90 tropovars.F90 surfacelayer.F90 GEOS_TopoGet.F90 GEOS_GmapMod.F90 GEOS_RemapMod.F90 G3_MPI_Util_Mod.F90 windfix.F90 - PWSSSP.F sphere.F calcdbz.F OVP.F90 + PWSSSP.F sphere.F calcdbz.F compute_radar_reflectivity.F90 OVP.F90 getco2.F90 atmOceanIntLayer.F90 ) diff --git a/GEOS_Shared/calcdbz.F b/GEOS_Shared/calcdbz.F index 8061bed3..0dddac6e 100644 --- a/GEOS_Shared/calcdbz.F +++ b/GEOS_Shared/calcdbz.F @@ -56,6 +56,8 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, REAL FACTOR_G,FACTOR_R,FACTOR_S REAL FACTORB_G,FACTORB_S REAL RHOAIR,Z_E + REAL ZE_R,ZE_S,ZE_G + c Constants used to calculate variable intercepts REAL R1,RON,RON2,SON,GON @@ -68,35 +70,38 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, REAL GAMMA_SEVEN,ALPHA REAL RHOWAT,CELKEL,PI,RD - c Constants used to calculate variable intercepts R1 = 1.D-15 RON = 8.D6 RON2 = 1.D10 SON = 2.D7 - GON = 5.D7 RON_MIN = 8.D6 RON_QR0 = 0.00010D0 RON_DELQR0 = 0.25D0*RON_QR0 RON_CONST1R = (RON2-RON_MIN)*0.5D0 RON_CONST2R = (RON2+RON_MIN)*0.5D0 -c Constant intercepts - RN0_R = 8.D6 - RN0_S = 2.D7 - RN0_G = 4.D6 - c Other constants GAMMA_SEVEN = 720.D0 RHOWAT = 1000.D0 RHO_R = RHOWAT RHO_S = 100.D0 - RHO_G = 400.D0 ALPHA = 0.224D0 CELKEL = 273.15D0 PI = 3.141592653589793D0 RD = 287.04D0 +c --- NEW: HAIL VS GRAUPEL CONSTANTS --- + IF (IVARINT .GE. 3) THEN + RHO_G = 917.D0 ! Hail density (kg/m^3) + GON = 4.D4 ! Hail max intercept (fewer, much larger stones) + RN0_G = 4.D4 ! Constant hail intercept + ELSE + RHO_G = 400.D0 ! Graupel density (kg/m^3) + GON = 5.D7 ! Graupel max intercept (many, smaller stones) + RN0_G = 4.D6 ! Constant graupel intercept + END IF + c Force all Q arrays to be 0.0 or greater. DO K = 1,BTDIM DO J = 1,SNDIM @@ -163,11 +168,19 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, c Calculate variable intercept parameters - IF (IVARINT.EQ.1) THEN + IF (IVARINT.GE.1) THEN TEMP_C = DMIN1(-0.001D0,TMK(I,J,K)-CELKEL) SONV = DMIN1(2.0D8,2.0D6*EXP(-0.12D0*TEMP_C)) +c simulate large snow aggregates by halving the intercept parameter. +c Apply this to both option 2 and 3! + IF (IVARINT.GE.2) THEN + IF (QSN(I,J,K).GT.R1) THEN + SONV = SONV * 0.5D0 + END IF + END IF + GONV = GON IF (QGR(I,J,K).GT.R1) THEN GONV = 2.38D0* (PI*RHO_G/ @@ -189,14 +202,31 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, END IF + c Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is c the sum of z_e for each hydrometeor species: - Z_E = FACTOR_R* (RHOAIR*QRA(I,J,K))**1.75D0/ - + RONV**.75D0 + FACTORB_S* - + (RHOAIR*QSN(I,J,K))**1.75D0/SONV**.75D0 + - + FACTORB_G* (RHOAIR*QGR(I,J,K))**1.75D0/ - + GONV**.75D0 +c Calculate individual equivalent reflectivity factors + ZE_R = (RHOAIR*QRA(I,J,K))**1.75D0 + ZE_R = FACTOR_R * ZE_R / RONV**.75D0 + + ZE_S = (RHOAIR*QSN(I,J,K))**1.75D0 + ZE_S = FACTORB_S * ZE_S / SONV**.75D0 + + ZE_G = (RHOAIR*QGR(I,J,K))**1.75D0 + ZE_G = FACTORB_G * ZE_G / GONV**.75D0 + +c --- PSEUDO-MIE PLATEAU FOR WET HAIL --- +c If using Hail (IVARINT>=3) with a liquid skin, Rayleigh math overestimates. +c Cap the hail's specific Ze contribution at 10^7 (approx 70 dBZ). + IF (IVARINT.GE.3 .AND. ILIQSKIN.EQ.2 .AND. + + TMK(I,J,K).GT.CELKEL) THEN + ZE_G = MIN(ZE_G, 1.0D7) + END IF + +c Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is +c the sum of z_e for each hydrometeor species: + Z_E = ZE_R + ZE_S + ZE_G c Adjust small values of Z_e so that dBZ is no lower than -30 Z_E = MAX(Z_E,.001D0) From af70d1ee6dc7842ac3a2954d6934bb109c346c32 Mon Sep 17 00:00:00 2001 From: William Putman Date: Tue, 12 May 2026 16:17:06 -0400 Subject: [PATCH 14/15] added new radar reflectivity module --- GEOS_Shared/compute_radar_reflectivity.F90 | 247 +++++++++++++++++++++ 1 file changed, 247 insertions(+) create mode 100644 GEOS_Shared/compute_radar_reflectivity.F90 diff --git a/GEOS_Shared/compute_radar_reflectivity.F90 b/GEOS_Shared/compute_radar_reflectivity.F90 new file mode 100644 index 00000000..b11f7391 --- /dev/null +++ b/GEOS_Shared/compute_radar_reflectivity.F90 @@ -0,0 +1,247 @@ +module GEOS_RadarMod + + implicit none + + private + public :: compute_radar_reflectivity + +contains + +FUNCTION compute_radar_reflectivity(PRS, TMK, QVP, QRAIN, QSNOW, & + QGRAUPEL, QHAIL, disable_variable_intercept_params, & + liqskin_snow, liqskin_graupel, liqskin_hail) RESULT(DBZ) + + IMPLICIT NONE + + ! State and Required Hydrometeor Arguments (1D column slices) + REAL, INTENT(IN) :: PRS(:) + REAL, INTENT(IN) :: TMK(:) + REAL, INTENT(INOUT) :: QVP(:) + REAL, INTENT(INOUT) :: QRAIN(:) + REAL, INTENT(INOUT) :: QSNOW(:) + + ! Optional Hydrometeor Arguments (1D column slices) + REAL, INTENT(INOUT), OPTIONAL :: QGRAUPEL(:) + REAL, INTENT(INOUT), OPTIONAL :: QHAIL(:) + + ! Optional Configuration Arguments + LOGICAL, INTENT(IN), OPTIONAL :: disable_variable_intercept_params + LOGICAL, INTENT(IN), OPTIONAL :: liqskin_snow + LOGICAL, INTENT(IN), OPTIONAL :: liqskin_graupel + LOGICAL, INTENT(IN), OPTIONAL :: liqskin_hail + + ! Return Array (Sized automatically based on the vertical column size) + REAL :: DBZ(SIZE(PRS, 1)) + + ! Local Configuration Variables + LOGICAL :: use_var_intercepts + LOGICAL :: has_graupel, has_hail + LOGICAL :: do_liq_snow, do_liq_graupel, do_liq_hail + + ! Local Variables + INTEGER :: K, LDIM + REAL :: tmk_val, prs_val, qv_val, qr_val, qs_val, qg_val, qh_val + REAL :: TEMP_C, VIRTUAL_T, RHOAIR, Z_E + REAL :: GONV, HONV, RONV, SONV + REAL :: FACTOR_G, FACTOR_H, FACTOR_R, FACTOR_S + REAL :: FACTOR_G_LIQ, FACTOR_H_LIQ, FACTOR_S_LIQ + REAL :: FACTORB_G, FACTORB_H, FACTORB_S + REAL :: ZE_R, ZE_S, ZE_G, ZE_H + REAL :: PI_RHO_G, PI_RHO_H, INV_RON_DELQR0 + + ! Constants used to calculate variable intercepts + REAL :: R1, RON, RON2, SON + REAL :: GON_CONST, HON_CONST + REAL :: RON_MIN, RON_QR0, RON_DELQR0 + REAL :: RON_CONST1R, RON_CONST2R + + ! Constant intercepts + REAL :: RN0_R, RN0_S, RN0_G, RN0_H + + ! Other constants + REAL :: RHO_R, RHO_S, RHO_G, RHO_H + REAL :: GAMMA_SEVEN, ALPHA, RHOWAT, CELKEL, PI, RD + + ! Extract vertical dimension dynamically + LDIM = SIZE(PRS, 1) + + ! Evaluate OPTIONAL arguments exactly once + has_graupel = PRESENT(QGRAUPEL) + has_hail = PRESENT(QHAIL) + + ! Variable intercepts default to ON + use_var_intercepts = .TRUE. + IF (PRESENT(disable_variable_intercept_params)) THEN + IF (disable_variable_intercept_params) use_var_intercepts = .FALSE. + END IF + + ! Liquid skin defaults to OFF + do_liq_snow = .FALSE. + do_liq_graupel = .FALSE. + do_liq_hail = .FALSE. + IF (PRESENT(liqskin_snow)) do_liq_snow = liqskin_snow + IF (PRESENT(liqskin_graupel)) do_liq_graupel = liqskin_graupel + IF (PRESENT(liqskin_hail)) do_liq_hail = liqskin_hail + + ! Constants + R1 = 1.D-15 + RON = 8.D6 + RON2 = 1.D10 + SON = 2.D7 + RON_MIN = 8.D6 + RON_QR0 = 0.00010D0 + RON_DELQR0 = 0.25D0 * RON_QR0 + RON_CONST1R = (RON2 - RON_MIN) * 0.5D0 + RON_CONST2R = (RON2 + RON_MIN) * 0.5D0 + + ! Precompute Division Constant + INV_RON_DELQR0 = 1.D0 / RON_DELQR0 + + GAMMA_SEVEN = 720.D0 + RHOWAT = 1000.D0 + RHO_R = RHOWAT + RHO_S = 100.D0 + ALPHA = 0.224D0 + CELKEL = 273.15D0 + PI = 3.141592653589793D0 + RD = 287.04D0 + + ! Graupel Constants + RHO_G = 400.D0 + GON_CONST = 5.D7 + RN0_G = 4.D6 + PI_RHO_G = PI * RHO_G + + ! Hail Constants + RHO_H = 917.D0 + HON_CONST = 4.D4 + RN0_H = 4.D4 + PI_RHO_H = PI * RHO_H + + ! Precalculate Base Factors (hoisted outside loops) + FACTOR_R = GAMMA_SEVEN * 1.D18 * (1.D0 / (PI * RHO_R))**1.75D0 + FACTOR_S = GAMMA_SEVEN * 1.D18 * (1.D0 / (PI * RHO_S))**1.75D0 * & + (RHO_S / RHOWAT)**2 * ALPHA + FACTOR_G = GAMMA_SEVEN * 1.D18 * (1.D0 / (PI_RHO_G))**1.75D0 * & + (RHO_G / RHOWAT)**2 * ALPHA + FACTOR_H = GAMMA_SEVEN * 1.D18 * (1.D0 / (PI_RHO_H))**1.75D0 * & + (RHO_H / RHOWAT)**2 * ALPHA + + ! Precalculate Liquid Skin Division Factors + FACTOR_S_LIQ = FACTOR_S / ALPHA + FACTOR_G_LIQ = FACTOR_G / ALPHA + FACTOR_H_LIQ = FACTOR_H / ALPHA + + ! Main Compute Loop - Pure Sequential over 1D Column + DO K = 1, LDIM + + ! 1. Apply negative bounds checks AND extract to local scalars + IF (QVP(K) < 0.0) QVP(K) = 0.0 + IF (QRAIN(K) < 0.0) QRAIN(K) = 0.0 + IF (QSNOW(K) < 0.0) QSNOW(K) = 0.0 + + qv_val = QVP(K) + qr_val = QRAIN(K) + qs_val = QSNOW(K) + tmk_val = TMK(K) + prs_val = PRS(K) + + ! Initialize optional species + ZE_G = 0.D0 + ZE_H = 0.D0 + + ! Thermodynamics + VIRTUAL_T = tmk_val * (0.622D0 + qv_val) / (0.622D0 * (1.D0 + qv_val)) + RHOAIR = prs_val / (RD * VIRTUAL_T) + + ! 2. Brightband Factor Assignment + FACTORB_S = FACTOR_S + FACTORB_G = FACTOR_G + FACTORB_H = FACTOR_H + + IF (tmk_val > CELKEL) THEN + IF (do_liq_snow) FACTORB_S = FACTOR_S_LIQ + IF (do_liq_graupel) FACTORB_G = FACTOR_G_LIQ + IF (do_liq_hail) FACTORB_H = FACTOR_H_LIQ + END IF + + ! 3. Variable Intercept Parameters + IF (use_var_intercepts) THEN + TEMP_C = DMIN1(-0.001D0, tmk_val - CELKEL) + SONV = DMIN1(2.0D8, 2.0D6 * EXP(-0.12D0 * TEMP_C)) + + ! Simulate large snow aggregates by halving the intercept parameter + IF (qs_val > R1) THEN + SONV = SONV * 0.5D0 + END IF + + RONV = RON2 + IF (qr_val > R1) THEN + RONV = RON_CONST1R * TANH((RON_QR0 - qr_val) * & + INV_RON_DELQR0) + RON_CONST2R + END IF + ELSE + RONV = RN0_R + SONV = RN0_S + END IF + + ! Calculate Required equivalent reflectivity factors + ZE_R = (RHOAIR * qr_val)**1.75D0 + ZE_R = FACTOR_R * ZE_R / RONV**.75D0 + + ZE_S = (RHOAIR * qs_val)**1.75D0 + ZE_S = FACTORB_S * ZE_S / SONV**.75D0 + + ! 4. Graupel Processing + IF (has_graupel) THEN + IF (QGRAUPEL(K) < 0.0) QGRAUPEL(K) = 0.0 + qg_val = QGRAUPEL(K) + + IF (use_var_intercepts) THEN + GONV = GON_CONST + IF (qg_val > R1) THEN + GONV = 2.38D0 * (PI_RHO_G / (RHOAIR * qg_val))**0.92D0 + GONV = MAX(1.D4, MIN(GONV, GON_CONST)) + END IF + ELSE + GONV = RN0_G + END IF + + ZE_G = (RHOAIR * qg_val)**1.75D0 + ZE_G = FACTORB_G * ZE_G / GONV**.75D0 + END IF + + ! 5. Hail Processing + IF (has_hail) THEN + IF (QHAIL(K) < 0.0) QHAIL(K) = 0.0 + qh_val = QHAIL(K) + + IF (use_var_intercepts) THEN + HONV = HON_CONST + IF (qh_val > R1) THEN + HONV = 2.38D0 * (PI_RHO_H / (RHOAIR * qh_val))**0.92D0 + HONV = MAX(1.D4, MIN(HONV, HON_CONST)) + END IF + ELSE + HONV = RN0_H + END IF + + ZE_H = (RHOAIR * qh_val)**1.75D0 + ZE_H = FACTORB_H * ZE_H / HONV**.75D0 + + ! --- PSEUDO-MIE PLATEAU FOR WET HAIL --- + IF (do_liq_hail .AND. tmk_val > CELKEL) THEN + ZE_H = MIN(ZE_H, 1.0D7) + END IF + END IF + + ! Total Z_e + Z_E = MAX(ZE_R + ZE_S + ZE_G + ZE_H, .001D0) + + ! Convert to dBZ + DBZ(K) = 10.D0 * LOG10(Z_E) + + END DO + +END FUNCTION compute_radar_reflectivity +end module GEOS_RadarMod From 5d77fa0388e39fecfe4e2b3371dc7616614a8896 Mon Sep 17 00:00:00 2001 From: William Putman Date: Wed, 13 May 2026 17:56:25 -0400 Subject: [PATCH 15/15] bug fixed in updated calcdbz when ivarint==0 --- GEOS_Shared/calcdbz.F | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/GEOS_Shared/calcdbz.F b/GEOS_Shared/calcdbz.F index 0dddac6e..048c0b17 100644 --- a/GEOS_Shared/calcdbz.F +++ b/GEOS_Shared/calcdbz.F @@ -81,6 +81,10 @@ SUBROUTINE CALCDBZ(DBZ,PRS,TMK,QVP,QRA,QSN,QGR,WEDIM,SNDIM,BTDIM, RON_CONST1R = (RON2-RON_MIN)*0.5D0 RON_CONST2R = (RON2+RON_MIN)*0.5D0 +c Constant intercepts + RN0_R = 8.D6 + RN0_S = 2.D7 + c Other constants GAMMA_SEVEN = 720.D0 RHOWAT = 1000.D0