!------------------------------------------------------------------------------------------- !Copyright (c) 2013-2016 by Wolfgang Kurtz, Guowei He and Mukund Pondkule (Forschungszentrum Juelich GmbH) ! !This file is part of TSMP-PDAF ! !TSMP-PDAF is free software: you can redistribute it and/or modify !it under the terms of the GNU Lesser General Public License as published by !the Free Software Foundation, either version 3 of the License, or !(at your option) any later version. ! !TSMP-PDAF is distributed in the hope that it will be useful, !but WITHOUT ANY WARRANTY; without even the implied warranty of !MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !GNU LesserGeneral Public License for more details. ! !You should have received a copy of the GNU Lesser General Public License !along with TSMP-PDAF. If not, see <http://www.gnu.org/licenses/>. !------------------------------------------------------------------------------------------- ! ! !------------------------------------------------------------------------------------------- !enkf_clm_mod.F90: Module for CLM !------------------------------------------------------------------------------------------- module enkf_clm_mod use iso_c_binding ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8, SHR_KIND_CL ! !ARGUMENTS: implicit none #if (defined CLMSA) integer :: COMM_model_clm integer :: clm_statevecsize integer :: clm_paramsize !hcp: Size of CLM parameter vector (f.e. LAI) integer :: clm_varsize integer :: clm_begg,clm_endg integer :: clm_begc,clm_endc integer :: clm_begp,clm_endp real(r8),allocatable :: clm_statevec(:) real(r8),allocatable :: clm_statevec_orig(:) integer,allocatable :: state_pdaf2clm_c_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) real(r8),allocatable :: clm_paramarr(:) !hcp CLM parameter vector (f.e. LAI) 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 integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol integer(c_int),bind(C,name="clmstatevec_colmean") :: clmstatevec_colmean integer(c_int),bind(C,name="clmstatevec_only_active") :: clmstatevec_only_active integer(c_int),bind(C,name="clmstatevec_max_layer") :: clmstatevec_max_layer integer(c_int),bind(C,name="clmt_printensemble") :: clmt_printensemble integer(c_int),bind(C,name="clmwatmin_switch") :: clmwatmin_switch real(c_double),bind(C,name="clmcrns_bd") :: clmcrns_bd integer :: nstep ! time step index real(r8) :: dtime ! time step increment (sec) integer :: ier ! error code character(kind=c_char),dimension(100),bind(C,name="outdir"),target :: outdir logical :: log_print ! true=> print diagnostics real(r8) :: eccf ! earth orbit eccentricity factor logical :: mpi_running ! true => MPI is initialized integer :: mpicom_glob ! MPI communicator character(len=SHR_KIND_CL) :: nlfilename = " " integer :: ierror, lengths_of_types, i logical :: flag integer(c_int),bind(C,name="clmprefixlen") :: clmprefixlen integer :: COMM_couple_clm ! CLM-version of COMM_couple ! (currently not used for eclm) logical :: newgridcell !only eclm contains #if defined CLMSA 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_varcon , only : ispval use ColumnType , only : col implicit none integer,intent(in) :: mype integer :: i integer :: j integer :: jj integer :: c integer :: g integer :: cg integer :: cc integer :: cccheck 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 ! Soil Moisture DA: State vector index arrays if(clmupdate_swc.eq.1) then ! 1) COL/GRC: CLM->PDAF 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 ! Default: inactive state_clm2pdaf_p = ispval end do end do ! All column variables in state vector if(clmstatevec_allcol.eq.1) then ! Only hydrologically active columns if(clmstatevec_only_active .eq. 1) then cc = 0 do i=1,nlevsoi ! Only take into account layers above input maximum layer if(i<=clmstatevec_max_layer) then do c=clm_begc,clm_endc ! 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_clm2pdaf_p(c,i) = cc end if end do end if end do ! All column variables in state vector simplifying the indexing else 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 end if ! Gridcell values or averages in state vector else ! Only hydrologically active columns if(clmstatevec_only_active.eq.1) then cc = 0 do i=1,nlevsoi ! Only layers above max_layer if(i<=clmstatevec_max_layer) then do g=clm_begg,clm_endg newgridcell = .true. do c=clm_begc,clm_endc if(col%gridcell(c) == g) then ! All hydrologically active columns above ! bedrock in a gridcell point to the state ! vector index of the gridcell if(col%hydrologically_active(c) .and. i<=col%nbedrock(c)) then if(newgridcell) then ! Update the index if first col found for grc, ! otherwise reproduce previous index cc = cc + 1 newgridcell = .false. end if state_clm2pdaf_p(c,i) = cc end if end if end do end do end if end do else 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 end if end if ! 2) COL/GRC: STATEVECSIZE if(clmstatevec_only_active.eq.1) then ! Use iterator cc for setting state vector size. ! ! Set `clm_varsize`, even though it is currently not used ! for `clmupdate_swc.eq.1` clm_varsize = cc clm_statevecsize = cc else if(clmstatevec_allcol.eq.1) then ! #cols * #levels clm_varsize = (endc-begc+1) * nlevsoi clm_statevecsize = (endc-begc+1) * nlevsoi else ! #grcs * #levels clm_varsize = (endg-begg+1) * nlevsoi clm_statevecsize = (endg-begg+1) * nlevsoi end if end if ! 3) COL/GRC: PDAF->CLM 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)) ! Defaults do cc=1,clm_statevecsize state_pdaf2clm_c_p(cc) = ispval state_pdaf2clm_j_p(cc) = ispval end do do cc=1,clm_statevecsize lay: do i=1,nlevsoi do c=clm_begc,clm_endc if (state_clm2pdaf_p(c,i) == cc) then ! Set column index and then exit loop state_pdaf2clm_c_p(cc) = c state_pdaf2clm_j_p(cc) = i exit lay end if end do end do lay #ifdef PDAF_DEBUG ! Check that all state vectors have been assigned c, i if(state_pdaf2clm_c_p(cc) == ispval) then write(*,*) 'cc: ', cc error stop "state_pdaf2clm_c_p not set at cc" end if if(state_pdaf2clm_j_p(cc) == ispval) then write(*,*) 'cc: ', cc error stop "state_pdaf2clm_j_p not set at cc" end if #endif end do endif if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" endif if(clmupdate_texture.eq.1) then clm_statevecsize = clm_statevecsize + 2*((endg-begg+1)*nlevsoi) endif if(clmupdate_texture.eq.2) then clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi) endif !hcp LST DA if(clmupdate_T.eq.1) then error stop "Not implemented: clmupdate_T.eq.1" endif !end hcp #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 #endif !write(*,*) 'clm_statevecsize is ',clm_statevecsize IF (allocated(clm_statevec)) deallocate(clm_statevec) if ((clmupdate_swc.ne.0) .or. (clmupdate_T.ne.0) .or. (clmupdate_texture.ne.0)) then !hcp added condition allocate(clm_statevec(clm_statevecsize)) end if ! Allocate statevector-duplicate for saving original column mean ! values used in computing increments during updating the state ! vector in column-mean-mode. IF (allocated(clm_statevec_orig)) deallocate(clm_statevec_orig) if (clmupdate_swc.ne.0 .and. clmstatevec_colmean.ne.0) then allocate(clm_statevec_orig(clm_statevecsize)) end if !write(*,*) 'clm_paramsize is ',clm_paramsize if (allocated(clm_paramarr)) deallocate(clm_paramarr) !hcp if ((clmupdate_T.ne.0)) then !hcp error stop "Not implemented clmupdate_T.NE.0" end if 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 ! use clm_varcon, only: nameg, namec ! use GetGlobalValuesMod, only: GetGlobalWrite use ColumnType , only : col use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer,intent(in) :: tstartcycle integer,intent(in) :: mype real(r8), pointer :: swc(:,:) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) integer :: i,j,jj,g,c,cc=0,offset=0 integer :: n_c character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output swc => waterstate_inst%h2osoi_vol_col psand => soilstate_inst%cellsand_col pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN IF(clmupdate_swc.NE.0) THEN ! TSMP-PDAF: Debug output of CLM swc WRITE(fn2, "(a,i5.5,a,i5.5,a)") "swcstate_", mype, ".integrate.", tstartcycle + 1, ".txt" OPEN(unit=71, file=fn2, action="write") WRITE (71,"(es22.15)") swc(:,:) CLOSE(71) END IF END IF #endif ! calculate shift when CRP data are assimilated if(clmupdate_swc.eq.2) then error stop "Not implemented clmupdate_swc.eq.2" endif if(clmupdate_swc.ne.0) then ! write swc values to state vector if (clmstatevec_colmean.eq.1) then do cc = 1, clm_statevecsize clm_statevec(cc) = 0.0 n_c = 0 ! Get gridcell and layer g = col%gridcell(state_pdaf2clm_c_p(cc)) j = state_pdaf2clm_j_p(cc) ! Loop over all columns do c=clm_begc,clm_endc ! Select columns in gridcell g if(col%gridcell(c).eq.g) then ! Select hydrologically active columns if(col%hydrologically_active(c)) then ! Add active column to swc-sum clm_statevec(cc) = clm_statevec(cc) + swc(c,j) n_c = n_c + 1 end if end if end do if(n_c == 0) then write(*,*) "WARNING: Gridcell g=", g write(*,*) "WARNING: Layer j=", j write(*,*) "Grid cell g at layer j without hydrologically active column! Setting SWC as in gridcell mode." clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) else ! Normalize sum to average clm_statevec(cc) = clm_statevec(cc) / real(n_c, r8) end if ! Save prior column mean state vector for computing ! increment in updating the state vector clm_statevec_orig(cc) = clm_statevec(cc) end do else do cc = 1, clm_statevecsize clm_statevec(cc) = swc(state_pdaf2clm_c_p(cc), state_pdaf2clm_j_p(cc)) end do end if endif !hcp LAI if(clmupdate_T.eq.1) then error stop "Not implemented: clmupdate_T.eq.1" endif !end hcp LAI ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" endif ! write texture values to state vector (if desired) if(clmupdate_texture.ne.0) then cc = 1 do i=1,nlevsoi do j=clm_begg,clm_endg clm_statevec(cc+1*clm_varsize+offset) = psand(j,i) clm_statevec(cc+2*clm_varsize+offset) = pclay(j,i) if(clmupdate_texture.eq.2) then !incl. organic matter values clm_statevec(cc+3*clm_varsize+offset) = porgm(j,i) end if cc = cc + 1 end do end do endif #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn, "(a,i5.5,a,i5.5,a)") "clmstate_", mype, ".integrate.", tstartcycle + 1, ".txt" OPEN(unit=71, file=fn, action="write") DO i = 1, clm_statevecsize WRITE (71,"(es22.15)") clm_statevec(i) END DO CLOSE(71) END IF #endif end subroutine set_clm_statevec subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") use clm_varpar , only : nlevsoi use clm_time_manager , only : update_DA_nstep use shr_kind_mod , only : r8 => shr_kind_r8 use ColumnType , only : col use clm_instMod, only : soilstate_inst, waterstate_inst use clm_varcon , only : denh2o, denice, watmin use clm_varcon , only : ispval use clm_varcon , only : spval implicit none integer,intent(in) :: tstartcycle integer,intent(in) :: mype real(r8), pointer :: swc(:,:) real(r8), pointer :: watsat(:,:) real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) real(r8), pointer :: dz(:,:) ! layer thickness depth (m) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) 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 character (len = 31) :: fn2 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn3 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn4 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn5 !TSMP-PDAF: function name for state vector outpu character (len = 32) :: fn6 !TSMP-PDAF: function name for state vector outpu logical :: swc_zero_before_update = .false. #ifdef PDAF_DEBUG IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn, "(a,i5.5,a,i5.5,a)") "clmstate_", mype, ".update.", tstartcycle, ".txt" OPEN(unit=71, file=fn, action="write") DO i = 1, clm_statevecsize WRITE (71,"(es22.15)") clm_statevec(i) END DO CLOSE(71) END IF #endif swc => waterstate_inst%h2osoi_vol_col watsat => soilstate_inst%watsat_col psand => soilstate_inst%cellsand_col pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col dz => col%dz h2osoi_liq => waterstate_inst%h2osoi_liq_col h2osoi_ice => waterstate_inst%h2osoi_ice_col #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(fn5, "(a,i5.5,a,i5.5,a)") "h2osoi_liq", mype, ".bef_up.", tstartcycle, ".txt" OPEN(unit=71, file=fn5, action="write") WRITE (71,"(es22.15)") h2osoi_liq(:,:) CLOSE(71) ! TSMP-PDAF: For debug runs, output the state vector in files WRITE(fn6, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".bef_up.", tstartcycle, ".txt" OPEN(unit=71, file=fn6, action="write") WRITE (71,"(es22.15)") h2osoi_ice(:,:) CLOSE(71) END IF END IF #endif ! calculate shift when CRP data are assimilated if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" endif ! CLM5: Update the Data Assimulation time-step to the current time ! step, since DA has been done. Used by CLM5 to skip BalanceChecks ! directly after the DA step. call update_DA_nstep() ! write updated swc back to CLM if(clmupdate_swc.ne.0) then ! Set minimum soil moisture for checking the state vector and ! for setting minimum swc for CLM if(clmwatmin_switch.eq.3) then ! CLM3.5 type watmin watmin_check = 0.00 watmin_set = 0.05 else if(clmwatmin_switch.eq.5) then ! CLM5.0 type watmin watmin_check = watmin watmin_set = watmin else ! Default watmin_check = 0.0 watmin_set = 0.0 end if ! cc = 0 do i=1,nlevsoi ! CLM3.5: iterate over grid cells ! CLM5.0: iterate over columns ! do j=clm_begg,clm_endg do j=clm_begc,clm_endc ! Update only those SWCs that are not excluded by ispval if(state_clm2pdaf_p(j,i) .ne. ispval) then if(swc(j,i).eq.0.0) then swc_zero_before_update = .true. ! Zero-SWC leads to zero denominator in computation of ! rliq/rice, therefore setting rliq/rice to special ! value rliq = spval rice = spval else swc_zero_before_update = .false. rliq = h2osoi_liq(j,i)/(dz(j,i)*denh2o*swc(j,i)) rice = h2osoi_ice(j,i)/(dz(j,i)*denice*swc(j,i)) !h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) end if if (clmstatevec_colmean.eq.1) then ! Update SWC column value with the increment-factor ! of the state vector update (state vector updates ! are means of cols in grc) swc_update = swc(j,i) * clm_statevec(state_clm2pdaf_p(j,i)) / clm_statevec_orig(state_clm2pdaf_p(j,i)) else ! Update SWC with updated state vector swc_update = clm_statevec(state_clm2pdaf_p(j,i)) end if if(swc_update.le.watmin_check) then swc(j,i) = watmin_set else if(swc_update.ge.watsat(j,i)) then swc(j,i) = watsat(j,i) else swc(j,i) = swc_update endif if (isnan(swc(j,i))) then swc(j,i) = watmin_set print *, "WARNING: swc at j,i is nan: ", j, i endif if(swc_zero_before_update) then ! This case should not appear for hydrologically ! active columns/layers, where always: swc > watmin ! ! If you want to make sure that no zero SWCs appear in ! the code, comment out the error stop #ifdef PDAF_DEBUG ! error stop "ERROR: Update of zero-swc" print *, "WARNING: Update of zero-swc" print *, "WARNING: Any new H2O added to h2osoi_liq(j,i) with j,i = ", j, i #endif h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o h2osoi_ice(j,i) = 0.0 else ! update liquid water content h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o*rliq ! update ice content h2osoi_ice(j,i) = swc(j,i) * dz(j,i)*denice*rice end if end if ! cc = cc + 1 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" endif ! end hcp TG, TV !! update liquid water content !do j=clm_begg,clm_endg ! do i=1,nlevsoi ! h2osoi_liq(j,i) = swc(j,i) * dz(j,i)*denh2o ! end do !end do ! write updated texture back to CLM if(clmupdate_texture.ne.0) then cc = 1 do i=1,nlevsoi do j=clm_begg,clm_endg psand(j,i) = clm_statevec(cc+1*clm_varsize+offset) pclay(j,i) = clm_statevec(cc+2*clm_varsize+offset) if(clmupdate_texture.eq.2) then ! incl. organic matter porgm(j,i) = clm_statevec(cc+3*clm_varsize+offset) end if cc = cc + 1 end do end do call clm_correct_texture call clm_texture_to_parameters endif end subroutine update_clm subroutine clm_correct_texture() use clm_varpar , only : nlevsoi use shr_kind_mod , only : r8 => shr_kind_r8 use clm_instMod, only : soilstate_inst use CNSharedParamsMod, only : CNParamsShareInst implicit none integer :: c,lev real(r8) :: clay,sand,ttot real(r8) :: orgm real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) psand => soilstate_inst%cellsand_col pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col do c = clm_begg,clm_endg do lev = 1,nlevsoi clay = pclay(c,lev) sand = psand(c,lev) if(sand.le.0.0) sand = 1.0 if(clay.le.0.0) clay = 1.0 ttot = sand + clay if(ttot.gt.100) then sand = sand/ttot * 100.0 clay = clay/ttot * 100.0 end if pclay(c,lev) = clay psand(c,lev) = sand if(clmupdate_texture.eq.2) then orgm = (porgm(c,lev) / CNParamsShareInst%organic_max) * 100.0 if(orgm.le.0.0) orgm = 0.0 if(orgm.ge.100.0) orgm = 100.0 porgm(c,lev) = (orgm / 100.0)* CNParamsShareInst%organic_max end if end do end do end subroutine clm_correct_texture subroutine clm_texture_to_parameters() use clm_varpar , only : nlevsoi use clm_varcon , only : zsoi, secspday use shr_kind_mod , only : r8 => shr_kind_r8 use clm_instMod, only : soilstate_inst use FuncPedotransferMod , only : pedotransf, get_ipedof use CNSharedParamsMod, only : CNParamsShareInst implicit none real(r8) :: om_tkm = 0.25_r8 ! thermal conductivity of organic soil (Farouki, 1986) [W/m/K] real(r8) :: om_watsat_lake = 0.9_r8 ! porosity of organic soil real(r8) :: om_hksat_lake = 0.1_r8 ! saturated hydraulic conductivity of organic soil [mm/s] real(r8) :: om_sucsat_lake = 10.3_r8 ! saturated suction for organic matter (Letts, 2000) real(r8) :: om_b_lake = 2.7_r8 ! Clapp Hornberger paramater for oragnic soil (Letts, 2000) (lake) real(r8) :: om_watsat ! porosity of organic soil real(r8) :: om_hksat ! saturated hydraulic conductivity of organic soil [mm/s] real(r8) :: om_sucsat ! saturated suction for organic matter (mm)(Letts, 2000) real(r8) :: om_csol = 2.5_r8 ! heat capacity of peat soil *10^6 (J/K m3) (Farouki, 1986) real(r8) :: om_tkd = 0.05_r8 ! thermal conductivity of dry organic soil (Farouki, 1981) real(r8) :: om_b ! Clapp Hornberger paramater for oragnic soil (Letts, 2000) real(r8) :: zsapric = 0.5_r8 ! depth (m) that organic matter takes on characteristics of sapric peat real(r8) :: pcalpha = 0.5_r8 ! percolation threshold real(r8) :: pcbeta = 0.139_r8 ! percolation exponent real(r8) :: perc_frac ! "percolating" fraction of organic soil real(r8) :: perc_norm ! normalize to 1 when 100% organic soil real(r8) :: uncon_hksat ! series conductivity of mineral/organic soil real(r8) :: uncon_frac ! fraction of "unconnected" soil real(r8) :: bd ! bulk density of dry soil material [kg/m^3] real(r8) :: tkm ! mineral conductivity real(r8) :: xksat ! maximum hydraulic conductivity of soil [mm/s] integer :: ipedof,c,lev real(r8) :: clay,sand,om_frac real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) real(r8), pointer :: porgm(:,:) psand => soilstate_inst%cellsand_col pclay => soilstate_inst%cellclay_col porgm => soilstate_inst%cellorg_col do c = clm_begg,clm_endg do lev = 1,nlevsoi clay = pclay(c,lev) sand = psand(c,lev) om_frac = porgm(c,lev) / CNParamsShareInst%organic_max ipedof=get_ipedof(0) call pedotransf(ipedof, sand, clay, & soilstate_inst%watsat_col(c,lev), & soilstate_inst%bsw_col(c,lev), & soilstate_inst%sucsat_col(c,lev), & xksat) om_watsat = max(0.93_r8 - 0.1_r8 *(zsoi(lev)/zsapric), 0.83_r8) om_b = min(2.7_r8 + 9.3_r8 *(zsoi(lev)/zsapric), 12.0_r8) om_sucsat = min(10.3_r8 - 0.2_r8 *(zsoi(lev)/zsapric), 10.1_r8) om_hksat = max(0.28_r8 - 0.2799_r8*(zsoi(lev)/zsapric), xksat) soilstate_inst%bd_col(c,lev) = & (1._r8 - soilstate_inst%watsat_col(c,lev))*2.7e3_r8 soilstate_inst%watsat_col(c,lev) = & (1._r8 - om_frac) * soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac tkm = (1._r8 - om_frac) * (8.80_r8*sand+2.92_r8*clay)/(sand+clay)+om_tkm*om_frac ! W/(m K) soilstate_inst%bsw_col(c,lev) = & (1._r8-om_frac) * (2.91_r8 + 0.159_r8*clay) + om_frac*om_b soilstate_inst%sucsat_col(c,lev) = & (1._r8-om_frac) * soilstate_inst%sucsat_col(c,lev) + om_sucsat*om_frac soilstate_inst%hksat_min_col(c,lev) = xksat ! perc_frac is zero unless perf_frac greater than percolation threshold if (om_frac > pcalpha) then perc_norm = (1._r8 - pcalpha)**(-pcbeta) perc_frac = perc_norm*(om_frac - pcalpha)**pcbeta else perc_frac = 0._r8 endif ! uncon_frac is fraction of mineral soil plus fraction of ! "nonpercolating" organic soil uncon_frac = (1._r8-om_frac)+(1._r8-perc_frac)*om_frac ! uncon_hksat is series addition of mineral/organic ! conductivites if (om_frac < 1._r8) then uncon_hksat = uncon_frac/((1._r8-om_frac)/xksat & +((1._r8-perc_frac)*om_frac)/om_hksat) else uncon_hksat = 0._r8 end if soilstate_inst%hksat_col(c,lev) = uncon_frac*uncon_hksat + & (perc_frac*om_frac)*om_hksat soilstate_inst%tkmg_col(c,lev) = tkm ** (1._r8 - & soilstate_inst%watsat_col(c,lev)) soilstate_inst%tksatu_col(c,lev) = & soilstate_inst%tkmg_col(c,lev)*0.57_r8**soilstate_inst%watsat_col(c,lev) soilstate_inst%tkdry_col(c,lev) = & ((0.135_r8*soilstate_inst%bd_col(c,lev) + 64.7_r8) / & (2.7e3_r8 - 0.947_r8*soilstate_inst%bd_col(c,lev)) & )*(1._r8-om_frac) + om_tkd*om_frac soilstate_inst%csol_col(c,lev) = & ((1._r8-om_frac)*(2.128_r8*sand+2.385_r8*clay) / (sand+clay) + & om_csol*om_frac)*1.e6_r8 soilstate_inst%watdry_col(c,lev) = & soilstate_inst%watsat_col(c,lev) * & (316230._r8/soilstate_inst%sucsat_col(c,lev) & ) ** (-1._r8/soilstate_inst%bsw_col(c,lev)) soilstate_inst%watopt_col(c,lev) = & soilstate_inst%watsat_col(c,lev) * & (158490._r8/soilstate_inst%sucsat_col(c,lev) & ) ** (-1._r8/soilstate_inst%bsw_col(c,lev)) ! secspday (day/sec) soilstate_inst%watfc_col(c,lev) = & soilstate_inst%watsat_col(c,lev) * & (0.1_r8 / (soilstate_inst%hksat_col(c,lev)*secspday) & )**(1._r8/(2._r8*soilstate_inst%bsw_col(c,lev)+3._r8)) end do end do end subroutine clm_texture_to_parameters ! subroutine average_swc_crp(profdat,profave) ! use clm_varcon , only : zsoi ! implicit none ! real(r8),intent(in) :: profdat(10) ! real(r8),intent(out) :: profave ! error stop "Not implemented average_swc_crp" ! end subroutine average_swc_crp #endif subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, & longxy, latixy, longxy_obs, latixy_obs) use spmdMod, only : npes, iam use domainMod, only : ldomain use decompMod, only : get_proc_total, get_proc_bounds, ldecomp implicit none real, intent(in) :: lon_clmobs(:) real, intent(in) :: lat_clmobs(:) integer, intent(in) :: dim_obs integer, allocatable, intent(inout) :: longxy(:) integer, allocatable, intent(inout) :: latixy(:) integer, allocatable, intent(inout) :: longxy_obs(:) integer, allocatable, intent(inout) :: latixy_obs(:) integer :: ni, nj, ii, jj, kk, cid, ier, ncells, nlunits integer :: ncols, counter integer :: npatches, ncohorts real :: minlon, minlat, maxlon, maxlat real(r8), pointer :: lon(:) real(r8), pointer :: lat(:) integer :: begg, endg ! per-proc gridcell ending gridcell indices lon => ldomain%lonc lat => ldomain%latc ni = ldomain%ni nj = ldomain%nj ! get total number of gridcells, landunits, ! columns, patches and cohorts on processor call get_proc_total(iam, ncells, nlunits, ncols, npatches, ncohorts) ! beg and end gridcell call get_proc_bounds(begg=begg, endg=endg) !print *,'ni, nj ', ni, nj !print *,'cells per processor ', ncells !print *,'begg, endg ', begg, endg ! allocate vector with size of elements in x directions * size of elements in y directions allocate(longxy(ncells), stat=ier) allocate(latixy(ncells), stat=ier) ! initialize vector with zero values longxy(:) = 0 latixy(:) = 0 ! fill vector with index values counter = 1 do ii = 1, nj do jj = 1, ni cid = (ii-1)*ni + jj do kk = begg, endg if(cid == ldecomp%gdc2glo(kk)) then latixy(counter) = ii longxy(counter) = jj counter = counter + 1 end if end do end do end do ! set intial values for max/min of lon/lat minlon = 999 minlat = 999 maxlon = -999 maxlat = -999 ! looping over all cell centers to get min/max longitude and latitude minlon = MINVAL(lon(:) + 180) maxlon = MAXVAL(lon(:) + 180) minlat = MINVAL(lat(:) + 90) maxlat = MAXVAL(lat(:) + 90) allocate(longxy_obs(dim_obs), stat=ier) allocate(latixy_obs(dim_obs), stat=ier) do i = 1, dim_obs if(((lon_clmobs(i) + 180) - minlon) /= 0 .and. & ((lat_clmobs(i) + 90) - minlat) /= 0) then longxy_obs(i) = ceiling(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon)) !+ 1 latixy_obs(i) = ceiling(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat)) !+ 1 !print *,'longxy_obs(i) , latixy_obs(i) ', longxy_obs(i) , latixy_obs(i) else if(((lon_clmobs(i) + 180) - minlon) == 0 .and. & ((lat_clmobs(i) + 90) - minlat) == 0) then longxy_obs(i) = 1 latixy_obs(i) = 1 else if(((lon_clmobs(i) + 180) - minlon) == 0) then longxy_obs(i) = 1 latixy_obs(i) = ceiling(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat)) else if(((lat_clmobs(i) + 90) - minlat) == 0) then longxy_obs(i) = ceiling(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon)) latixy_obs(i) = 1 endif end do ! deallocate temporary arrays !deallocate(longxy) !deallocate(latixy) !deallocate(longxy_obs) !deallocate(latixy_obs) end subroutine domain_def_clm !> @author Mukund Pondkule, Johannes Keller !> @date 27.03.2023 !> @brief Set indices of grid cells with lon/lat smaller than observation locations !> @details !> This routine sets the indices of grid cells with lon/lat !> smaller than observation locations. subroutine get_interp_idx(lon_clmobs, lat_clmobs, dim_obs, longxy_obs_floor, latixy_obs_floor) USE domainMod, ONLY: ldomain ! USE decompMod, ONLY: get_proc_total, get_proc_bounds_atm, adecomp ! USE spmdMod, ONLY: npes, iam ! number of processors for clm and processor number ! USE mod_read_obs, ONLY: lon_clmobs, lat_clmobs ! USE clmtype, ONLY : clm3 ! USE mod_assimilation, ONLY: dim_obs !#endif implicit none real, intent(in) :: lon_clmobs(:) real, intent(in) :: lat_clmobs(:) integer, intent(in) :: dim_obs integer, allocatable, intent(inout) :: longxy_obs_floor(:) integer, allocatable, intent(inout) :: latixy_obs_floor(:) integer :: i integer :: ni, nj ! integer :: ii, jj, kk, cid integer :: ier ! integer :: ncells, nlunits,ncols, npfts ! integer :: counter real :: minlon, minlat, maxlon, maxlat real(r8), pointer :: lon(:) real(r8), pointer :: lat(:) ! integer :: begg, endg ! per-proc gridcell ending gridcell indices lon => ldomain%lonc lat => ldomain%latc ni = ldomain%ni nj = ldomain%nj ! ! get total number of gridcells, landunits, ! ! columns and pfts for any processor ! call get_proc_total(iam, ncells, nlunits, ncols, npfts) ! ! beg and end gridcell for atm ! call get_proc_bounds_atm(begg, endg) ! set intial values for max/min of lon/lat minlon = 999 minlat = 999 maxlon = -999 maxlat = -999 ! looping over all cell centers to get min/max longitude and latitude minlon = MINVAL(lon(:) + 180) maxlon = MAXVAL(lon(:) + 180) minlat = MINVAL(lat(:) + 90) maxlat = MAXVAL(lat(:) + 90) if(allocated(longxy_obs_floor)) deallocate(longxy_obs_floor) allocate(longxy_obs_floor(dim_obs), stat=ier) if(allocated(latixy_obs_floor)) deallocate(latixy_obs_floor) allocate(latixy_obs_floor(dim_obs), stat=ier) do i = 1, dim_obs if(((lon_clmobs(i) + 180) - minlon) /= 0 .and. ((lat_clmobs(i) + 90) - minlat) /= 0) then longxy_obs_floor(i) = floor(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon)) !+ 1 latixy_obs_floor(i) = floor(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat)) !+ 1 !print *,'longxy_obs(i) , latixy_obs(i) ', longxy_obs(i) , latixy_obs(i) else if(((lon_clmobs(i) + 180) - minlon) == 0 .and. ((lat_clmobs(i) + 90) - minlat) == 0) then longxy_obs_floor(i) = 1 latixy_obs_floor(i) = 1 else if(((lon_clmobs(i) + 180) - minlon) == 0) then longxy_obs_floor(i) = 1 latixy_obs_floor(i) = floor(((lat_clmobs(i) + 90) - minlat) * nj / (maxlat - minlat)) else if(((lat_clmobs(i) + 90) - minlat) == 0) then longxy_obs_floor(i) = floor(((lon_clmobs(i) + 180) - minlon) * ni / (maxlon - minlon)) latixy_obs_floor(i) = 1 endif end do end subroutine get_interp_idx #if defined CLMSA subroutine init_clm_l_size(dim_l) use clm_varpar , only : nlevsoi implicit none integer, intent(out) :: dim_l integer :: nshift if(clmupdate_swc.eq.1) then dim_l = nlevsoi nshift = nlevsoi endif if(clmupdate_swc.eq.2) then error stop "Not implemented: clmupdate_swc.eq.2" ! dim_l = nlevsoi + 1 ! nshift = nlevsoi + 1 endif if(clmupdate_texture.eq.1) then dim_l = 2*nlevsoi + nshift endif if(clmupdate_texture.eq.2) then dim_l = 3*nlevsoi + nshift endif end subroutine init_clm_l_size #endif end module enkf_clm_mod