!------------------------------------------------------------------------------------------- !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 use shr_orb_mod ! shr_orb_params, SHR_ORB_UNDEF_REAL use clm_varorb , only : eccen, mvelpp, lambm0, obliqr, obliq, & iyear_AD, nmvelp use clm_comp , only : clm_init0, clm_init1, clm_init2, clm_run1, clm_run2 use clm_time_manager, only : is_last_step, advance_timestep, get_nstep use atmdrvMod , only : atmdrv, atmdrv_init use abortutils , only : endrun use controlMod , only : control_setNL use clm_mct_mod ! mct_world_init use spmdMod ! mpicom, comp_id, masterproc, spmd_init use ESMF_Mod ! ESMF_Initialize() use perf_mod #if ((defined COUP_OAS_COS || defined COUP_OAS_PFL) && (!defined CLMSA)) use oas_clm_vardef , only : kl_comm #endif ! ! !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(:) 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(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="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 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 implicit none integer,intent(in) :: mype 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_swc.eq.1) then if(clmstatevec_allcol.eq.1) then error stop "Not implemented: clmstatevec_allcol.ne.0" else ! One value per grid-cell clm_varsize = (endg-begg+1) * nlevsoi clm_statevecsize = (endg-begg+1) * nlevsoi end if endif if(clmupdate_swc.eq.2) then clm_varsize = (endg-begg+1) * nlevsoi clm_statevecsize = (endg-begg+1) * (nlevsoi+1) endif if(clmupdate_texture.eq.1) then clm_statevecsize = clm_statevecsize + 2*((endg-begg+1)*nlevsoi) endif if(clmupdate_texture.eq.2) then error stop "Not implemented: clmupdate_texture.eq.2" endif !hcp LST DA if(clmupdate_T.eq.1) then clm_varsize = endg-begg+1 clm_paramsize = endg-begg+1 !LAI clm_statevecsize = (endg-begg+1)*2 !TG, then TV 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)) allocate(state_pdaf2clm_c_p(clm_statevecsize)) allocate(state_pdaf2clm_j_p(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 allocate(clm_paramarr(clm_paramsize)) end if end subroutine define_clm_statevec subroutine set_clm_statevec(tstartcycle, mype) USE clmtype , only : clm3 USE clm_varpar , only : nlevsoi 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 :: tvege(:) !hcp: CLM vegetation temperature (Kelvin) real(r8), pointer :: tgrou(:) !hcp: CLM ground temperature (Kelvin) real(r8), pointer :: ttlai(:) !hcp: CLM one-sided leaf area index, no burying by snow integer :: i,j,cc=1,offset=0 real(r8) :: swcave real(r8) :: swctmp(10) character (len = 34) :: fn !TSMP-PDAF: function name for state vector output character (len = 34) :: fn2 !TSMP-PDAF: function name for swc output swc => clm3%g%l%c%cws%h2osoi_vol psand => clm3%g%l%c%cps%psand pclay => clm3%g%l%c%cps%pclay tvege => clm3%g%l%c%p%pes%t_veg !hcp tgrou => clm3%g%l%c%ces%t_grnd !hcp ttlai => clm3%g%l%c%p%pps%tlai !hcp #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 offset = clm_endg-clm_begg+1 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 error stop "Not implemented: clmstatevec_allcol.ne.0" else do j=clm_begg,clm_endg ! SWC from the first column of each gridcell clm_statevec(cc+offset) = swc(j,i) state_pdaf2clm_c_p(cc+offset) = j state_pdaf2clm_j_p(cc+offset) = i cc = cc + 1 end do end if end do endif !hcp LAI if(clmupdate_T.eq.1) then cc = 1 do j=clm_begg,clm_endg clm_statevec(cc) = tgrou(j) clm_statevec(cc+clm_varsize) = tvege(j) clm_paramarr(cc) = ttlai(j) cc = cc + 1 end do write(*,*) 'b4 update, tgrou(beg) tvege(beg) ttlai(beg)=',tgrou(clm_begg), tvege(clm_begg), ttlai(clm_begg) endif !end hcp LAI ! write average swc to state vector (CRP assimilation) if(clmupdate_swc.eq.2) then cc = 1 do j=clm_begg,clm_endg do i=1,nlevsoi swctmp(i) = swc(j,i) end do call average_swc_crp(swctmp,swcave) clm_statevec(cc) = swcave cc = cc + 1 end do 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 error stop "Not implemented: clmupdate_texture.eq.2" 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 clmtype , only : clm3 use clm_varpar , only : nlevsoi use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varcon , only : denh2o,denice 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 :: tvege(:) !hcp real(r8), pointer :: tgrou(:) !hcp 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 = 0.01_r8 ! minimum soil moisture as in CLM5.0 (mm) real(r8) :: watmin_check ! minimum soil moisture for checking clm_statevec (mm) real(r8) :: watmin_set ! minimum soil moisture for setting swc (mm) integer :: i,j,cc=1,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 => clm3%g%l%c%cws%h2osoi_vol watsat => clm3%g%l%c%cps%watsat psand => clm3%g%l%c%cps%psand pclay => clm3%g%l%c%cps%pclay tvege => clm3%g%l%c%p%pes%t_veg !hcp tgrou => clm3%g%l%c%ces%t_grnd !hcp dz => clm3%g%l%c%cps%dz h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice #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 offset = clm_endg-clm_begg+1 endif ! 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 = 1 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 ! 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 error stop "Not implemented: clmstatevec_allcol.ne.0" else cc = (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. ! 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(clm_statevec(cc+offset).le.watmin_check) then swc(j,i) = watmin_set else if(clm_statevec(cc+offset).ge.watsat(j,i)) then swc(j,i) = watsat(j,i) else swc(j,i) = clm_statevec(cc+offset) 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 ! 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 cc = 1 do j=clm_begg,clm_endg tgrou(j) = clm_statevec(cc) tvege(j) = clm_statevec(cc+clm_varsize) cc = cc + 1 end do write(*,*) 'After update, tgrou(beg) tvege(beg)=',tgrou(clm_begg), tvege(clm_begg) 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.eq.1) 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) 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 clmtype use clm_varpar , only : nlevsoi use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer :: c,lev real(r8) :: clay,sand,ttot real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) psand => clm3%g%l%c%cps%psand pclay => clm3%g%l%c%cps%pclay 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 end do end do end subroutine clm_correct_texture subroutine clm_texture_to_parameters() use clmtype use clm_varpar , only : nlevsoi use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer :: c,lev real(r8) :: clay,sand ! temporaries real(r8) :: bd ! bulk density of dry soil material [kg/m^3] real(r8) :: xksat ! maximum hydraulic conductivity of soil [mm/s] real(r8) :: tkm ! mineral conductivity real(r8), pointer :: watsat(:,:) ! volumetric soil water at saturation (porosity) (nlevsoi) real(r8), pointer :: bsw(:,:) ! Clapp and Hornberger "b" (nlevsoi) real(r8), pointer :: bsw2(:,:) ! Clapp and Hornberger "b" for CN code real(r8), pointer :: psisat(:,:) ! soil water potential at saturation for CN code (MPa) real(r8), pointer :: vwcsat(:,:) ! volumetric water content at saturation for CN code (m3/m3) real(r8), pointer :: hksat(:,:) ! hydraulic conductivity at saturation (mm H2O /s) (nlevsoi) real(r8), pointer :: sucsat(:,:) ! minimum soil suction (mm) (nlevsoi) real(r8), pointer :: tkmg(:,:) ! thermal conductivity, soil minerals [W/m-K] (new) (nlevsoi) real(r8), pointer :: tksatu(:,:) ! thermal conductivity, saturated soil [W/m-K] (new) (nlevsoi) real(r8), pointer :: tkdry(:,:) ! thermal conductivity, dry soil (W/m/Kelvin) (nlevsoi) real(r8), pointer :: csol(:,:) ! heat capacity, soil solids (J/m**3/Kelvin) (nlevsoi) real(r8), pointer :: watdry(:,:) ! btran parameter for btran=0 real(r8), pointer :: watopt(:,:) ! btran parameter for btran = 1 real(r8), pointer :: psand(:,:) real(r8), pointer :: pclay(:,:) watsat => clm3%g%l%c%cps%watsat tkmg => clm3%g%l%c%cps%tkmg bsw => clm3%g%l%c%cps%bsw bsw2 => clm3%g%l%c%cps%bsw2 psisat => clm3%g%l%c%cps%psisat vwcsat => clm3%g%l%c%cps%vwcsat hksat => clm3%g%l%c%cps%hksat sucsat => clm3%g%l%c%cps%sucsat tkmg => clm3%g%l%c%cps%tkmg tksatu => clm3%g%l%c%cps%tksatu tkdry => clm3%g%l%c%cps%tkdry csol => clm3%g%l%c%cps%csol watdry => clm3%g%l%c%cps%watdry watopt => clm3%g%l%c%cps%watopt psand => clm3%g%l%c%cps%psand pclay => clm3%g%l%c%cps%pclay do c = clm_begg,clm_endg do lev = 1,nlevsoi clay = pclay(c,lev) sand = psand(c,lev) watsat(c,lev) = 0.489_r8 - 0.00126_r8*sand bd = (1._r8-watsat(c,lev))*2.7e3_r8 xksat = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand) ) ! mm/s tkm = (8.80_r8*sand+2.92_r8*clay)/(sand+clay) ! W/(m K) bsw(c,lev) = 2.91_r8 + 0.159_r8*clay bsw2(c,lev) = -(3.10_r8 + 0.157_r8*clay - 0.003_r8*sand) psisat(c,lev) = -(exp((1.54_r8 - 0.0095_r8*sand + 0.0063_r8*(100.0_r8-sand-clay))*log(10.0_r8))*9.8e-5_r8) vwcsat(c,lev) = (50.5_r8 - 0.142_r8*sand - 0.037_r8*clay)/100.0_r8 hksat(c,lev) = xksat sucsat(c,lev) = 10._r8 * ( 10._r8**(1.88_r8-0.0131_r8*sand) ) tkmg(c,lev) = tkm ** (1._r8- watsat(c,lev)) tksatu(c,lev) = tkmg(c,lev)*0.57_r8**watsat(c,lev) tkdry(c,lev) = (0.135_r8*bd + 64.7_r8) / (2.7e3_r8 - 0.947_r8*bd) csol(c,lev) = (2.128_r8*sand+2.385_r8*clay) / (sand+clay)*1.e6_r8 ! J/(m3 K) watdry(c,lev) = watsat(c,lev) * (316230._r8/sucsat(c,lev)) ** (-1._r8/bsw(c,lev)) watopt(c,lev) = watsat(c,lev) * (158490._r8/sucsat(c,lev)) ** (-1._r8/bsw(c,lev)) 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 real(r8) :: w(10) real(r8) :: mold=0,mnew=0,delta=1,g,wtmp=0 integer :: iter=0,i,j w = 1/10 do while(delta>0.0001) ! calculate initial weighted mean mold = 0 do i=1,10 mold = mold + w(i)*profdat(i) end do ! calculate g factor g = -5.8/(log(0.14)*(mold+0.0829)) ! calculate new weights w(1) = 1-exp(-zsoi(1)*100/g) do i=2,10 ! calculate sum of previous weights wtmp = 0 do j=1,(i-1) wtmp = wtmp + w(j) end do w(i) = 1-exp(-zsoi(i)*100/g)-wtmp end do ! calculate new weighted mean mnew = 0 do i=1,10 mnew = mnew + w(i)*profdat(i) end do ! compare old and new weighted mean iter=iter+1 delta = abs(mnew-mold) end do profave = mnew end subroutine average_swc_crp #endif subroutine domain_def_clm(lon_clmobs, lat_clmobs, dim_obs, & longxy, latixy, longxy_obs, latixy_obs) !#if defined CLMSA USE domainMod, ONLY: alatlon 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(:) 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 :: npfts real :: minlon, minlat, maxlon, maxlat real(r8), pointer :: lon(:) real(r8), pointer :: lat(:) integer :: begg, endg ! per-proc gridcell ending gridcell indices lon => alatlon%lonc lat => alatlon%latc ni = alatlon%ni nj = alatlon%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) !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 if(allocated(longxy)) deallocate(longxy) allocate(longxy(ncells), stat=ier) if(allocated(latixy)) deallocate(latixy) 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 == adecomp%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) if(allocated(longxy_obs)) deallocate(longxy_obs) allocate(longxy_obs(dim_obs), stat=ier) if(allocated(latixy_obs)) deallocate(latixy_obs) 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: alatlon ! 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 => alatlon%lonc lat => alatlon%latc ni = alatlon%ni nj = alatlon%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 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 error stop "Not implemented: clmupdate_texture.eq.2" endif end subroutine init_clm_l_size #endif end module enkf_clm_mod