LST-DA for eCLM-PDAF

interface/framework/init_dim_obs_f_pdaf.F90

@@ -129,6 +129,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
 #ifdef CLMFIVE
   use GridcellType, only: grc
   use ColumnType, only : col
+  use PatchType, only : patch
   ! use GetGlobalValuesMod, only: GetGlobalWrite
   ! use clm_varcon, only: nameg
   use enkf_clm_mod, only: state_clm2pdaf_p
@@ -145,6 +146,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
   USE enkf_clm_mod, only: domain_def_clm
   USE enkf_clm_mod, only: get_interp_idx
   use enkf_clm_mod, only: clmstatevec_allcol
+  use enkf_clm_mod, only: clmupdate_swc
+  use enkf_clm_mod, only: clmupdate_T
   !hcp end
 #endif
 #endif
@@ -164,9 +167,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
   integer :: ierror
   INTEGER :: max_var_id         ! Multi-scale DA
   INTEGER :: sum_dim_obs_p
+  INTEGER :: p                ! CLM Patch index
   INTEGER :: c                ! CLM Column index
   INTEGER :: g                ! CLM Gridcell index
   INTEGER :: cg
+  INTEGER :: pg
+  INTEGER :: pc
   INTEGER :: i,j,k        ! Counters
   INTEGER :: cnt          ! Counters
   INTEGER :: cnt_interp   ! Counter for interpolation grid cells
@@ -935,6 +941,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
      do i = 1, dim_obs
         obs(i) = clm_obs(i)
 
+       if(clmupdate_swc==1) then
+
        do g = begg,endg
          newgridcell = .true.
 
@@ -1004,9 +1012,101 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
              end if
 
            end if
-
          end do
        end do
+
+       else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then
+#ifdef CLMFIVE
+         ! patch loop
+         do g = begg,endg
+           newgridcell = .true.
+
+           do p = begp,endp
+
+             pg = patch%gridcell(p)
+             pc = patch%column(p)
+
+             if(pg == g) then
+               if(newgridcell) then
+                 ! Sets first patch/column in a gridcell. TODO: Make
+                 ! patch / column information part of the observation
+                 ! file
+
+                 if(is_use_dr) then
+                   deltax = abs(lon(g)-clmobs_lon(i))
+                   deltay = abs(lat(g)-clmobs_lat(i))
+                 end if
+
+                 if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. &
+                   ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then
+
+                   ! Set index in state vector, LST will be computed
+                   ! for first patch appearing here
+                   obs_index_p(cnt) = state_clm2pdaf_p(p,1)
+
+                   !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
+                   obs_p(cnt) = clm_obs(i)
+                   if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i)
+                   cnt = cnt + 1
+
+                 end if
+
+                 newgridcell = .false.
+
+               end if
+             end if
+
+           end do
+         end do
+#else
+         ! gridcell loop
+         do g = begg,endg
+
+           if(is_use_dr) then
+             deltax = abs(lon(g)-clmobs_lon(i))
+             deltay = abs(lat(g)-clmobs_lat(i))
+           end if
+
+           if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. &
+             ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then
+             obs_index_p(cnt) = g-begg+1
+           end if
+
+           !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
+           obs_p(cnt) = clm_obs(i)
+           if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i)
+           cnt = cnt + 1
+
+         end do
+#endif
+       else
+
+         print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p."
+         print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop."
+         ! call abort_parallel()
+
+         ! gridcell loop with layer information as default!
+         do g = begg,endg
+
+           if(is_use_dr) then
+             deltax = abs(lon(g)-clmobs_lon(i))
+             deltay = abs(lat(g)-clmobs_lat(i))
+           end if
+
+           if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. &
+             ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then
+             obs_index_p(cnt) = g-begg+1 +  ((endg-begg+1) * (clmobs_layer(i)-1))
+           end if
+
+           !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
+           obs_p(cnt) = clm_obs(i)
+           if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i)
+           cnt = cnt + 1
+
+         end do
+
+       end if
+
      end do
 
      if(obs_interp_switch==1) then

interface/framework/init_dim_obs_pdaf.F90

@@ -125,6 +125,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
 #ifdef CLMFIVE
   use GridcellType, only: grc
   use ColumnType, only : col
+  use PatchType, only : patch
   ! use GetGlobalValuesMod, only: GetGlobalWrite
   ! use clm_varcon, only: nameg
   use enkf_clm_mod, only: state_clm2pdaf_p
