From 71df8678e98fb6db38fda5f68348ebf9077ac902 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 17 Aug 2022 18:01:47 -0600 Subject: [PATCH 1/7] Fix memory leak in deallocate_microphysics for recloud_p, reice_p, and resnow_p In the deallocate_microphysics routine, the logic for deallocating the recloud_p, reice_p, and resnow_p fields was incorrect, leading to a memory leak at the end of a simulation. This commit corrects that logic, thereby avoiding a memory leak. --- .../physics/mpas_atmphys_driver_microphysics.F | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index 1d6e3235e6..bcc1e155a4 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -214,9 +214,9 @@ subroutine deallocate_microphysics(configs) if(allocated(graupelncv_p) ) deallocate(graupelncv_p ) !cloud water,cloud ice,and snow effective radii: - if(.not.allocated(recloud_p) ) allocate(recloud_p(ims:ime,kms:kme,jms:jme) ) - if(.not.allocated(reice_p) ) allocate(reice_p(ims:ime,kms:kme,jms:jme) ) - if(.not.allocated(resnow_p) ) allocate(resnow_p(ims:ime,kms:kme,jms:jme) ) + if(allocated(recloud_p) ) deallocate(recloud_p ) + if(allocated(reice_p) ) deallocate(reice_p ) + if(allocated(resnow_p) ) deallocate(resnow_p ) microp2_select: select case(microp_scheme) From 9e015db1ed9018780596b12956e997119c8d6ded Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 23 Jun 2023 13:43:05 -0600 Subject: [PATCH 2/7] Update version number to 8.0.1 --- README.md | 2 +- src/core_atmosphere/Registry.xml | 2 +- src/core_init_atmosphere/Registry.xml | 2 +- src/core_landice/Registry.xml | 2 +- src/core_ocean/Registry.xml | 2 +- src/core_seaice/Registry.xml | 2 +- src/core_sw/Registry.xml | 2 +- src/core_test/Registry.xml | 2 +- 8 files changed, 8 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 4b67663698..5765206560 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -MPAS-v8.0.0 +MPAS-v8.0.1 ==== The Model for Prediction Across Scales (MPAS) is a collaborative project for diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index c01e9e3db7..173b483774 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1,5 +1,5 @@ - + diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index 7ac677d7d0..e40f5c3722 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -1,5 +1,5 @@ - + diff --git a/src/core_landice/Registry.xml b/src/core_landice/Registry.xml index e57490356f..42e2a99254 100644 --- a/src/core_landice/Registry.xml +++ b/src/core_landice/Registry.xml @@ -1,5 +1,5 @@ - + diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml index c64049a5ea..c084373916 100644 --- a/src/core_ocean/Registry.xml +++ b/src/core_ocean/Registry.xml @@ -1,5 +1,5 @@ - + - + - + diff --git a/src/core_test/Registry.xml b/src/core_test/Registry.xml index 9d2aeeced3..f2ffbbf18d 100644 --- a/src/core_test/Registry.xml +++ b/src/core_test/Registry.xml @@ -1,5 +1,5 @@ - + From 10f63eed1158a45e144094ec514ab02c11a864df Mon Sep 17 00:00:00 2001 From: Manda Chasteen Date: Wed, 10 Aug 2022 15:02:06 -0600 Subject: [PATCH 3/7] Fix PV registry descriptions In the Registry, the pv_vertex, pv_edge, and pv_cell variables are described as "absolute vorticity/rho_zz", but the density component was previously removed from the model code. This request updates the descriptions and units of these variables to be consistent with the absolute vertical vorticity. In Line 4816 of mpas_atm_time_integration.F (v7.3): "the original definition of pv_edge had a factor of 1/density. We have removed that factor given that it was not integral to any conservation property of the system" Relative vertical vorticity is calculated beginning at Line 5606 (in v7.3) as: do iVertex=vertexStart,vertexEnd vorticity(1:nVertLevels,iVertex) = 0.0 do i=1,vertexDegree iEdge = edgesOnVertex(i,iVertex) s = edgesOnVertex_sign(i,iVertex) * dcEdge(iEdge) !DIR$ IVDEP do k=1,nVertLevels vorticity(k,iVertex) = vorticity(k,iVertex) + s * u(k,iEdge) end do end do !DIR$ IVDEP do k=1,nVertLevels vorticity(k,iVertex) = vorticity(k,iVertex) * invAreaTriangle(iVertex) end do end do (with 'u' as the uncoupled normal component of horizontal velocity), and planetary vorticity is added to get absolute vorticity at Line 5754: pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) Thus, pv_vertex is the absolute vertical vorticity at a vertex. --- src/core_atmosphere/Registry.xml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index e4cf270af4..72072c058a 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1565,8 +1565,8 @@ - + @@ -1574,11 +1574,11 @@ - + - + Date: Wed, 28 Jun 2023 09:57:13 -0600 Subject: [PATCH 4/7] Fix units and description for rt_diabatic_tend in atmosphere core Registry.xml The rt_diabatic_tend field represents the tendency of the modified potential temperature due to cloud microphysics in K/s; accordingly, this commit corrects the units and description attributes for rt_diabatic_tend in the atmosphere core's Registry.xml file. --- src/core_atmosphere/Registry.xml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index c01e9e3db7..1707e60c4e 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1791,8 +1791,8 @@ - + From 4a2f1bfb4e586bd173712b2166861cc03d130dc8 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Tue, 27 Jun 2023 17:33:01 -0600 Subject: [PATCH 5/7] Allow real-valued attributes to be read with different precision using SMIOL Previously, the mpas_io_get_att_real0d() routine relied on its 'precision' optional argument to determine whether to attempt to read real-valued attributes in a precision different from MPAS's native precision. However, in some parts of the MPAS framework (specifically, in mpas_bootstrap_framework_phase1()), the MPAS_io_get_att() routine is called without the 'precision' argument, leading to problems in reading, e.g., 'sphere_radius' when the precision of the input file is not the MPAS native precision. Now, the mpas_io_get_att_real0d() routine first tries to read an attribute in MPAS's native precision; if that fails, a second attempt is made to read the attribute using the opposite precision (single vs. double). If both attempts fail, the attribute argument is set to a fill value and an error code is returned. This allows MPAS to read real-valued attributes using SMIOL during the bootstrapping process regardless of the precision of the input file. --- src/framework/mpas_io.F | 89 +++++++++++++++++++++++++---------------- 1 file changed, 54 insertions(+), 35 deletions(-) diff --git a/src/framework/mpas_io.F b/src/framework/mpas_io.F index 01ab167243..f699eae7fc 100644 --- a/src/framework/mpas_io.F +++ b/src/framework/mpas_io.F @@ -4603,77 +4603,45 @@ subroutine MPAS_io_get_att_real0d(handle, attName, attValue, fieldname, precisio local_precision = MPAS_IO_NATIVE_PRECISION end if - ! Query attribute value #ifdef MPAS_PIO_SUPPORT + ! Query attribute value pio_ierr = PIO_inq_att(handle % pio_file, varid, attName, xtype, len) if (pio_ierr /= PIO_noerr) then if (present(ierr)) ierr = MPAS_IO_ERR_PIO return end if -#endif if ((local_precision == MPAS_IO_SINGLE_PRECISION) .and. & (MPAS_IO_NATIVE_PRECISION /= MPAS_IO_SINGLE_PRECISION)) then -#ifdef MPAS_PIO_SUPPORT if (xtype /= PIO_REAL) then if (present(ierr)) ierr=MPAS_IO_ERR_WRONG_ATT_TYPE return end if pio_ierr = PIO_get_att(handle % pio_file, varid, attName, singleVal) -#endif - -#ifdef MPAS_SMIOL_SUPPORT - if (present(fieldname)) then - local_ierr = SMIOLf_inquire_att(handle % smiol_file, trim(fieldname), trim(attName), singleVal) - else - local_ierr = SMIOLf_inquire_att(handle % smiol_file, '', trim(attName), singleVal) - end if -#endif attValue = real(singleVal,RKIND) else if ((local_precision == MPAS_IO_DOUBLE_PRECISION) .and. & (MPAS_IO_NATIVE_PRECISION /= MPAS_IO_DOUBLE_PRECISION)) then -#ifdef MPAS_PIO_SUPPORT if (xtype /= PIO_DOUBLE) then if (present(ierr)) ierr=MPAS_IO_ERR_WRONG_ATT_TYPE return end if pio_ierr = PIO_get_att(handle % pio_file, varid, attName, doubleVal) -#endif - -#ifdef MPAS_SMIOL_SUPPORT - if (present(fieldname)) then - local_ierr = SMIOLf_inquire_att(handle % smiol_file, trim(fieldname), trim(attName), doubleVal) - else - local_ierr = SMIOLf_inquire_att(handle % smiol_file, '', trim(attName), doubleVal) - end if -#endif attValue = real(doubleVal,RKIND) else -#ifdef MPAS_PIO_SUPPORT if (xtype /= PIO_DOUBLE .and. xtype /= PIO_REAL) then if (present(ierr)) ierr=MPAS_IO_ERR_WRONG_ATT_TYPE return end if pio_ierr = PIO_get_att(handle % pio_file, varid, attName, attValue) -#endif - -#ifdef MPAS_SMIOL_SUPPORT - if (present(fieldname)) then - local_ierr = SMIOLf_inquire_att(handle % smiol_file, trim(fieldname), trim(attName), attValue) - else - local_ierr = SMIOLf_inquire_att(handle % smiol_file, '', trim(attName), attValue) - end if -#endif end if -#ifdef MPAS_PIO_SUPPORT if (pio_ierr /= PIO_noerr) then if (present(ierr)) ierr = MPAS_IO_ERR_PIO return @@ -4681,14 +4649,65 @@ subroutine MPAS_io_get_att_real0d(handle, attName, attValue, fieldname, precisio #endif #ifdef MPAS_SMIOL_SUPPORT + ! + ! Try to read the attribute in the MPAS native precision + ! + if (present(fieldname)) then + local_ierr = SMIOLf_inquire_att(handle % smiol_file, trim(fieldname), & + trim(attName), attValue) + else + local_ierr = SMIOLf_inquire_att(handle % smiol_file, '', & + trim(attName), attValue) + end if + + ! + ! If that fails, perhaps the attribute is in a different precision from + ! the native MPAS precision + ! + if (local_ierr == SMIOL_WRONG_ARG_TYPE) then + if (MPAS_IO_NATIVE_PRECISION == MPAS_IO_DOUBLE_PRECISION) then + + ! + ! Try again, but read a single-precision value + ! + if (present(fieldname)) then + local_ierr = SMIOLf_inquire_att(handle % smiol_file, trim(fieldname), & + trim(attName), singleVal) + else + local_ierr = SMIOLf_inquire_att(handle % smiol_file, '', & + trim(attName), singleVal) + end if + attValue = real(singleVal,RKIND) + + else + + ! + ! Try again, but read a double-precision value + ! + if (present(fieldname)) then + local_ierr = SMIOLf_inquire_att(handle % smiol_file, trim(fieldname), & + trim(attName), doubleVal) + else + local_ierr = SMIOLf_inquire_att(handle % smiol_file, '', & + trim(attName), doubleVal) + end if + attValue = real(doubleVal,RKIND) + + end if + end if + + ! + ! If all of the above were unsuccessful, set attValue to a fill value + ! and return an error + ! if (local_ierr /= SMIOL_SUCCESS) then + attValue = MPAS_REAL_FILLVAL if (local_ierr == SMIOL_WRONG_ARG_TYPE) then if (present(ierr)) ierr = MPAS_IO_ERR_WRONG_ATT_TYPE - return else if (present(ierr)) ierr = MPAS_IO_ERR_PIO - return end if + return end if #endif From cf002350d265886d4fd254dc85ab1f8ef688ec75 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Thu, 29 Jun 2023 13:17:30 -0600 Subject: [PATCH 6/7] Add missing physics_mmm module search path in atmosphere Makefiles Some compilers (e.g., nvfortran) require the paths of modules that are used indirectly to be present in the module search path provided with the -I... option. This commit adds -I./physics/physics_mmm to the set of compiler flags in src/core_atmosphere/Makefile and -I../physics/physics_mmm to the set of compiler flags in src/core_atmosphere/dynamics/Makefile to fix compilation errors with certain Fortran compilers. --- src/core_atmosphere/Makefile | 4 ++-- src/core_atmosphere/dynamics/Makefile | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core_atmosphere/Makefile b/src/core_atmosphere/Makefile index c03e2fc3d0..60ce1b2703 100644 --- a/src/core_atmosphere/Makefile +++ b/src/core_atmosphere/Makefile @@ -79,7 +79,7 @@ clean: $(RM) $@ $*.mod ifeq "$(GEN_F90)" "true" $(CPP) $(CPPFLAGS) $(PHYSICS) $(CPPINCLUDES) -I./inc $< > $*.f90 - $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators -I./physics -I./dynamics -I./diagnostics -I./physics/physics_wrf -I../external/esmf_time_f90 + $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators -I./physics -I./dynamics -I./diagnostics -I./physics/physics_wrf -I./physics/physics_mmm -I../external/esmf_time_f90 else - $(FC) $(CPPFLAGS) $(PHYSICS) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I./inc -I../framework -I../operators -I./physics -I./dynamics -I./diagnostics -I./physics/physics_wrf -I../external/esmf_time_f90 + $(FC) $(CPPFLAGS) $(PHYSICS) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I./inc -I../framework -I../operators -I./physics -I./dynamics -I./diagnostics -I./physics/physics_wrf -I./physics/physics_mmm -I../external/esmf_time_f90 endif diff --git a/src/core_atmosphere/dynamics/Makefile b/src/core_atmosphere/dynamics/Makefile index 732c12926e..6892633c68 100644 --- a/src/core_atmosphere/dynamics/Makefile +++ b/src/core_atmosphere/dynamics/Makefile @@ -20,7 +20,7 @@ clean: $(RM) $@ $*.mod ifeq "$(GEN_F90)" "true" $(CPP) $(CPPFLAGS) $(PHYSICS) $(CPPINCLUDES) $< > $*.f90 - $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I.. -I../../framework -I../../operators -I../physics -I../physics/physics_wrf -I../../external/esmf_time_f90 + $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I.. -I../../framework -I../../operators -I../physics -I../physics/physics_wrf -I../physics/physics_mmm -I../../external/esmf_time_f90 else - $(FC) $(CPPFLAGS) $(PHYSICS) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I.. -I../../framework -I../../operators -I../physics -I../physics/physics_wrf -I../../external/esmf_time_f90 + $(FC) $(CPPFLAGS) $(PHYSICS) $(FFLAGS) -c $*.F $(CPPINCLUDES) $(FCINCLUDES) -I.. -I../../framework -I../../operators -I../physics -I../physics/physics_wrf -I../physics/physics_mmm -I../../external/esmf_time_f90 endif From b8f7757d998c9f83b40cc29f4a94b55d81956b69 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 23 Jun 2023 16:49:44 -0600 Subject: [PATCH 7/7] Fix OpenMP error in deallocation of rthdynten without certain cumulus schemes When neither the Grell-Freitas nor the Tiedtke/nTiedtke cumulus schemes were used, the rthdynten array was allocated locally in the compute_dyn_tend(...) routine, since it is a required argument to the compute_dyn_tend_work(...) routine. The local rthdynten pointer was initialized and therefore implicitly had the SAVE attribute, which, according to 2.15.1.2 of the OpenMP 4.5 standard, made the rthdynten pointer shared: Local variables declared in called routines in the region and that have the save attribute, or that are data initialized, are shared unless they appear in a threadprivate directive. The allocation and deallocation of rthdynten happened inside of a threaded region, leading to potential model failures when multiple threads attempted to deallocate the same shared rthdynten pointer. The simple solution is to remove packages from the rthdynten field in the atmosphere core's Registry.xml file, and to remove the additional code that handled selective local allocation of this array in the compute_dyn_tend(...) routine. Since rthdynten is allocated by the compute_dyn_tend(...) routine if it is not otherwise allocated due to packages, simply removing packages from the rthdynten field in the atmosphere core's Registry.xml file should have little impact on the memory usage of the dynamics, only increasing memory use outside of the dynamics in cases where neither the Grell-Freitas nor the Tiedtke/nTiedtke schemes are used. However, this potentially larger memory use comes with the benefit of significantly simplified logic in the dynamics. --- src/core_atmosphere/Registry.xml | 7 ++--- .../dynamics/mpas_atm_time_integration.F | 28 ++----------------- 2 files changed, 6 insertions(+), 29 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 173b483774..8157bf7e9d 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -2962,6 +2962,9 @@ description="Total dry air tendency from physics" persistence="scratch" /> + + - - diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index 7b80f1e087..63b2ddd25e 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -1142,9 +1142,7 @@ subroutine atm_srk3(domain, dt, itimestep, exchange_halo_group) call mpas_pool_get_array(tend_physics, 'rqvdynten', rqvdynten) call mpas_pool_get_array(state, 'theta_m', theta_m, 2) - if (associated(tend_physics)) then - call mpas_pool_get_array(tend_physics, 'rthdynten', rthdynten) - end if + call mpas_pool_get_array(tend_physics, 'rthdynten', rthdynten) !NOTE: The calculation of the tendency due to horizontal and vertical advection for the water vapor mixing ratio !requires that the subroutine atm_advance_scalars_mono was called on the third Runge Kutta step, so that a halo @@ -3888,7 +3886,7 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, real (kind=RKIND), dimension(:,:), pointer :: rr_save - real (kind=RKIND), dimension(:,:), pointer :: rthdynten => null() + real (kind=RKIND), dimension(:,:), pointer :: rthdynten real (kind=RKIND), dimension(:,:,:), pointer :: scalars @@ -3931,8 +3929,6 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, real (kind=RKIND), pointer :: config_rayleigh_damp_u_timescale_days integer, pointer :: config_number_rayleigh_damp_u_levels, config_number_cam_damping_levels - logical :: inactive_rthdynten - call mpas_pool_get_config(mesh, 'sphere_radius', r_earth) call mpas_pool_get_config(configs, 'config_coef_3rd_order', coef_3rd_order) @@ -3982,9 +3978,7 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, call mpas_pool_get_array(diag, 'h_divergence', h_divergence) call mpas_pool_get_array(diag, 'exner', exner) - if (associated(tend_physics)) then - call mpas_pool_get_array(tend_physics, 'rthdynten', rthdynten) - end if + call mpas_pool_get_array(tend_physics, 'rthdynten', rthdynten) call mpas_pool_get_array(mesh, 'weightsOnEdge', weightsOnEdge) call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge) @@ -4061,18 +4055,6 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, call mpas_pool_get_array(mesh, 'cf2', cf2) call mpas_pool_get_array(mesh, 'cf3', cf3) - ! - ! rthdynten is currently associated with packages, and if those packages - ! are not active at run-time, we need to produce an rthdynten array for - ! use in the atm_compute_dyn_tend_work routine - ! - inactive_rthdynten = .false. - if (.not. associated(rthdynten)) then - allocate(rthdynten(nVertLevels,nCells+1)) - rthdynten(:,nCells+1) = 0.0_RKIND - inactive_rthdynten = .true. - end if - call atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels, & nCellsSolve, nEdgesSolve, vertexDegree, maxEdges, maxEdges2, num_scalars, moist_start, moist_end, & fEdge, dvEdge, dcEdge, invDcEdge, invDvEdge, invAreaCell, invAreaTriangle, meshScalingDel2, meshScalingDel4, & @@ -4094,10 +4076,6 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) - if (inactive_rthdynten) then - deallocate(rthdynten) - end if - end subroutine atm_compute_dyn_tend