From 7f0e71c766c760022133b6252caab4c0ce6ef894 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 16 Jan 2017 14:27:12 +0100 Subject: [PATCH 01/24] Bug fix to propagate GEN_F90 to src/external/esmf_time_f90 for compilation on IBM Bluegene --- src/external/Makefile | 2 +- src/external/esmf_time_f90/Makefile | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/external/Makefile b/src/external/Makefile index a64539e270..4409d9c704 100644 --- a/src/external/Makefile +++ b/src/external/Makefile @@ -3,7 +3,7 @@ all: esmf_time ezxml-lib esmf_time: - ( cd esmf_time_f90; $(MAKE) FC="$(FC) $(FFLAGS)" CPP="$(CPP)" CPPFLAGS="$(CPPFLAGS) -DHIDE_MPI" ) + ( cd esmf_time_f90; $(MAKE) FC="$(FC) $(FFLAGS)" CPP="$(CPP)" CPPFLAGS="$(CPPFLAGS) -DHIDE_MPI" GEN_F90=$(GEN_F90) ) ezxml-lib: ( cd ezxml; $(MAKE) ) diff --git a/src/external/esmf_time_f90/Makefile b/src/external/esmf_time_f90/Makefile index 192a52c2b6..0b520c0853 100644 --- a/src/external/esmf_time_f90/Makefile +++ b/src/external/esmf_time_f90/Makefile @@ -57,4 +57,10 @@ clean: .F90.o: $(RM) $@ $*.mod +ifeq "$(GEN_F90)" "true" + $(CPP) $(CPPFLAGS) $< > $*.f90 $(CPPINCLUDES) $(FCINCLUDES) -I. + $(FC) $(FFLAGS) -c $*.f90 $(CPPINCLUDES) $(FCINCLUDES) -I. +else $(FC) $(CPPFLAGS) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I. +endif +# $(FC) $(CPPFLAGS) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I. From f3956e27e14fa8a7ef9f324f87a9f846e24d6a5c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 16 Jan 2017 18:39:19 +0100 Subject: [PATCH 02/24] Fix to compile Thompson MP scheme on IBM Bluegene; compiler cannot handle conversions of data types as arguments of max() --- src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F b/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F index f027f56d14..20190e3fd3 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F +++ b/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F @@ -3019,7 +3019,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ! -tpc_wev(idx_d, idx_c, idx_n)*orho*odt) prw_vcd(k) = MAX(DBLE(-rc(k)*0.99*orho*odt), prw_vcd(k)) pnc_wcd(k) = MAX(DBLE(-nc(k)*0.99*orho*odt), & - -tnc_wev(idx_d, idx_c, idx_n)*orho*odt) + DBLE(-tnc_wev(idx_d, idx_c, idx_n)*orho*odt)) endif else From 2423776a9e067751cc2171ccbe98aa1eac697c3f Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 17 Jan 2017 07:24:26 +0100 Subject: [PATCH 03/24] Bug fixes for landice core for compilation on IBM Bluegene --- src/core_landice/mpas_li_diagnostic_vars.F | 4 ++-- src/core_landice/mpas_li_setup.F | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/core_landice/mpas_li_diagnostic_vars.F b/src/core_landice/mpas_li_diagnostic_vars.F index 8186e846b1..2402c97e5d 100644 --- a/src/core_landice/mpas_li_diagnostic_vars.F +++ b/src/core_landice/mpas_li_diagnostic_vars.F @@ -625,7 +625,7 @@ subroutine vertical_remap_cism_loops(layerThickness, thickness, tracers, meshPoo do iCell = 1, nCells zhi = min (layerInterfaceSigma_Input(k1+1,iCell), layerInterfaceSigma(k2+1)) zlo = max (layerInterfaceSigma_Input(k1,iCell), layerInterfaceSigma(k2)) - hOverlap = max (zhi-zlo, 0.d0) * thickness(iCell) + hOverlap = max (DBLE(zhi-zlo), 0.d0) * thickness(iCell) hTsum(iCell,nt,k2) = htsum(iCell,nt,k2) & + tracers(nt,k1,iCell) * hOverlap enddo ! iCell @@ -781,7 +781,7 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers do nt = 1, nTracers zhi = min (layerInterfaceSigma_Input(k1+1), layerInterfaceSigma(k2+1)) zlo = max (layerInterfaceSigma_Input(k1), layerInterfaceSigma(k2)) - hOverlap = max (zhi-zlo, 0.d0) * thisThk + hOverlap = max (DBLE(zhi-zlo), 0.d0) * thisThk hTsum(nt,k2) = htsum(nt,k2) + tracers(nt,k1,iCell) * hOverlap enddo ! nt enddo ! k1 diff --git a/src/core_landice/mpas_li_setup.F b/src/core_landice/mpas_li_setup.F index 6e06c0ebdf..507b417458 100644 --- a/src/core_landice/mpas_li_setup.F +++ b/src/core_landice/mpas_li_setup.F @@ -178,7 +178,7 @@ subroutine li_setup_vertical_grid(meshPool, err) write(0,*) 'Error: The sum of layerThicknessFractions is different from 1.0 by more than 0.001.' err = 1 end if - write (6,*), 'Adjusting upper layerThicknessFrac by small amount because sum of layerThicknessFractions is slightly different from 1.0.' + write (6,*) 'Adjusting upper layerThicknessFrac by small amount because sum of layerThicknessFractions is slightly different from 1.0.' ! TODO - distribute the residual amongst all layers (and then put the residual of that in a single layer layerThicknessFractions(1) = layerThicknessFractions(1) - (fractionTotal - 1.0_RKIND) endif From 8ac643105ea43ebaa8cf30990fbd6602feedc670 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 18 Jan 2017 17:12:43 -0700 Subject: [PATCH 04/24] Update version number to v5.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_sw/Registry.xml | 2 +- src/core_test/Registry.xml | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 6484abe1e6..02f7bb0966 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -MPAS-v5.0 +MPAS-v5.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 253970d1b0..e91767a5d7 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 6a6cccb384..b96a3a3e15 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 b623ac65a7..47daaff5c3 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 23d39b7844..ae880c0004 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 9e969559ec..6ab6c0780b 100644 --- a/src/core_test/Registry.xml +++ b/src/core_test/Registry.xml @@ -1,5 +1,5 @@ - + From 27b3101a6c297b0c50c649fa33082db5451c7770 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 27 Jan 2017 15:46:16 -0700 Subject: [PATCH 05/24] Avoid setting invalid pointers to scalar constituents in mpas_atmphys_interface module Several scalar constituents, e.g., 'ni' and 'nr', are only allocated by specific microphysics schemes through the use of packages. In the routines microphysics_to_MPAS and microphysics_from_MPAS, we should only set pointers to sub-arrays of the scalars var_array for constituents that are guaranteed to be allocated based on the choice of microphysics scheme. Thanks to Wei Huang for finding this problem in MPAS v5.0. --- .../physics/mpas_atmphys_interface.F | 24 +++++++++++-------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_interface.F b/src/core_atmosphere/physics/mpas_atmphys_interface.F index 06b0e972ef..fe4e2e1377 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_interface.F +++ b/src/core_atmosphere/physics/mpas_atmphys_interface.F @@ -508,11 +508,6 @@ subroutine microphysics_from_MPAS(mesh,state,time_lev,diag,diag_physics,its,ite) qv => scalars(index_qv,:,:) qc => scalars(index_qc,:,:) qr => scalars(index_qr,:,:) - qi => scalars(index_qi,:,:) - qs => scalars(index_qs,:,:) - qg => scalars(index_qg,:,:) - ni => scalars(index_ni,:,:) - nr => scalars(index_nr,:,:) !initialize variables needed in the cloud microphysics schemes: do j = jts, jte @@ -539,6 +534,10 @@ subroutine microphysics_from_MPAS(mesh,state,time_lev,diag,diag_physics,its,ite) microp_select_init: select case(microp_scheme) case ("mp_thompson","mp_wsm6") + qi => scalars(index_qi,:,:) + qs => scalars(index_qs,:,:) + qg => scalars(index_qg,:,:) + do j = jts, jte do k = kts, kte do i = its, ite @@ -552,6 +551,9 @@ subroutine microphysics_from_MPAS(mesh,state,time_lev,diag,diag_physics,its,ite) microp2_select: select case(microp_scheme) case("mp_thompson") + ni => scalars(index_ni,:,:) + nr => scalars(index_nr,:,:) + do j = jts,jte do i = its,ite muc_p(i,j) = mu_c(i) @@ -654,11 +656,6 @@ subroutine microphysics_to_MPAS(mesh,state,time_lev,diag,diag_physics,tend,itime qv => scalars(index_qv,:,:) qc => scalars(index_qc,:,:) qr => scalars(index_qr,:,:) - qi => scalars(index_qi,:,:) - qs => scalars(index_qs,:,:) - qg => scalars(index_qg,:,:) - ni => scalars(index_ni,:,:) - nr => scalars(index_nr,:,:) call mpas_pool_get_array(tend,'rt_diabatic_tend',rt_diabatic_tend) @@ -711,6 +708,10 @@ subroutine microphysics_to_MPAS(mesh,state,time_lev,diag,diag_physics,tend,itime microp_select_init: select case(microp_scheme) case ("mp_thompson","mp_wsm6") + qi => scalars(index_qi,:,:) + qs => scalars(index_qs,:,:) + qg => scalars(index_qg,:,:) + do j = jts, jte do k = kts, kte do i = its, ite @@ -724,6 +725,9 @@ subroutine microphysics_to_MPAS(mesh,state,time_lev,diag,diag_physics,tend,itime microp2_select: select case(microp_scheme) case("mp_thompson") + ni => scalars(index_ni,:,:) + nr => scalars(index_nr,:,:) + do j = jts, jte do k = kts, kte do i = its, ite From dc3c19a7f3d17b84d96c6496a8c7b8d6587fdfc8 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 24 Jan 2017 18:51:41 +0100 Subject: [PATCH 06/24] Removed CPP directives -DSINGLE_PRECISION from Fortran compiler flags as this breaks builds on systems where GEN_F90 is required (e.g. Bluegene) --- Makefile | 1 - 1 file changed, 1 deletion(-) diff --git a/Makefile b/Makefile index 97ae93f709..bf80bc87b1 100644 --- a/Makefile +++ b/Makefile @@ -509,7 +509,6 @@ ifeq "$(OPENMP)" "true" endif #OPENMP IF ifeq "$(PRECISION)" "single" - FFLAGS += "-DSINGLE_PRECISION" CFLAGS += "-DSINGLE_PRECISION" CXXFLAGS += "-DSINGLE_PRECISION" override CPPFLAGS += "-DSINGLE_PRECISION" From 7f82b610771f7ca7244ee287db54bdc9aa86dbe8 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sat, 28 Jan 2017 18:46:48 +0100 Subject: [PATCH 07/24] Changes to real/double conversion as requested in amendment of 2423776a9e067751cc2171ccbe98aa1eac697c3f, however 0.0d0_RKIND does not work, need to use 0.0_RKIND --- src/core_landice/mpas_li_diagnostic_vars.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_landice/mpas_li_diagnostic_vars.F b/src/core_landice/mpas_li_diagnostic_vars.F index 2402c97e5d..53fd7f24c2 100644 --- a/src/core_landice/mpas_li_diagnostic_vars.F +++ b/src/core_landice/mpas_li_diagnostic_vars.F @@ -625,7 +625,7 @@ subroutine vertical_remap_cism_loops(layerThickness, thickness, tracers, meshPoo do iCell = 1, nCells zhi = min (layerInterfaceSigma_Input(k1+1,iCell), layerInterfaceSigma(k2+1)) zlo = max (layerInterfaceSigma_Input(k1,iCell), layerInterfaceSigma(k2)) - hOverlap = max (DBLE(zhi-zlo), 0.d0) * thickness(iCell) + hOverlap = max (zhi-zlo, 0.0_RKIND) * thickness(iCell) hTsum(iCell,nt,k2) = htsum(iCell,nt,k2) & + tracers(nt,k1,iCell) * hOverlap enddo ! iCell @@ -781,7 +781,7 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers do nt = 1, nTracers zhi = min (layerInterfaceSigma_Input(k1+1), layerInterfaceSigma(k2+1)) zlo = max (layerInterfaceSigma_Input(k1), layerInterfaceSigma(k2)) - hOverlap = max (DBLE(zhi-zlo), 0.d0) * thisThk + hOverlap = max (zhi-zlo, 0.0_RKIND) * thisThk hTsum(nt,k2) = htsum(nt,k2) + tracers(nt,k1,iCell) * hOverlap enddo ! nt enddo ! k1 From 48a58fdab2c367a14c5dbe37b05c404b6a76e9c4 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 1 Feb 2017 09:00:17 +0100 Subject: [PATCH 08/24] Updated INSTALL documentation with known problems on IBM Bluegene and Mac OSX --- INSTALL | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/INSTALL b/INSTALL index 1bdd0d89f1..f16c6d8a14 100644 --- a/INSTALL +++ b/INSTALL @@ -12,7 +12,11 @@ and the gfortran compiler. The gfortran compiler can be installed using homebrew (http://brew.sh), or using pre-compiled binaries from the MacOSX HPC website (http://hpc.soureforge.net), or it can be compiled by the user from the GNU sources. -The standard clang compiler does not support OpenMP. Users wanting to compile MPAS +Caveats: + +(1) MPAS cannot be compiled with gfortran-clang if GEN_F90=true. + +(2) The standard clang compiler does not support OpenMP. Users wanting to compile MPAS with OpenMP support on MacOSX will have to install the LLVM clang compiler, which is accomplished easiest with homebrew. Since this alternative clang compiler is not in the standard search/library path, the user will have to modify the call to the clang @@ -57,3 +61,10 @@ would become ... assuming that the LLVM clang compiler is installed in /usr/local/opt/llvm. + + +bluegene: Compiling MPAS on IBM Bluegene using the xl compilers +---------- +All MPAS cores except the ocean compile on IBM Bluegene using the xl compilers. The ocean +core currently does not work on IBM Bluegene. Known limitations: OPENMP must be disabled +(OPENMP=false) for compiling, since the xl compilers do not support nested OpenMP directives. From aca6d9438314f49c7671dbcc0d83432ca616608b Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 7 Feb 2017 08:14:54 +0100 Subject: [PATCH 09/24] Bug fix to re-enable the read and interpolation of pressure levels as output by calc_ecmpwf_p.exe from WPS --- src/core_init_atmosphere/mpas_init_atm_cases.F | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index a92b1a9cd9..247b07cea2 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -3346,6 +3346,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state trim(field % field) == 'ST028100' .or. & trim(field % field) == 'ST100255' .or. & trim(field % field) == 'PRES' .or. & + trim(field % field) == 'PRESSURE' .or. & trim(field % field) == 'SNOW' .or. & trim(field % field) == 'SEAICE' .or. & trim(field % field) == 'SKINTEMP') then @@ -3510,6 +3511,13 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state lonPoints => lonCell call mpas_pool_get_array(fg, 'p', destField2d) ndims = 2 + else if (trim(field % field) == 'PRESSURE') then +write(0,*) 'Interpolating PRESSURE at ', k, vert_level(k) + nInterpPoints = nCells + latPoints => latCell + lonPoints => lonCell + call mpas_pool_get_array(fg, 'p', destField2d) + ndims = 2 else if (trim(field % field) == 'PMSL') then write(0,*) 'Interpolating PMSL' nInterpPoints = nCells From 57f7cb7d8207b8bdd751ab755ea113a8825a4897 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 8 Feb 2017 06:36:31 +0100 Subject: [PATCH 10/24] Tiny formatting corrections in mpas_init_atm_cases.F --- src/core_init_atmosphere/mpas_init_atm_cases.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index 247b07cea2..361e0d4974 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -3374,8 +3374,8 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state trim(field % field) == 'SKINTEMP') then k = 1 else if (trim(field % field) /= 'PMSL' .and. & - trim(field % field) /='PSFC' .and. & - trim(field % field) /= 'SOILHGT') then + trim(field % field) /= 'PSFC' .and. & + trim(field % field) /= 'SOILHGT') then ! Since the hash table can only store integers, transfer the bit pattern from ! the real-valued xlvl into an integer; that the result is not an integer version From 9e0d753eb8ff01b18514c2446194214bb2ded3c7 Mon Sep 17 00:00:00 2001 From: Laura Fowler Date: Fri, 10 Feb 2017 17:31:06 -0700 Subject: [PATCH 11/24] * In ./src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F and mpas_atmphys_driver_radiation_sw.F, revised the initialization of the variables re_cloud, re_ice, re_snow, and rre_cloud, rre_ice, and rre_snow, when config_microp_scheme is set to off. --- .../mpas_atmphys_driver_radiation_lw.F | 117 +++++++++++++----- .../mpas_atmphys_driver_radiation_sw.F | 66 ++++++---- 2 files changed, 125 insertions(+), 58 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F index 52e841d4c4..0db3900c2d 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F @@ -82,6 +82,10 @@ module mpas_atmphys_driver_radiation_lw ! * in the call to rrtmg_lwrad, substituted the variables qv_p, qc_p, qi_p, and qs_p with qvrad_p, qcrad_p, ! qirad_p, and qsrad_p initialized in subroutine cloudiness_from_MPAS. ! Laura D. Fowler (laura@ucar.edu) / 2016-07-09. +! * in subroutines radiation_lw_from_MPAS and radiation_lw_to_MPAS, revised the initialization of re_cloud, +! re_ice, re_snow, and rre_cloud, rre_ice, and rre_snow to handle the case when the cloud microphysics +! parameterization is turned off, i.e. config_microp_scheme='off'. +! Laura D. Fowler (laura@ucar.edu) / 2017-02-10. contains @@ -292,6 +296,7 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi !local pointers: logical,pointer:: config_o3climatology + logical,pointer:: config_microp_re real(kind=RKIND),dimension(:),pointer :: latCell,lonCell real(kind=RKIND),dimension(:),pointer :: skintemp,snow,xice,xland @@ -309,6 +314,7 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi !----------------------------------------------------------------------------------------------------------------- call mpas_pool_get_config(configs,'config_o3climatology',config_o3climatology) + call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re ) call mpas_pool_get_array(mesh,'latCell',latCell) call mpas_pool_get_array(mesh,'lonCell',lonCell) @@ -383,22 +389,48 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi radiation_lw_select: select case (trim(radt_lw_scheme)) case("rrtmg_lw") - call mpas_pool_get_array(diag_physics,'re_cloud',re_cloud) - call mpas_pool_get_array(diag_physics,'re_ice' ,re_ice ) - call mpas_pool_get_array(diag_physics,'re_snow' ,re_snow ) + microp_select: select case(microp_scheme) + case("mp_thompson","mp_wsm6") + if(config_microp_re) then + call mpas_pool_get_array(diag_physics,'re_cloud',re_cloud) + call mpas_pool_get_array(diag_physics,'re_ice' ,re_ice ) + call mpas_pool_get_array(diag_physics,'re_snow' ,re_snow ) + + do j = jts,jte + do k = kts,kte + do i = its,ite + recloud_p(i,k,j) = re_cloud(k,i) + reice_p(i,k,j) = re_ice(k,i) + resnow_p(i,k,j) = re_snow(k,i) + enddo + enddo + enddo + else + has_reqc = 0 + has_reqi = 0 + has_reqs = 0 + do j = jts,jte + do k = kts,kte + do i = its,ite + recloud_p(i,k,j) = 0._RKIND + reice_p(i,k,j) = 0._RKIND + resnow_p(i,k,j) = 0._RKIND + enddo + enddo + enddo + endif + do j = jts,jte + do k = kts,kte + do i = its,ite + rrecloud_p(i,k,j) = 0._RKIND + rreice_p(i,k,j) = 0._RKIND + rresnow_p(i,k,j) = 0._RKIND + enddo + enddo + enddo - do j = jts,jte - do k = kts,kte - do i = its,ite - recloud_p(i,k,j) = re_cloud(k,i) - reice_p(i,k,j) = re_ice(k,i) - resnow_p(i,k,j) = re_snow(k,i) - rrecloud_p(i,k,j) = 0._RKIND - rreice_p(i,k,j) = 0._RKIND - rresnow_p(i,k,j) = 0._RKIND - enddo - enddo - enddo + case default + end select microp_select do j = jts,jte do k = kts,kte+2 @@ -559,16 +591,21 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi end subroutine radiation_lw_from_MPAS !================================================================================================================= - subroutine radiation_lw_to_MPAS(diag_physics,tend_physics,its,ite) + subroutine radiation_lw_to_MPAS(configs,diag_physics,tend_physics,its,ite) !================================================================================================================= !input arguments: + type(mpas_pool_type),intent(in):: configs + +!inout arguments: type(mpas_pool_type),intent(inout):: diag_physics type(mpas_pool_type),intent(inout):: tend_physics integer,intent(in):: its,ite !local pointers: + logical,pointer:: config_microp_re + real(kind=RKIND),dimension(:),pointer :: glw,lwcf,lwdnb,lwdnbc,lwdnt,lwdntc,lwupb,lwupbc, & lwupt,lwuptc,olrtoa real(kind=RKIND),dimension(:,:),pointer:: rthratenlw @@ -581,6 +618,8 @@ subroutine radiation_lw_to_MPAS(diag_physics,tend_physics,its,ite) !----------------------------------------------------------------------------------------------------------------- + call mpas_pool_get_config(configs,'config_microp_re',config_microp_re) + call mpas_pool_get_array(diag_physics,'glw' ,glw ) call mpas_pool_get_array(diag_physics,'lwcf' ,lwcf ) call mpas_pool_get_array(diag_physics,'lwdnb' ,lwdnb ) @@ -620,21 +659,36 @@ subroutine radiation_lw_to_MPAS(diag_physics,tend_physics,its,ite) radiation_lw_select: select case (trim(radt_lw_scheme)) case("rrtmg_lw") - call mpas_pool_get_array(diag_physics,'rre_cloud',rre_cloud) - call mpas_pool_get_array(diag_physics,'rre_ice' ,rre_ice ) - call mpas_pool_get_array(diag_physics,'rre_snow' ,rre_snow ) - - do j = jts,jte - do k = kts,kte - do i = its,ite - rre_cloud(k,i) = rrecloud_p(i,k,j) - rre_ice(k,i) = rreice_p(i,k,j) - rre_snow(k,i) = rresnow_p(i,k,j) - enddo - enddo - enddo - case default + microp_select: select case(microp_scheme) + case("mp_thompson","mp_wsm6") + call mpas_pool_get_array(diag_physics,'rre_cloud',rre_cloud) + call mpas_pool_get_array(diag_physics,'rre_ice' ,rre_ice ) + call mpas_pool_get_array(diag_physics,'rre_snow' ,rre_snow ) + if(config_microp_re) then + do j = jts,jte + do k = kts,kte + do i = its,ite + rre_cloud(k,i) = rrecloud_p(i,k,j) + rre_ice(k,i) = rreice_p(i,k,j) + rre_snow(k,i) = rresnow_p(i,k,j) + enddo + enddo + enddo + else + do j = jts,jte + do k = kts,kte + do i = its,ite + rre_cloud(k,i) = 0._RKIND + rre_ice(k,i) = 0._RKIND + rre_snow(k,i) = 0._RKIND + enddo + enddo + enddo + endif + + case default + end select microp_select end select radiation_lw_select @@ -785,6 +839,7 @@ subroutine driver_radiation_lw(xtime_s,configs,mesh,state,time_lev,diag_physics, ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , & its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte & ) + case ("cam_lw") xtime_m = xtime_s/60. @@ -849,7 +904,7 @@ subroutine driver_radiation_lw(xtime_s,configs,mesh,state,time_lev,diag_physics, end select radiation_lw_select !copy local arrays to MPAS grid: - call radiation_lw_to_MPAS(diag_physics,tend_physics,its,ite) + call radiation_lw_to_MPAS(configs,diag_physics,tend_physics,its,ite) !write(0,*) '--- end subroutine driver_radiation_lw.' diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F index d04c90f701..ba99644efc 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F @@ -77,6 +77,9 @@ module mpas_atmphys_driver_radiation_sw ! * in the call to rrtmg_swrad, substituted the variables qv_p, qc_p, qi_p, and qs_p with qvrad_p, qcrad_p, ! qirad_p, and qsrad_p initialized in subroutine cloudiness_from_MPAS. ! Laura D. Fowler (laura@ucar.edu) / 2016-07-09. +! * in subroutines radiation_sw_from_MPAS, revised the initialization of re_cloud, re_ice, re_snow, to +! handle the case when the cloud microphysics parameterization is turned off, i.e. config_microp_scheme='off'. +! Laura D. Fowler (laura@ucar.edu) / 2017-02-10. contains @@ -296,6 +299,7 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i !local pointers: logical,pointer:: config_o3climatology + logical,pointer:: config_microp_re real(kind=RKIND),dimension(:),pointer :: latCell,lonCell real(kind=RKIND),dimension(:),pointer :: skintemp,snow,xice,xland @@ -308,6 +312,7 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i !----------------------------------------------------------------------------------------------------------------- call mpas_pool_get_config(configs,'config_o3climatology',config_o3climatology) + call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re ) call mpas_pool_get_array(mesh,'latCell',latCell) call mpas_pool_get_array(mesh,'lonCell',lonCell) @@ -384,19 +389,40 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i radiation_sw_select: select case (trim(radt_sw_scheme)) case("rrtmg_sw") - call mpas_pool_get_array(diag_physics,'re_cloud',re_cloud) - call mpas_pool_get_array(diag_physics,'re_ice' ,re_ice ) - call mpas_pool_get_array(diag_physics,'re_snow' ,re_snow ) - do j = jts,jte - do k = kts,kte - do i = its,ite - recloud_p(i,k,j) = re_cloud(k,i) - reice_p(i,k,j) = re_ice(k,i) - resnow_p(i,k,j) = re_snow(k,i) - enddo - enddo - enddo + microp_select: select case(microp_scheme) + case("mp_thompson","mp_wsm6") + if(config_microp_re) then + call mpas_pool_get_array(diag_physics,'re_cloud',re_cloud) + call mpas_pool_get_array(diag_physics,'re_ice' ,re_ice ) + call mpas_pool_get_array(diag_physics,'re_snow' ,re_snow ) + + do j = jts,jte + do k = kts,kte + do i = its,ite + recloud_p(i,k,j) = re_cloud(k,i) + reice_p(i,k,j) = re_ice(k,i) + resnow_p(i,k,j) = re_snow(k,i) + enddo + enddo + enddo + else + has_reqc = 0 + has_reqi = 0 + has_reqs = 0 + do j = jts,jte + do k = kts,kte + do i = its,ite + recloud_p(i,k,j) = 0._RKIND + reice_p(i,k,j) = 0._RKIND + resnow_p(i,k,j) = 0._RKIND + enddo + enddo + enddo + endif + + case default + end select microp_select do j = jts,jte do k = kts,kte+2 @@ -409,7 +435,7 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i enddo enddo - !ozone volum mixing ratio: + !ozone volume mixing ratio: if(config_o3climatology) then do k = 1, num_oznLevels pin_p(k) = pin(k) @@ -544,7 +570,6 @@ subroutine radiation_sw_to_MPAS(diag_physics,tend_physics,its,ite) !local pointers: real(kind=RKIND),dimension(:),pointer :: coszr,gsw,swcf,swdnb,swdnbc,swdnt,swdntc, & swupb,swupbc,swupt,swuptc -!real(kind=RKIND),dimension(:,:),pointer:: swdnflx,swdnflxc,swupflx,swupflxc real(kind=RKIND),dimension(:,:),pointer:: rthratensw !----------------------------------------------------------------------------------------------------------------- @@ -560,10 +585,6 @@ subroutine radiation_sw_to_MPAS(diag_physics,tend_physics,its,ite) call mpas_pool_get_array(diag_physics,'swupbc' , swupbc ) call mpas_pool_get_array(diag_physics,'swupt' ,swupt ) call mpas_pool_get_array(diag_physics,'swuptc' ,swuptc ) -!call mpas_pool_get_array(diag_physics,'swdnflx' ,swdnflx ) -!call mpas_pool_get_array(diag_physics,'swdnflxc' ,swdnflxc ) -!call mpas_pool_get_array(diag_physics,'swupflx' ,swupflx ) -!call mpas_pool_get_array(diag_physics,'swupflxc' ,swupflxc ) call mpas_pool_get_array(tend_physics,'rthratensw',rthratensw) do j = jts,jte @@ -581,15 +602,6 @@ subroutine radiation_sw_to_MPAS(diag_physics,tend_physics,its,ite) swupt(i) = swupt_p(i,j) swuptc(i) = swuptc_p(i,j) enddo -!not needed: -!do k = kts,kte+2 -!do i = its,ite -! swdnflx(k,i) = swdnflx_p(i,k,j) -! swdnflxc(k,i) = swdnflxc_p(i,k,j) -! swupflx(k,i) = swupflx_p(i,k,j) -! swupflxc(k,i) = swupflxc_p(i,k,j) -!enddo -!enddo do k = kts,kte do i = its,ite From 570a032cdab38810c84c036c70543d74081e102a Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 8 Feb 2017 14:21:18 +0100 Subject: [PATCH 12/24] Provide compatibility with ECMWF 100-289 cm soil layer This commit introduces changes to allow the init_atmosphere core to properly handle the ECMWF 100-289 cm soil layer (for both soil temperature and soil moisture). In older versions of the WPS, the 100-289 cm layer was decoded as spanning the depths 100-255 cm; however, the actual bottom depth of this layer in ECMWF is 289 cm, and this was later corrected in the WPS. In order to be compatible with this new, correct layer, additional block of code are added in the init_atm_case_gfs( ) routine. At present, compatibility with the 100-255 cm layer is maintained, and the 100-289 cm layer is treated in a completely independent way. --- .../mpas_init_atm_cases.F | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index a92b1a9cd9..a08ae6fceb 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -3336,6 +3336,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state trim(field % field) == 'SM007028' .or. & trim(field % field) == 'SM028100' .or. & trim(field % field) == 'SM100255' .or. & + trim(field % field) == 'SM100289' .or. & trim(field % field) == 'ST000010' .or. & trim(field % field) == 'ST010040' .or. & trim(field % field) == 'ST040100' .or. & @@ -3345,6 +3346,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state trim(field % field) == 'ST007028' .or. & trim(field % field) == 'ST028100' .or. & trim(field % field) == 'ST100255' .or. & + trim(field % field) == 'ST100289' .or. & trim(field % field) == 'PRES' .or. & trim(field % field) == 'SNOW' .or. & trim(field % field) == 'SEAICE' .or. & @@ -3359,6 +3361,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state trim(field % field) == 'SM007028' .or. & trim(field % field) == 'SM028100' .or. & trim(field % field) == 'SM100255' .or. & + trim(field % field) == 'SM100289' .or. & trim(field % field) == 'ST000010' .or. & trim(field % field) == 'ST010040' .or. & trim(field % field) == 'ST040100' .or. & @@ -3368,6 +3371,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state trim(field % field) == 'ST007028' .or. & trim(field % field) == 'ST028100' .or. & trim(field % field) == 'ST100255' .or. & + trim(field % field) == 'ST100289' .or. & trim(field % field) == 'SNOW' .or. & trim(field % field) == 'SEAICE' .or. & trim(field % field) == 'SKINTEMP') then @@ -3711,6 +3715,26 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state ndims = 2 dzs_fg(k,:) = 255.-100. zs_fg(k,:) = 255. + else if (trim(field % field) == 'SM100289') then +write(0,*) 'Interpolating SM100289' + + interp_list(1) = FOUR_POINT + interp_list(2) = W_AVERAGE4 + interp_list(3) = SEARCH + interp_list(4) = 0 + + maskval = 0.0 + masked = 0 + fillval = 1.0 + + nInterpPoints = nCells + latPoints => latCell + lonPoints => lonCell + call mpas_pool_get_array(fg, 'sm_fg', destField2d) + k = 4 + ndims = 2 + dzs_fg(k,:) = 289.-100. + zs_fg(k,:) = 289. else if (trim(field % field) == 'ST000010') then write(0,*) 'Interpolating ST000010' @@ -3900,6 +3924,27 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state ndims = 2 dzs_fg(k,:) = 255.-100. zs_fg(k,:) = 255. + else if (trim(field % field) == 'ST100289') then +write(0,*) 'Interpolating ST100289' + + interp_list(1) = SIXTEEN_POINT + interp_list(2) = FOUR_POINT + interp_list(3) = W_AVERAGE4 + interp_list(4) = SEARCH + interp_list(5) = 0 + + maskval = 0.0 + masked = 0 + fillval = 285.0 + + nInterpPoints = nCells + latPoints => latCell + lonPoints => lonCell + call mpas_pool_get_array(fg, 'st_fg', destField2d) + k = 4 + ndims = 2 + dzs_fg(k,:) = 289.-100. + zs_fg(k,:) = 289. else if (trim(field % field) == 'SNOW') then write(0,*) 'Interpolating SNOW' From 3e93d8e7ad0c4c921690b1240eb7d6a930b6cc90 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 8 Feb 2017 10:17:38 +0100 Subject: [PATCH 13/24] Change name of preprocessed Fortran files from .f90 to .tmp.f90 in esmf_time_f90 The default macOS filesystem is case-insensitve and therefore does not distinguish between .F90 and .f90 files, which creates a problem when compiling the ESMF time manager code with GEN_F90=true, since the .F90 files will be overwritten with the preprocessed .f90 files without the .tmp infix --- src/external/esmf_time_f90/Makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/external/esmf_time_f90/Makefile b/src/external/esmf_time_f90/Makefile index 0b520c0853..79da5c5c04 100644 --- a/src/external/esmf_time_f90/Makefile +++ b/src/external/esmf_time_f90/Makefile @@ -58,8 +58,8 @@ clean: .F90.o: $(RM) $@ $*.mod ifeq "$(GEN_F90)" "true" - $(CPP) $(CPPFLAGS) $< > $*.f90 $(CPPINCLUDES) $(FCINCLUDES) -I. - $(FC) $(FFLAGS) -c $*.f90 $(CPPINCLUDES) $(FCINCLUDES) -I. + $(CPP) $(CPPFLAGS) $< > $*.tmp.f90 $(CPPINCLUDES) $(FCINCLUDES) -I. + $(FC) $(FFLAGS) -c $*.tmp.f90 $(CPPINCLUDES) $(FCINCLUDES) -I. else $(FC) $(CPPFLAGS) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I. endif From 1a2c0c25e36d513be30ad649db2c37e91dcd2068 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 8 Feb 2017 10:22:44 +0100 Subject: [PATCH 14/24] Replace angle brackets with quotes in include statements in esmf_time_f90 For non-standard include files in src/external/esmf_time_f90, #include <...> is replaced with #include "...". --- src/external/esmf_time_f90/ESMF.F90 | 2 +- src/external/esmf_time_f90/ESMF_AlarmClockMod.F90 | 2 +- src/external/esmf_time_f90/ESMF_AlarmMod.F90 | 2 +- src/external/esmf_time_f90/ESMF_BaseTimeMod.F90 | 2 +- src/external/esmf_time_f90/ESMF_CalendarMod.F90 | 2 +- src/external/esmf_time_f90/ESMF_ClockMod.F90 | 2 +- src/external/esmf_time_f90/ESMF_TimeIntervalMod.F90 | 2 +- src/external/esmf_time_f90/ESMF_TimeMgr.inc | 2 +- src/external/esmf_time_f90/ESMF_TimeMod.F90 | 2 +- src/external/esmf_time_f90/MeatMod.F90 | 2 +- 10 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/external/esmf_time_f90/ESMF.F90 b/src/external/esmf_time_f90/ESMF.F90 index 8eb5b7a181..f970f71ed5 100644 --- a/src/external/esmf_time_f90/ESMF.F90 +++ b/src/external/esmf_time_f90/ESMF.F90 @@ -13,7 +13,7 @@ MODULE ESMF USE ESMF_AlarmClockMod USE ESMF_Stubs ! add new dummy interfaces and typedefs here as needed USE MeatMod -#include +#include "ESMF_TimeMgr.inc" INTEGER, PARAMETER :: ESMF_MAX_ALARMS=MAX_ALARMS ! END MODULE ESMF diff --git a/src/external/esmf_time_f90/ESMF_AlarmClockMod.F90 b/src/external/esmf_time_f90/ESMF_AlarmClockMod.F90 index 63932f91f7..3e3e52ed27 100644 --- a/src/external/esmf_time_f90/ESMF_AlarmClockMod.F90 +++ b/src/external/esmf_time_f90/ESMF_AlarmClockMod.F90 @@ -19,7 +19,7 @@ module ESMF_AlarmClockMod ! !------------------------------------------------------------------------------ ! INCLUDES -#include +#include "ESMF_TimeMgr.inc" !=============================================================================== !BOPI diff --git a/src/external/esmf_time_f90/ESMF_AlarmMod.F90 b/src/external/esmf_time_f90/ESMF_AlarmMod.F90 index 67400ae7e6..fe0f0b405d 100644 --- a/src/external/esmf_time_f90/ESMF_AlarmMod.F90 +++ b/src/external/esmf_time_f90/ESMF_AlarmMod.F90 @@ -20,7 +20,7 @@ module ESMF_AlarmMod ! !------------------------------------------------------------------------------ ! INCLUDES -#include +#include "ESMF_TimeMgr.inc" !=============================================================================== !BOPI diff --git a/src/external/esmf_time_f90/ESMF_BaseTimeMod.F90 b/src/external/esmf_time_f90/ESMF_BaseTimeMod.F90 index 6eb4573afe..ad9b561de1 100644 --- a/src/external/esmf_time_f90/ESMF_BaseTimeMod.F90 +++ b/src/external/esmf_time_f90/ESMF_BaseTimeMod.F90 @@ -21,7 +21,7 @@ module ESMF_BaseTimeMod !------------------------------------------------------------------------------ ! INCLUDES -#include +#include "ESMF_TimeMgr.inc" ! !=============================================================================== !BOPI diff --git a/src/external/esmf_time_f90/ESMF_CalendarMod.F90 b/src/external/esmf_time_f90/ESMF_CalendarMod.F90 index dc874bdb9c..3a22c4e5f6 100644 --- a/src/external/esmf_time_f90/ESMF_CalendarMod.F90 +++ b/src/external/esmf_time_f90/ESMF_CalendarMod.F90 @@ -20,7 +20,7 @@ module ESMF_CalendarMod ! !------------------------------------------------------------------------------ ! INCLUDES -#include +#include "ESMF_TimeMgr.inc" !============================================================================== !BOPI diff --git a/src/external/esmf_time_f90/ESMF_ClockMod.F90 b/src/external/esmf_time_f90/ESMF_ClockMod.F90 index d634362e8d..627c4793e6 100644 --- a/src/external/esmf_time_f90/ESMF_ClockMod.F90 +++ b/src/external/esmf_time_f90/ESMF_ClockMod.F90 @@ -19,7 +19,7 @@ module ESMF_ClockMod ! !------------------------------------------------------------------------------ ! INCLUDES -#include +#include "ESMF_TimeMgr.inc" !============================================================================== !BOPI diff --git a/src/external/esmf_time_f90/ESMF_TimeIntervalMod.F90 b/src/external/esmf_time_f90/ESMF_TimeIntervalMod.F90 index 5d8be4e738..ccc8ccc53d 100644 --- a/src/external/esmf_time_f90/ESMF_TimeIntervalMod.F90 +++ b/src/external/esmf_time_f90/ESMF_TimeIntervalMod.F90 @@ -22,7 +22,7 @@ module ESMF_TimeIntervalMod ! !------------------------------------------------------------------------------ ! INCLUDES -#include +#include "ESMF_TimeMgr.inc" ! !=============================================================================== !BOPI diff --git a/src/external/esmf_time_f90/ESMF_TimeMgr.inc b/src/external/esmf_time_f90/ESMF_TimeMgr.inc index e41a1f8514..140b335b36 100644 --- a/src/external/esmf_time_f90/ESMF_TimeMgr.inc +++ b/src/external/esmf_time_f90/ESMF_TimeMgr.inc @@ -29,7 +29,7 @@ by both C++ and F90 compilers. !EOP #endif -#include +#include "ESMF_Macros.inc" #define SECONDS_PER_DAY 86400_ESMF_KIND_I8 #define SECONDS_PER_HOUR 3600_ESMF_KIND_I8 diff --git a/src/external/esmf_time_f90/ESMF_TimeMod.F90 b/src/external/esmf_time_f90/ESMF_TimeMod.F90 index 587b21a1d6..60b17fb01b 100644 --- a/src/external/esmf_time_f90/ESMF_TimeMod.F90 +++ b/src/external/esmf_time_f90/ESMF_TimeMod.F90 @@ -19,7 +19,7 @@ module ESMF_TimeMod ! !------------------------------------------------------------------------------ ! INCLUDES -#include +#include "ESMF_TimeMgr.inc" !============================================================================== !BOPI diff --git a/src/external/esmf_time_f90/MeatMod.F90 b/src/external/esmf_time_f90/MeatMod.F90 index c9ab2ec145..c6fedb2f59 100644 --- a/src/external/esmf_time_f90/MeatMod.F90 +++ b/src/external/esmf_time_f90/MeatMod.F90 @@ -1,7 +1,7 @@ module MeatMod -#include +#include "ESMF_TimeMgr.inc" use ESMF_BaseMod From b5d13e11c9dcfa583c75ab0cc1f74682d017bee2 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 8 Feb 2017 10:42:51 +0100 Subject: [PATCH 15/24] Remove .tmp.f90 files from esmf_time_f90 directory when 'make clean' is run --- src/external/esmf_time_f90/Makefile | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/external/esmf_time_f90/Makefile b/src/external/esmf_time_f90/Makefile index 79da5c5c04..790852e6d7 100644 --- a/src/external/esmf_time_f90/Makefile +++ b/src/external/esmf_time_f90/Makefile @@ -53,7 +53,7 @@ wrf_error_fatal.o: wrf_message.o: clean: - rm -rf *.o *.mod *.a + rm -rf *.tmp.f90 *.o *.mod *.a .F90.o: $(RM) $@ $*.mod @@ -63,4 +63,3 @@ ifeq "$(GEN_F90)" "true" else $(FC) $(CPPFLAGS) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I. endif -# $(FC) $(CPPFLAGS) $(FFLAGS) -c $*.F90 $(CPPINCLUDES) $(FCINCLUDES) -I. From f68596df2dc246ce5c4585ddb9bfd99d6707670c Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Tue, 7 Mar 2017 13:09:59 -0700 Subject: [PATCH 16/24] Ensure that temperature is computed if MSLP is needed in isobaric_diagnostics In the isobaric_diagnostics module, we have logic that attempts to compute fields only if they appear in an output stream, or if they are needed by another diagnostic that appears in an output stream. The temperature field was previously only computed if any of the isobaric temperature, dewpoint, or relative humidity fields appeared in an output stream. This caused problems if only MSLP was being written to a stream, since a valid temperature field was not available for the computation of MSLP. Now, the computation of temperature in the isobaric_diagnostics module also depends on whether MSLP appears in an output stream. --- src/core_atmosphere/diagnostics/isobaric_diagnostics.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_atmosphere/diagnostics/isobaric_diagnostics.F b/src/core_atmosphere/diagnostics/isobaric_diagnostics.F index 7ce0ba7183..4a9e4de650 100644 --- a/src/core_atmosphere/diagnostics/isobaric_diagnostics.F +++ b/src/core_atmosphere/diagnostics/isobaric_diagnostics.F @@ -553,7 +553,7 @@ subroutine interp_diagnostics(mesh, state, time_lev, diag) enddo enddo - if (NEED_TEMP .or. NEED_RELHUM .or. NEED_DEWPOINT) then + if (NEED_TEMP .or. NEED_RELHUM .or. NEED_DEWPOINT .or. need_mslp) then !calculation of temperature at cell centers: do iCell = 1,nCells do k = 1,nVertLevels From dbfcaf96014abfea660be96a566ff5fb6c3cadf3 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 17 Feb 2017 14:58:11 -0700 Subject: [PATCH 17/24] Add namelist option config_extrap_airtemp for extrapolation above/below first-guess levels This commit adds a new namelist option, config_extrap_airtemp, in the new namelist record &interpolation_control to specify the method by which air temperature is extrapolated beyond the first-guess levels. Note that this commit only defines the namelist option in the Registry.xml file, but does not actually add code to make use of this new option. --- src/core_init_atmosphere/Registry.xml | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index 6a6cccb384..415aa1f386 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -174,6 +174,15 @@ + + + + + + Date: Wed, 22 Feb 2017 14:13:47 -0700 Subject: [PATCH 18/24] Add check for lapse-rate extrapolation above first-guess data in vertical_interp This commit adds logic in the vertical_interp function to ensure that we don't use the lapse-rate extrapolation method (extrap=2) above the top of the first-guess dataset; this extrapolation type is only intended to be used to extrapolate downward. --- src/core_init_atmosphere/mpas_init_atm_cases.F | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index a92b1a9cd9..4afe75e9ad 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -4819,7 +4819,7 @@ integer function nearest_edge(target_lat, target_lon, & end function nearest_edge - real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val) + real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val, ierr) implicit none @@ -4827,10 +4827,11 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf integer, intent(in) :: nz real (kind=RKIND), dimension(2,nz), intent(in) :: zf ! zf(1,:) is column of vertical coordinate values, zf(2,:) is column of field values integer, intent(in), optional :: order - integer, intent(in), optional :: extrap + integer, intent(in), optional :: extrap ! can take values 0 = constant, 1 = linear (default), 2 = lapse-rate real (kind=RKIND), intent(in), optional :: surface_val real (kind=RKIND), intent(in), optional :: sealev_val - + integer, intent(out), optional :: ierr + integer :: k, lm, lp real (kind=RKIND) :: wm, wp real (kind=RKIND) :: slope @@ -4838,7 +4839,8 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf integer :: interp_order, extrap_type real (kind=RKIND) :: surface, sealevel - + if (present(ierr)) ierr = 0 + if (present(order)) then interp_order = order else @@ -4883,6 +4885,10 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf else if (extrap_type == 1) then slope = (zf(2,nz) - zf(2,nz-1)) / (zf(1,nz) - zf(1,nz-1)) vertical_interp = zf(2,nz) + slope * (target_z - zf(1,nz)) + else if (extrap_type == 2) then + write(0,*) 'ERROR: extrap_type == 2 not implemented for target_z >= zf(1,nz)' + if (present(ierr)) ierr = 1 + return end if return end if From c6fb303bc52f15e16738b518ff3c3cdbf44ed84f Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 22 Feb 2017 14:18:17 -0700 Subject: [PATCH 19/24] Use 'config_extrap_airtemp' to control vertical extrapolation of temperature When vertically interpolating air temperature, use the extrapolation method specified by the config_extrap_airtemp namelist option to control the behavior when extrapolation is necessary. --- .../mpas_init_atm_cases.F | 28 +++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index 4afe75e9ad..3b06e7b18b 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -2599,6 +2599,9 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state integer, pointer :: config_theta_adv_order real (kind=RKIND), pointer :: config_coef_3rd_order + character (len=StrKIND), pointer :: config_extrap_airtemp + integer :: extrap_airtemp + real (kind=RKIND), dimension(:), pointer :: latCell, lonCell real (kind=RKIND), dimension(:), pointer :: latEdge, lonEdge real (kind=RKIND), dimension(:), pointer :: angleEdge @@ -2670,6 +2673,20 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state call mpas_pool_get_config(configs, 'config_smooth_surfaces', config_smooth_surfaces) call mpas_pool_get_config(configs, 'config_theta_adv_order', config_theta_adv_order) call mpas_pool_get_config(configs, 'config_coef_3rd_order', config_coef_3rd_order) + + call mpas_pool_get_config(configs, 'config_extrap_airtemp', config_extrap_airtemp) + if (trim(config_extrap_airtemp) == 'constant') then + extrap_airtemp = 0 + else if (trim(config_extrap_airtemp) == 'linear') then + extrap_airtemp = 1 + else if (trim(config_extrap_airtemp) == 'lapse-rate') then + extrap_airtemp = 2 + else + call mpas_dmpar_global_abort('*************************************************************', deferredAbort=.true.) + call mpas_dmpar_global_abort('* Invalid value for namelist variable config_extrap_airtemp *', deferredAbort=.true.) + call mpas_dmpar_global_abort('*************************************************************') + end if + write(0,*) "Using option '" // trim(config_extrap_airtemp) // "' for vertical extrapolation of temperature" parinfo => block % parinfo dminfo => block % domain % dminfo @@ -3991,7 +4008,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state ! ! In addition to interpolating wind fields to cell edges, we should - ! also intperolate to cell centers at the surface in order to + ! also interpolate to cell centers at the surface in order to ! produce U10 and V10 fields ! is_sfc_field = .false. @@ -4407,7 +4424,14 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state do k = 1, nVertLevels target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell)) t(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, & - sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=1) + sorted_arr(:,1:nfglevels_actual-1), order=1, & + extrap=extrap_airtemp, ierr=istatus) + if (istatus /= 0) then + write(errstring,'(a,i4,a,i10)') 'Error in interpolation of t(k,iCell) for k=', k, ', iCell=', iCell + call mpas_dmpar_global_abort('*****************************************************************', deferredAbort=.true.) + call mpas_dmpar_global_abort(trim(errstring), deferredAbort=.true.) + call mpas_dmpar_global_abort('*****************************************************************') + end if end do From 24d4a1dfb06f18f307b17221289e9aabe7b5827a Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 31 Mar 2017 20:19:37 +0000 Subject: [PATCH 20/24] Initialize horizontal fluxes for ghost-cell faces in monotonic transport In the atm_advance_scalars_mono_work( ) routine, the flux_arr field holds the horizontal fluxes through cell faces. This array was previously only set for edges that bordered owned cells, leaving the fluxes between ghost cells uninitialized. Later in the routine, however, values from the flux_arr field are used for all edges, and this can lead to floating-point exceptions. The fix in this commit simply sets the flux_arr field to zero for edges that separate two ghost cells in the loop over edges. --- src/core_atmosphere/dynamics/mpas_atm_time_integration.F | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index b2861fe7dc..e44eeeca1f 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -3568,7 +3568,7 @@ subroutine atm_advance_scalars_mono_work(block, state, nCells, nEdges, nVertLeve if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve) then ! only for owned cells - ! speclal treatment of calculations involving edges between hexagonal cells + ! special treatment of calculations involving edges between hexagonal cells ! original code retained in select "default" case ! be sure to see additional declarations near top of subroutine select case(nAdvCellsForEdge(iEdge)) @@ -3602,6 +3602,8 @@ subroutine atm_advance_scalars_mono_work(block, state, nCells, nEdges, nVertLeve end do end select + else + flux_arr(:,iEdge) = 0.0_RKIND end if end do From 1fc9585d26804bf6119167ecf126cdcd4e2c639b Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 19 Apr 2017 15:35:54 -0600 Subject: [PATCH 21/24] Comment-out code in MYNN surface-layer scheme to bound 2-m theta --- src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F b/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F index 554d1f05ec..1584f3c2e1 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F +++ b/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F @@ -1054,10 +1054,10 @@ SUBROUTINE SFCLAY1D_mynn(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,rho1d, & TH2(I)=THGB(I)+DTG*PSIT2/PSIT !*** BE CERTAIN THAT THE 2-M THETA IS BRACKETED BY !*** THE VALUES AT THE SURFACE AND LOWEST MODEL LEVEL. - IF ((TH1D(I)>THGB(I) .AND. (TH2(I)TH1D(I))) .OR. & - (TH1D(I)THGB(I) .OR. TH2(I)THGB(I) .AND. (TH2(I)TH1D(I))) .OR. & +! (TH1D(I)THGB(I) .OR. TH2(I) Date: Thu, 9 Mar 2017 14:09:09 -0700 Subject: [PATCH 22/24] Introduced scaling of the gravity-wave absorbing layer coefficient This commit modifies the formulation of the upper gravity-wave absorbing layer to be scale-aware by multiplying the damping coefficient (config_xnutr) by dx/dx_fine (dx_fine is the mesh spacing within the highest-resolution region of the variable resolution mesh) when setting the damping layer configuration at the beginning of an MPAS integration. The mesh density function used to generate the mesh is also used to determine this ratio. Analysis of the gravity wave absorbing layer characteristics and test results indicate that this scaling is needed for the absorbing layer to effectively damp vertically propagating waves in the coarse-resolution region of variable-resolution meshes. --- src/core_atmosphere/mpas_atm_core.F | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F index 66a2d90242..823692dd41 100644 --- a/src/core_atmosphere/mpas_atm_core.F +++ b/src/core_atmosphere/mpas_atm_core.F @@ -1062,6 +1062,8 @@ subroutine atm_compute_damping_coefs(mesh, configs) real (kind=RKIND), pointer :: config_xnutr, config_zd real (kind=RKIND) :: z, zt, m1, pii real (kind=RKIND), dimension(:,:), pointer :: dss, zgrid + real (kind=RKIND), dimension(:), pointer :: meshDensity + real (kind=RKIND) :: dx_scale_power m1 = -1.0 pii = acos(m1) @@ -1069,12 +1071,14 @@ subroutine atm_compute_damping_coefs(mesh, configs) call mpas_pool_get_dimension(mesh, 'nCells', nCells) call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels) + call mpas_pool_get_array(mesh, 'meshDensity', meshDensity) call mpas_pool_get_array(mesh, 'dss', dss) call mpas_pool_get_array(mesh, 'zgrid', zgrid) call mpas_pool_get_config(configs, 'config_zd', config_zd) call mpas_pool_get_config(configs, 'config_xnutr', config_xnutr) + dx_scale_power = 1.0 dss(:,:) = 0.0 do iCell=1,nCells zt = zgrid(nVertLevels+1,iCell) @@ -1082,6 +1086,7 @@ subroutine atm_compute_damping_coefs(mesh, configs) z = 0.5*(zgrid(k,iCell) + zgrid(k+1,iCell)) if (z > config_zd) then dss(k,iCell) = config_xnutr*sin(0.5*pii*(z-config_zd)/(zt-config_zd))**2.0 + dss(k,iCell) = dss(k,iCell) / meshDensity(iCell)**(0.25*dx_scale_power) end if end do end do From 099032412ed8504839dc895b95af22228b928aae Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Tue, 17 Jan 2017 18:23:14 -0700 Subject: [PATCH 23/24] Introduce a new, scale-aware 3-d divergence damping scheme for acoustic modes This commit replaces the old divergence damping mechanism with a re-engineered filter that is scale-aware. Analysis of the linear behavior of the filter as it was configured for the v5.0 release suggested that it should be stable for all typical MPAS variable-mesh configurations, but certain instances were discovered in which the filter caused the model to become unstable. To stabilize the filter, it has been configured such that, at any place on a variable-resolution mesh, the divergence damping term has the same effective physical eddy viscosity as it would have for a global uniform mesh with the same mesh spacing. Thus, the damping rate in the 3D divergence damping filter is treated in the same scale-aware manner as the other horizontal filters. Also, the application of the 3D divergence damping has been moved immediately after the acoustic step to allow for the application of the filter on the final acoustic step in the Runge-Kutta timesplit integration. Finally, this commit changes the default value of the divergence damping filter coefficient config_smdiv = 0.1. --- src/core_atmosphere/Registry.xml | 10 +- .../dynamics/mpas_atm_time_integration.F | 159 ++++++++++++------ 2 files changed, 112 insertions(+), 57 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 253970d1b0..493127c2c8 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -224,16 +224,11 @@ description="Off-centering parameter for the vertically implicit acoustic integration" possible_values="Positive real values"/> - - - - - diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index b2861fe7dc..4733d39d39 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -184,7 +184,6 @@ subroutine atm_srk3(domain, dt, itimestep) type (field2DReal), pointer :: pressure_p_field type (field2DReal), pointer :: rtheta_p_field type (field2DReal), pointer :: rtheta_pp_field - type (field2DReal), pointer :: divergence_3d_field type (field2DReal), pointer :: tend_u_field type (field2DReal), pointer :: u_field type (field2DReal), pointer :: w_field @@ -729,8 +728,29 @@ subroutine atm_srk3(domain, dt, itimestep) call mpas_pool_get_field(diag, 'rtheta_pp', rtheta_pp_field) call mpas_dmpar_exch_halo_field(rtheta_pp_field, (/ 1 /)) - call mpas_pool_get_field(diag, 'divergence_3d', divergence_3d_field) - call mpas_dmpar_exch_halo_field(divergence_3d_field, (/ 1 /)) +! complete update of horizontal momentum by including 3d divergence damping at the end of the acoustic step + + call mpas_timer_start('atm_divergence_damping_3d') + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', mesh) + call mpas_pool_get_subpool(block % structs, 'state', state) + call mpas_pool_get_subpool(block % structs, 'diag', diag) + + call mpas_pool_get_dimension(block % dimensions, 'nThreads', nThreads) + call mpas_pool_get_dimension(block % dimensions, 'edgeThreadStart', edgeThreadStart) + call mpas_pool_get_dimension(block % dimensions, 'edgeThreadEnd', edgeThreadEnd) + +!$OMP PARALLEL DO + do thread=1,nThreads + call atm_divergence_damping_3d( state, diag, mesh, block % configs, rk_sub_timestep(rk_step), & + edgeThreadStart(thread), edgeThreadEnd(thread) ) + end do +!$OMP END PARALLEL DO + + block => block % next + end do + call mpas_timer_stop('atm_divergence_damping_3d') end do ! end of acoustic steps loop @@ -1834,7 +1854,7 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, wwAvg, rho_pp, cofwt, coftz, zxu, & a_tri, alpha_tri, gamma_tri, dss, & tend_ru, tend_rho, tend_rt, tend_rw, & - zgrid, cofwr, cofwz, w, divergence_3d + zgrid, cofwr, cofwz, w ! redefine ru_p to be perturbation from time t, change 3a ! temporary real (kind=RKIND), dimension(:,:), pointer :: ru @@ -1849,7 +1869,7 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell real (kind=RKIND), dimension(:,:), pointer :: edgesOnCell_sign - real (kind=RKIND), pointer :: epssm, smdiv, smdiv_p_forward + real (kind=RKIND), pointer :: epssm real (kind=RKIND), pointer :: cf1, cf2, cf3 @@ -1876,7 +1896,6 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, call mpas_pool_get_array(diag, 'ruAvg', ruAvg) call mpas_pool_get_array(diag, 'wwAvg', wwAvg) call mpas_pool_get_array(diag, 'rho_pp', rho_pp) - call mpas_pool_get_array(diag, 'divergence_3d', divergence_3d) call mpas_pool_get_array(diag, 'cofwt', cofwt) call mpas_pool_get_array(diag, 'coftz', coftz) call mpas_pool_get_array(diag, 'cofrz', cofrz) @@ -1919,18 +1938,15 @@ subroutine atm_advance_acoustic_step( state, diag, tend, mesh, configs, nCells, call mpas_pool_get_array(diag, 'rw_save', rw_save) ! epssm is the offcentering coefficient for the vertically implicit integration. - ! smdiv is the 3D divergence-damping coefficients. call mpas_pool_get_config(configs, 'config_epssm', epssm) - call mpas_pool_get_config(configs, 'config_smdiv', smdiv) - call mpas_pool_get_config(configs, 'config_smdiv_p_forward', smdiv_p_forward) call atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd, & rho_zz, theta_m, ru_p, rw_p, rtheta_pp, rtheta_pp_old, zz, exner, cqu, ruAvg, wwAvg, & rho_pp, cofwt, coftz, zxu, a_tri, alpha_tri, gamma_tri, dss, tend_ru, tend_rho, tend_rt, & - tend_rw, zgrid, cofwr, cofwz, w, ru, ru_save, rw, rw_save, divergence_3d, fzm, fzp, rdzw, dcEdge, invDcEdge, & + tend_rw, zgrid, cofwr, cofwz, w, ru, ru_save, rw, rw_save, fzm, fzp, rdzw, dcEdge, invDcEdge, & invAreaCell, cofrz, dvEdge, nEdgesOnCell, cellsOnEdge, edgesOnCell, edgesOnCell_sign, & - dts, small_step, epssm, smdiv, smdiv_p_forward, cf1, cf2, cf3 & + dts, small_step, epssm, cf1, cf2, cf3 & ) end subroutine atm_advance_acoustic_step @@ -1940,9 +1956,9 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd, & rho_zz, theta_m, ru_p, rw_p, rtheta_pp, rtheta_pp_old, zz, exner, cqu, ruAvg, wwAvg, & rho_pp, cofwt, coftz, zxu, a_tri, alpha_tri, gamma_tri, dss, tend_ru, tend_rho, tend_rt, & - tend_rw, zgrid, cofwr, cofwz, w, ru, ru_save, rw, rw_save, divergence_3d, fzm, fzp, rdzw, dcEdge, invDcEdge, & + tend_rw, zgrid, cofwr, cofwz, w, ru, ru_save, rw, rw_save, fzm, fzp, rdzw, dcEdge, invDcEdge, & invAreaCell, cofrz, dvEdge, nEdgesOnCell, cellsOnEdge, edgesOnCell, edgesOnCell_sign, & - dts, small_step, epssm, smdiv, smdiv_p_forward, cf1, cf2, cf3 & + dts, small_step, epssm, cf1, cf2, cf3 & ) use mpas_atm_dimensions @@ -1964,7 +1980,6 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart real (kind=RKIND), dimension(nVertLevels,nCells+1) :: rtheta_pp real (kind=RKIND), dimension(nVertLevels,nCells+1) :: rtheta_pp_old - real (kind=RKIND), dimension(nVertLevels,nCells+1) :: divergence_3d real (kind=RKIND), dimension(nVertLevels,nCells+1) :: zz real (kind=RKIND), dimension(nVertLevels,nCells+1) :: exner real (kind=RKIND), dimension(nVertLevels,nEdges+1) :: cqu @@ -2008,7 +2023,7 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart real (kind=RKIND), dimension(maxEdges,nCells+1) :: edgesOnCell_sign integer, intent(in) :: small_step - real (kind=RKIND), intent(in) :: dts, epssm, smdiv, smdiv_p_forward, cf1, cf2, cf3 + real (kind=RKIND), intent(in) :: dts, epssm,cf1, cf2, cf3 real (kind=RKIND), dimension(nVertLevels) :: ts, rs @@ -2027,17 +2042,6 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart if(small_step /= 1) then ! not needed on first small step - do iCell = cellStart,cellEnd - ! acoustic step divergence damping - forward weight rtheta_pp - see Klemp et al MWR 2007 - do k = 1,nVertLevels - rtheta_pp_tmp = rtheta_pp(k,iCell) - rtheta_pp(k,iCell) = (rtheta_pp(k,iCell) + smdiv_p_forward * (rtheta_pp(k,iCell)-rtheta_pp_old(k,iCell)))*zz(k,iCell) - rtheta_pp_old(k,iCell) = rtheta_pp_tmp - end do - end do - -!$OMP BARRIER - ! forward-backward acoustic step integration. ! begin by updating the horizontal velocity u, ! and accumulating the contribution from the updated u to the other tendencies. @@ -2058,12 +2062,10 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP do k=1,nVertLevels -!! pgrad = ((zz_rtheta_pp(k,cell2)-rtheta_pp_old(k,cell1))*invDcEdge(iEdge) )/(.5*(zz(k,cell2)+zz(k,cell1))) pgrad = ((rtheta_pp(k,cell2)-rtheta_pp(k,cell1))*invDcEdge(iEdge) )/(.5*(zz(k,cell2)+zz(k,cell1))) pgrad = cqu(k,iEdge)*0.5*c2*(exner(k,cell1)+exner(k,cell2))*pgrad pgrad = pgrad + 0.5*zxu(k,iEdge)*gravity*(rho_pp(k,cell1)+rho_pp(k,cell2)) - ru_p(k,iEdge) = ru_p(k,iEdge) + dts*(tend_ru(k,iEdge) - pgrad) & - - smdiv*dcEdge(iEdge)*(divergence_3d(k,cell2)-divergence_3d(k,cell1)) + ru_p(k,iEdge) = ru_p(k,iEdge) + dts*(tend_ru(k,iEdge) - pgrad) end do ! accumulate ru_p for use later in scalar transport @@ -2088,11 +2090,10 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart !DIR$ IVDEP do k=1,nVertLevels - ru_p(k,iEdge) = dts*tend_ru(k,iEdge) - smdiv*dcEdge(iEdge)*(tend_rho(k,cell2)-tend_rho(k,cell1)) + ru_p(k,iEdge) = dts*tend_ru(k,iEdge) end do !DIR$ IVDEP do k=1,nVertLevels -!! ruAvg(k,iEdge) = ruAvg(k,iEdge) + ru_p(k,iEdge) ruAvg(k,iEdge) = ru_p(k,iEdge) end do @@ -2101,15 +2102,17 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart end if ! test for first acoustic step -!$OMP BARRIER - if (small_step == 1) then ! initialize here on first small timestep. do iCell=cellStart,cellEnd rtheta_pp_old(1:nVertLevels,iCell) = 0.0 end do + else + do iCell=cellStart,cellEnd + rtheta_pp_old(1:nVertLevels,iCell) = rtheta_pp(1:nVertLevels,iCell) + end do end if -!!!OMP BARRIER -- not needed, since rtheta_pp_old not used below when small_step == 1 +!$OMP BARRIER do iCell=cellSolveStart,cellSolveEnd ! loop over all owned cells to solve @@ -2120,15 +2123,7 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart wwAvg(1:nVertLevels+1,iCell) = 0.0 rho_pp(1:nVertLevels,iCell) = 0.0 rtheta_pp(1:nVertLevels,iCell) = 0.0 -!MGD moved to loop above over all cells -! rtheta_pp_old(1:nVertLevels,iCell) = 0.0 rw_p(:,iCell) = 0.0 - divergence_3d(1:nVertLevels,iCell) = 0. - else ! reset rtheta_pp to input value; - ! rtheta_pp_old stores input value for use in div damping on next acoustic step. - ! Save rho_pp to compute d_rho_pp/dt to get divergence for next acoustic filter application. - rtheta_pp(1:nVertLevels,iCell) = rtheta_pp_old(1:nVertLevels,iCell) - divergence_3d(1:nVertLevels,iCell) = rho_pp(1:nVertLevels,iCell) end if do i=1,nEdgesOnCell(iCell) @@ -2210,7 +2205,6 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart rho_pp(k,iCell) = rs(k) - cofrz(k) *(rw_p(k+1,iCell)-rw_p(k ,iCell)) rtheta_pp(k,iCell) = ts(k) - rdzw(k)*(coftz(k+1,iCell)*rw_p(k+1,iCell) & -coftz(k ,iCell)*rw_p(k ,iCell)) - divergence_3d(k,iCell) = (rho_pp(k,iCell) - divergence_3d(k,iCell))*rdts end do end do ! end of loop over cells @@ -2218,6 +2212,76 @@ subroutine atm_advance_acoustic_step_work(nCells, nEdges, nCellsSolve, cellStart end subroutine atm_advance_acoustic_step_work + subroutine atm_divergence_damping_3d( state, diag, mesh, configs, dts, edgeStart, edgeEnd ) + + ! This subroutine updates the horizontal momentum with the 3d divergence damping component. + + implicit none + + type (mpas_pool_type), intent(inout) :: state + type (mpas_pool_type), intent(inout) :: diag + type (mpas_pool_type), intent(inout) :: mesh + type (mpas_pool_type), intent(in) :: configs + real (kind=RKIND), intent(in) :: dts + integer, intent(in) :: edgeStart, edgeEnd + + real (kind=RKIND), dimension(:,:), pointer :: theta_m, ru_p, rtheta_pp, rtheta_pp_old +! real (kind=RKIND), dimension(:), pointer :: dcEdge + real (kind=RKIND), pointer :: smdiv, config_len_disp + + integer, dimension(:,:), pointer :: cellsOnEdge + integer, pointer :: nCellsSolve + integer, pointer :: nVertLevels + + real (kind=RKIND) :: divCell1, divCell2, rdts, coef_divdamp + integer :: cell1, cell2, iEdge, k + + call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge) +! call mpas_pool_get_array(mesh, 'dcEdge', dcEdge) + call mpas_pool_get_array(state, 'theta_m', theta_m, 1) + call mpas_pool_get_array(diag, 'rtheta_pp', rtheta_pp) + call mpas_pool_get_array(diag, 'rtheta_pp_old', rtheta_pp_old) + call mpas_pool_get_array(diag, 'ru_p', ru_p) + + call mpas_pool_get_dimension(mesh, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels) + + call mpas_pool_get_config(configs, 'config_smdiv', smdiv) + call mpas_pool_get_config(configs, 'config_len_disp', config_len_disp) + + rdts = 1.0_RKIND / dts + coef_divdamp = 2.0_RKIND * smdiv * config_len_disp * rdts + + do iEdge=edgeStart,edgeEnd ! MGD do we really just need edges touching owned cells? + + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + ! update edges for block-owned cells + if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then + +!DIR$ IVDEP + do k=1,nVertLevels + +!! unscaled 3d divergence damping +!! divCell1 = -(rtheta_pp(k,cell1)-rtheta_pp_old(k,cell1))*rdts +!! divCell2 = -(rtheta_pp(k,cell2)-rtheta_pp_old(k,cell2))*rdts +!! ru_p(k,iEdge) = ru_p(k,iEdge) + 2.*smdiv*dcEdge(iEdge)*(divCell2-divCell1) & +!! /(theta_m(k,cell1)+theta_m(k,cell2)) + +!! scaled 3d divergence damping + divCell1 = -(rtheta_pp(k,cell1)-rtheta_pp_old(k,cell1)) + divCell2 = -(rtheta_pp(k,cell2)-rtheta_pp_old(k,cell2)) + ru_p(k,iEdge) = ru_p(k,iEdge) + coef_divdamp*(divCell2-divCell1) & + /(theta_m(k,cell1)+theta_m(k,cell2)) + + end do + end if ! edges for block-owned cells + end do ! end loop over edges + + end subroutine atm_divergence_damping_3d + + subroutine atm_recover_large_step_variables( state, diag, tend, mesh, configs, dt, ns, rk_step, & cellStart, cellEnd, vertexStart, vertexEnd, edgeStart, edgeEnd, & cellSolveStart, cellSolveEnd, vertexSolveStart, vertexSolveEnd, edgeSolveStart, edgeSolveEnd) @@ -3928,7 +3992,7 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, real(kind=RKIND), dimension(:,:), pointer :: tend_w_pgf, tend_w_buoy - real (kind=RKIND), pointer :: coef_3rd_order, c_s, smdiv + real (kind=RKIND), pointer :: coef_3rd_order, c_s logical, pointer :: config_mix_full character (len=StrKIND), pointer :: config_horiz_mixing real (kind=RKIND), pointer :: config_del4u_div_factor @@ -3956,7 +4020,6 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, call mpas_pool_get_config(configs, 'config_visc4_2dsmag', config_visc4_2dsmag) call mpas_pool_get_config(configs, 'config_len_disp', config_len_disp) call mpas_pool_get_config(configs, 'config_smagorinsky_coef', c_s) - call mpas_pool_get_config(configs, 'config_smdiv', smdiv) call mpas_pool_get_array(state, 'rho_zz', rho_zz, 2) call mpas_pool_get_array(state, 'u', u, 2) @@ -4086,7 +4149,7 @@ subroutine atm_compute_dyn_tend(tend, tend_physics, state, diag, mesh, configs, cellsOnEdge, verticesOnEdge, edgesOnCell, edgesOnEdge, cellsOnCell, edgesOnVertex, nEdgesOnCell, nEdgesOnEdge, & latCell, latEdge, angleEdge, u_init, advCellsForEdge, nAdvCellsForEdge, adv_coefs, adv_coefs_3rd, & rdzu, rdzw, fzm, fzp, qv_init, t_init, cf1, cf2, cf3, r_earth, ur_cell, vr_cell, defc_a, defc_b, & - tend_w_pgf, tend_w_buoy, coef_3rd_order, c_s, smdiv, config_mix_full, config_horiz_mixing, config_del4u_div_factor, & + tend_w_pgf, tend_w_buoy, coef_3rd_order, c_s, config_mix_full, config_horiz_mixing, config_del4u_div_factor, & config_h_mom_eddy_visc2, config_v_mom_eddy_visc2, config_h_theta_eddy_visc2, config_v_theta_eddy_visc2, & config_h_theta_eddy_visc4, config_h_mom_eddy_visc4, config_visc4_2dsmag, config_len_disp, rk_step, dt, & tend_rtheta_adv, rthdynten, & @@ -4111,7 +4174,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm cellsOnEdge, verticesOnEdge, edgesOnCell, edgesOnEdge, cellsOnCell, edgesOnVertex, nEdgesOnCell, nEdgesOnEdge, & latCell, latEdge, angleEdge, u_init, advCellsForEdge, nAdvCellsForEdge, adv_coefs, adv_coefs_3rd, & rdzu, rdzw, fzm, fzp, qv_init, t_init, cf1, cf2, cf3, r_earth, ur_cell, vr_cell, defc_a, defc_b, & - tend_w_pgf, tend_w_buoy, coef_3rd_order, c_s, smdiv, config_mix_full, config_horiz_mixing, config_del4u_div_factor, & + tend_w_pgf, tend_w_buoy, coef_3rd_order, c_s, config_mix_full, config_horiz_mixing, config_del4u_div_factor, & config_h_mom_eddy_visc2, config_v_mom_eddy_visc2, config_h_theta_eddy_visc2, config_v_theta_eddy_visc2, & config_h_theta_eddy_visc4, config_h_mom_eddy_visc4, config_visc4_2dsmag, config_len_disp, rk_step, dt, & tend_rtheta_adv, rthdynten, & @@ -4220,7 +4283,7 @@ subroutine atm_compute_dyn_tend_work(nCells, nEdges, nVertices, nVertLevels_dumm real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: tend_w_pgf real (kind=RKIND), dimension(nVertLevels+1,nCells+1) :: tend_w_buoy - real (kind=RKIND) :: coef_3rd_order, c_s, smdiv + real (kind=RKIND) :: coef_3rd_order, c_s logical :: config_mix_full character (len=StrKIND) :: config_horiz_mixing real (kind=RKIND) :: config_del4u_div_factor From 10d5a12975619e1585f8a8f15c72603d37073389 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 12 May 2017 13:38:38 -0600 Subject: [PATCH 24/24] Note in description of config_len_disp that it is also used by 3-d div damping Beginning with the changes merged in 80a1fa0, the config_len_disp namelist option in MPAS-Atmosphere is also used by the 3-d divergence damping mechanism. --- src/core_atmosphere/Registry.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 1bcc69e91a..f89d56c973 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -141,7 +141,7 @@