@@ -141,6 +142,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
   USE enkf_clm_mod, only: domain_def_clm
   USE enkf_clm_mod, only: get_interp_idx
   use enkf_clm_mod, only: clmstatevec_allcol
+  use enkf_clm_mod, only: clmupdate_swc
+  use enkf_clm_mod, only: clmupdate_T
   !hcp end
 #endif
 #endif
@@ -160,9 +163,12 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
   integer :: ierror
   INTEGER :: max_var_id         ! Multi-scale DA
   INTEGER :: sum_dim_obs_p
+  INTEGER :: p                ! CLM Patch index
   INTEGER :: c                ! CLM Column index
   INTEGER :: g                ! CLM Gridcell index
   INTEGER :: cg
+  INTEGER :: pg
+  INTEGER :: pc
   INTEGER :: i,j,k        ! Counters
   INTEGER :: cnt          ! Counters
   INTEGER :: cnt_interp   ! Counter for interpolation grid cells
@@ -928,6 +934,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
      do i = 1, dim_obs
         obs(i) = clm_obs(i)
 
+       if(clmupdate_swc==1) then
+
        do g = begg,endg
          newgridcell = .true.
 
@@ -997,9 +1005,101 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
              end if
 
            end if
-
          end do
        end do
+
+       else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then
+#ifdef CLMFIVE
+         ! patch loop
+         do g = begg,endg
+           newgridcell = .true.
+
+           do p = begp,endp
+
+             pg = patch%gridcell(p)
+             pc = patch%column(p)
+
+             if(pg == g) then
+               if(newgridcell) then
+                 ! Sets first patch/column in a gridcell. TODO: Make
+                 ! patch / column information part of the observation
+                 ! file
+
+                 if(is_use_dr) then
+                   deltax = abs(lon(g)-clmobs_lon(i))
+                   deltay = abs(lat(g)-clmobs_lat(i))
+                 end if
+
+                 if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. &
+                   ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then
+
+                   ! Set index in state vector, LST will be computed
+                   ! for first patch appearing here
+                   obs_index_p(cnt) = state_clm2pdaf_p(p,1)
+
+                   !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
+                   obs_p(cnt) = clm_obs(i)
+                   if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i)
+                   cnt = cnt + 1
+
+                 end if
+
+                 newgridcell = .false.
+
+               end if
+             end if
+
+           end do
+         end do
+#else
+         ! gridcell loop
+         do g = begg,endg
+
+           if(is_use_dr) then
+             deltax = abs(lon(g)-clmobs_lon(i))
+             deltay = abs(lat(g)-clmobs_lat(i))
+           end if
+
+           if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. &
+             ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then
+             obs_index_p(cnt) = g-begg+1
+           end if
+
+           !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
+           obs_p(cnt) = clm_obs(i)
+           if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i)
+           cnt = cnt + 1
+
+         end do
+#endif
+       else
+
+         print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p."
+         print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop."
+         ! call abort_parallel()
+
+         ! gridcell loop with layer information as default!
+         do g = begg,endg
+
+           if(is_use_dr) then
+             deltax = abs(lon(g)-clmobs_lon(i))
+             deltay = abs(lat(g)-clmobs_lat(i))
+           end if
+
+           if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. &
+             ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then
+             obs_index_p(cnt) = g-begg+1 +  ((endg-begg+1) * (clmobs_layer(i)-1))
+           end if
+
+           !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
+           obs_p(cnt) = clm_obs(i)
+           if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i)
+           cnt = cnt + 1
+
+         end do
+
+       end if
+
      end do
 
      if(obs_interp_switch==1) then

interface/framework/obs_op_pdaf.F90

@@ -133,13 +133,14 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 
   lpointobs = .false.
 
+  ! 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`
+
   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 &
@@ -151,6 +152,18 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
 !  write(*,*) 'model LST', m_state_p(:)
 !  write(*,*) 'TG', state_p(obs_index_p(:))
 !  write(*,*) 'TV', state_p(clm_varsize+obs_index_p(:))
+
+endif
+
+if (clmupdate_T==2 .OR. clmupdate_T==3 .OR. clmupdate_T==4 .OR. clmupdate_T==5) then
+
+  lpointobs = .false.
+
+  DO i = 1, dim_obs_p
+    ! first implementation: simulated LST equals TSKIN
+    m_state_p(i) = state_p(obs_index_p(i))
+  END DO
+
 endif
 #endif