Skip to content

Commit f5c4af8

Browse files
partial update_bc. Does not compile.
commiting before I rip out the lbc domain see #753 (comment) I think lbc is in the same place as anl_domain Need to run a regional mpas case to get this info
1 parent 70e60cb commit f5c4af8

File tree

3 files changed

+315
-151
lines changed

3 files changed

+315
-151
lines changed

models/mpas_atm/model_mod.f90

Lines changed: 89 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -157,12 +157,18 @@ module model_mod
157157
get_bdy_mask, &
158158
find_closest_cell_center, &
159159
cell_ok_to_interpolate, &
160-
uv_cell_to_edges
160+
uv_increments_cell_to_edges, &
161+
uv_field_cell_to_edges, &
162+
on_boundary_cell, &
163+
on_boundary_edge, &
164+
get_analysis_weight, &
165+
cell_next_to_boundary_edge
161166

162167
public :: update_u_from_reconstruct, &
163168
set_lbc_variables, &
164169
force_u_into_state, &
165-
use_increments_for_u_update
170+
use_increments_for_u_update, &
171+
anl_domid, lbc_domid ! HK todo accessor functions?
166172

167173
! set_lbc_variables sets the lbc_variables string array
168174
! force_u_into_state sets a logical add_u_to_state_list that forces u to be in state
@@ -3794,6 +3800,23 @@ function on_boundary_edge(edgeid)
37943800

37953801
end function on_boundary_edge
37963802

3803+
!------------------------------------------------------------
3804+
!> Determine if this edge is on the boundary
3805+
3806+
function cell_next_to_boundary_edge(edgeid)
3807+
3808+
integer, intent(in) :: edgeid
3809+
logical :: cell_next_to_boundary_edge
3810+
3811+
if (global_grid .or. .not. allocated(bdyMaskEdge)) return
3812+
3813+
cell_next_to_boundary_edge = .false.
3814+
cell_next_to_boundary_edge = on_boundary_cell(cellsOnEdge(1,edgeid))
3815+
cell_next_to_boundary_edge = on_boundary_cell(cellsOnEdge(2,edgeid))
3816+
3817+
end function cell_next_to_boundary_edge
3818+
3819+
37973820
!------------------------------------------------------------
37983821
!> Determine if this edge is the outermost edge (bdyMaskEdge = 7)
37993822

@@ -4354,56 +4377,91 @@ subroutine inside_triangle(t1, t2, t3, r, lat, lon, inside, weights)
43544377
end subroutine inside_triangle
43554378

43564379
!------------------------------------------------------------
4380+
function uv_increments_cell_to_edges(zonal_wind, meridional_wind) result(u_inc)
4381+
4382+
! Project u, v wind increments at cell centers onto the edges.
4383+
! Increments at the outermost edge are set to 0.0 because we do not update/compute
4384+
! uedge in the outermost edge in the regional MPAS.
4385+
4386+
real(r8), intent(in) :: zonal_wind(:,:) ! u wind increments from filter
4387+
real(r8), intent(in) :: meridional_wind(:,:) ! v wind increments from filter
4388+
real(r8) :: u_inc(size(zonal_wind, 1), size(zonal_wind, 2)) ! normal velocity increment on the edges
4389+
4390+
real(r8) :: uedge(size(zonal_wind, 1), size(zonal_wind, 2))
4391+
4392+
if ( .not. module_initialized ) call static_init_model
4393+
4394+
uedge(:,:) = 0.0_r8 ! initialize increments to 0 for outermost edge
4395+
call project_uv_cell_to_edges(zonal_wind, meridional_wind, uedge)
4396+
u_inc = uedge
4397+
4398+
end function uv_increments_cell_to_edges
4399+
4400+
!------------------------------------------------------------------
4401+
4402+
subroutine uv_field_cell_to_edges(zonal_wind, meridional_wind, uedge)
43574403

4358-
subroutine uv_cell_to_edges(zonal_wind, meridional_wind, uedge, full_u)
4404+
! Replaces a edge normal velocity field, by projecting u, v wind fields at cell
4405+
! centers onto the edges.
43594406

