TSMP-PDAF Update

bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90

@@ -94,7 +94,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
        dampfac_state_time_dependent_in, dampfac_param_time_dependent_in
   use mod_tsmp, &
       only: idx_map_subvec2state_fortran, tag_model_parflow, enkf_subvecsize, &
-      nx_glob, ny_glob, nz_glob, crns_flag, &
+      nx_glob, ny_glob, nz_glob, crns_flag, da_print_obs_index, &
 #ifndef CLMSA
 #ifndef OBS_ONLY_CLM
       xcoord, ycoord, zcoord, xcoord_fortran, ycoord_fortran, &
@@ -170,6 +170,8 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
 #endif
 #endif
 
+  character (len = 27) :: fn    !TSMP-PDAF: function name for obs_index_p output
+
   ! ****************************************
   ! *** Initialize observation dimension ***
   ! ****************************************
@@ -866,6 +868,18 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
 
      end if
 
+#ifdef PDAF_DEBUG
+     IF (da_print_obs_index > 0) THEN
+       ! TSMP-PDAF: For debug runs, output the state vector in files
+       WRITE(fn, "(a,i5.5,a,i5.5,a)") "obs_index_p_", mype_world, ".", step, ".txt"
+       OPEN(unit=71, file=fn, action="write")
+       DO i = 1, dim_obs_p
+         WRITE (71,"(i10)") obs_index_p(i)
+       END DO
+       CLOSE(71)
+     END IF
+#endif
+
   end if
   end if
   call mpi_allreduce(MPI_IN_PLACE,obs_nc2pdaf,dim_obs,MPI_INTEGER,MPI_SUM,comm_filter,ierror)

bldsva/intf_DA/pdaf/framework/mod_tsmp.F90

@@ -35,6 +35,7 @@ module mod_tsmp
     integer(c_int), bind(c)  :: tag_model_parflow = 1
     integer(c_int), bind(c)  :: tag_model_cosmo   = 2
     integer(c_int), bind(c)  :: crns_flag
+    integer(c_int), bind(c)  :: da_print_obs_index
     type(c_ptr), bind(c)     :: pf_statevec
     type(c_ptr), bind(c)     :: xcoord, ycoord, zcoord
     real(c_double), pointer  :: xcoord_fortran(:), ycoord_fortran(:), zcoord_fortran(:)

bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90

@@ -115,6 +115,8 @@ SUBROUTINE obs_op_pdaf(step, dim_p, dim_obs_p, state_p, m_state_p)
   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
      m_state_p(i) &
     = (exp(-0.5*clm_paramarr(obs_index_p(i))) &
                      *state_p(obs_index_p(i))**4 & 

bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm.F90

@@ -32,7 +32,7 @@
 !> @details
 !> Initialization routines for the land model analogous to
 !> clm3_5/src/main/program_off.F90
-subroutine clm_init(finname) bind(C,name="clm_init")
+subroutine clm_init(finname, pdaf_id, pdaf_max, mype) bind(C,name="clm_init")
   use iso_C_binding
   use enkf_clm_mod, only: &
 #if (defined COUP_OAS_COS || defined COUP_OAS_PFL)
@@ -61,6 +61,9 @@ subroutine clm_init(finname) bind(C,name="clm_init")
 
 !  character(c_char),target   :: finname
   character(kind=c_char,len=1),dimension(100),intent(in) :: finname
+  integer(c_int), intent(in) :: pdaf_id !unused: only for clm5_0
+  integer(c_int), intent(in) :: pdaf_max !unused: only for clm5_0
+  integer(c_int), intent(in) :: mype
   integer(c_int) :: counter
   !character(100),pointer :: pchar
 
@@ -143,17 +146,19 @@ subroutine clm_init(finname) bind(C,name="clm_init")
 
 
 #if defined CLMSA
-  call define_clm_statevec
+  call define_clm_statevec(mype)
 #endif
 
 end subroutine clm_init
 
 
-subroutine clm_advance(ntstep) bind(C,name="clm_advance")
+subroutine clm_advance(ntstep, tstartcycle, mype) bind(C,name="clm_advance")
   use enkf_clm_mod
   use iso_c_binding
   use mod_clm_statistics
   integer(c_int),intent(in) :: ntstep
+  integer(c_int),intent(in) :: tstartcycle
+  integer(c_int),intent(in) :: mype
   integer :: counter=0
 
    do counter=1,ntstep
@@ -192,7 +197,8 @@ subroutine clm_advance(ntstep) bind(C,name="clm_advance")
   end do
 
 #if defined CLMSA
-  call set_clm_statevec()
+  ! Calling PDAF Function to set state vector before assimiliation
+  call set_clm_statevec(tstartcycle, mype)
 #endif
 
 end subroutine clm_advance