!------------------------------------------------------------------------------------------- !Copyright (c) 2013-2016 by Wolfgang Kurtz and Guowei He (Forschungszentrum Juelich GmbH) ! !This file is part of TSMP-PDAF ! !TSMP-PDAF is free software: you can redistribute it and/or modify !it under the terms of the GNU Lesser General Public License as published by !the Free Software Foundation, either version 3 of the License, or !(at your option) any later version. ! !TSMP-PDAF is distributed in the hope that it will be useful, !but WITHOUT ANY WARRANTY; without even the implied warranty of !MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !GNU LesserGeneral Public License for more details. ! !You should have received a copy of the GNU Lesser General Public License !along with TSMP-PDAF. If not, see <http://www.gnu.org/licenses/>. !------------------------------------------------------------------------------------------- ! ! !------------------------------------------------------------------------------------------- !obs_op_pdaf.F90: TSMP-PDAF implementation of routine ! 'obs_op_pdaf' (PDAF online coupling) !------------------------------------------------------------------------------------------- !$Id: obs_op_pdaf.F90 1441 2013-10-04 10:33:42Z lnerger $ !BOP ! ! !ROUTINE: obs_op_pdaf --- Implementation of observation operator ! ! !INTERFACE: SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p) ! !DESCRIPTION: ! User-supplied routine for PDAF. ! Used in the filters: SEIK/EnKF/ETKF/ESTKF ! ! The routine is called during the analysis step. ! It has to perform the operation of the ! observation operator acting on a state vector. ! For domain decomposition, the action is on the ! PE-local sub-domain of the state and has to ! provide the observed sub-state for the PE-local ! domain. ! ! !REVISION HISTORY: ! 2013-02 - Lars Nerger - Initial code ! Later revisions - see svn log ! ! !USES: USE mod_assimilation, & ONLY: obs_index_p, obs_p, & #ifndef CLMSA #ifndef OBS_ONLY_CLM sc_p, & #endif #endif obs_interp_indices_p, & obs_interp_weights_p use mod_tsmp, & only: obs_interp_switch, & soilay, & soilay_fortran, & nz_glob, & da_crns_depth_tol, & crns_flag ! tcycle USE, INTRINSIC :: iso_c_binding #if defined CLMSA USE enkf_clm_mod, & ONLY : clm_varsize, clm_paramarr, clmupdate_swc, clmupdate_T, clmcrns_bd #ifdef CLMFIVE USE clm_instMod, & ONLY : soilstate_inst #endif #endif IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: step ! Currrent time step INTEGER, INTENT(in) :: dim_p ! PE-local dimension of state INTEGER, INTENT(in) :: dim_obs_p ! Dimension of observed state REAL, INTENT(in) :: state_p(dim_p) ! PE-local model state REAL, INTENT(out) :: m_state_p(dim_obs_p) ! PE-local observed state integer :: i, j, k, z, n integer :: icorner logical :: lpointobs !If true: no special observation; use point observation ! !CALLING SEQUENCE: ! Called by: PDAF_seik_analysis, PDAF_seik_analysis_newT (as U_obs_op) ! Called by: PDAF_enkf_analysis_rlm, PDAF_enkf_analysis_rsm !EOP ! *** local variables *** ! hcp test with hardcoding variable declaration real(8), dimension(:), allocatable :: soide !soil depth !real(8), dimension(0:12), parameter :: & ! soide=(/0.d0, 0.02d0, 0.05d0, 0.1d0, 0.17d0, 0.3d0, 0.5d0, & ! 0.8d0, 1.3d0, 2.d0, 3.d0, 5.d0, 12.d0/) !soil depth real(8) :: tot, avesm, avesm_temp, Dp integer :: nsc ! end of hcp #ifndef PARFLOW_STAND_ALONE #ifndef OBS_ONLY_PARFLOW ! Variables used in crns version 2 REAL :: weights_r1(920), weights_r2(920), weights_r3(920) Real :: weights_layer(8) Integer :: nweights(8) Real :: d86_r1, d86_r2, d86_r3 REAL :: r1 = 1.0 REAL :: r2 = 20.0 REAL :: r3 = 85.0 REAL :: bd, y REAL :: sum_r1, sum_r2, sum_r3, totw #endif #endif ! ********************************************* ! *** Perform application of measurement *** ! *** operator H on vector or matrix column *** ! ********************************************* ! If no special observation operator is compiled, use point observations lpointobs = .true. #if defined CLMSA if (clmupdate_T.EQ.1) then lpointobs = .false. DO i = 1, dim_obs_p ! Equation for LST computation from Kustas2009, Eq(7) ! http://dx.doi.org/10.1016/j.agrformet.2009.05.016 ! ! Comment: Fractional vegetation cover (Eq(8) from Kustas2009) ! currently implemented with simplified settings: Vegetation ! clumping parameter `Omega=1`; radiometer view angle `phi=0` m_state_p(i) & = (exp(-0.5*clm_paramarr(obs_index_p(i))) & *state_p(obs_index_p(i))**4 & +(1.-exp(-0.5*clm_paramarr(obs_index_p(i)))) & *state_p(clm_varsize+obs_index_p(i))**4)**0.25 END DO ! write(*,*) 'now is cycle ', tcycle ! write(*,*) 'obs_p =', obs_p(:) ! write(*,*) 'model LST', m_state_p(:) ! write(*,*) 'TG', state_p(obs_index_p(:)) ! write(*,*) 'TV', state_p(clm_varsize+obs_index_p(:)) endif #endif #ifndef CLMSA #ifndef OBS_ONLY_CLM if (crns_flag.EQ.1) then !Schroen et al HESS 2017 modelled CRNS averaging lpointobs = .false. call C_F_POINTER(soilay,soilay_fortran,[nz_glob]) Allocate(soide(0:nz_glob)) soide(0)=0.d0 do i=1,nz_glob soide(i)=soide(i-1)+soilay_fortran(nz_glob-i+1) enddo do i = 1, dim_obs_p !Initial average soil moisture for 1st iteration avesm=0.d0 do j=1,nz_glob avesm=avesm+(soide(j)-soide(j-1))*state_p(sc_p(j,i))/soide(nz_glob) enddo avesm_temp=0.d0 !iteration do while (abs(avesm-avesm_temp)/avesm .GE. da_crns_depth_tol) !Averaging, conventional profile, Schroen et al HESS 2017 Eq. (3) avesm_temp=avesm Dp=0.058d0/(avesm+0.0829d0) !Sum weight*soil_moisture avesm=0.d0; nsc=nz_glob do j=1,nz_glob if ((soide(j-1).LT.Dp).AND.(Dp.LE.soide(j))) then nsc=j endif enddo do j=1, nsc-1 avesm=avesm+(1.d0-0.5d0*(soide(j)+soide(j-1))/Dp)*(soide(j)-soide(j-1)) & *state_p(sc_p(j,i))/Dp enddo avesm=avesm+(1.d0-0.5d0*(Dp+soide(nsc-1))/Dp)*(Dp-soide(nsc-1)) & *state_p(sc_p(nsc,i))/Dp !Sum weight tot=0.d0 do j=1, nsc-1 tot = tot+(1.d0-0.5d0*(soide(j)+soide(j-1))/Dp)*(soide(j)-soide(j-1)) & /Dp enddo tot = tot+(1.d0-0.5d0*(Dp+soide(nsc-1))/Dp)*(Dp-soide(nsc-1)) & /Dp avesm=avesm/tot enddo m_state_p(i)=avesm enddo deallocate(soide) end if #endif #endif #ifndef PARFLOW_STAND_ALONE #ifndef OBS_ONLY_PARFLOW #ifdef CLMFIVE if (crns_flag.EQ.2) then lpointobs = .false. ! CRNS implementation based on Schrön et al. 2017 using ! d86 for 3 different radius values DO i = 1, dim_obs_p ! Bulk density average for the 8 considered layers IF (clmcrns_bd > 0.0) THEN bd = clmcrns_bd ELSE bd = 0.0 DO j = 1, 8 bd = bd + soilstate_inst%bd_col(obs_index_p(i),j) ! bulk density END DO bd = bd / 8.0 * 0.001 ! average and convert from kg/m^3 to g/cm^3 ENDIF ! CRNS observed value y = obs_p(i) ! CRNS observation ! Penetration depth calculations D86(bd, r, y) d86_r1 = (1/bd*(8.321+0.14249*(0.96655+exp(-0.01*r1))*(20+y)/(0.0429+y))) d86_r2 = (1/bd*(8.321+0.14249*(0.96655+exp(-0.01*r2))*(20+y)/(0.0429+y))) d86_r3 = (1/bd*(8.321+0.14249*(0.96655+exp(-0.01*r3))*(20+y)/(0.0429+y))) ! Then calculate the weights for thin (1mm) slices of the layers for 85cm sum_r1 = 0.0 sum_r2 = 0.0 sum_r3 = 0.0 DO j = 1, 920 ! depth in mm but in calculation used in cm: weights_r1(j) = exp(-2*(j/10)/d86_r1) weights_r2(j) = exp(-2*(j/10)/d86_r2) weights_r3(j) = exp(-2*(j/10)/d86_r3) sum_r1 = sum_r1 + weights_r1(j) sum_r2 = sum_r2 + weights_r2(j) sum_r3 = sum_r3 + weights_r3(j) END DO ! Normalize the weights: weights_r1(:) = weights_r1(:) / sum_r1 weights_r2(:) = weights_r2(:) / sum_r2 weights_r3(:) = weights_r3(:) / sum_r3 ! assign average weights to each layer nweights(:) = 0 weights_layer(:) = 0.0 ! z index for different layers, manually here, could be done better with model layer depth DO j = 1, 920 IF (j > 680) then z = 8 ELSEIF (j < 680 .and. j > 480) then z = 7 ELSEIF (j < 480 .and. j > 320) then z = 6 ELSEIF (j < 320 .and. j > 200) then z = 5 ELSEIF (j < 200 .and. j > 120) then z = 4 ELSEIF (j < 120 .and. j > 60) then z = 3 ELSEIF (j < 60 .and. j > 20) then z = 2 ELSEIF (j < 20) then z = 1 ENDIF weights_layer(z) = weights_layer(z) + weights_r1(j) + weights_r2(j) + weights_r3(j) nweights(z) = nweights(z) + 1 END DO ! Normalize the weights totw = 0.0 DO j = 1, 8 weights_layer(j) = weights_layer(j) / (nweights(j) * 3) totw = totw + weights_layer(j) END DO weights_layer(:) = weights_layer(:) / totw ! Finally use the weights to calculate the weighted average of the state variable avesm = 0.0 DO j = 1, 8 avesm = avesm + weights_layer(j) * state_p(obs_index_p(i) + (j-1)) ! This assumes that obs_index_p(i) for obs i is the index of ! the first layer of the gridcell where obs i is END DO ! Assign new average as the state variable m_state_p(i) = avesm ! end loop over observations END DO end if #endif #endif #endif if(obs_interp_switch == 1) then lpointobs = .false. do i = 1, dim_obs_p m_state_p(i) = 0 do icorner = 1, 4 m_state_p(i) = m_state_p(i) + state_p(obs_interp_indices_p(i,icorner)) * obs_interp_weights_p(i,icorner) enddo enddo end if if(lpointobs) then DO i = 1, dim_obs_p m_state_p(i) = state_p(obs_index_p(i)) END DO end if END SUBROUTINE obs_op_pdaf