Tsmp pdaf patched lstda cap

docs/users_guide/running_tsmp_pdaf/input_enkfpf.md

@@ -63,6 +63,10 @@ statevec_max_layer =
 t_printensemble =
 watmin_switch =
 swc_mask_snow =
+T_mask_snow =
+increment_type =
+T_mask_T =
+T_max_increment =
 
 [COSMO]
 nprocs      =
@@ -484,12 +488,14 @@ CLM (standalone only).
 `CLM:update_T`: (integer) Flag for updating of ground and vegetation
 temperature.
 
-Currently only CLM3.5
-
 -  0: No update of ground and vegetation temperature
 
 -  1: Update of ground and vegetation temperature
 
+-  2: Update of soil and vegetation temperature (prognostic variables)
+   based on comparison of satellite LST with skin temperature
+   (diagnostic variable)
+
 ### CLM:print_swc ###
 
 `CLM:print_swc`: (integer) If set to `1`, the updated soil moisture
@@ -601,6 +607,50 @@ Only takes effect if `CLM:update_swc``is switched on.
 
 Default setting is `0`: No masking of columns with snow cover.
 
+### CLM:T_mask_snow ###
+
+`CLM:T_mask_snow`: (integer) Switch for masking columns with snow
+cover from T updates.
+
+Snow covers larger than 1mm are switched off for the update.
+
+Only takes effect if `CLM:update_T``is switched on.
+
+Default setting is `0`: No masking of columns with snow cover.
+
+### CLM:increment_type ###
+
+`CLM:increment_type`: (integer) Switch for changing increment type in
+T-update.
+
+- `0`: Multiplicative increment
+- `1`: Additive increment
+
+Only takes effect if `CLM:update_T` is switched on.
+
+Default setting is `0`: Multiplicative increment.
+
+### CLM:T_max_increment ###
+
+`CLM:T_max_increment`: (double) Maximum T-increment to update
+(additively).
+
+Only takes effect if `CLM:update_T` is switched on and
+`CLM:increment_type` is set to `1`.
+
+Default setting is `5.0`: Updates larger than 5K are not applied.
+
+### CLM:T_mask_T ###
+
+`CLM:T_mask_T`: (double) Temperature difference to add to freezing
+temperature resulting in a threshold temperature for masking out T
+updates.  Whenever the soil temperature of the surface layer falls
+below this threshold temperature, no T is not updated.
+
+Only takes effect if `CLM:update_T``is switched on.
+
+Default setting is `0.`: Masking updates below freezing temperaturs.
+
 (enkfpf:cosmo)=
 ## [COSMO] ##
 
@@ -900,6 +950,8 @@ Default: 0, output turned off.
  |           | `statevec_max_layer`    | 25            |
  |           | `t_printensemble`       | -2            |
  |           | `watmin_switch`         | 0             |
+ |           | `T_mask_snow`           | 0             |
+ |           | `T_mask_T`              | 0.0           |
  | `[COSMO]` |                         |               |
  |           | `nprocs`                | 0             |
  |           | `dtmult`                | 0             |

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,11 @@ 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 :: i,j,k        ! Counters
   INTEGER :: cnt          ! Counters
   INTEGER :: cnt_interp   ! Counter for interpolation grid cells
@@ -935,6 +940,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.
 
@@ -1005,9 +1012,100 @@ 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) then
+#ifdef CLMFIVE
+         ! patch loop
+         do g = begg,endg
+           newgridcell = .true.
+
+           do p = begp,endp
+
+             pg = patch%gridcell(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.
 
@@ -998,9 +1006,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) 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/localize_covar_pdaf.F90

@@ -40,7 +40,6 @@ SUBROUTINE localize_covar_pdaf(dim_p, dim_obs, HP, HPH)
    USE shr_kind_mod , only : r8 => shr_kind_r8
    USE mod_read_obs, ONLY: clmobs_lon
    USE mod_read_obs, ONLY: clmobs_lat
-   USE enkf_clm_mod, ONLY: clmupdate_T
    USE enkf_clm_mod, ONLY: clm_begc
    USE enkf_clm_mod, ONLY: clm_endc
    USE enkf_clm_mod, ONLY: state_pdaf2clm_c_p