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 3045f721..8b46f851 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 @@ -120,7 +120,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) use ColumnType, only : col ! use GetGlobalValuesMod, only: GetGlobalWrite ! use clm_varcon, only: nameg - use enkf_clm_mod, only: col_index_hydr_act + use enkf_clm_mod, only: state_clm2pdaf_p use enkf_clm_mod, only: clmstatevec_only_active use enkf_clm_mod, only: clmstatevec_max_layer #else @@ -951,7 +951,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f) print *, "clmstatevec_max_layer=", clmstatevec_max_layer call abort_parallel() end if - obs_index_p(cnt) = col_index_hydr_act(c,clmobs_layer(i)) + 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)) 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 6c3c1ed8..e93d110e 100755 --- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 @@ -116,7 +116,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) use ColumnType, only : col ! use GetGlobalValuesMod, only: GetGlobalWrite ! use clm_varcon, only: nameg - use enkf_clm_mod, only: col_index_hydr_act + use enkf_clm_mod, only: state_clm2pdaf_p use enkf_clm_mod, only: clmstatevec_only_active use enkf_clm_mod, only: clmstatevec_max_layer #else @@ -944,7 +944,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p) print *, "clmstatevec_max_layer=", clmstatevec_max_layer call abort_parallel() end if - obs_index_p(cnt) = col_index_hydr_act(c,clmobs_layer(i)) + 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)) diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_5.F90 index 7cce0355..ab11532e 100644 --- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_5.F90 +++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_5.F90 @@ -224,11 +224,17 @@ subroutine clm_finalize() bind(C,name="clm_finalize") ! use ESMF, only : ESMF_Initialize, ESMF_Finalize use cime_comp_mod, only : cime_final + use enkf_clm_mod, only : cleanup_clm_statevec implicit none call cime_final() +#if defined CLMSA + ! TSMP-PDAF: Deallocate arrays from `define_clm_statevec` + call cleanup_clm_statevec() +#endif + !-------------------------------------------------------------------------- ! Clean-up !-------------------------------------------------------------------------- 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 45b85993..0a5cf6a8 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 @@ -46,7 +46,7 @@ module enkf_clm_mod ! 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) - integer, allocatable :: col_index_hydr_act(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1) + integer, allocatable :: state_clm2pdaf_p(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1) integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture @@ -93,7 +93,10 @@ subroutine define_clm_statevec(mype) integer,intent(in) :: mype integer :: i + integer :: j + integer :: jj integer :: c + integer :: g integer :: cc integer :: cccheck @@ -123,10 +126,8 @@ subroutine define_clm_statevec(mype) if(clmstatevec_only_active .eq. 1) then - IF (allocated(col_index_hydr_act)) deallocate(col_index_hydr_act) - allocate(col_index_hydr_act(begc:endc,min(nlevsoi,clmstatevec_max_layer))) - - clm_statevecsize = 0 + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,min(nlevsoi,clmstatevec_max_layer))) cc = 0 @@ -138,47 +139,119 @@ subroutine define_clm_statevec(mype) ! and layers above bedrock if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then cc = cc + 1 - col_index_hydr_act(c,i) = cc + state_clm2pdaf_p(c,i) = cc else - col_index_hydr_act(c,i) = ispval + state_clm2pdaf_p(c,i) = ispval end if end if end do end do - clm_statevecsize = cc - - ! ! Check against other method of computation - ! cccheck = 0 - ! do c=clm_begc,clm_endc - ! if(col%hydrologically_active(c)) then - ! ! Only count non-bedrock layers - ! cccheck = cccheck + col%nbedrock(c) - ! ! Possible -1 to leave out layer that is partly bedrock - ! end if - ! end do - - ! if(clm_statevecsize .ne. cccheck) then - ! print *, "clm_statevecsize", clm_statevecsize - ! print *, "cccheck", cccheck - ! error stop "clm_statevecsize not equal to cccheck" - ! end if - ! Set `clm_varsize`, even though it is currently not used ! for `clmupdate_swc.eq.1` clm_varsize = clm_statevecsize + clm_statevecsize = cc + + 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 i=1,nlevsoi + do c=clm_begc,clm_endc + ! Only take into account layers above input maximum layer + if(i<=clmstatevec_max_layer) then + ! Only take into account hydrologically active columns + ! and layers above bedrock + if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then + cc = cc + 1 + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + end if + end if + end do + end do else + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + + do i=1,nlevsoi + do c=clm_begc,clm_endc + state_clm2pdaf_p(c,i) = (c - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1) + end do + end do + ! #cols values per grid-cell clm_varsize = (endc-begc+1) * nlevsoi clm_statevecsize = (endc-begc+1) * nlevsoi + + 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 i=1,nlevsoi + do c=clm_begc,clm_endc + cc = cc + 1 + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + end do + end do + end if else + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + + do i=1,nlevsoi + do c=clm_begc,clm_endc + ! All columns in a gridcell are assigned the updated + ! gridcell-SWC + state_clm2pdaf_p(c,i) = (col%gridcell(c) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) + end do + end do + ! One value per grid-cell clm_varsize = (endg-begg+1) * nlevsoi clm_statevecsize = (endg-begg+1) * nlevsoi + 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 i=1,nlevsoi + do j=clm_begg,clm_endg + + ! SWC from the first column of each gridcell + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + ! Possibliy: Add state_pdaf2clm_g_p + state_pdaf2clm_c_p(cc) = jj + state_pdaf2clm_j_p(cc) = i + end if + end if + end do + + cc = cc + 1 + end do + end do + + + end if endif @@ -207,13 +280,9 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_statevecsize is ',clm_statevecsize IF (allocated(clm_statevec)) deallocate(clm_statevec) - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) if ((clmupdate_swc.ne.0) .or. (clmupdate_T.ne.0) .or. (clmupdate_texture.ne.0)) then !hcp added condition allocate(clm_statevec(clm_statevecsize)) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) end if !write(*,*) 'clm_paramsize is ',clm_paramsize @@ -224,6 +293,18 @@ subroutine define_clm_statevec(mype) end subroutine define_clm_statevec + subroutine cleanup_clm_statevec() + + implicit none + + ! Deallocate arrays from `define_clm_statevec` + IF (allocated(clm_statevec)) deallocate(clm_statevec) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + + end subroutine cleanup_clm_statevec + subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : soilstate_inst, waterstate_inst use clm_varpar , only : nlevsoi @@ -263,52 +344,10 @@ subroutine set_clm_statevec(tstartcycle, mype) endif if(clmupdate_swc.ne.0) then - ! write swc values to state vector - cc = 1 - do i=1,nlevsoi - - if(clmstatevec_allcol.eq.1) then - - do jj=clm_begc,clm_endc - ! SWC from all columns of each gridcell - if(clmstatevec_only_active.eq.1) then - if(i<=clmstatevec_max_layer .and. col%hydrologically_active(jj) .and. i<=col%nbedrock(jj) ) then - clm_statevec(cc+offset) = swc(jj,i) - state_pdaf2clm_c_p(cc+offset) = jj - state_pdaf2clm_j_p(cc+offset) = i - cc = cc + 1 - end if - else - clm_statevec(cc+offset) = swc(jj,i) - state_pdaf2clm_c_p(cc+offset) = jj - state_pdaf2clm_j_p(cc+offset) = i - cc = cc + 1 - end if - - end do - - else - - do j=clm_begg,clm_endg - ! SWC from the first column of each gridcell - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - clm_statevec(cc+offset) = swc(jj,i) - state_pdaf2clm_c_p(cc+offset) = jj - state_pdaf2clm_j_p(cc+offset) = i - endif - endif - end do - cc = cc + 1 - end do - - end if - - end do + ! write swc values to state vector + do cc = 1, clm_statevecsize + clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do endif !hcp LAI @@ -378,6 +417,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: rliq,rice real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) 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 character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu @@ -466,23 +506,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! do j=clm_begg,clm_endg do j=clm_begc,clm_endc - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.1) then - if(clmstatevec_only_active.eq.1) then - if(i<=clmstatevec_max_layer .and. col%hydrologically_active(j) .and. i<=col%nbedrock(j) ) then - cc = col_index_hydr_act(j,i) - else - cycle - end if - else - cc = (j - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1) - end if - else - cc = (col%gridcell(j) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) - end if - if(swc(j,i).eq.0.0) then swc_zero_before_update = .true. @@ -499,12 +522,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) end if - if(clm_statevec(cc+offset).le.watmin_check) then + swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + + if(swc_update.le.watmin_check) then swc(j,i) = watmin_set - else if(clm_statevec(cc+offset).ge.watsat(j,i)) then + else if(swc_update.ge.watsat(j,i)) then swc(j,i) = watsat(j,i) else - swc(j,i) = clm_statevec(cc+offset) + swc(j,i) = swc_update endif if (isnan(swc(j,i))) then diff --git a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_5.F90 index 7cce0355..ab11532e 100644 --- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_5.F90 +++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_5.F90 @@ -224,11 +224,17 @@ subroutine clm_finalize() bind(C,name="clm_finalize") ! use ESMF, only : ESMF_Initialize, ESMF_Finalize use cime_comp_mod, only : cime_final + use enkf_clm_mod, only : cleanup_clm_statevec implicit none call cime_final() +#if defined CLMSA + ! TSMP-PDAF: Deallocate arrays from `define_clm_statevec` + call cleanup_clm_statevec() +#endif + !-------------------------------------------------------------------------- ! Clean-up !-------------------------------------------------------------------------- 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 45b85993..0a5cf6a8 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 @@ -46,7 +46,7 @@ module enkf_clm_mod ! 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) - integer, allocatable :: col_index_hydr_act(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1) + integer, allocatable :: state_clm2pdaf_p(:,:) !Index of column in hydraulic active state vector (nlevsoi,endc-begc+1) integer(c_int),bind(C,name="clmupdate_swc") :: clmupdate_swc integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture @@ -93,7 +93,10 @@ subroutine define_clm_statevec(mype) integer,intent(in) :: mype integer :: i + integer :: j + integer :: jj integer :: c + integer :: g integer :: cc integer :: cccheck @@ -123,10 +126,8 @@ subroutine define_clm_statevec(mype) if(clmstatevec_only_active .eq. 1) then - IF (allocated(col_index_hydr_act)) deallocate(col_index_hydr_act) - allocate(col_index_hydr_act(begc:endc,min(nlevsoi,clmstatevec_max_layer))) - - clm_statevecsize = 0 + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,min(nlevsoi,clmstatevec_max_layer))) cc = 0 @@ -138,47 +139,119 @@ subroutine define_clm_statevec(mype) ! and layers above bedrock if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then cc = cc + 1 - col_index_hydr_act(c,i) = cc + state_clm2pdaf_p(c,i) = cc else - col_index_hydr_act(c,i) = ispval + state_clm2pdaf_p(c,i) = ispval end if end if end do end do - clm_statevecsize = cc - - ! ! Check against other method of computation - ! cccheck = 0 - ! do c=clm_begc,clm_endc - ! if(col%hydrologically_active(c)) then - ! ! Only count non-bedrock layers - ! cccheck = cccheck + col%nbedrock(c) - ! ! Possible -1 to leave out layer that is partly bedrock - ! end if - ! end do - - ! if(clm_statevecsize .ne. cccheck) then - ! print *, "clm_statevecsize", clm_statevecsize - ! print *, "cccheck", cccheck - ! error stop "clm_statevecsize not equal to cccheck" - ! end if - ! Set `clm_varsize`, even though it is currently not used ! for `clmupdate_swc.eq.1` clm_varsize = clm_statevecsize + clm_statevecsize = cc + + 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 i=1,nlevsoi + do c=clm_begc,clm_endc + ! Only take into account layers above input maximum layer + if(i<=clmstatevec_max_layer) then + ! Only take into account hydrologically active columns + ! and layers above bedrock + if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then + cc = cc + 1 + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + end if + end if + end do + end do else + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + + do i=1,nlevsoi + do c=clm_begc,clm_endc + state_clm2pdaf_p(c,i) = (c - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1) + end do + end do + ! #cols values per grid-cell clm_varsize = (endc-begc+1) * nlevsoi clm_statevecsize = (endc-begc+1) * nlevsoi + + 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 i=1,nlevsoi + do c=clm_begc,clm_endc + cc = cc + 1 + state_pdaf2clm_c_p(cc) = c + state_pdaf2clm_j_p(cc) = i + end do + end do + end if else + + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + allocate(state_clm2pdaf_p(begc:endc,nlevsoi)) + + do i=1,nlevsoi + do c=clm_begc,clm_endc + ! All columns in a gridcell are assigned the updated + ! gridcell-SWC + state_clm2pdaf_p(c,i) = (col%gridcell(c) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) + end do + end do + ! One value per grid-cell clm_varsize = (endg-begg+1) * nlevsoi clm_statevecsize = (endg-begg+1) * nlevsoi + 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 i=1,nlevsoi + do j=clm_begg,clm_endg + + ! SWC from the first column of each gridcell + newgridcell = .true. + do jj=clm_begc,clm_endc + g = col%gridcell(jj) + if (g .eq. j) then + if (newgridcell) then + newgridcell = .false. + ! Possibliy: Add state_pdaf2clm_g_p + state_pdaf2clm_c_p(cc) = jj + state_pdaf2clm_j_p(cc) = i + end if + end if + end do + + cc = cc + 1 + end do + end do + + + end if endif @@ -207,13 +280,9 @@ subroutine define_clm_statevec(mype) !write(*,*) 'clm_statevecsize is ',clm_statevecsize IF (allocated(clm_statevec)) deallocate(clm_statevec) - IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) - IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) if ((clmupdate_swc.ne.0) .or. (clmupdate_T.ne.0) .or. (clmupdate_texture.ne.0)) then !hcp added condition allocate(clm_statevec(clm_statevecsize)) - allocate(state_pdaf2clm_c_p(clm_statevecsize)) - allocate(state_pdaf2clm_j_p(clm_statevecsize)) end if !write(*,*) 'clm_paramsize is ',clm_paramsize @@ -224,6 +293,18 @@ subroutine define_clm_statevec(mype) end subroutine define_clm_statevec + subroutine cleanup_clm_statevec() + + implicit none + + ! Deallocate arrays from `define_clm_statevec` + IF (allocated(clm_statevec)) deallocate(clm_statevec) + IF (allocated(state_pdaf2clm_c_p)) deallocate(state_pdaf2clm_c_p) + IF (allocated(state_pdaf2clm_j_p)) deallocate(state_pdaf2clm_j_p) + IF (allocated(state_clm2pdaf_p)) deallocate(state_clm2pdaf_p) + + end subroutine cleanup_clm_statevec + subroutine set_clm_statevec(tstartcycle, mype) use clm_instMod, only : soilstate_inst, waterstate_inst use clm_varpar , only : nlevsoi @@ -263,52 +344,10 @@ subroutine set_clm_statevec(tstartcycle, mype) endif if(clmupdate_swc.ne.0) then - ! write swc values to state vector - cc = 1 - do i=1,nlevsoi - - if(clmstatevec_allcol.eq.1) then - - do jj=clm_begc,clm_endc - ! SWC from all columns of each gridcell - if(clmstatevec_only_active.eq.1) then - if(i<=clmstatevec_max_layer .and. col%hydrologically_active(jj) .and. i<=col%nbedrock(jj) ) then - clm_statevec(cc+offset) = swc(jj,i) - state_pdaf2clm_c_p(cc+offset) = jj - state_pdaf2clm_j_p(cc+offset) = i - cc = cc + 1 - end if - else - clm_statevec(cc+offset) = swc(jj,i) - state_pdaf2clm_c_p(cc+offset) = jj - state_pdaf2clm_j_p(cc+offset) = i - cc = cc + 1 - end if - - end do - - else - - do j=clm_begg,clm_endg - ! SWC from the first column of each gridcell - newgridcell = .true. - do jj=clm_begc,clm_endc - g = col%gridcell(jj) - if (g .eq. j) then - if (newgridcell) then - newgridcell = .false. - clm_statevec(cc+offset) = swc(jj,i) - state_pdaf2clm_c_p(cc+offset) = jj - state_pdaf2clm_j_p(cc+offset) = i - endif - endif - end do - cc = cc + 1 - end do - - end if - - end do + ! write swc values to state vector + do cc = 1, clm_statevecsize + clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) + end do endif !hcp LAI @@ -378,6 +417,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8) :: rliq,rice real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) 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 character (len = 31) :: fn !TSMP-PDAF: function name for state vector outpu @@ -466,23 +506,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! do j=clm_begg,clm_endg do j=clm_begc,clm_endc - ! Set cc (the state vector index) from the - ! CLM5-grid-index and the `CLM5-layer-index times - ! num_gridcells` - if(clmstatevec_allcol.eq.1) then - if(clmstatevec_only_active.eq.1) then - if(i<=clmstatevec_max_layer .and. col%hydrologically_active(j) .and. i<=col%nbedrock(j) ) then - cc = col_index_hydr_act(j,i) - else - cycle - end if - else - cc = (j - clm_begc + 1) + (i - 1)*(clm_endc - clm_begc + 1) - end if - else - cc = (col%gridcell(j) - clm_begg + 1) + (i - 1)*(clm_endg - clm_begg + 1) - end if - if(swc(j,i).eq.0.0) then swc_zero_before_update = .true. @@ -499,12 +522,14 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) end if - if(clm_statevec(cc+offset).le.watmin_check) then + swc_update = clm_statevec(state_clm2pdaf_p(j,i)) + + if(swc_update.le.watmin_check) then swc(j,i) = watmin_set - else if(clm_statevec(cc+offset).ge.watsat(j,i)) then + else if(swc_update.ge.watsat(j,i)) then swc(j,i) = watsat(j,i) else - swc(j,i) = clm_statevec(cc+offset) + swc(j,i) = swc_update endif if (isnan(swc(j,i))) then