diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F index a146cb3d00..5757e16baa 100644 --- a/src/core_init_atmosphere/mpas_init_atm_cases.F +++ b/src/core_init_atmosphere/mpas_init_atm_cases.F @@ -16,6 +16,7 @@ module init_atm_cases use atm_advection use mpas_RBF_interpolation use mpas_vector_reconstruction + use init_atm_vinterp, only: vertical_interp use mpas_timer use mpas_init_atm_static use mpas_init_atm_surface @@ -7595,101 +7596,6 @@ 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, ierr) - - implicit none - - real (kind=RKIND), intent(in) :: target_z - 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 ! 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 - - integer :: interp_order, extrap_type - real (kind=RKIND) :: surface, sealevel - - if (present(ierr)) ierr = 0 - - if (present(order)) then - interp_order = order - else - interp_order = 2 - end if - - if (present(extrap)) then - extrap_type = extrap - else - extrap_type = 1 - end if - - if (present(surface_val)) then - surface = surface_val - else - surface = 200100.0 - end if - - if (present(sealev_val)) then - sealevel = sealev_val - else - sealevel = 201300.0 - end if - - ! - ! Extrapolation required - ! - if (target_z < zf(1,1)) then - if (extrap_type == 0) then - vertical_interp = zf(2,1) - else if (extrap_type == 1) then - slope = (zf(2,2) - zf(2,1)) / (zf(1,2) - zf(1,1)) - vertical_interp = zf(2,1) + slope * (target_z - zf(1,1)) - else if (extrap_type == 2) then - vertical_interp = zf(2,1) - (target_z - zf(1,1))*0.0065 - end if - return - end if - if (target_z >= zf(1,nz)) then - if (extrap_type == 0) then - vertical_interp = zf(2,nz) - 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 - call mpas_log_write('extrap_type == 2 not implemented for target_z >= zf(1,nz)', messageType=MPAS_LOG_ERR) - if (present(ierr)) ierr = 1 - return - end if - return - end if - - - ! - ! No extrapolation required - ! - do k=1,nz-1 - if (target_z >= zf(1,k) .and. target_z < zf(1,k+1)) then - lm = k - lp = k+1 - wm = (zf(1,k+1) - target_z) / (zf(1,k+1) - zf(1,k)) - wp = (target_z - zf(1,k)) / (zf(1,k+1) - zf(1,k)) - exit - end if - end do - - vertical_interp = wm*zf(2,lm) + wp*zf(2,lp) - - return - - end function vertical_interp - - !---------------------------------------------------------------------------------------------------------- real (kind=RKIND) function env_qv( temperature, pressure, rh_max ) diff --git a/src/core_init_atmosphere/mpas_init_atm_vinterp.F b/src/core_init_atmosphere/mpas_init_atm_vinterp.F index c5163d2708..164cfac546 100644 --- a/src/core_init_atmosphere/mpas_init_atm_vinterp.F +++ b/src/core_init_atmosphere/mpas_init_atm_vinterp.F @@ -22,7 +22,10 @@ module init_atm_vinterp ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - 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) + + use mpas_log, only : mpas_log_write + use mpas_derived_types, only : MPAS_LOG_ERR implicit none @@ -30,10 +33,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 @@ -41,7 +45,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 @@ -75,6 +80,8 @@ real (kind=RKIND) function vertical_interp(target_z, nz, zf, order, extrap, surf else if (extrap_type == 1) then slope = (zf(2,2) - zf(2,1)) / (zf(1,2) - zf(1,1)) vertical_interp = zf(2,1) + slope * (target_z - zf(1,1)) + else if (extrap_type == 2) then + vertical_interp = zf(2,1) - (target_z - zf(1,1))*0.0065 end if return end if @@ -84,6 +91,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 + call mpas_log_write('extrap_type == 2 not implemented for target_z >= zf(1,nz)', messageType=MPAS_LOG_ERR) + if (present(ierr)) ierr = 1 + return end if return end if