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..68fe74c8cd 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -1,5 +1,5 @@ - + @@ -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_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 da3d8fa605..e217323947 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 - 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 + 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 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 @@ -1457,7 +1484,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, & diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F b/src/core_atmosphere/physics/mpas_atmphys_driver_microphysics.F index 1d6e3235e6..27fc07e768 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 @@ -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 @@ -286,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 @@ -300,6 +302,28 @@ subroutine driver_microphysics(configs,mesh,state,time_lev,diag,diag_physics,ten 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), 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. + 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('') !call mpas_log_write('---enter subroutine driver_microphysics:') @@ -320,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( & @@ -345,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 , & @@ -385,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. & @@ -414,6 +438,14 @@ 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 + 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) + end do + end do + !... deallocation of all microphysics arrays: !$OMP BARRIER !$OMP MASTER @@ -439,7 +471,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 @@ -483,7 +515,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 @@ -531,7 +563,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: @@ -546,7 +578,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) @@ -555,7 +587,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 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 diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index 815804d8b9..8d80216aaa 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -1,5 +1,5 @@ - + @@ -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"/> - + 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 @@ - +