diff --git a/interface/framework/init_dim_obs_f_pdaf.F90 b/interface/framework/init_dim_obs_f_pdaf.F90 index a0280f9c..252bcc35 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,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) integer :: ierror INTEGER :: max_var_id ! Multi-scale DA INTEGER :: sum_dim_obs_p + INTEGER :: p ! CLM Patch index INTEGER :: c ! CLM Column index INTEGER :: g ! CLM Gridcell index INTEGER :: cg + INTEGER :: pg + INTEGER :: pc INTEGER :: i,j,k ! Counters INTEGER :: cnt ! Counters INTEGER :: cnt_interp ! Counter for interpolation grid cells @@ -935,6 +941,8 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) do i = 1, dim_obs obs(i) = clm_obs(i) + if(clmupdate_swc==1) then + do g = begg,endg newgridcell = .true. @@ -1004,9 +1012,101 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) end if end if - end do end do + + else if(clmupdate_T==1 .or. clmupdate_T==2 .or. clmupdate_T==3 .or. clmupdate_T==4) then +#ifdef CLMFIVE + ! patch loop + do g = begg,endg + newgridcell = .true. + + do p = begp,endp + + pg = patch%gridcell(p) + pc = patch%column(p) + + if(pg == g) then + if(newgridcell) then + ! Sets first patch/column in a gridcell. TODO: Make + ! patch / column information part of the observation + ! file + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + + ! Set index in state vector, LST will be computed + ! for first patch appearing here + obs_index_p(cnt) = state_clm2pdaf_p(p,1) + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end if + + newgridcell = .false. + + end if + end if + + end do + end do +#else + ! gridcell loop + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do +#endif + else + + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p." + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop." + ! call abort_parallel() + + ! gridcell loop with layer information as default! + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do + + end if + end do if(obs_interp_switch==1) then diff --git a/interface/framework/init_dim_obs_pdaf.F90 b/interface/framework/init_dim_obs_pdaf.F90 index 6ce590e6..4855fe47 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 .or. clmupdate_T==3 .or. clmupdate_T==4 .or. clmupdate_T==5) then +#ifdef CLMFIVE + ! patch loop + do g = begg,endg + newgridcell = .true. + + do p = begp,endp + + pg = patch%gridcell(p) + pc = patch%column(p) + + if(pg == g) then + if(newgridcell) then + ! Sets first patch/column in a gridcell. TODO: Make + ! patch / column information part of the observation + ! file + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + + ! Set index in state vector, LST will be computed + ! for first patch appearing here + obs_index_p(cnt) = state_clm2pdaf_p(p,1) + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end if + + newgridcell = .false. + + end if + end if + + end do + end do +#else + ! gridcell loop + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do +#endif + else + + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING unsupported update in setting obs_index_p." + print *, "TSMP-PDAF mype(w)=", mype_world, ": WARNING using default gridcell loop." + ! call abort_parallel() + + ! gridcell loop with layer information as default! + do g = begg,endg + + if(is_use_dr) then + deltax = abs(lon(g)-clmobs_lon(i)) + deltay = abs(lat(g)-clmobs_lat(i)) + end if + + if(((is_use_dr).and.(deltax<=clmobs_dr(1)).and.(deltay<=clmobs_dr(2))).or. & + ((.not. is_use_dr).and.(longxy_obs(i) == longxy(g-begg+1)) .and. (latixy_obs(i) == latixy(g-begg+1)))) then + obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1)) + end if + + !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt) + obs_p(cnt) = clm_obs(i) + if(multierr==1) clm_obserr_p(cnt) = clm_obserr(i) + cnt = cnt + 1 + + end do + + end if + end do if(obs_interp_switch==1) then diff --git a/interface/framework/obs_op_pdaf.F90 b/interface/framework/obs_op_pdaf.F90 index 1ce5db81..c5859515 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 .OR. clmupdate_T==3 .OR. clmupdate_T==4 .OR. clmupdate_T==5) then + + lpointobs = .false. + + DO i = 1, dim_obs_p + ! first implementation: simulated LST equals TSKIN + m_state_p(i) = state_p(obs_index_p(i)) + END DO + endif #endif diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 7462a6f1..690c1f3a 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -216,6 +216,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 +483,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 +518,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..9d59d484 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -46,6 +46,7 @@ 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(:) ! clm_paramarr: Contains LAI used in obs_op_pdaf for computing model @@ -139,8 +140,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 @@ -166,8 +167,8 @@ subroutine define_clm_statevec(mype) !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 +368,229 @@ 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 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)) + + 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==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==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==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 + + if(clmupdate_T==5) 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 = 2 * (endp-begp+1) !TSKIN 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 !TV + state_pdaf2clm_c_p(cc+clm_varsize) = patch%column(p) !TV + state_pdaf2clm_j_p(cc+clm_varsize) = 1 + + ! Check that cc is state_clm2pdaf_p + if (cc /= state_clm2pdaf_p(p,1)) then + print *, "ERROR: Index problem in clmupdate_T==5" + print *, "ERROR: Differences cc,state_clm2pdaf_p(p,1) = ", cc, state_clm2pdaf_p(p,1) + end if + + end do + + endif + + end subroutine define_clm_statevec_T + subroutine cleanup_clm_statevec() implicit none @@ -431,8 +655,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 +760,108 @@ 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 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 + integer :: n_c + 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 state vector + if(clmupdate_T==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 updating state vector with skin, soil and vegetation temperature + if(clmupdate_T==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 updating state vector with soil temperature + if(clmupdate_T==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 + + ! skin temperature updating state vector with skin and vegetation temperature + if(clmupdate_T==5) 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_veg(state_pdaf2clm_p_p(cc+clm_varsize)) + 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 +936,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 +1125,96 @@ 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 + + implicit none + + integer,intent(in) :: tstartcycle + integer,intent(in) :: mype + + integer :: i + integer :: j + integer :: c + integer :: p + integer :: lev + + real(r8), pointer :: t_grnd(:) + real(r8), pointer :: t_soisno(:,:) + real(r8), pointer :: t_veg(:) + real(r8), pointer :: t_skin(:) + + 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 state vector + if(clmupdate_T==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 updating skin, soil and vegetation temperature + if(clmupdate_T==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 updating soil temperature + if(clmupdate_T==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 + + ! skin temperature updating skin and vegetation temperature + if(clmupdate_T==5) then + do p = clm_begp, clm_endp + t_skin(p) = clm_statevec(state_clm2pdaf_p(p,1)) + t_veg(p) = clm_statevec(state_clm2pdaf_p(p,1) + clm_varsize) + 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