diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md index 790b4aa0..6ba12242 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md +++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md @@ -484,12 +484,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 diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index a0280f9c..ea2d7531 100755 --- a/interface/framework/init_dim_obs_f_pdaf.F90 +++ b/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. @@ -1004,9 +1011,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) 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 diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 6ce590e6..f25490c7 100755 --- a/interface/framework/init_dim_obs_pdaf.F90 +++ b/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) 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 diff --git a/interface/framework/localize_covar_pdaf.F90 b/interface/framework/localize_covar_pdaf.F90 index b0865fb8..eb93a453 100644 --- a/interface/framework/localize_covar_pdaf.F90 +++ b/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 diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 1ce5db81..78d3e05c 100644 --- a/interface/framework/obs_op_pdaf.F90 +++ b/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) 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 diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 7462a6f1..8e3a3f3d 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -61,6 +61,7 @@ module enkf_clm_mod integer,allocatable :: state_pdaf2clm_c_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) integer,allocatable :: state_loc2clm_c_p(:) + integer,allocatable :: state_loc2clm_p_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) @@ -216,6 +217,14 @@ subroutine set_clm_statevec(tstartcycle, mype) CLOSE(71) END IF + IF(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + END IF #endif @@ -475,32 +484,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end do -#ifdef PDAF_DEBUG - IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN - - IF(clmupdate_swc.NE.0) THEN - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn3, action="write") - WRITE (71,"(es22.15)") h2osoi_liq(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn4, action="write") - WRITE (71,"(es22.15)") h2osoi_ice(:,:) - CLOSE(71) - - ! TSMP-PDAF: For debug runs, output the state vector in files - WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" - OPEN(unit=71, file=fn2, action="write") - WRITE (71,"(es22.15)") swc(:,:) - CLOSE(71) - END IF - - END IF -#endif - endif !hcp: TG, TV @@ -536,6 +519,39 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") call clm_texture_to_parameters endif +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN + + IF(clmupdate_swc.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn3, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn3, action="write") + WRITE (71,"(es22.15)") h2osoi_liq(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn4, action="write") + WRITE (71,"(es22.15)") h2osoi_ice(:,:) + CLOSE(71) + + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") swc(:,:) + CLOSE(71) + END IF + + IF(clmupdate_T.NE.0) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn2, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn2, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + END IF + + END IF +#endif + end subroutine update_clm subroutine clm_correct_texture() diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index 81ce1fd2..195d7168 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -46,8 +46,10 @@ module enkf_clm_mod real(r8),allocatable :: clm_statevec(:) real(r8),allocatable :: clm_statevec_orig(:) integer,allocatable :: state_pdaf2clm_c_p(:) + integer,allocatable :: state_pdaf2clm_p_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) integer,allocatable :: state_loc2clm_c_p(:) + integer,allocatable :: state_loc2clm_p_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) @@ -139,8 +141,8 @@ subroutine define_clm_statevec(mype) end if !hcp LST DA - if(clmupdate_T==1) then - error stop "Not implemented: clmupdate_T.eq.1" + if(clmupdate_T/=0) then + call define_clm_statevec_T(mype) end if !end hcp @@ -158,16 +160,16 @@ subroutine define_clm_statevec(mype) ! Allocate statevector-duplicate for saving original column mean ! values used in computing increments during updating the state - ! vector in column-mean-mode. + ! vector in column-mean-mode (SWC) or gridcell-mean-mode (T). IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) - if (clmupdate_swc/=0 .and. clmstatevec_colmean/=0) then + if ((clmupdate_swc/=0 .and. clmstatevec_colmean/=0) .or. clmupdate_T==2) then allocate(clm_statevec_orig(clm_statevecsize)) end if !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp - if ((clmupdate_T/=0)) then !hcp - error stop "Not implemented clmupdate_T.NE.0" + if ((clmupdate_T==1)) then !hcp + allocate(clm_paramarr(clm_paramsize)) end if end subroutine define_clm_statevec @@ -367,6 +369,191 @@ subroutine define_clm_statevec_swc(mype) end subroutine define_clm_statevec_swc + subroutine define_clm_statevec_T(mype) + use decompMod , only : get_proc_bounds + use clm_varpar , only : nlevgrnd + use clm_varcon , only : ispval + use PatchType , only : patch + + implicit none + + integer,intent(in) :: mype + + integer :: j + integer :: jj + integer :: lev + integer :: p + integer :: cc + + integer :: begp, endp ! per-proc beginning and ending pft indices + integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begl, endl ! per-proc beginning and ending landunit indices + integer :: begg, endg ! per-proc gridcell ending gridcell indices + + + call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp) + +#ifdef PDAF_DEBUG + WRITE(*,"(a,i5,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a,i10,a)") & + "TSMP-PDAF mype(w)=", mype, " define_clm_statevec, CLM5-bounds (g,l,c,p)----",& + begg,",",endg,",",begl,",",endl,",",begc,",",endc,",",begp,",",endp," -------" +#endif + + clm_begg = begg + clm_endg = endg + clm_begc = begc + clm_endc = endc + clm_begp = begp + clm_endp = endp + + if(clmupdate_T==1) then + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1)) + + do p=clm_begp,clm_endp + state_clm2pdaf_p(p,1) = (p - clm_begp + 1) + end do + + clm_varsize = endp-begp+1 + clm_paramsize = endp-begp+1 !LAI + clm_statevecsize = (endp-begp+1)*2 !TG, then TV + + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TG + state_pdaf2clm_c_p(cc) = patch%column(p) !TG + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + end do + + endif + !end hcp + + if(clmupdate_T==2) then + + ! 1) PATCH/GRC: CLM->PDAF + ! Allocate with full dimension for all variables/layers + ! + ! Here, the second index of STATE_CLM2PDAF_P is NOT simply the + ! layer index of CLM, but rather a variable index of the state + ! vector, in order tomake the index mapping 1:1. + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begp:endp,1:(2+nlevgrnd))) + ! ^ + ! dimension layout: + ! 1: TSKIN + ! 2:(1+nlevgrnd): TSOIL layers + ! (2+nlevgrnd): TVEG + do lev=1,(2+nlevgrnd) + do p=begp,endp + ! Default: inactive + state_clm2pdaf_p(p,lev) = ispval + end do + end do + + do p=clm_begp,clm_endp + ! All patches in a gridcell are assigned the index of + ! the gridcell-average in the state vector + cc = (patch%gridcell(p) - clm_begg + 1) + + ! TSKIN (variable index 1) + state_clm2pdaf_p(p,1) = cc + + ! TSOIL layers (variable indices 2 to 1+nlevgrnd) + do lev=1,nlevgrnd + state_clm2pdaf_p(p, 1+lev) = cc + lev*(clm_endg - clm_begg + 1) + end do + + ! TVEG (variable index 2+nlevgrnd) + state_clm2pdaf_p(p, 2+nlevgrnd) = cc + (1+nlevgrnd)*(clm_endg - clm_begg + 1) + end do + + ! 2) PATCH/GRC: STATEVECSIZE + clm_varsize = endg-begg+1 + clm_statevecsize = (1 + nlevgrnd + 1)* (endg-begg+1) !TSKIN, then nlevgrnd times TSOIL and TV + + ! 3) PATCH/GRC: PDAF->CLM + ! + ! Inverse mapping from state vector index to CLM indices. + ! STATE_PDAF2CLM_*_P is the inverse of STATE_CLM2PDAF_P: + ! state_pdaf2clm_p_p(cc) = p (patch index) + ! state_pdaf2clm_c_p(cc) = c (column index) + ! state_pdaf2clm_j_p(cc) = j (variable index, NOT CLM level) + ! + ! Variable index j layout (same as STATE_CLM2PDAF_P second dimension): + ! j = 1: TSKIN + ! j = 2:(1+nlevgrnd): TSOIL layers (CLM level = j - 1) + ! j = (2+nlevgrnd): TVEG + ! + ! NOTE: To get the CLM level from j for TSOIL, use: clm_level = j - 1 + ! + IF (allocated(state_pdaf2clm_p_p)) deallocate(state_pdaf2clm_p_p) + allocate(state_pdaf2clm_p_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + allocate(state_pdaf2clm_c_p(clm_statevecsize)) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + allocate(state_pdaf2clm_j_p(clm_statevecsize)) + + ! Default: inactive + do cc=1,clm_statevecsize + state_pdaf2clm_p_p(cc) = ispval + state_pdaf2clm_c_p(cc) = ispval + state_pdaf2clm_j_p(cc) = ispval + end do + + cc = 0 + + do cc=1,clm_statevecsize + + do p=clm_begp,clm_endp + if(state_clm2pdaf_p(p, 1) == cc) then + ! Set indices for the three temperature variables and then + ! exit loops + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 ! variable index 1: TSKIN + do lev=1,nlevgrnd + ! variable index 2 to 1+nlevgrnd: TSOIL layers + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p !TSOIL + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc + lev*clm_varsize) = 1 + lev ! variable index TSOIL + end do + state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+(1+nlevgrnd)*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+(1+nlevgrnd)*clm_varsize) = 2 + nlevgrnd ! variable index: TVEG + exit + end if + end do + + end do + + endif + + if(clmupdate_T < 1 .or. clmupdate_T > 2) then + error stop "LST-DA only implemented for clmupdate_T==1 and clmupdate_T==2." + end if + + end subroutine define_clm_statevec_T + subroutine cleanup_clm_statevec() implicit none @@ -431,8 +618,8 @@ subroutine set_clm_statevec(tstartcycle, mype) endif !hcp LAI - if(clmupdate_T==1) then - error stop "Not implemented: clmupdate_T.eq.1" + if(clmupdate_T/=0) then + call set_clm_statevec_T(tstartcycle,mype) endif !end hcp LAI @@ -536,6 +723,140 @@ subroutine set_clm_statevec_swc() end subroutine set_clm_statevec_swc + subroutine set_clm_statevec_T(tstartcycle,mype) + use clm_instMod, only : temperature_inst + use clm_instMod, only : canopystate_inst + use clm_varpar , only : nlevgrnd + use PatchType , only : patch + use ColumnType , only : col + use shr_kind_mod, only: r8 => shr_kind_r8 + + implicit none + + integer,intent(in) :: tstartcycle + integer,intent(in) :: mype + + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + real(r8), pointer :: tlai(:) + integer :: j,g,cc,c,p + integer :: n_c,n_p + integer :: lev + + character (len = 34) :: fn !TSMP-PDAF: function name for swc output + + ! LST variables + t_grnd => temperature_inst%t_grnd_col + t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch + t_soisno => temperature_inst%t_soisno_col + tlai => canopystate_inst%tlai_patch + + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble == -1) THEN + + IF(clmupdate_T/=0) THEN + ! TSMP-PDAF: Debug output of CLM t_soisno, first layer + WRITE(fn, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".integrate.", tstartcycle + 1, ".txt" + OPEN(unit=71, file=fn, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF + + END IF +#endif + + !hcp LAI + if(clmupdate_T==1) then + do cc = 1, clm_varsize + ! t_grnd iterated over cols + ! t_veg iterated over patches + clm_statevec(cc) = t_grnd(state_pdaf2clm_c_p(cc)) + clm_statevec(cc+clm_varsize) = t_veg( state_pdaf2clm_p_p(cc+clm_varsize)) + end do + + do cc = 1, clm_paramsize + ! Works only if clm_paramsize corresponds to clm_varsize (also + ! the order) + clm_paramarr(cc) = tlai(state_pdaf2clm_p_p(cc)) + end do + endif + !end hcp LAI + + ! Skin temperature updating state vector with skin, soil and vegetation temperature. + ! Uses gridcell mean: averages over all patches (for TSKIN, TVEG) or + ! columns (for TSOIL) within each gridcell. + if(clmupdate_T==2) then + + do cc = 1, clm_varsize + + ! Get gridcell from the reference patch + g = patch%gridcell(state_pdaf2clm_p_p(cc)) + + ! --- TSKIN: average over patches in gridcell --- + clm_statevec(cc) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc) = clm_statevec(cc) + t_skin(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc) = clm_statevec(cc) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TSKIN averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TSKIN)" + end if + + ! --- TSOIL: average over columns in gridcell (per layer) --- + do lev=1,nlevgrnd + clm_statevec(cc+lev*clm_varsize) = 0.0 + n_c = 0 + do c=clm_begc,clm_endc + if(col%gridcell(c)==g) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) + t_soisno(c,lev) + n_c = n_c + 1 + end if + end do + if(n_c > 0) then + clm_statevec(cc+lev*clm_varsize) = clm_statevec(cc+lev*clm_varsize) / real(n_c, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " layer=", lev, " has no columns for TSOIL averaging" + error stop "Gridcell without columns in set_clm_statevec_T (TSOIL)" + end if + end do + + ! --- TVEG: average over patches in gridcell --- + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = 0.0 + n_p = 0 + do p=clm_begp,clm_endp + if(patch%gridcell(p)==g) then + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) + t_veg(p) + n_p = n_p + 1 + end if + end do + if(n_p > 0) then + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = clm_statevec(cc+(1+nlevgrnd)*clm_varsize) / real(n_p, r8) + else + write(*,*) "ERROR: Gridcell g=", g, " has no patches for TVEG averaging" + error stop "Gridcell without patches in set_clm_statevec_T (TVEG)" + end if + + end do + + ! Save prior gridcell mean state vector for computing + ! increment in updating the state vector + do cc = 1, clm_statevecsize + clm_statevec_orig(cc) = clm_statevec(cc) + end do + + endif + + end subroutine set_clm_statevec_T subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_time_manager , only : update_DA_nstep @@ -610,8 +931,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") endif !hcp: TG, TV - if(clmupdate_T==1) then - error stop "Not implemented: clmupdate_T.eq.1" + if(clmupdate_T/=0) then + call update_clm_T(tstartcycle, mype) endif ! end hcp TG, TV @@ -799,6 +1120,114 @@ subroutine update_clm_swc(tstartcycle, mype) end subroutine update_clm_swc + subroutine update_clm_T(tstartcycle,mype) + + use PatchType , only : patch + use clm_varpar , only : nlevgrnd + use clm_instMod, only : temperature_inst + use IEEE_ARITHMETIC, only: ieee_is_nan + + implicit none + + integer,intent(in) :: tstartcycle + integer,intent(in) :: mype + + integer :: i + integer :: j + integer :: c + integer :: p + integer :: cc + integer :: lev + + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + + real(r8) :: increment_factor + real(r8) :: t_update + + character (len = 31) :: fn !TSMP-PDAF: function name for swc output + + ! LST + t_grnd => temperature_inst%t_grnd_col + t_soisno => temperature_inst%t_soisno_col + t_veg => temperature_inst%t_veg_patch + t_skin => temperature_inst%t_skin_patch + ! tlai => canopystate_inst%tlai_patch + + !hcp: TG, TV + if(clmupdate_T==1) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1)) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + end do + endif + ! end hcp TG, TV + + ! Skin temperature updating skin, soil and vegetation temperature. + ! Uses gridcell mean increment factor: applies the ratio of + ! (new gridcell mean / old gridcell mean) to each patch/column value. + if(clmupdate_T==2) then + do p = clm_begp, clm_endp + c = patch%column(p) + + ! --- TSKIN: update with increment factor --- + cc = state_clm2pdaf_p(p,1) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_skin(p) * increment_factor + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_skin update is NaN at p=", p + else + t_skin(p) = t_update + end if + end if + + ! --- TSOIL: update with increment factor for each layer --- + do lev=1,nlevgrnd + cc = state_clm2pdaf_p(p, 1+lev) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_soisno(c,lev) * increment_factor + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_soisno update is NaN at c=", c, " lev=", lev + else + t_soisno(c,lev) = t_update + end if + end if + end do + + ! --- TVEG: update with increment factor --- + cc = state_clm2pdaf_p(p, 2+nlevgrnd) + ! Skip if no significant change in gridcell mean + if(abs(clm_statevec(cc) - clm_statevec_orig(cc)) > 1.0e-7) then + increment_factor = clm_statevec(cc) / clm_statevec_orig(cc) + t_update = t_veg(p) * increment_factor + if (ieee_is_nan(t_update)) then + print *, "WARNING: t_veg update is NaN at p=", p + else + t_veg(p) = t_update + end if + end if + + end do + endif + +#ifdef PDAF_DEBUG + IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble == -1) THEN + ! TSMP-PDAF: For debug runs, output the state vector in files + WRITE(fn, "(a,i5.5,a,i5.5,a)") "t_soisno_", mype, ".update.", tstartcycle, ".txt" + OPEN(unit=71, file=fn, action="write") + WRITE (71,"(es22.15)") t_soisno(:,1) + CLOSE(71) + END IF +#endif + + end subroutine update_clm_T subroutine update_clm_texture(tstartcycle, mype) use clm_varpar , only : nlevsoi @@ -1249,6 +1678,7 @@ subroutine init_n_domains_clm(n_domains_p) use decompMod, only : get_proc_bounds use clm_varcon , only : ispval use ColumnType , only : col + use PatchType , only : patch implicit none @@ -1256,6 +1686,7 @@ subroutine init_n_domains_clm(n_domains_p) integer :: domain_p integer :: begg, endg ! per-proc gridcell ending gridcell indices integer :: begc, endc ! per-proc beginning and ending column indices + integer :: begp, endp ! per-proc beginning and ending pft indices integer :: g integer :: c @@ -1263,7 +1694,7 @@ subroutine init_n_domains_clm(n_domains_p) ! TODO: remove unnecessary calls of get_proc_bounds (use clm_begg, ! clm_endg, etc) - call get_proc_bounds(begg=begg, endg=endg, begc=begc, endc=endc) + call get_proc_bounds(begg=begg, endg=endg, begc=begc, endc=endc, begp=begp, endp=endp) if(clmupdate_swc==1) then if(clmstatevec_allcol==1) then @@ -1275,6 +1706,12 @@ subroutine init_n_domains_clm(n_domains_p) ! -> DIM_L: number of layers in gridcell n_domains_p = endg - begg + 1 end if + else if (clmupdate_T/=0) then + ! For LSTDA: gridcells are domains + ! Each gridcell is a local domain + ! -> DIM_L: number of temperature variables (each with gridcell + ! -> averages) + n_domains_p = endg - begg + 1 else ! Process-local number of gridcells Default, possibly not tested ! for other updates except SWC @@ -1289,80 +1726,111 @@ subroutine init_n_domains_clm(n_domains_p) ! local domain domain_p (from the column, the gridcell can be ! derived) - ! Allocate state_loc2clm_c_p with preliminary n_domains_p - IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) - allocate(state_loc2clm_c_p(n_domains_p)) - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = ispval - end do - - if(clmstatevec_only_active == 1) then - - ! Reset n_domains_p - n_domains_p = 0 - domain_p = 0 - - if(clmstatevec_allcol == 1) then - ! COLUMNS - - ! Each hydrologically active layer is a local domain - ! -> DIM_L: number of layers in hydrologically active column - do c=clm_begc,clm_endc - ! Skip state vector loop and directly check if column is - ! hydrologically active - if(col%hydrologically_active(c)) then - domain_p = domain_p + 1 - n_domains_p = n_domains_p + 1 - state_loc2clm_c_p(domain_p) = c - end if - end do + if(clmupdate_swc==1) then + ! Allocate state_loc2clm_c_p with preliminary n_domains_p + IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) + allocate(state_loc2clm_c_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = ispval + end do - else - ! GRIDCELLS + if(clmstatevec_only_active == 1) then - ! For gridcells - do g = clm_begg,clm_endg + ! Reset n_domains_p + n_domains_p = 0 + domain_p = 0 - ! Search the state vector for col in grc - do cc = 1,clm_statevecsize + if(clmstatevec_allcol == 1) then + ! COLUMNS - if (col%gridcell(state_pdaf2clm_c_p(cc)) == g) then - ! Set local domain index + ! Each hydrologically active column is a local domain + ! -> DIM_L: number of layers in hydrologically active column + do c=clm_begc,clm_endc + ! Skip state vector loop and directly check if column is + ! hydrologically active + if(col%hydrologically_active(c)) then domain_p = domain_p + 1 - ! Set new number of local domains n_domains_p = n_domains_p + 1 - ! Set CLM-column-index corresponding to local domain - state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(cc) - ! Exit state vector loop, when fitting column is found - exit + state_loc2clm_c_p(domain_p) = c end if end do - end do + else + ! GRIDCELLS + + ! For gridcells + do g = clm_begg,clm_endg + + ! Search the state vector for col in grc + do cc = 1,clm_statevecsize + + if (col%gridcell(state_pdaf2clm_c_p(cc)) == g) then + ! Set local domain index + domain_p = domain_p + 1 + ! Set new number of local domains + n_domains_p = n_domains_p + 1 + ! Set CLM-column-index corresponding to local domain + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(cc) + ! Exit state vector loop, when fitting column is found + exit + end if + end do - end if + end do - else + end if - ! Set state_loc2clm_c_p for non-excluding hydrologically - ! inactive cols/grcs - if(clmstatevec_allcol == 1) then - ! COLUMNS - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = clm_begc + domain_p - 1 - end do else - ! GRIDCELLS - do domain_p=1,n_domains_p - state_loc2clm_c_p(domain_p) = clm_begg + domain_p - 1 - end do - end if + ! Set state_loc2clm_c_p for non-excluding hydrologically + ! inactive cols/grcs + if(clmstatevec_allcol == 1) then + ! COLUMNS + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = clm_begc + domain_p - 1 + end do + else + ! GRIDCELLS + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) + end do + end if + + end if end if ! Possibly: Warning when final n_domains_p actually excludes ! hydrologically inactive gridcells + if(clmupdate_T/=0) then + ! For LSTDA: Gridcells are domains + + ! Allocate state_loc2clm_c_p with preliminary n_domains_p + IF (allocated(state_loc2clm_c_p)) deallocate(state_loc2clm_c_p) + allocate(state_loc2clm_c_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = ispval + end do + IF (allocated(state_loc2clm_p_p)) deallocate(state_loc2clm_p_p) + allocate(state_loc2clm_p_p(n_domains_p)) + do domain_p=1,n_domains_p + state_loc2clm_p_p(domain_p) = ispval + end do + + ! Columns from gridcells + ! Patches from gridcells + ! + ! domain_p is equal to cc for the first variable TSKIN + ! Therefore, state_pdaf2clm_c_p / state_pdaf2clm_p_p can be used + ! to derive the corresponding column / patch + do domain_p=1,n_domains_p + state_loc2clm_c_p(domain_p) = state_pdaf2clm_c_p(domain_p) + state_loc2clm_p_p(domain_p) = state_pdaf2clm_p_p(domain_p) + end do + + end if + + end subroutine init_n_domains_clm @@ -1373,6 +1841,7 @@ end subroutine init_n_domains_clm !> This routine sets DIM_L, the local state vector dimension. subroutine init_dim_l_clm(domain_p, dim_l) use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use ColumnType , only : col implicit none @@ -1407,6 +1876,16 @@ subroutine init_dim_l_clm(domain_p, dim_l) dim_l = 3*nlevsoi + nshift endif + if(clmupdate_T==1) then + ! TG + TV: 2 temperatures per patch + dim_l = 2 + endif + + if(clmupdate_T==2) then + ! TSKIN + TSOIL(nlevgrnd layers) + TV + dim_l = 1 + nlevgrnd + 1 + endif + end subroutine init_dim_l_clm !> @author Wolfgang Kurtz, Johannes Keller @@ -1438,11 +1917,21 @@ subroutine g2l_state_clm(domain_p, dim_p, state_p, dim_l, state_l) ! ENDDO ! Column index inside gridcell index domain_p + if(clmupdate_swc==1) then DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index: i state_l(i) = state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) END DO + end if + + if(clmupdate_T==1 .or. clmupdate_T==2) then + DO i = 1, dim_l + ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P + ! Variable index: i + state_l(i) = state_p(state_clm2pdaf_p(state_loc2clm_p_p(domain_p),i)) + END DO + end if end subroutine g2l_state_clm @@ -1476,11 +1965,21 @@ subroutine l2g_state_clm(domain_p, dim_l, state_l, dim_p, state_p) ! ENDDO ! Column index inside gridcell index domain_p + if(clmupdate_swc==1) then DO i = 1, dim_l ! Column index from DOMAIN_P via STATE_LOC2CLM_C_P ! Layer index i state_p(state_clm2pdaf_p(state_loc2clm_c_p(domain_p),i)) = state_l(i) END DO + end if + + if(clmupdate_T==1 .or. clmupdate_T==2) then + DO i = 1, dim_l + ! Patch index from DOMAIN_P via STATE_LOC2CLM_P_P + ! Variable index: i + state_p(state_clm2pdaf_p(state_loc2clm_p_p(domain_p),i)) = state_l(i) + END DO + end if end subroutine l2g_state_clm #endif