diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 index d3f1484d3..98b44597c 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 @@ -119,6 +119,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 @@ -135,6 +136,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 @@ -156,9 +159,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 @@ -921,68 +927,161 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) cnt = 1 do i = 1, dim_obs - obs(i) = clm_obs(i) - do g = begg,endg - newgridcell = .true. + obs(i) = clm_obs(i) - do c = begc,endc + if(clmupdate_swc.eq.1) then - cg = mycgridcell(c) + do g = begg,endg + newgridcell = .true. - if(cg .eq. g) then + do c = begc,endc - if(newgridcell) then + cg = mycgridcell(c) - if(is_use_dr) then - deltax = abs(lon(g)-clmobs_lon(i)) - deltay = abs(lat(g)-clmobs_lat(i)) - end if + if(cg .eq. g) then + + if(newgridcell) then + + 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.le.clmobs_dr(1)).and.(deltay.le.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 + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.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 - ! Different settings of observation-location-index in - ! state vector depending on the method of state - ! vector assembling. - if(clmstatevec_allcol.eq.1) then + ! Different settings of observation-location-index in + ! state vector depending on the method of state + ! vector assembling. + if(clmstatevec_allcol.eq.1) then #ifdef CLMFIVE - if(clmstatevec_only_active.eq.1) then - - ! Error if observation deeper than clmstatevec_max_layer - if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." - print *, "i=", i - print *, "c=", c - print *, "clmobs_layer(i)=", clmobs_layer(i) - print *, "col%nbedrock(c)=", col%nbedrock(c) - print *, "clmstatevec_max_layer=", clmstatevec_max_layer - call abort_parallel() - end if - obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) - else + if(clmstatevec_only_active.eq.1) then + + ! Error if observation deeper than clmstatevec_max_layer + if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." + print *, "i=", i + print *, "c=", c + print *, "clmobs_layer(i)=", clmobs_layer(i) + print *, "col%nbedrock(c)=", col%nbedrock(c) + print *, "clmstatevec_max_layer=", clmstatevec_max_layer + call abort_parallel() + end if + obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) + else #endif - obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) + obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) #ifdef CLMFIVE - end if + end if #endif - else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + else + 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.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 end if - !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) - obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) - cnt = cnt + 1 + newgridcell = .false. + end if - newgridcell = .false. + end if + + end do + end do + + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3 .or. clmupdate_T.eq.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 .eq. 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.le.clmobs_dr(1)).and.(deltay.le.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.eq.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.le.clmobs_dr(1)).and.(deltay.le.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.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + end do - 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.le.clmobs_dr(1)).and.(deltay.le.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.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do + + end if + end do if(obs_interp_switch.eq.1) then diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 index e1a6ee4e0..c17dc28f8 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 @@ -115,6 +115,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 @@ -131,6 +132,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 @@ -152,9 +155,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 @@ -914,68 +920,161 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) cnt = 1 do i = 1, dim_obs - obs(i) = clm_obs(i) - do g = begg,endg - newgridcell = .true. + obs(i) = clm_obs(i) - do c = begc,endc + if(clmupdate_swc.eq.1) then - cg = mycgridcell(c) + do g = begg,endg + newgridcell = .true. - if(cg .eq. g) then + do c = begc,endc - if(newgridcell) then + cg = mycgridcell(c) - if(is_use_dr) then - deltax = abs(lon(g)-clmobs_lon(i)) - deltay = abs(lat(g)-clmobs_lat(i)) - end if + if(cg .eq. g) then + + if(newgridcell) then + + 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.le.clmobs_dr(1)).and.(deltay.le.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 + if(((is_use_dr).and.(deltax.le.clmobs_dr(1)).and.(deltay.le.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 - ! Different settings of observation-location-index in - ! state vector depending on the method of state - ! vector assembling. - if(clmstatevec_allcol.eq.1) then + ! Different settings of observation-location-index in + ! state vector depending on the method of state + ! vector assembling. + if(clmstatevec_allcol.eq.1) then #ifdef CLMFIVE - if(clmstatevec_only_active.eq.1) then - - ! Error if observation deeper than clmstatevec_max_layer - if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then - print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." - print *, "i=", i - print *, "c=", c - print *, "clmobs_layer(i)=", clmobs_layer(i) - print *, "col%nbedrock(c)=", col%nbedrock(c) - print *, "clmstatevec_max_layer=", clmstatevec_max_layer - call abort_parallel() - end if - obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) - else + if(clmstatevec_only_active.eq.1) then + + ! Error if observation deeper than clmstatevec_max_layer + if(clmobs_layer(i) > min(clmstatevec_max_layer, col%nbedrock(c))) then + print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR observation layer deeper than clmstatevec_max_layer or bedrock." + print *, "i=", i + print *, "c=", c + print *, "clmobs_layer(i)=", clmobs_layer(i) + print *, "col%nbedrock(c)=", col%nbedrock(c) + print *, "clmstatevec_max_layer=", clmstatevec_max_layer + call abort_parallel() + end if + obs_index_p(cnt) = state_clm2pdaf_p(c,clmobs_layer(i)) + else #endif - obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) + obs_index_p(cnt) = c-begc+1 + ((endc-begc+1) * (clmobs_layer(i)-1)) #ifdef CLMFIVE - end if + end if #endif - else - obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + else + 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.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 end if - !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) - obs_p(cnt) = clm_obs(i) - if(multierr.eq.1) clm_obserr_p(cnt) = clm_obserr(i) - cnt = cnt + 1 + newgridcell = .false. + end if - newgridcell = .false. + end if + + end do + end do + + else if(clmupdate_T.eq.1 .or. clmupdate_T.eq.2 .or. clmupdate_T.eq.3 .or. clmupdate_T.eq.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 .eq. 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.le.clmobs_dr(1)).and.(deltay.le.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.eq.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.le.clmobs_dr(1)).and.(deltay.le.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.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + end do - 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.le.clmobs_dr(1)).and.(deltay.le.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.eq.1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do + + end if + end do if(obs_interp_switch.eq.1) then diff --git a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 index 97becfc66..62ef12e5a 100644 --- a/bldsva/intf_DA/pdaf/framework/obs_op_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/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.EQ.2 .OR. clmupdate_T.EQ.3 .OR. clmupdate_T.EQ.4) 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/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 b/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 index 5a8114de0..58aa5e9f9 100755 --- a/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 +++ b/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 @@ -214,6 +214,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 @@ -473,32 +481,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 < 0) 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 @@ -534,6 +516,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 < 0) 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/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 index 45b849411..bdb7df336 100755 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 @@ -42,6 +42,7 @@ module enkf_clm_mod integer :: clm_begp,clm_endp real(r8),allocatable :: clm_statevec(:) integer,allocatable :: state_pdaf2clm_c_p(:) + integer,allocatable :: state_pdaf2clm_p_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) @@ -86,8 +87,10 @@ subroutine define_clm_statevec(mype) use shr_kind_mod, only: r8 => shr_kind_r8 use decompMod , only : get_proc_bounds use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use clm_varcon , only : ispval use ColumnType , only : col + use PatchType , only : patch implicit none @@ -96,6 +99,8 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: lev + integer :: p integer :: c integer :: g integer :: cc @@ -272,10 +277,150 @@ subroutine define_clm_statevec(mype) !hcp LST DA if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + + 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)) + + 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 + ! skin temperature state vector + if(clmupdate_T.eq.2) 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 = 3* (endp-begp+1) !TSKIN, then TG and 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)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TG + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TG + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 + end do + + endif + + if(clmupdate_T.eq.3) 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 = (1 + nlevgrnd + 1)* (endp-begp+1) !TSKIN, then nlevgrnd times TSOIL and 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)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 + do lev=1,nlevgrnd + ! ivar = 2-26: TSOIL + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) + state_pdaf2clm_j_p(cc + lev*clm_varsize) = lev + 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) = 1 + end do + + endif + + if(clmupdate_T.eq.4) 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 = 1* (endp-begp+1) !TSOIL + + 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)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSOIL + state_pdaf2clm_c_p(cc) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc) = 1 + end do + + endif + + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize WRITE(*, '(a,x,a,i5,x,a,i10)') "TSMP-PDAF-debug", "mype(w)=", mype, "define_clm_statevec: clm_statevecsize=", clm_statevecsize @@ -290,8 +435,8 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp - if ((clmupdate_T.ne.0)) then !hcp - error stop "Not implemented clmupdate_T.NE.0" + if ((clmupdate_T.eq.1)) then !hcp + allocate(clm_paramarr(clm_paramsize)) end if end subroutine define_clm_statevec @@ -309,11 +454,16 @@ subroutine cleanup_clm_statevec() end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst + use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col + use PatchType , only : patch use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer,intent(in) :: tstartcycle @@ -322,7 +472,13 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + 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 :: i,j,jj,g,cc=0,offset=0 + integer :: lev character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -331,6 +487,14 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + ! 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 < 0) THEN @@ -342,6 +506,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 @@ -359,10 +531,50 @@ subroutine set_clm_statevec(tstartcycle, mype) !hcp LAI if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + 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 state vector + if(clmupdate_T.eq.2) then + do cc = 1, clm_varsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc+clm_varsize) = t_grnd(state_pdaf2clm_c_p(cc+clm_varsize)) + clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) + end do + endif + + ! skin temperature state vector updating soil temperature + if(clmupdate_T.eq.3) then + do cc = 1, clm_varsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + do lev=1,nlevgrnd + clm_statevec(cc+lev*clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+lev*clm_varsize), state_pdaf2clm_j_p(cc+lev*clm_varsize)) + end do + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize)) + end do + endif + + ! soil temperature state vector updating soil temperature + if(clmupdate_T.eq.4) then + do cc = 1, clm_varsize + clm_statevec(cc) = t_soisno(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -400,10 +612,14 @@ end subroutine set_clm_statevec subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 + use PatchType , only : patch use ColumnType , only : col - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst use clm_varcon , only : denh2o, denice, watmin use clm_varcon , only : ispval use clm_varcon , only : spval @@ -419,6 +635,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) @@ -427,7 +648,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop - integer :: i,j,jj,g,cc=0,offset=0 + integer :: i,j,jj,g,c,p,cc=0,offset=0 + integer :: lev character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn3 !TSMP-PDAF: function name for state vector outpu @@ -459,6 +681,13 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + ! 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 + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -575,40 +804,48 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end do -#ifdef PDAF_DEBUG - IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) 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 if(clmupdate_T.EQ.1) then - error stop "Not implemented: clmupdate_T.eq.1" + 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 state vector + if(clmupdate_T.EQ.2) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) + end do + endif + + ! skin temperature state vector updating soil temperature + if(clmupdate_T.EQ.3) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + do lev=1,nlevgrnd + t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1) + (lev)*clm_varsize) + end do + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + (1+nlevgrnd)*clm_varsize) + end do + endif + + ! soil temperature state vector updating soil temperature + if(clmupdate_T.EQ.4) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1)) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi @@ -634,6 +871,40 @@ 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 < 0) 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) + CLOSE(71) + END IF + + END IF +#endif + end subroutine update_clm subroutine clm_correct_texture() diff --git a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 index 45b849411..bdb7df336 100755 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 @@ -42,6 +42,7 @@ module enkf_clm_mod integer :: clm_begp,clm_endp real(r8),allocatable :: clm_statevec(:) integer,allocatable :: state_pdaf2clm_c_p(:) + integer,allocatable :: state_pdaf2clm_p_p(:) integer,allocatable :: state_pdaf2clm_j_p(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model ! LST in LST assimilation (clmupdate_T) @@ -86,8 +87,10 @@ subroutine define_clm_statevec(mype) use shr_kind_mod, only: r8 => shr_kind_r8 use decompMod , only : get_proc_bounds use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use clm_varcon , only : ispval use ColumnType , only : col + use PatchType , only : patch implicit none @@ -96,6 +99,8 @@ subroutine define_clm_statevec(mype) integer :: i integer :: j integer :: jj + integer :: lev + integer :: p integer :: c integer :: g integer :: cc @@ -272,10 +277,150 @@ subroutine define_clm_statevec(mype) !hcp LST DA if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + + 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)) + + 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 + ! skin temperature state vector + if(clmupdate_T.eq.2) 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 = 3* (endp-begp+1) !TSKIN, then TG and 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)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 + state_pdaf2clm_p_p(cc+clm_varsize) = p !TG + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TG + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + state_pdaf2clm_p_p(cc+2*clm_varsize) = p !TV + state_pdaf2clm_c_p(cc+2*clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+2*clm_varsize) = 1 + end do + + endif + + if(clmupdate_T.eq.3) 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 = (1 + nlevgrnd + 1)* (endp-begp+1) !TSKIN, then nlevgrnd times TSOIL and 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)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSKIN + state_pdaf2clm_c_p(cc) = patch%column(p) !TSKIN + state_pdaf2clm_j_p(cc) = 1 + do lev=1,nlevgrnd + ! ivar = 2-26: TSOIL + state_pdaf2clm_p_p(cc + lev*clm_varsize) = p + state_pdaf2clm_c_p(cc + lev*clm_varsize) = patch%column(p) + state_pdaf2clm_j_p(cc + lev*clm_varsize) = lev + 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) = 1 + end do + + endif + + if(clmupdate_T.eq.4) 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 = 1* (endp-begp+1) !TSOIL + + 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)) + + cc = 0 + + do p=clm_begp,clm_endp + cc = cc + 1 + state_pdaf2clm_p_p(cc) = p !TSOIL + state_pdaf2clm_c_p(cc) = patch%column(p) !TSOIL + state_pdaf2clm_j_p(cc) = 1 + end do + + endif + + #ifdef PDAF_DEBUG ! Debug output of clm_statevecsize WRITE(*, '(a,x,a,i5,x,a,i10)') "TSMP-PDAF-debug", "mype(w)=", mype, "define_clm_statevec: clm_statevecsize=", clm_statevecsize @@ -290,8 +435,8 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp - if ((clmupdate_T.ne.0)) then !hcp - error stop "Not implemented clmupdate_T.NE.0" + if ((clmupdate_T.eq.1)) then !hcp + allocate(clm_paramarr(clm_paramsize)) end if end subroutine define_clm_statevec @@ -309,11 +454,16 @@ subroutine cleanup_clm_statevec() end subroutine cleanup_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst + use clm_instMod, only : canopystate_inst use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col + use PatchType , only : patch use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer,intent(in) :: tstartcycle @@ -322,7 +472,13 @@ subroutine set_clm_statevec(tstartcycle, mype) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + 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 :: i,j,jj,g,cc=0,offset=0 + integer :: lev character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output @@ -331,6 +487,14 @@ subroutine set_clm_statevec(tstartcycle, mype) pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col + ! 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 < 0) THEN @@ -342,6 +506,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 @@ -359,10 +531,50 @@ subroutine set_clm_statevec(tstartcycle, mype) !hcp LAI if(clmupdate_T.eq.1) then - error stop "Not implemented: clmupdate_T.eq.1" + 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 state vector + if(clmupdate_T.eq.2) then + do cc = 1, clm_varsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + clm_statevec(cc+clm_varsize) = t_grnd(state_pdaf2clm_c_p(cc+clm_varsize)) + clm_statevec(cc+2*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+2*clm_varsize)) + end do + endif + + ! skin temperature state vector updating soil temperature + if(clmupdate_T.eq.3) then + do cc = 1, clm_varsize + ! t_skin iterated over patches + clm_statevec(cc) = t_skin(state_pdaf2clm_p_p(cc)) + do lev=1,nlevgrnd + clm_statevec(cc+lev*clm_varsize) = t_soisno(state_pdaf2clm_c_p(cc+lev*clm_varsize), state_pdaf2clm_j_p(cc+lev*clm_varsize)) + end do + clm_statevec(cc+(1+nlevgrnd)*clm_varsize) = t_veg(state_pdaf2clm_p_p(cc+(1+nlevgrnd)*clm_varsize)) + end do + endif + + ! soil temperature state vector updating soil temperature + if(clmupdate_T.eq.4) then + do cc = 1, clm_varsize + clm_statevec(cc) = t_soisno(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do + endif + ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" @@ -400,10 +612,14 @@ end subroutine set_clm_statevec subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_varpar , only : nlevsoi + use clm_varpar , only : nlevgrnd use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 + use PatchType , only : patch use ColumnType , only : col - use clm_instMod, only : soilstate_inst, waterstate_inst + use clm_instMod, only : soilstate_inst + use clm_instMod, only : waterstate_inst + use clm_instMod, only : temperature_inst use clm_varcon , only : denh2o, denice, watmin use clm_varcon , only : ispval use clm_varcon , only : spval @@ -419,6 +635,11 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) @@ -427,7 +648,8 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) real(r8) :: swc_update ! updated SWC in loop - integer :: i,j,jj,g,cc=0,offset=0 + integer :: i,j,jj,g,c,p,cc=0,offset=0 + integer :: lev character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn3 !TSMP-PDAF: function name for state vector outpu @@ -459,6 +681,13 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col + ! 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 + #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN @@ -575,40 +804,48 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") end do end do -#ifdef PDAF_DEBUG - IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) 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 if(clmupdate_T.EQ.1) then - error stop "Not implemented: clmupdate_T.eq.1" + 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 state vector + if(clmupdate_T.EQ.2) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_grnd(c) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + 2*clm_varsize) + end do + endif + + ! skin temperature state vector updating soil temperature + if(clmupdate_T.EQ.3) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + do lev=1,nlevgrnd + t_soisno(c,lev) = clm_statevec(state_clm2pdaf_p(p,1) + (lev)*clm_varsize) + end do + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + (1+nlevgrnd)*clm_varsize) + end do + endif + + ! soil temperature state vector updating soil temperature + if(clmupdate_T.EQ.4) then + do p = clm_begp, clm_endp + c = patch%column(p) + t_soisno(c,1) = clm_statevec(state_clm2pdaf_p(p,1)) + end do + endif + !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi @@ -634,6 +871,40 @@ 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 < 0) 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) + CLOSE(71) + END IF + + END IF +#endif + end subroutine update_clm subroutine clm_correct_texture() diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index 4d1082a1a..0831ac59e 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -469,14 +469,49 @@ CLM (standalone only). ### CLM:update_T ### -`CLM:update_T`: (integer) Flag for updating of ground and vegetation -temperature. +`CLM:update_T`: (integer) Flag for updating of CLM's ground and +vegetation temperature. + +`CLM:update_T=0` +: No update of ground and vegetation temperature + +`CLM:update_T=1` +: Update of ground and vegetation temperature, based on + Kustas2009-observation operator + +`CLM:update_T=2` +: Update of ground and vegetation temperature based on TSKIN in + observation operator + +`CLM:update_T=3` +: Update of first-layer-soil-temperature and vegetation temperature + based on TSKIN in observation operator + +#### Kustas2009 observation operator + +Variable names in this subsection reflect the conventions from +Kustas2009, not from CLM model. + +LST is compared to the surface temperature $T_R$ defined by +Kustas2009: + +{attribution="Kustas2009 (), Eq(7)"} +>\begin{align*} +>T_R(\phi) \approx \left[ f_C(\phi) T_{C}^{4} + (1 - f_C(\phi)) T_{S}^{4} \right]^{\frac{1}{4}} +>\end{align*} +> +> where $T_C$ is canopy temperature, $T_{S}$ is soil temperature, and +> $f_C(\phi)$ is the fractional vegetation cover observed at the +> radiometer view angle $\phi$. For a canopy with a spherical leaf +> angle distribution and leaf area index $LAI$, +> +>\begin{align*} +> f_C(\phi) = 1 - \exp\left(\frac{-0.5 \, \Omega \, LAI}{\cos\phi}\right) +>\end{align*} +> +> where the factor $\Omega$ indicates the degree to which vegetation +> is clumped as in row crops or sparsely vegetated shrubland canopies. -Currently only CLM3.5 - -- 0: No update of ground and vegetation temperature - -- 1: Update of ground and vegetation temperature ### CLM:print_swc ###