Skip to content

Commit

Permalink
chore: move non-public routines after public routines
Browse files Browse the repository at this point in the history
  • Loading branch information
@hkershaw-brown
hkershaw-brown committed Sep 12, 2023
1 parent 2e8d4b3 commit 216a3f4
Showing 1 changed file with 151 additions and 153 deletions.
304 changes: 151 additions & 153 deletions models/wrf/model_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -444,154 +444,6 @@ subroutine model_interpolate(state_handle, ens_size, location, qty_in, expected_

end subroutine model_interpolate


!------------------------------------------------------------------
function force_non_negative_if_required(n, qty, fld)

integer, intent(in) :: n
integer, intent(in) :: qty
real(r8), intent(in) :: fld(n)

real(r8) :: force_non_negative_if_required(n)

select case (qty)
case (QTY_RAINWATER_MIXING_RATIO, &
QTY_GRAUPEL_MIXING_RATIO, &
QTY_HAIL_MIXING_RATIO, &
QTY_SNOW_MIXING_RATIO, &
QTY_ICE_MIXING_RATIO, &
QTY_CLOUDWATER_MIXING_RATIO, &
QTY_DROPLET_NUMBER_CONCENTR, &
QTY_ICE_NUMBER_CONCENTRATION, &
QTY_SNOW_NUMBER_CONCENTR, &
QTY_RAIN_NUMBER_CONCENTR, &
QTY_GRAUPEL_NUMBER_CONCENTR, &
QTY_HAIL_NUMBER_CONCENTR )
force_non_negative_if_required = max(0.0_r8, fld)
case default
force_non_negative_if_required = fld
end select

end function force_non_negative_if_required

!------------------------------------------------------------------
function surface_qty(qty)

integer, intent(in) :: qty
logical :: surface_qty

select case (qty)

case (QTY_2M_TEMPERATURE, &
QTY_2M_SPECIFIC_HUMIDITY, &
QTY_10M_U_WIND_COMPONENT, &
QTY_10M_V_WIND_COMPONENT, &
QTY_SURFACE_TYPE, &
QTY_SKIN_TEMPERATURE)
surface_qty = .true.
case default
surface_qty = .false.

end select

end function surface_qty
!------------------------------------------------------------------
! WRF has separate model qtys for surface variables
function update_qty_if_location_is_surface(qty_in, location) result(qty)

integer, intent(in) :: qty_in
type(location_type), intent(in) :: location
integer :: qty

if (.not. is_vertical(location,"SURFACE")) then
qty = qty_in
return
endif

select case (qty)

case (QTY_U_WIND_COMPONENT); qty = QTY_10M_U_WIND_COMPONENT
case (QTY_V_WIND_COMPONENT); qty = QTY_10M_V_WIND_COMPONENT
case (QTY_POTENTIAL_TEMPERATURE); qty = QTY_2M_TEMPERATURE
case (QTY_SPECIFIC_HUMIDITY); qty = QTY_2M_SPECIFIC_HUMIDITY
case (QTY_VAPOR_MIXING_RATIO); qty = QTY_2M_SPECIFIC_HUMIDITY ! Vapor Mixing Ratio (QV, Q2)

end select

end function update_qty_if_location_is_surface

!------------------------------------------------------------------
function simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym)

integer, intent(in) :: ens_size
type(ensemble_type), intent(in) :: state_handle
integer, intent(in) :: qty
integer, intent(in) :: id
integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners
integer, intent(in) :: k(ens_size) ! k may be different across the ensemble
real(r8), intent(in) :: dxm, dx, dy, dym

real(r8) :: simple_interpolation(ens_size)

integer :: e ! loop variable
! lower left, upper left, lower right, upper right
integer(i8), dimension(ens_size) :: ill, iul, ilr, iur ! dart index at four corners
real(r8), dimension(ens_size) :: x_ill, x_iul, x_ilr, x_iur ! state value at four corners
integer :: var_id


var_id = get_varid_from_kind(wrf_dom(id), qty)

do e = 1, ens_size
! x, y, z, domain, variable
ill(e) = get_dart_vector_index(ll(1), ll(2), k(e), wrf_dom(id), var_id)
iul(e) = get_dart_vector_index(ul(1), ul(2), k(e), wrf_dom(id), var_id)
ilr(e) = get_dart_vector_index(lr(1), lr(2), k(e), wrf_dom(id), var_id)
iur(e) = get_dart_vector_index(ur(1), ur(2), k(e), wrf_dom(id), var_id)

