From d1df097ca3795ef76235e937849555807d711d85 Mon Sep 17 00:00:00 2001 From: Nick Szapiro Date: Tue, 24 Apr 2018 19:08:57 -0500 Subject: [PATCH 01/12] Calculate dtheta_dt_mp by finite difference around microphysics call. Implementing in microphysics_to_MPAS like rt_diabatic_tend doesn't look possible since we don't have the previous qv to recover the previous theta. --- .../mpas_atmphys_driver_microphysics.F | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index 1d6e3235e6..b60ea3b89a 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -299,6 +299,30 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !local variables and arrays: logical:: log_microphysics integer:: i,icell,icount,istep,j,k,kk + +!Calculate dtheta_dt_mp (NS: 2018-04-24) +!Store theta, call microphysics->updates theta, calculate tendency + integer, pointer :: nCells, nVertLevels, index_qv + real(kind=RKIND) :: rvord + real(kind=RKIND), dimension(:,:), pointer :: dtheta_dt_mp, theta_m + real (kind=RKIND), dimension(:,:,:), pointer :: scalars + + call mpas_pool_get_dimension(mesh, 'nCells', nCells) + call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels) + call mpas_pool_get_dimension(state, 'index_qv', index_qv) + + call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev) + call mpas_pool_get_array(diag, 'dtheta_dt_mp', dtheta_dt_mp) + call mpas_pool_get_array(state, 'scalars', scalars, time_lev) + + !initialize to theta before microphysics call. + !I think the physics call actually updates variables, not generate a tendency. + rvord = R_v/R_d + do icell=1,nCells + do k=1,nVertLevels + dtheta_dt_mp(k,icell) = theta_m(k,icell) / (1._RKIND + rvord * scalars(index_qv,k,icell)) + end do + end do !----------------------------------------------------------------------------------------------------------------- !call mpas_log_write('') @@ -414,6 +438,18 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten ! dynamics grid: call microphysics_to_MPAS(configs,mesh,state,time_lev,diag,diag_physics,tend,itimestep,its,ite) +!Calculate dtheta_dt_mp (NS: 2018-04-24) +!From stored dtheta_dt_mp, calculate changeInTheta/timeStep +!TODO: Figure out dt_microp vs dt_dyn + call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev) + call mpas_pool_get_array(diag, 'dtheta_dt_mp', dtheta_dt_mp) + call mpas_pool_get_array(state, 'scalars', scalars, time_lev) + do icell=1,nCells + do k=1,nVertLevels + dtheta_dt_mp(k,icell) = ( theta_m(k,icell)/(1._RKIND+rvord*scalars(index_qv,k,icell)) - dtheta_dt_mp(k,icell) )/(dt_microp) + end do + end do + !... deallocation of all microphysics arrays: !$OMP BARRIER !$OMP MASTER From 87163d90fc8da5959e7293205df2e3f906ff2968 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 22 Sep 2021 19:25:50 +0000 Subject: [PATCH 02/12] Increment version number to 7.2 --- 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 1d28ea2a6d..3884a9cb0f 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -MPAS-v7.1 +MPAS-v7.2 ==== 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 5bce0c36ac..8119a7c0ea 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 815804d8b9..502d56ddb1 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 8e7911e231..ff1b26a23a 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 4a9710e10f..d261ded66f 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 9f889b4dd9..f9944c7d50 100644 --- a/src/core_test/Registry.xml +++ b/src/core_test/Registry.xml @@ -1,5 +1,5 @@ - + From 88333e948e93ba7443c84d83be6710d57ceb6a45 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 6 Oct 2021 17:37:42 -0600 Subject: [PATCH 03/12] Obtain 'rvord' from mpas_constants for computation of dtheta_dt_mp The computation of dtheta_dt_mp in the driver_microphysics routine previously relied on a locally computed value of R_v / R_d, but this constant is available from the mpas_constants module. --- .../physics/mpas_atmphys_driver_microphysics.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index b60ea3b89a..cbb316c37b 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -279,6 +279,8 @@ end subroutine microphysics_init subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,tend,itimestep,its,ite) !================================================================================================================= + use mpas_constants, only : rvord + !input arguments: type(mpas_pool_type),intent(in):: configs type(mpas_pool_type),intent(in):: mesh @@ -303,7 +305,6 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !Calculate dtheta_dt_mp (NS: 2018-04-24) !Store theta, call microphysics->updates theta, calculate tendency integer, pointer :: nCells, nVertLevels, index_qv - real(kind=RKIND) :: rvord real(kind=RKIND), dimension(:,:), pointer :: dtheta_dt_mp, theta_m real (kind=RKIND), dimension(:,:,:), pointer :: scalars @@ -317,7 +318,6 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !initialize to theta before microphysics call. !I think the physics call actually updates variables, not generate a tendency. - rvord = R_v/R_d do icell=1,nCells do k=1,nVertLevels dtheta_dt_mp(k,icell) = theta_m(k,icell) / (1._RKIND + rvord * scalars(index_qv,k,icell)) From 38dfb1f4d4bec9663287f6de623847382bc089c0 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 6 Oct 2021 17:41:49 -0600 Subject: [PATCH 04/12] Remove trailing whitespace and adjust indentation of dtheta_dt_mp computation This commit has no impact on results and simply removes trailing whitespace from several lines in the mpas_atmphys_driver_microphysics module and adjusts the indentation of the code that computes dtheta_dt_microphysics. --- .../mpas_atmphys_driver_microphysics.F | 52 +++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index cbb316c37b..b7406f83c2 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -10,7 +10,7 @@ module mpas_atmphys_driver_microphysics use mpas_kind_types use mpas_pool_routines use mpas_timer, only : mpas_timer_start, mpas_timer_stop - + use mpas_atmphys_constants use mpas_atmphys_init_microphysics use mpas_atmphys_interface @@ -47,7 +47,7 @@ module mpas_atmphys_driver_microphysics ! WRF physics called from microphysics_driver: ! -------------------------------------------- ! * module_mp_kessler : Kessler cloud microphysics. -! * module_mp_thompson: Thompson cloud microphysics. +! * module_mp_thompson: Thompson cloud microphysics. ! * module_mp_wsm6 : WSM6 cloud microphysics. ! ! comments: @@ -58,12 +58,12 @@ module mpas_atmphys_driver_microphysics ! ---------------------------------------- ! * removed call to the Thompson cloud microphysics scheme until scheme is updated to that in WRF revision 3.5. ! Laura D. Fowler (laura@ucar.edu) / 2013-05-29. -! * added subroutine compute_relhum to calculate the relative humidity using the functions rslf and rsif from +! * added subroutine compute_relhum to calculate the relative humidity using the functions rslf and rsif from ! the Thompson cloud microphysics scheme. -! Laura D. Fowler (laura@ucar.edu) / 2013-07-12. +! Laura D. Fowler (laura@ucar.edu) / 2013-07-12. ! * removed the argument tend from the call to microphysics_from_MPAS (not needed). ! Laura D. Fowler (laura@ucar.edu) / 2013-11-07. -! * in call to subroutine wsm6, replaced the variable g (that originally pointed to gravity) with gravity, +! * in call to subroutine wsm6, replaced the variable g (that originally pointed to gravity) with gravity, ! for simplicity. ! Laura D. Fowler (laura@ucar.edu) / 2014-03-21. ! * throughout the sourcecode, replaced all "var_struct" defined arrays by local pointers. @@ -270,7 +270,7 @@ subroutine microphysics_init(dminfo,configs,mesh,sfc_input,diag_physics) call wsm6init(rho_a,rho_r,rho_s,cliq,cpv,hail_opt,.false.) case default - + end select microp_select end subroutine microphysics_init @@ -288,7 +288,7 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten integer,intent(in):: time_lev integer,intent(in):: itimestep integer,intent(in):: its,ite - + !inout arguments: type(mpas_pool_type),intent(inout):: state type(mpas_pool_type),intent(inout):: diag @@ -301,7 +301,7 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !local variables and arrays: logical:: log_microphysics integer:: i,icell,icount,istep,j,k,kk - + !Calculate dtheta_dt_mp (NS: 2018-04-24) !Store theta, call microphysics->updates theta, calculate tendency integer, pointer :: nCells, nVertLevels, index_qv @@ -311,17 +311,17 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten call mpas_pool_get_dimension(mesh, 'nCells', nCells) call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(state, 'index_qv', index_qv) - + call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev) call mpas_pool_get_array(diag, 'dtheta_dt_mp', dtheta_dt_mp) call mpas_pool_get_array(state, 'scalars', scalars, time_lev) - - !initialize to theta before microphysics call. + + !initialize to theta before microphysics call. !I think the physics call actually updates variables, not generate a tendency. do icell=1,nCells - do k=1,nVertLevels - dtheta_dt_mp(k,icell) = theta_m(k,icell) / (1._RKIND + rvord * scalars(index_qv,k,icell)) - end do + do k=1,nVertLevels + dtheta_dt_mp(k,icell) = theta_m(k,icell) / (1._RKIND + rvord * scalars(index_qv,k,icell)) + end do end do !----------------------------------------------------------------------------------------------------------------- @@ -344,7 +344,7 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !... call to different cloud microphysics schemes: microp_select: select case(microp_scheme) - + case ("mp_kessler") call mpas_timer_start('Kessler') call kessler( & @@ -369,7 +369,7 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten th = th_p , qv = qv_p , qc = qc_p , & qr = qr_p , qi = qi_p , qs = qs_p , & qg = qg_p , ni = ni_p , nr = nr_p , & - pii = pi_p , p = pres_p , dz = dz_p , & + pii = pi_p , p = pres_p , dz = dz_p , & w = w_p , dt_in = dt_microp , itimestep = itimestep , & rainnc = rainnc_p , rainncv = rainncv_p , snownc = snownc_p , & snowncv = snowncv_p , graupelnc = graupelnc_p , graupelncv = graupelncv_p , & @@ -409,12 +409,12 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten call mpas_timer_stop('WSM6') case default - + end select microp_select !... calculate the 10cm radar reflectivity and relative humidity, if needed: if (l_diags) then - + !ensure that we only call compute_radar_reflectivity() if we are using an MPS that supports !the computation of simulated radar reflectivity: if(trim(microp_scheme) == "mp_wsm6" .or. & @@ -445,9 +445,9 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten call mpas_pool_get_array(diag, 'dtheta_dt_mp', dtheta_dt_mp) call mpas_pool_get_array(state, 'scalars', scalars, time_lev) do icell=1,nCells - do k=1,nVertLevels - dtheta_dt_mp(k,icell) = ( theta_m(k,icell)/(1._RKIND+rvord*scalars(index_qv,k,icell)) - dtheta_dt_mp(k,icell) )/(dt_microp) - end do + do k=1,nVertLevels + dtheta_dt_mp(k,icell) = ( theta_m(k,icell)/(1._RKIND+rvord*scalars(index_qv,k,icell)) - dtheta_dt_mp(k,icell) )/(dt_microp) + end do end do !... deallocation of all microphysics arrays: @@ -475,7 +475,7 @@ subroutine precip_from_MPAS(configs,diag_physics,its,ite) !local pointers: character(len=StrKIND),pointer:: microp_scheme integer,pointer:: nCellsSolve - real,dimension(:),pointer:: graupelncv,rainncv,snowncv,sr + real,dimension(:),pointer:: graupelncv,rainncv,snowncv,sr !local variables and arrays: integer:: i,j @@ -519,7 +519,7 @@ subroutine precip_from_MPAS(configs,diag_physics,its,ite) snowncv(i) = 0._RKIND graupelncv(i) = 0._RKIND sr(i) = 0._RKIND - enddo + enddo case default @@ -567,7 +567,7 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) call mpas_pool_get_array(diag_physics,'sr' ,sr ) do i = its,ite - precipw(i) = 0._RKIND + precipw(i) = 0._RKIND enddo !variables common to all cloud microphysics schemes: @@ -582,7 +582,7 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) !time-step precipitation: rainncv(i) = rainnc_p(i,j) - + !accumulated precipitation: rainnc(i) = rainnc(i) + rainncv(i) @@ -591,7 +591,7 @@ subroutine precip_to_MPAS(configs,diag_physics,its,ite) i_rainnc(i) = i_rainnc(i) + 1 rainnc(i) = rainnc(i) - config_bucket_rainnc endif - + enddo enddo From db4f0d020449f69281fdc4eb381ae11fee97e436 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 6 Oct 2021 17:51:32 -0600 Subject: [PATCH 05/12] Use dt_dyn rather than dt_microp when computing dtheta_dt_mp Because the finite-difference computation of dtheta_dt_mp happens outside of the microphysics subcycling loop, dt_dyn, rather than dt_microp, should be used when computing dtheta_dt_mp. --- src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index b7406f83c2..01a8a3252e 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -440,13 +440,12 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !Calculate dtheta_dt_mp (NS: 2018-04-24) !From stored dtheta_dt_mp, calculate changeInTheta/timeStep -!TODO: Figure out dt_microp vs dt_dyn call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev) call mpas_pool_get_array(diag, 'dtheta_dt_mp', dtheta_dt_mp) call mpas_pool_get_array(state, 'scalars', scalars, time_lev) do icell=1,nCells do k=1,nVertLevels - dtheta_dt_mp(k,icell) = ( theta_m(k,icell)/(1._RKIND+rvord*scalars(index_qv,k,icell)) - dtheta_dt_mp(k,icell) )/(dt_microp) + dtheta_dt_mp(k,icell) = (theta_m(k,icell)/(1._RKIND+rvord*scalars(index_qv,k,icell)) - dtheta_dt_mp(k,icell))/(dt_dyn) end do end do From 3e06ce92b02b7ddb813668e7ca0c4df5a18f32e4 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 6 Oct 2021 17:58:19 -0600 Subject: [PATCH 06/12] Remove redundant mpas_pool_get_array calls for dtheta_dt_mp computation The calls to get array pointers for theta_m, dtheta_dt_mp, and scalars arrays immediately before the computation of dtheta_dt_mp in driver_microphysics are unnecessary since time_lev is invariant and similar calls are already made earlier in the routine. --- src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F | 3 --- 1 file changed, 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 01a8a3252e..27fc07e768 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F @@ -440,9 +440,6 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten !Calculate dtheta_dt_mp (NS: 2018-04-24) !From stored dtheta_dt_mp, calculate changeInTheta/timeStep - call mpas_pool_get_array(state, 'theta_m', theta_m, time_lev) - call mpas_pool_get_array(diag, 'dtheta_dt_mp', dtheta_dt_mp) - call mpas_pool_get_array(state, 'scalars', scalars, time_lev) do icell=1,nCells do k=1,nVertLevels dtheta_dt_mp(k,icell) = (theta_m(k,icell)/(1._RKIND+rvord*scalars(index_qv,k,icell)) - dtheta_dt_mp(k,icell))/(dt_dyn) From 883ab4f0ab62d6b47ccbfb4b0286c35e290fd862 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 13 Oct 2021 20:46:28 +0000 Subject: [PATCH 07/12] Initialize qvb to zero before its first use in init_atm_case_squall_line When computing the reference sounding in the squall line/supercell initialization, the qvb array was used before it has been given an initial value. A simple zero initial value for qvb is appropriate (and it's almost certainly what was assumed by the code). --- src/core_init_atmosphere/mpas_init_atm_cases.F | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index e08d55e48e..a31f0829fe 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -1652,6 +1652,8 @@ subroutine init_atm_case_squall_line(dminfo, mesh, nCells, nVertLevels, state, d ! ! for reference sounding ! + qvb(:) = 0.0_RKIND + do itr=1,30 pitop = 1.-.5*dzw(1)*gravity*(1.+scalars(index_qv,1,1))/(cp*t(1,1)*zz(1,1)) From f232665e5b7bb8a5484d22cca601a2b26527d311 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 13 Oct 2021 16:28:21 -0600 Subject: [PATCH 08/12] Fix bug in vertical interp of r.h. and s.h. when levels are given top-to-bottom The code in the real-data initialization routine for vertically interpolating relative humidity and specific humidity assumed that first-guess levels would be given in bottom-to-top order when attempting to vertically extrapolate to model levels below the lowest first-guess level. The relevant code for relative humidity read as follows -- the code for specific humidity is similar. if (target_z < z_fg(1,iCell) .and. k < nVertLevels) relhum(k,iCell) = relhum(k+1,iCell) If first-guess levels are not given in bottom-to-top order, then z_fg(1,iCell) does not necessarily contain the height of the surface in the first-guess data, resulting in a copy of vertically interpolated relative humidity level k+1 to level k. One possible fix for this issue might be to compare target_z with sorted_arr(1,1). Since sorted_arr is always sorted in ascending order by its first index, the check to decide when to copy/extrapolate relative humidity would be independent of the order in which first-guess levels were provided. If the comparison were to be made against sorted_arr(1,1) rather than against z_fg(1,iCell), the effect would be to copy the lowest *interpolated* level downward to model levels below the first-guess ground. An alternative fix adopted in this commit is to simply delete the logic to explicitly copy/extrapolate downward, since both relative humidity and specific humidity are vertically interpolated with extrap=0, which specifies constant- value extrapolation. In this case, the lowest *first-guess* level -- rather than the lowest *interpolated* level -- is copied downward. Note: The issue being addressed by this commit is a dependence on the order in which levels are given in the input intermediate file, and not on the direction in which the first-guess model indexes its levels. --- src/core_init_atmosphere/mpas_init_atm_cases.F | 2 -- 1 file changed, 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 e08d55e48e..cfd3f3d27b 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -4563,7 +4563,6 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell)) relhum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, & sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0) - if (target_z < z_fg(1,iCell) .and. k < nVertLevels) relhum(k,iCell) = relhum(k+1,iCell) end do @@ -4581,7 +4580,6 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state target_z = 0.5 * (zgrid(k,iCell) + zgrid(k+1,iCell)) spechum(k,iCell) = vertical_interp(target_z, nfglevels_actual-1, & sorted_arr(:,1:nfglevels_actual-1), order=1, extrap=0) - if (target_z < z_fg(1,iCell) .and. k < nVertLevels) spechum(k,iCell) = spechum(k+1,iCell) end do From 7e514a379276e464040126065bb1137e1a5ead5c Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Wed, 13 Oct 2021 18:39:23 -0600 Subject: [PATCH 09/12] Correct check on availability of dtheta_dt_mix when computing depv_dt_mix The code to compute depv_dt_mix in the calc_pvBudget routine previously checked that the dtheta_dt_mp variable had been successfully been retrieved from the diag pool, when instead the check should have been on the dtheta_dt_mix variable. Note that the fix in this commit isn't expected to have any impact in practice, since both dtheta_dt_mp and dtheta_dt_mix should always be found in the diag pool. --- src/core_atmosphere/diagnostics/pv_diagnostics.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_atmosphere/diagnostics/pv_diagnostics.F b/src/core_atmosphere/diagnostics/pv_diagnostics.F index da3d8fa605..1db3429c6e 100644 --- a/src/core_atmosphere/diagnostics/pv_diagnostics.F +++ b/src/core_atmosphere/diagnostics/pv_diagnostics.F @@ -1457,7 +1457,7 @@ subroutine calc_pvBudget(state, time_lev, diag, mesh, tend, tend_physics) depv_dt_mp(k,iCell) = 0.0_RKIND end if - if (associated(dtheta_dt_mp)) then + if (associated(dtheta_dt_mix)) then call calc_grad_cell(gradtheta, & iCell, k, nVertLevels, nEdgesOnCell(iCell), verticesOnCell, kiteAreasOnVertex, & cellsOnCell, edgesOnCell, cellsOnEdge, dvEdge, edgeNormalVectors, & From 6d4e344e278a026489bf77c3066fdaee3b8e07d7 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Tue, 19 Oct 2021 23:38:54 +0000 Subject: [PATCH 10/12] Fix parallel reproducibility error in 'smstav' field This commit fixes reproducibility issues in several fields within the Noah LSM over land-ice points when running with different MPI task counts; however, only one of these fields -- smstav, the surface moisture availability field -- persists outside of the physics driver and is written to MPAS-Atmosphere restart files. Within the lsm and lsm_mosaic routines, the soilw, flx4, fvb, fbur, and fgsn variables are not set for land-ice points in the i-loop over grid points, and at the bottom of the i-loop, where 2-d array elements are set based on these variables, the array elements then take on values from the last non-land-ice grid point. For different horizontal domain decompositions, the last non-land-ice point will in general vary, leading to different patterns over land-ice points. The fix for this issue involves initializing soilw, flx4, fvb, fbur, and fgsn at the beginning of each i-loop iteration. Since the smstav field is a diagnostic, correcting this reproducibility error has no other effect on model simulations. Note that, although the lsm_mosaic routine has been fixed, this routine is not actually compiled by MPAS-Atmosphere. The changes in this commit have also been committed to the WRF model repository. --- .../physics/physics_wrf/module_sf_noahdrv.F | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F b/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F index 3dde0130d4..24ec87ea97 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F +++ b/src/core_atmosphere/physics/physics_wrf/module_sf_noahdrv.F @@ -666,12 +666,6 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & ! ! END FASDAS ! - FLX4 = 0.0 !BSINGH - Initialized to 0.0 - FVB = 0.0 !BSINGH - Initialized to 0.0 - FBUR = 0.0 !BSINGH - Initialized to 0.0 - FGSN = 0.0 !BSINGH - Initialized to 0.0 - SOILW = 0.0 !BSINGH - Initialized to 0.0 - sigma_sb=5.67e-08 ! MEK MAY 2007 @@ -743,6 +737,12 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5 ! convert from mixing ratio to specific humidity Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j)) +! initializing local variables + SOILW=0. + FLX4=0. + FVB=0. + FBUR=0. + FGSN=0. ! ! Q2SAT=QGH(I,j) Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity @@ -2962,6 +2962,13 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & !----------------------------------------------------------------------- ILOOP : DO I=its,ite +! initializing local variables + SOILW = 0. + FLX4 = 0. + FVB = 0. + FBUR = 0. + FGSN = 0. + IF (((XLAND(I,J)-1.5).LT.0.) .AND. (XICE(I,J) < XICE_THRESHOLD) ) THEN IVGTYP_dominant(I,J)=IVGTYP(I,J) ! save this From 3d2086e27a7c009b00ece1073f3c297e69cf6145 Mon Sep 17 00:00:00 2001 From: Nick Szapiro Date: Sun, 26 Feb 2017 23:59:18 -0600 Subject: [PATCH 11/12] Halo exchange inTropo in flood fill to find DT. Originally, it was assumed that each (MPI) domain would have >0 cells with "right" DT found by flood filling w/in domain, and that would be sufficient to connect each domain's troposphere. However, for "small" domains over the Arctic say during winter w/ strong inversions, the entire surface can be capped by high PV. So, we need to communicate between domains during the flood fill or else we can find a wrong DT near the surface. Only in very rare edge cases should more than two halo exchanges occur. --- .../diagnostics/Registry_pv.xml | 3 + .../diagnostics/pv_diagnostics.F | 125 +++++++++++------- 2 files changed, 79 insertions(+), 49 deletions(-) diff --git a/src/core_atmosphere/diagnostics/Registry_pv.xml b/src/core_atmosphere/diagnostics/Registry_pv.xml index fdf5d3b674..28dc7f5138 100644 --- a/src/core_atmosphere/diagnostics/Registry_pv.xml +++ b/src/core_atmosphere/diagnostics/Registry_pv.xml @@ -62,5 +62,8 @@ + + diff --git a/src/core_atmosphere/diagnostics/pv_diagnostics.F b/src/core_atmosphere/diagnostics/pv_diagnostics.F index 0d48c036ea..a0a8365a78 100644 --- a/src/core_atmosphere/diagnostics/pv_diagnostics.F +++ b/src/core_atmosphere/diagnostics/pv_diagnostics.F @@ -591,29 +591,41 @@ subroutine floodFill_tropo(mesh, diag, pvuVal) ! (2) flood fill troposphere (<2pvu) from troposphere seeds near surface. !Somewhat paradoxically, the bottom of the stratosphere is lower than the top of the troposphere. + !Originally, it was assumed that each (MPI) domain would have >0 cells with "right" DT found by flood filling. + !However, for "small" domains over the Arctic say during winter, the entire surface can be capped by high PV. + !So, we need to communicate between domains during the flood fill or else we find the DT at the surface. + !The extreme limiting case is if we had every cell as its own domain; then, it's clear that there has to be communication. + !The "output" is iLev_DT, which is the vertical index for the level >= pvuVal. If >nVertLevels, pvuVal above column. If <2, pvuVal below column. !Communication between blocks during the flood fill may be needed to treat some edge cases appropriately. - use mpas_pool_routines, only : mpas_pool_get_dimension, mpas_pool_get_array + use mpas_pool_routines, only : mpas_pool_get_dimension, mpas_pool_get_array, mpas_pool_get_field + use mpas_dmpar, only : mpas_dmpar_max_int,mpas_dmpar_exch_halo_field + use mpas_derived_types, only : dm_info, field2DInteger implicit none type (mpas_pool_type), intent(in) :: mesh type (mpas_pool_type), intent(inout) :: diag real(kind=RKIND), intent(in) :: pvuVal - - integer :: iCell, k, nChanged, iNbr, iCellNbr, levInd - integer, pointer :: nCells, nVertLevels + + integer :: iCell, k, nChanged, iNbr, iCellNbr, levInd, haloChanged, global_haloChanged + integer, pointer :: nCells, nVertLevels, nCellsSolve integer, dimension(:), pointer :: nEdgesOnCell, iLev_DT - integer, dimension(:,:), pointer :: cellsOnCell - + integer, dimension(:,:), pointer :: cellsOnCell, inTropo + + type (field2DInteger), pointer :: inTropo_f + real(kind=RKIND) :: sgnHemi, sgn real(kind=RKIND),dimension(:),pointer:: latCell real(kind=RKIND), dimension(:,:), pointer :: ertel_pv - integer, dimension(:,:), allocatable :: candInTropo, inTropo !whether in troposphere + type (dm_info), pointer :: dminfo + + integer, dimension(:,:), allocatable :: candInTropo !whether in troposphere call mpas_pool_get_dimension(mesh, 'nCells', nCells) + call mpas_pool_get_dimension(mesh, 'nCellsSolve', nCellsSolve) call mpas_pool_get_dimension(mesh, 'nVertLevels', nVertLevels) call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell) @@ -622,9 +634,9 @@ subroutine floodFill_tropo(mesh, diag, pvuVal) call mpas_pool_get_array(diag, 'ertel_pv', ertel_pv) !call mpas_pool_get_array(diag, 'iLev_DT_trop', iLev_DT) call mpas_pool_get_array(diag, 'iLev_DT', iLev_DT) + call mpas_pool_get_array(diag, 'inTropo', inTropo) allocate(candInTropo(nVertLevels, nCells+1)) - allocate(inTropo(nVertLevels, nCells+1)) candInTropo(:,:) = 0 inTropo(:,:) = 0 !store whether each level above DT to avoid repeating logic. we'll use cand as a isVisited marker further below. @@ -645,7 +657,7 @@ subroutine floodFill_tropo(mesh, diag, pvuVal) do k=1,levInd if (candInTropo(k,iCell) .GT. 0) then inTropo(k,iCell) = 1 - candInTropo(k,iCell) = 0 + !candInTropo(k,iCell) = 0 nChanged = nChanged+1 end if end do @@ -653,48 +665,63 @@ subroutine floodFill_tropo(mesh, diag, pvuVal) !flood fill from the given seeds. since I don't know enough fortran, !we'll just brute force a continuing loop rather than queue. - do while(nChanged .GT. 0) - nChanged = 0 - do iCell=1,nCells - do k=1,nVertLevels - !update if candidate and neighbor in troposphere - if (candInTropo(k,iCell) .GT. 0) then - !nbr below - if (k .GT. 1) then - if (inTropo(k-1,iCell) .GT. 0) then - inTropo(k,iCell) = 1 - candInTropo(k,iCell) = 0 - nChanged = nChanged+1 - cycle + call mpas_pool_get_field(diag, 'inTropo', inTropo_f) + dminfo => inTropo_f % block % domain % dminfo + global_haloChanged = 1 + do while(global_haloChanged .GT. 0) !any cell in a halo has changed, to propagate to other domains + global_haloChanged = 0 !aggregate the number of changed cells w/in the loop below + do while(nChanged .GT. 0) + nChanged = 0 + do iCell=1,nCells !should we look for neighbors of hallo cells? + !do iCell=1,nCellsSolve !should we look for neighbors of hallo cells? + do k=1,nVertLevels + !update if candidate and neighbor in troposphere + if ((candInTropo(k,iCell) .GT. 0) .AND. (inTropo(k,iCell).LT.1) ) then + !nbr below + if (k .GT. 1) then + if (inTropo(k-1,iCell) .GT. 0) then + inTropo(k,iCell) = 1 + !candInTropo(k,iCell) = 0 + nChanged = nChanged+1 + cycle + end if end if - end if - - !side nbrs - do iNbr = 1, nEdgesOnCell(iCell) - iCellNbr = cellsOnCell(iNbr,iCell) - if (inTropo(k,iCellNbr) .GT. 0) then - inTropo(k,iCell) = 1 - candInTropo(k,iCell) = 0 - nChanged = nChanged+1 - cycle - end if - end do - - !nbr above - if (k .LT. nVertLevels) then - if (inTropo(k+1,iCell) .GT. 0) then - inTropo(k,iCell) = 1 - candInTropo(k,iCell) = 0 - nChanged = nChanged+1 - cycle + + !side nbrs + do iNbr = 1, nEdgesOnCell(iCell) + iCellNbr = cellsOnCell(iNbr,iCell) + if (inTropo(k,iCellNbr) .GT. 0) then + inTropo(k,iCell) = 1 + !candInTropo(k,iCell) = 0 + nChanged = nChanged+1 + exit + end if + end do + + !nbr above + if (k .LT. nVertLevels) then + if (inTropo(k+1,iCell) .GT. 0) then + inTropo(k,iCell) = 1 + !candInTropo(k,iCell) = 0 + nChanged = nChanged+1 + cycle + end if end if - end if - - end if !candIn - end do !levels - end do !cells - !here's where a communication would be needed for edge cases !!! - end do !while + + end if !candIn + end do !levels + end do !cells + global_haloChanged = global_haloChanged+nChanged + end do !while w/in domain + !communicate to other domains for edge case where a chunk of a block hasn't gotten to fill + nChanged = global_haloChanged + call mpas_dmpar_max_int(dminfo, nChanged, global_haloChanged) + if (global_haloChanged .GT. 0) then !communicate inTropo everywhere + call mpas_dmpar_exch_halo_field(inTropo_f) + end if + nChanged = global_haloChanged !so each block will iterate again if anything changed + end do !while haloChanged + deallocate(candInTropo) !Fill iLev_DT with the lowest level above the tropopause (If DT above column, iLev>nVertLevels. If DT below column, iLev=0. do iCell=1,nCells From 6065bd472d8caa012afe09321bd32ff79ec13cd5 Mon Sep 17 00:00:00 2001 From: Michael Duda Date: Fri, 14 Jan 2022 22:04:46 +0000 Subject: [PATCH 12/12] Correct units and descriptions for GWDO fields var2d, con, oa[1-4], ol[1-4] This commit corrects the units and description attributes for various GWDO fields in both the init_atmosphere and atmosphere core Registry.xml files. --- src/core_atmosphere/Registry.xml | 24 ++++++++++++------------ src/core_init_atmosphere/Registry.xml | 26 +++++++++++++------------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index e4cf270af4..9b393e4146 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -3067,35 +3067,35 @@ - + - + + description="asymmetry of subgrid-scale orography for westerly flow"/> + description="asymmetry of subgrid-scale orography for southerly flow"/> + description="asymmetry of subgrid-scale orography for south-westerly flow"/> + description="asymmetry of subgrid-scale orography for north-westerly flow"/> + description="effective orographic length for westerly flow"/> + description="effective orographic length for southerly flow"/> + description="effective orographic length for south-westerly flow"/> + description="effective orographic length for north-westerly flow"/> diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index 515b881d01..e74592abb9 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -245,7 +245,7 @@ - + - + + description="asymmetry of subgrid-scale orography for westerly flow"/> + description="asymmetry of subgrid-scale orography for southerly flow"/> + description="asymmetry of subgrid-scale orography for south-westerly flow"/> + description="asymmetry of subgrid-scale orography for north-westerly flow"/> + description="effective orographic length for westerly flow"/> + description="effective orographic length for southerly flow"/> + description="effective orographic length for south-westerly flow"/> + description="effective orographic length for north-westerly flow"/>