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. 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" 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..f89d56c973 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1,5 +1,5 @@ - + @@ -141,7 +141,7 @@ - - - - - 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 diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index b2861fe7dc..25ba38e1da 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) @@ -3568,7 +3632,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 +3666,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 @@ -3928,7 +3994,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 +4022,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 +4151,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 +4176,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 +4285,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 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 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 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 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 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) - + @@ -174,6 +174,15 @@ + + + + + + block % parinfo dminfo => block % domain % dminfo @@ -3336,6 +3353,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,7 +3363,9 @@ 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) == 'PRESSURE' .or. & trim(field % field) == 'SNOW' .or. & trim(field % field) == 'SEAICE' .or. & trim(field % field) == 'SKINTEMP') then @@ -3359,6 +3379,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,13 +3389,14 @@ 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 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 @@ -3510,6 +3532,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 @@ -3711,6 +3740,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 +3949,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' @@ -3991,7 +4061,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 +4477,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 @@ -4819,7 +4896,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 +4904,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 +4916,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 +4962,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 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_landice/mpas_li_diagnostic_vars.F b/src/core_landice/mpas_li_diagnostic_vars.F index 8186e846b1..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 (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 (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 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 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 @@ - + 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/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/Makefile b/src/external/esmf_time_f90/Makefile index 192a52c2b6..790852e6d7 100644 --- a/src/external/esmf_time_f90/Makefile +++ b/src/external/esmf_time_f90/Makefile @@ -53,8 +53,13 @@ wrf_error_fatal.o: wrf_message.o: clean: - rm -rf *.o *.mod *.a + rm -rf *.tmp.f90 *.o *.mod *.a .F90.o: $(RM) $@ $*.mod +ifeq "$(GEN_F90)" "true" + $(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 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