enddo

call get_state_array(x_ill, ill, state_handle)
call get_state_array(x_iul, iul, state_handle)
call get_state_array(x_ilr, ilr, state_handle)
call get_state_array(x_iur, iur, state_handle)

simple_interpolation = dym*( dxm*x_ill(:) + dx*x_ilr(:) ) + dy*( dxm*x_iul(:) + dx*x_iur(:) )

end function simple_interpolation

!------------------------------------------------------------------
function vertical_interpolation(ens_size, k, zloc, fld1, fld2)

integer, intent(in) :: ens_size
integer, intent(in) :: k(ens_size)
real(r8), intent(in) :: zloc(ens_size)
real(r8), intent(in) :: fld1(ens_size), fld2(ens_size)

real(r8) :: vertical_interpolation(ens_size)

real(r8) :: dz, dzm
integer :: z ! level
integer :: e

do e = 1, ens_size
call toGrid(zloc(e), z, dz, dzm)
! comment from original code:
! If you get here and zloc < 1.0, then z will be 0, and
! we should extrapolate. fld(1,:) and fld(2,:) where computed
! at levels 1 and 2.

if (z >= 1) then
! Linearly interpolate between levels
vertical_interpolation(e) = dzm*fld1(e) + dz*fld2(e)
else
! Extrapolate below first level.
vertical_interpolation(e) = fld1(e) - (fld2(e)-fld1(e))*dzm
endif
enddo


end function vertical_interpolation
!------------------------------------------------------------------
! Returns the smallest increment in time that the model is capable
! of advancing the state in a given implementation, or the shortest
Expand All @@ -606,8 +458,6 @@ function shortest_time_between_assimilations()

end function shortest_time_between_assimilations



!------------------------------------------------------------------
! Given an integer index into the state vector, returns the
! associated location and optionally the physical quantity.
Expand Down Expand Up @@ -808,7 +658,156 @@ subroutine end_model()
end subroutine end_model

!------------------------------------------------------------------
!----------------------------------------------------------------------
! end of public routines
!------------------------------------------------------------------
function force_non_negative_if_required(n, qty, fld)

integer, intent(in) :: n
integer, intent(in) :: qty
real(r8), intent(in) :: fld(n)

real(r8) :: force_non_negative_if_required(n)

select case (qty)
case (QTY_RAINWATER_MIXING_RATIO, &
QTY_GRAUPEL_MIXING_RATIO, &
QTY_HAIL_MIXING_RATIO, &
QTY_SNOW_MIXING_RATIO, &
QTY_ICE_MIXING_RATIO, &
QTY_CLOUDWATER_MIXING_RATIO, &
QTY_DROPLET_NUMBER_CONCENTR, &
QTY_ICE_NUMBER_CONCENTRATION, &
QTY_SNOW_NUMBER_CONCENTR, &
QTY_RAIN_NUMBER_CONCENTR, &
QTY_GRAUPEL_NUMBER_CONCENTR, &
QTY_HAIL_NUMBER_CONCENTR )
force_non_negative_if_required = max(0.0_r8, fld)
case default
force_non_negative_if_required = fld
end select

end function force_non_negative_if_required

!------------------------------------------------------------------
function surface_qty(qty)

integer, intent(in) :: qty
logical :: surface_qty

select case (qty)

case (QTY_2M_TEMPERATURE, &
QTY_2M_SPECIFIC_HUMIDITY, &
QTY_10M_U_WIND_COMPONENT, &
QTY_10M_V_WIND_COMPONENT, &
QTY_SURFACE_TYPE, &
QTY_SKIN_TEMPERATURE)
surface_qty = .true.
case default
surface_qty = .false.

end select

end function surface_qty

!------------------------------------------------------------------
! WRF has separate model qtys for surface variables
function update_qty_if_location_is_surface(qty_in, location) result(qty)

integer, intent(in) :: qty_in
type(location_type), intent(in) :: location
integer :: qty

if (.not. is_vertical(location,"SURFACE")) then
qty = qty_in
return
endif

select case (qty)