4360-
! Project u, v wind (increments) at cell centers onto the edges.
4407+
real(r8), intent(in) :: zonal_wind(:,:) ! u wind from filter
4408+
real(r8), intent(in) :: meridional_wind(:,:) ! v wind from filter
4409+
real(r8), intent(inout):: uedge(:,:) ! normal velocity on the edges
4410+
4411+
integer, parameter :: R3 = 3
4412+
real(r8) :: east(R3,nCells), north(R3,nCells)
4413+
real(r8) :: lonCell_rad(nCells), latCell_rad(nCells)
4414+
integer :: iCell, iEdge, cell1, cell2
4415+
4416+
if ( .not. module_initialized ) call static_init_model
4417+
4418+
call project_uv_cell_to_edges(zonal_wind, meridional_wind, uedge)
4419+
4420+
end subroutine uv_field_cell_to_edges
4421+
4422+
!------------------------------------------------------------------
4423+
4424+
subroutine project_uv_cell_to_edges(zonal_wind, meridional_wind, uedge)
4425+
4426+
! Replaces a edge normal velocity , by projecting u, v wind at cell
4427+
! centers onto the edges.
4428+
! The original field on the outermost edge is left unchanged, because
4429+
! we do not update/compute uedge in the outermost edge in the regional MPAS.
4430+
!
43614431
! FIXME:
43624432
! we can hard-code R3 here since it comes from the (3d) x/y/z cartesian coordinate.
43634433
! We read cellsOnEdge and edgeNormalVectors in get_grid.
4364-
! Here "uedge" is an edge normal wind which can be either a full field or an increment
4365-
! depending on the optional input (full_u).
4366-
! if full_u = .true., then the edge wind is replaced by averaging zonal and meridional
4367-
! winds at cell centers.
4368-
! if full_u does not exist, uedge is the analysis increment from the updated cell-center winds.
4369-
! We do not update/compute uedge in the outermost edge in the regional MPAS.
4434+
!
43704435
! This routine followed the updating part in tend_toEdges in
43714436
! MPAS/src/core_atmosphere/physics/mpas_atmphys_todynamics.F.
43724437

4373-
real(r8), intent(in) :: zonal_wind(:,:) ! u wind updated from filter
4374-
real(r8), intent(in) :: meridional_wind(:,:) ! v wind updated from filter
4375-
real(r8), intent(inout):: uedge(:,:) ! normal velocity (increment) on the edges
4376-
logical, intent(in), optional :: full_u ! compute a full field, not an increment
4438+
real(r8), intent(in) :: zonal_wind(:,:) ! u wind from filter
4439+
real(r8), intent(in) :: meridional_wind(:,:) ! v wind from filter
4440+
real(r8), intent(inout):: uedge(:,:) ! normal velocity on the edges
43774441

4378-
! Local variables
43794442
integer, parameter :: R3 = 3
43804443
real(r8) :: east(R3,nCells), north(R3,nCells)
43814444
real(r8) :: lonCell_rad(nCells), latCell_rad(nCells)
43824445
integer :: iCell, iEdge, cell1, cell2
43834446

43844447
if ( .not. module_initialized ) call static_init_model
43854448

4386-
! Initialization for increments
4387-
if ( .not. present(full_u)) then
4388-
uedge(:,:) = 0.0_r8
4389-
endif
4390-
43914449
! Back to radians (locally)
43924450
lonCell_rad = lonCell*deg2rad
43934451
latCell_rad = latCell*deg2rad
43944452

43954453
! Compute unit vectors in east and north directions for each cell:
43964454
do iCell = 1, nCells
43974455

4398-
east(1,iCell) = -sin(lonCell_rad(iCell))
4399-
east(2,iCell) = cos(lonCell_rad(iCell))
4400-
east(3,iCell) = 0.0_r8
4401-
call r3_normalize(east(1,iCell), east(2,iCell), east(3,iCell))
4456+
east(1,iCell) = -sin(lonCell_rad(iCell))
4457+
east(2,iCell) = cos(lonCell_rad(iCell))
4458+
east(3,iCell) = 0.0_r8
4459+
call r3_normalize(east(1,iCell), east(2,iCell), east(3,iCell))
44024460

4403-
north(1,iCell) = -cos(lonCell_rad(iCell))*sin(latCell_rad(iCell))
4404-
north(2,iCell) = -sin(lonCell_rad(iCell))*sin(latCell_rad(iCell))
4405-
north(3,iCell) = cos(latCell_rad(iCell))
4406-
call r3_normalize(north(1,iCell), north(2,iCell), north(3,iCell))
4461+
north(1,iCell) = -cos(lonCell_rad(iCell))*sin(latCell_rad(iCell))
4462+
north(2,iCell) = -sin(lonCell_rad(iCell))*sin(latCell_rad(iCell))
4463+
north(3,iCell) = cos(latCell_rad(iCell))
4464+
call r3_normalize(north(1,iCell), north(2,iCell), north(3,iCell))
44074465

44084466
enddo
44094467

@@ -4429,11 +4487,11 @@ subroutine uv_cell_to_edges(zonal_wind, meridional_wind, uedge, full_u)
44294487
endif ! if(.not.on_outermost_edge(iEdge)) then
44304488
enddo ! do iEdge = 1, nEdges
44314489

4432-
end subroutine uv_cell_to_edges
4433-
4434-
4490+
end subroutine project_uv_cell_to_edges
4491+
44354492
!------------------------------------------------------------------
44364493

4494+
44374495
!==================================================================
44384496
! The following (private) routines were borrowed from the MPAS code
44394497
!==================================================================

0 commit comments

Comments
 (0)