!------------------------------------------------------------------------------- ! Description: ! ! Interpolate all model fields to observation positions. ! ! Copyright: ! ! This software was developed within the context of the EUMETSAT Satellite ! Application Facility on Numerical Weather Prediction (NWP SAF), under the ! Cooperation Agreement dated 29 June 2011, between EUMETSAT and the ! Met Office, UK, by one or more partners within the NWP SAF. The partners ! in the NWP SAF are the Met Office, ECMWF, KNMI and MeteoFrance. ! ! Copyright 2017, EUMETSAT, All Rights Reserved. ! !------------------------------------------------------------------------------- subroutine radsim_interp(obs, model_in, model_out, keep_input) use radsim_mod_constants, only : & real64, & status_error use radsim_mod_types, only : & model_type, & obs_type implicit none include 'radsim_error_report.interface' include 'radsim_grid_init.interface' include 'radsim_interp_horiz.interface' include 'radsim_interp_index.interface' type(obs_type), intent(in) :: obs ! Observation data type(model_type), intent(inout) :: model_in ! Input model data type(model_type), intent(inout) :: model_out ! Output model data logical, optional, intent(in) :: keep_input ! Keep input model data integer :: i, j integer :: nobs, nlevels integer :: arrayshape(2) real(real64), allocatable :: input_array(:,:) real(real64), pointer :: output_array1d(:) real(real64), pointer :: output_array2d(:,:) logical :: keep integer, allocatable :: x_index(:) integer, allocatable :: y_index(:) !------------------------------------------------------------------------------- nobs = size(obs%lon) if ( present(keep_input) ) then keep = keep_input else keep = .true. end if ! Create indices locating position of observation points in the model arrays call radsim_interp_index( & model_in % grid, & obs % lon, & obs % lat, & x_index, & y_index, & mask = obs % qcflags /= 0 ) ! Allocate space for nlevels+1 because half-level pressures may have an ! additional value. allocate(input_array(model_in%nprofs, model_in%nlevels+1)) do i = 1, 50 nlevels = 0 select case(i) !-------------------- ! Single-level fields !-------------------- case(1) if ( associated(model_in % pstar) ) then nlevels = 1 input_array(:,1) = model_in % pstar if (.not. keep) deallocate(model_in % pstar) if ( associated(model_out % pstar) .and. model_out % nprofs /= nobs ) then deallocate(model_out % pstar) end if if ( .not. associated(model_out % pstar) ) allocate(model_out % pstar(nobs)) output_array1d => model_out % pstar end if case(2) if ( associated(model_in % t2) ) then nlevels = 1 input_array(:,1) = model_in % t2 if (.not. keep) deallocate(model_in % t2) if ( associated(model_out % t2) .and. model_out % nprofs /= nobs ) then deallocate(model_out % t2) end if if ( .not. associated(model_out % t2) ) allocate(model_out % t2(nobs)) output_array1d => model_out % t2 end if case(3) if ( associated(model_in % td2) ) then nlevels = 1 input_array(:,1) = model_in % td2 if (.not. keep) deallocate(model_in % td2) if ( associated(model_out % td2) .and. model_out % nprofs /= nobs ) then deallocate(model_out % td2) end if if ( .not. associated(model_out % td2) ) allocate(model_out % td2(nobs)) output_array1d => model_out % td2 end if case(4) if ( associated(model_in % q2) ) then nlevels = 1 input_array(:,1) = model_in % q2 if (.not. keep) deallocate(model_in % q2) if ( associated(model_out % q2) .and. model_out % nprofs /= nobs ) then deallocate(model_out % q2) end if if ( .not. associated(model_out % q2) ) allocate(model_out % q2(nobs)) output_array1d => model_out % q2 end if case(5) if ( associated(model_in % rh2) ) then nlevels = 1 input_array(:,1) = model_in % rh2 if (.not. keep) deallocate(model_in % rh2) if ( associated(model_out % rh2) .and. model_out % nprofs /= nobs ) then deallocate(model_out % rh2) end if if ( .not. associated(model_out % rh2) ) allocate(model_out % rh2(nobs)) output_array1d => model_out % rh2 end if case(6) if ( associated(model_in % tskin) ) then nlevels = 1 input_array(:,1) = model_in % tskin if (.not. keep) deallocate(model_in % tskin) if ( associated(model_out % tskin) .and. model_out % nprofs /= nobs ) then deallocate(model_out % tskin) end if if ( .not. associated(model_out % tskin) ) allocate(model_out % tskin(nobs)) output_array1d => model_out % tskin end if case(7) if ( associated(model_in % u10) ) then nlevels = 1 input_array(:,1) = model_in % u10 if (.not. keep) deallocate(model_in % u10) if ( associated(model_out % u10) .and. model_out % nprofs /= nobs ) then deallocate(model_out % u10) end if if ( .not. associated(model_out % u10) ) allocate(model_out % u10(nobs)) output_array1d => model_out % u10 end if case(8) if ( associated(model_in % v10) ) then nlevels = 1 input_array(:,1) = model_in % v10 if (.not. keep) deallocate(model_in % v10) if ( associated(model_out % v10) .and. model_out % nprofs /= nobs ) then deallocate(model_out % v10) end if if ( .not. associated(model_out % v10) ) allocate(model_out % v10(nobs)) output_array1d => model_out % v10 end if case(9) if ( associated(model_in % seaice) ) then nlevels = 1 input_array(:,1) = model_in % seaice where(input_array(:,1) < 0.0 .or. input_array(:,1) > 1.0) input_array(:,1) = 0.0 ! remove missing data if (.not. keep) deallocate(model_in % seaice) if ( associated(model_out % seaice) .and. model_out % nprofs /= nobs ) then deallocate(model_out % seaice) end if if ( .not. associated(model_out % seaice) ) allocate(model_out % seaice(nobs)) output_array1d => model_out % seaice end if case(10) if ( associated(model_in % zsurf) ) then nlevels = 1 input_array(:,1) = model_in % zsurf where(input_array(:,1) < 0.0) input_array(:,1) = 0.0 ! remove missing data if (.not. keep) deallocate(model_in % zsurf) if ( associated(model_out % zsurf) .and. model_out % nprofs /= nobs ) then deallocate(model_out % zsurf) end if if ( .not. associated(model_out % zsurf) ) allocate(model_out % zsurf(nobs)) output_array1d => model_out % zsurf end if case(11) if ( associated(model_in % lsm) ) then nlevels = 1 input_array(:,1) = model_in % lsm if (.not. keep) deallocate(model_in % lsm) if ( associated(model_out % lsm) .and. model_out % nprofs /= nobs ) then deallocate(model_out % lsm) end if if ( .not. associated(model_out % lsm) ) allocate(model_out % lsm(nobs)) output_array1d => model_out % lsm end if !------------------- ! Multi-level fields !------------------- case(21) if ( associated(model_in % p) ) then nlevels = size(model_in % p, dim=2) input_array(:,1:nlevels) = model_in % p if (.not. keep) deallocate(model_in % p) if ( associated(model_out % p) .and. model_out % nprofs /= nobs ) then deallocate(model_out % p) end if if ( .not. associated(model_out % p) ) allocate(model_out % p(nobs,nlevels)) output_array2d => model_out % p end if case(22) if ( associated(model_in % ph) ) then nlevels = size(model_in % ph, dim=2) input_array(:,1:nlevels) = model_in % ph if (.not. keep) deallocate(model_in % ph) if ( associated(model_out % ph) .and. model_out % nprofs /= nobs ) then deallocate(model_out % ph) end if if ( .not. associated(model_out % ph) ) allocate(model_out % ph(nobs,nlevels)) output_array2d => model_out % ph end if case(23) if ( associated(model_in % t) ) then nlevels = size(model_in % t, dim=2) input_array(:,1:nlevels) = model_in % t if (.not. keep) deallocate(model_in % t) if ( associated(model_out % t) .and. model_out % nprofs /= nobs ) then deallocate(model_out % t) end if if ( .not. associated(model_out % t) ) allocate(model_out % t(nobs,nlevels)) output_array2d => model_out % t end if case(24) if ( associated(model_in % theta) ) then nlevels = size(model_in % theta, dim=2) input_array(:,1:nlevels) = model_in % theta if (.not. keep) deallocate(model_in % theta) if ( associated(model_out % theta) .and. model_out % nprofs /= nobs ) then deallocate(model_out % theta) end if if ( .not. associated(model_out % theta) ) allocate(model_out % theta(nobs,nlevels)) output_array2d => model_out % theta end if case(25) if ( associated(model_in % q) ) then nlevels = size(model_in % q, dim=2) input_array(:,1:nlevels) = model_in % q if (.not. keep) deallocate(model_in % q) if ( associated(model_out % q) .and. model_out % nprofs /= nobs ) then deallocate(model_out % q) end if if ( .not. associated(model_out % q) ) allocate(model_out % q(nobs,nlevels)) output_array2d => model_out % q end if case(26) if ( associated(model_in % rh) ) then nlevels = size(model_in % rh, dim=2) input_array(:,1:nlevels) = model_in % rh if (.not. keep) deallocate(model_in % rh) if ( associated(model_out % rh) .and. model_out % nprofs /= nobs ) then deallocate(model_out % rh) end if if ( .not. associated(model_out % rh) ) allocate(model_out % rh(nobs,nlevels)) output_array2d => model_out % rh end if case(27) if ( associated(model_in % clw) ) then nlevels = size(model_in % clw, dim=2) input_array(:,1:nlevels) = model_in % clw if (.not. keep) deallocate(model_in % clw) if ( associated(model_out % clw) .and. model_out % nprofs /= nobs ) then deallocate(model_out % clw) end if if ( .not. associated(model_out % clw) ) allocate(model_out % clw(nobs,nlevels)) output_array2d => model_out % clw end if case(28) if ( associated(model_in % ciw) ) then nlevels = size(model_in % ciw, dim=2) input_array(:,1:nlevels) = model_in % ciw if (.not. keep) deallocate(model_in % ciw) if ( associated(model_out % ciw) .and. model_out % nprofs /= nobs ) then deallocate(model_out % ciw) end if if ( .not. associated(model_out % ciw) ) allocate(model_out % ciw(nobs,nlevels)) output_array2d => model_out % ciw end if case(29) if ( associated(model_in % rain) ) then nlevels = size(model_in % rain, dim=2) input_array(:,1:nlevels) = model_in % rain if (.not. keep) deallocate(model_in % rain) if ( associated(model_out % rain) .and. model_out % nprofs /= nobs ) then deallocate(model_out % rain) end if if ( .not. associated(model_out % rain) ) allocate(model_out % rain(nobs,nlevels)) output_array2d => model_out % rain end if case(30) if ( associated(model_in % cfrac) ) then nlevels = size(model_in % cfrac, dim=2) input_array(:,1:nlevels) = model_in % cfrac if (.not. keep) deallocate(model_in % cfrac) if ( associated(model_out % cfrac) .and. model_out % nprofs /= nobs ) then deallocate(model_out % cfrac) end if if ( .not. associated(model_out % cfrac) ) allocate(model_out % cfrac(nobs,nlevels)) output_array2d => model_out % cfrac end if case(31) if ( associated(model_in % o3) ) then nlevels = size(model_in % o3, dim=2) input_array(:,1:nlevels) = model_in % o3 if (.not. keep) deallocate(model_in % o3) if ( associated(model_out % o3) .and. model_out % nprofs /= nobs ) then deallocate(model_out % o3) end if if ( .not. associated(model_out % o3) ) allocate(model_out % o3(nobs,nlevels)) output_array2d => model_out % o3 end if case(32) if ( associated(model_in % co2) ) then nlevels = size(model_in % co2, dim=2) input_array(:,1:nlevels) = model_in % co2 if (.not. keep) deallocate(model_in % co2) if ( associated(model_out % co2) .and. model_out % nprofs /= nobs ) then deallocate(model_out % co2) end if if ( .not. associated(model_out % co2) ) allocate(model_out % co2(nobs,nlevels)) output_array2d => model_out % co2 end if case(33) if ( associated(model_in % n2o) ) then nlevels = size(model_in % n2o, dim=2) input_array(:,1:nlevels) = model_in % n2o if (.not. keep) deallocate(model_in % n2o) if ( associated(model_out % n2o) .and. model_out % nprofs /= nobs ) then deallocate(model_out % n2o) end if if ( .not. associated(model_out % n2o) ) allocate(model_out % n2o(nobs,nlevels)) output_array2d => model_out % n2o end if case(34) if ( associated(model_in % co) ) then nlevels = size(model_in % co, dim=2) input_array(:,1:nlevels) = model_in % co if (.not. keep) deallocate(model_in % co) if ( associated(model_out % co) .and. model_out % nprofs /= nobs ) then deallocate(model_out % co) end if if ( .not. associated(model_out % co) ) allocate(model_out % co(nobs,nlevels)) output_array2d => model_out % co end if case(35) if ( associated(model_in % ch4) ) then nlevels = size(model_in % ch4, dim=2) input_array(:,1:nlevels) = model_in % ch4 if (.not. keep) deallocate(model_in % ch4) if ( associated(model_out % ch4) .and. model_out % nprofs /= nobs ) then deallocate(model_out % ch4) end if if ( .not. associated(model_out % ch4) ) allocate(model_out % ch4(nobs,nlevels)) output_array2d => model_out % ch4 end if case(36) if ( associated(model_in % z) ) then nlevels = size(model_in % z, dim=2) input_array(:,1:nlevels) = model_in % z if (.not. keep) deallocate(model_in % z) if ( associated(model_out % z) .and. model_out % nprofs /= nobs ) then deallocate(model_out % z) end if if ( .not. associated(model_out % z) ) allocate(model_out % z(nobs,nlevels)) output_array2d => model_out % z end if case(37) if ( associated(model_in % cfrac_liq) ) then nlevels = size(model_in % cfrac_liq, dim=2) input_array(:,1:nlevels) = model_in % cfrac_liq if (.not. keep) deallocate(model_in % cfrac_liq) if ( associated(model_out % cfrac_liq) .and. model_out % nprofs /= nobs ) then deallocate(model_out % cfrac_liq) end if if ( .not. associated(model_out % cfrac_liq) ) allocate(model_out % cfrac_liq(nobs,nlevels)) output_array2d => model_out % cfrac_liq end if case(38) if ( associated(model_in % cfrac_ice) ) then nlevels = size(model_in % cfrac_ice, dim=2) input_array(:,1:nlevels) = model_in % cfrac_ice if (.not. keep) deallocate(model_in % cfrac_ice) if ( associated(model_out % cfrac_ice) .and. model_out % nprofs /= nobs ) then deallocate(model_out % cfrac_ice) end if if ( .not. associated(model_out % cfrac_ice) ) allocate(model_out % cfrac_ice(nobs,nlevels)) output_array2d => model_out % cfrac_ice end if case(39) if ( associated(model_in % cfrac_conv) ) then nlevels = size(model_in % cfrac_conv, dim=2) input_array(:,1:nlevels) = model_in % cfrac_conv if (.not. keep) deallocate(model_in % cfrac_conv) if ( associated(model_out % cfrac_conv) .and. model_out % nprofs /= nobs ) then deallocate(model_out % cfrac_conv) end if if ( .not. associated(model_out % cfrac_conv) ) allocate(model_out % cfrac_conv(nobs,nlevels)) output_array2d => model_out % cfrac_conv end if case(40) if ( associated(model_in % conv_cloud) ) then nlevels = size(model_in % conv_cloud, dim=2) input_array(:,1:nlevels) = model_in % conv_cloud if (.not. keep) deallocate(model_in % conv_cloud) if ( associated(model_out % conv_cloud) .and. model_out % nprofs /= nobs ) then deallocate(model_out % conv_cloud) end if if ( .not. associated(model_out % conv_cloud) ) allocate(model_out % conv_cloud(nobs,nlevels)) output_array2d => model_out % conv_cloud end if case(41) if ( associated(model_in % snow) ) then nlevels = size(model_in % snow, dim=2) input_array(:,1:nlevels) = model_in % snow if (.not. keep) deallocate(model_in % snow) if ( associated(model_out % snow) .and. model_out % nprofs /= nobs ) then deallocate(model_out % snow) end if if ( .not. associated(model_out % snow) ) allocate(model_out % snow(nobs,nlevels)) output_array2d => model_out % snow end if case default cycle end select if ( nlevels > 1 ) then do j = 1, nlevels call radsim_interp_horiz( & model_in % rmdi, & model_in % grid, & input_array(:,j), & obs % lon, & obs % lat, & x_index, & y_index, & output_array2d(:,j) ) end do else if ( nlevels == 1 ) then call radsim_interp_horiz( & model_in % rmdi, & model_in % grid, & input_array(:,1), & obs % lon, & obs % lat, & x_index, & y_index, & output_array1d ) end if end do deallocate(x_index) deallocate(y_index) ! Copy meta-data model_out % nprofs = nobs model_out % nlevels = model_in % nlevels model_out % rmdi = model_in % rmdi model_out % imdi = model_in % imdi if ( .not. associated(model_out % validity_time) ) then arrayshape = shape(model_in % validity_time) allocate(model_out % validity_time(arrayshape(1),arrayshape(2))) end if model_out % validity_time = model_in % validity_time ! Model data is now on the observation grid. Reset and assign the correct ! values. call radsim_grid_init(model_out % grid) allocate(model_out % grid % lat(nobs)) allocate(model_out % grid % lon(nobs)) model_out % grid % lat = obs % lat model_out % grid % lon = obs % lon end subroutine radsim_interp