case (QTY_U_WIND_COMPONENT); qty = QTY_10M_U_WIND_COMPONENT
case (QTY_V_WIND_COMPONENT); qty = QTY_10M_V_WIND_COMPONENT
case (QTY_POTENTIAL_TEMPERATURE); qty = QTY_2M_TEMPERATURE
case (QTY_SPECIFIC_HUMIDITY); qty = QTY_2M_SPECIFIC_HUMIDITY
case (QTY_VAPOR_MIXING_RATIO); qty = QTY_2M_SPECIFIC_HUMIDITY ! Vapor Mixing Ratio (QV, Q2)

end select

end function update_qty_if_location_is_surface

!------------------------------------------------------------------
function simple_interpolation(ens_size, state_handle, qty, id, ll, ul, lr, ur, k, dxm, dx, dy, dym)

integer, intent(in) :: ens_size
type(ensemble_type), intent(in) :: state_handle
integer, intent(in) :: qty
integer, intent(in) :: id
integer, intent(in) :: ll(2), ul(2), lr(2), ur(2) ! (x,y) at four corners
integer, intent(in) :: k(ens_size) ! k may be different across the ensemble
real(r8), intent(in) :: dxm, dx, dy, dym

real(r8) :: simple_interpolation(ens_size)

integer :: e ! loop variable
! lower left, upper left, lower right, upper right
integer(i8), dimension(ens_size) :: ill, iul, ilr, iur ! dart index at four corners
real(r8), dimension(ens_size) :: x_ill, x_iul, x_ilr, x_iur ! state value at four corners
integer :: var_id


var_id = get_varid_from_kind(wrf_dom(id), qty)

do e = 1, ens_size
! x, y, z, domain, variable
ill(e) = get_dart_vector_index(ll(1), ll(2), k(e), wrf_dom(id), var_id)
iul(e) = get_dart_vector_index(ul(1), ul(2), k(e), wrf_dom(id), var_id)
ilr(e) = get_dart_vector_index(lr(1), lr(2), k(e), wrf_dom(id), var_id)
iur(e) = get_dart_vector_index(ur(1), ur(2), k(e), wrf_dom(id), var_id)

enddo

call get_state_array(x_ill, ill, state_handle)
call get_state_array(x_iul, iul, state_handle)
call get_state_array(x_ilr, ilr, state_handle)
call get_state_array(x_iur, iur, state_handle)

simple_interpolation = dym*( dxm*x_ill(:) + dx*x_ilr(:) ) + dy*( dxm*x_iul(:) + dx*x_iur(:) )

end function simple_interpolation

!------------------------------------------------------------------
function vertical_interpolation(ens_size, k, zloc, fld1, fld2)

integer, intent(in) :: ens_size
integer, intent(in) :: k(ens_size)
real(r8), intent(in) :: zloc(ens_size)
real(r8), intent(in) :: fld1(ens_size), fld2(ens_size)

real(r8) :: vertical_interpolation(ens_size)

real(r8) :: dz, dzm
integer :: z ! level
integer :: e

do e = 1, ens_size
call toGrid(zloc(e), z, dz, dzm)
! comment from original code:
! If you get here and zloc < 1.0, then z will be 0, and
! we should extrapolate. fld(1,:) and fld(2,:) where computed
! at levels 1 and 2.

if (z >= 1) then
! Linearly interpolate between levels
vertical_interpolation(e) = dzm*fld1(e) + dz*fld2(e)
else
! Extrapolate below first level.
vertical_interpolation(e) = fld1(e) - (fld2(e)-fld1(e))*dzm
endif
enddo

end function vertical_interpolation

!------------------------------------------------------------------
subroutine get_wrf_date(tstring, year, month, day, hour, minute, second)

character(len=19), intent(in) :: tstring ! YYYY-MM-DD_hh:mm:ss
Expand Down Expand Up @@ -1191,7 +1190,7 @@ subroutine get_model_pressure_profile(id, ll, ul, lr, ur, dx, dy, dxm, dym, &
v_p(0,:) = extrap_4pressure(lev2_pres1(:), lev2_pres2(:), lev2_pres3(:), lev2_pres4(:), dx, dxm, dy, dym, ens_size, &
edgep=v_p(1,:))

endif
endif

end subroutine get_model_pressure_profile

Expand Down Expand Up @@ -3033,4 +3032,3 @@ subroutine vortex()
end subroutine vortex

end module model_mod

0 comments on commit 216a3f4

Please sign in to comment.