module SoilMoistureStreamMod ! ********************************************************************** ! --------------------------- IMPORTANT NOTE --------------------------- ! ! This file is here temporarily in order to exercise the CDEPS stream code for this 3-d ! stream. In cases using the NUOPC driver/mediator, this version is used instead of the ! version in src/biogeophys. Changes to the science here should also be made in the ! similar file in src/biogeophys. Once we start using CDEPS by default, we can remove ! the version in src/biogeophys and move this version into there. ! ********************************************************************** #include "shr_assert.h" !----------------------------------------------------------------------- ! !DESCRIPTION: ! Read in soil moisture from data stream ! ! !USES: use ESMF use dshr_strdata_mod , only : shr_strdata_type, shr_strdata_print use dshr_strdata_mod , only : shr_strdata_init_from_inline, shr_strdata_advance use dshr_methods_mod , only : dshr_fldbun_getfldptr use dshr_stream_mod , only : shr_stream_file_null use shr_kind_mod , only : r8 => shr_kind_r8, cl => shr_kind_CL, cxx => shr_kind_CXX use shr_log_mod , only : errMsg => shr_log_errMsg use shr_mpi_mod , only : shr_mpi_bcast use decompMod , only : bounds_type use abortutils , only : endrun use clm_varctl , only : iulog, use_soil_moisture_streams use controlMod , only : NLFilename use LandunitType , only : lun use ColumnType , only : col use SoilStateType , only : soilstate_type use WaterStateBulkType , only : waterstatebulk_type use perf_mod , only : t_startf, t_stopf use spmdMod , only : masterproc, mpicom, iam use clm_time_manager , only : get_calendar, get_curr_date use clm_nlUtilsMod , only : find_nlgroup_name use clm_varpar , only : nlevsoi use clm_varcon , only : denh2o, denice, watmin, spval use landunit_varcon , only : istsoil, istcrop use lnd_comp_shr , only : mesh, model_meshfile, model_clock ! ! !PUBLIC TYPES: implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: public :: PrescribedSoilMoistureInit ! position datasets for soil moisture public :: PrescribedSoilMoistureAdvance ! Advance the soil moisture stream (outside of Open-MP loops) public :: PrescribedSoilMoistureInterp ! interpolates between two periods of soil moisture data ! !PRIVATE MEMBER DATA: type(shr_strdata_type) :: sdat_soilm ! soil moisture input data stream character(*), parameter :: stream_var_name = "H2OSOI" ! base string for field string character(len=CL) :: stream_lev_dimname = 'levsoi' ! name of vertical layer dimension integer, allocatable :: g_to_ig(:) ! Array matching gridcell index to data index logical :: soilm_ignore_data_if_missing ! If should ignore overridding a point with soil moisture data ! from the streams file, if the streams file shows that point ! as missing (namelist item) ! !PRIVATE TYPES: character(len=*), parameter, private :: sourcefile = & __FILE__ character(*), parameter :: u_FILE_u = & __FILE__ !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- ! ! soil_moisture_init ! !----------------------------------------------------------------------- subroutine PrescribedSoilMoistureInit(bounds) ! ! Initialize data stream information for soil moisture. ! ! !USES: ! ! !ARGUMENTS: implicit none type(bounds_type), intent(in) :: bounds ! bounds ! ! !LOCAL VARIABLES: integer :: rc ! error coe integer :: i ! index integer :: stream_year_first_soilm ! first year in Ustar stream to use integer :: stream_year_last_soilm ! last year in Ustar stream to use integer :: model_year_align_soilm ! align stream_year_first_soilm with integer :: nu_nml ! unit for namelist file integer :: nml_error ! namelist i/o error flag integer :: soilm_offset ! Offset in time for dataset (sec) character(len=CL) :: stream_fldfilename_soilm ! ustar stream filename to read character(len=CL) :: soilm_tintalgo = 'linear' ! Time interpolation alogrithm character(len=CL) :: stream_mapalgo = 'bilinear' real(r8) :: stream_dtlimit = 15._r8 character(len=CL) :: stream_taxmode = 'cycle' character(*), parameter :: subName = "('PrescribedSoilMoistureInit')" character(*), parameter :: F00 = "('(PrescribedSoilMoistureInit) ',4a)" !----------------------------------------------------------------------- ! ! deal with namelist variables here in init ! namelist /soil_moisture_streams/ & stream_year_first_soilm, & stream_year_last_soilm, & model_year_align_soilm, & soilm_tintalgo, & soilm_offset, & soilm_ignore_data_if_missing, & stream_fldfilename_soilm ! Default values for namelist stream_year_first_soilm = 1 ! first year in stream to use stream_year_last_soilm = 1 ! last year in stream to use model_year_align_soilm = 1 ! align stream_year_first_soilm with this model year stream_fldfilename_soilm = shr_stream_file_null soilm_offset = 0 soilm_ignore_data_if_missing = .false. ! Read soilm_streams namelist if (masterproc) then open( newunit=nu_nml, file=trim(NLFilename), status='old', iostat=nml_error ) call find_nlgroup_name(nu_nml, 'soil_moisture_streams', status=nml_error) if (nml_error == 0) then read(nu_nml, nml=soil_moisture_streams,iostat=nml_error) if (nml_error /= 0) then call endrun(subname // ':: ERROR reading soil_moisture_streams namelist') end if else call endrun(subname // ':: ERROR finding soilm_streams namelist') end if close(nu_nml) endif call shr_mpi_bcast(stream_year_first_soilm, mpicom) call shr_mpi_bcast(stream_year_last_soilm, mpicom) call shr_mpi_bcast(model_year_align_soilm, mpicom) call shr_mpi_bcast(stream_fldfilename_soilm, mpicom) call shr_mpi_bcast(soilm_tintalgo, mpicom) call shr_mpi_bcast(soilm_offset, mpicom) call shr_mpi_bcast(soilm_ignore_data_if_missing, mpicom) if (masterproc) then write(iulog,*) ' ' write(iulog,*) 'soil_moisture_stream settings:' write(iulog,*) ' stream_year_first_soilm = ',stream_year_first_soilm write(iulog,*) ' stream_year_last_soilm = ',stream_year_last_soilm write(iulog,*) ' model_year_align_soilm = ',model_year_align_soilm write(iulog,*) ' stream_fldfilename_soilm = ',trim(stream_fldfilename_soilm) write(iulog,*) ' soilm_tintalgo = ',trim(soilm_tintalgo) write(iulog,*) ' soilm_offset = ',soilm_offset if ( soilm_ignore_data_if_missing ) then write(iulog,*) ' Do NOT override a point with streams data if the streams data is missing' else write(iulog,*) ' Abort, if you find a model point where the input streams data is set to missing value' end if endif ! Initialize the cdeps data type sdat_soilm ! TODO: for now stream_meshfile is the same as the model meshfile - must generalize this if want to have ! stream be at a different resolution call shr_strdata_init_from_inline(sdat_soilm, & my_task = iam, & logunit = iulog, & compname = 'LND', & model_clock = model_clock, & model_mesh = mesh, & stream_meshfile = model_meshfile, & stream_lev_dimname = trim(stream_lev_dimname), & stream_mapalgo = trim(stream_mapalgo), & stream_filenames = (/trim(stream_fldfilename_soilm)/), & stream_fldlistFile = (/trim(stream_var_name)/), & stream_fldListModel = (/trim(stream_var_name)/), & stream_yearFirst = stream_year_first_soilm, & stream_yearLast = stream_year_last_soilm, & stream_yearAlign = model_year_align_soilm, & stream_offset = soilm_offset, & stream_taxmode = stream_taxmode, & stream_dtlimit = stream_dtlimit, & stream_tintalgo = soilm_tintalgo, & rc = rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then call ESMF_Finalize(endflag=ESMF_END_ABORT) end if if (masterproc) then call shr_strdata_print(sdat_soilm,'soil moisture data') endif end subroutine PrescribedSoilMoistureInit !----------------------------------------------------------------------- ! ! PrescribedSoilMoistureAdvance ! !----------------------------------------------------------------------- subroutine PrescribedSoilMoistureAdvance( bounds ) ! ! Advanace the prescribed soil moisture stream ! ! !USES: ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds ! ! !LOCAL VARIABLES: integer :: rc integer :: g, ig integer :: ier ! error code integer :: year ! year (0, ...) for nstep+1 integer :: mon ! month (1, ..., 12) for nstep+1 integer :: day ! day of month (1, ..., 31) for nstep+1 integer :: sec ! seconds into current date for nstep+1 integer :: mcdate ! Current model date (yyyymmdd) !----------------------------------------------------------------------- rc = ESMF_SUCCESS call get_curr_date(year, mon, day, sec) mcdate = year*10000 + mon*100 + day ! Advance sdat stream call shr_strdata_advance(sdat_soilm, ymd=mcdate, tod=sec, logunit=iulog, istr='soil_moisture', rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then call ESMF_Finalize(endflag=ESMF_END_ABORT) end if ! Map gridcell to 1->local_size (g_to_ig is a module variable) ier = 0 if ( .not. allocated(g_to_ig) )then allocate (g_to_ig(bounds%begg:bounds%endg), stat=ier) if (ier /= 0) then write(iulog,*) 'Prescribed soil moisture allocation error' call endrun(msg=errMsg(sourcefile, __LINE__)) end if ig = 0 do g = bounds%begg,bounds%endg ig = ig+1 g_to_ig(g) = ig end do end if end subroutine PrescribedSoilMoistureAdvance !----------------------------------------------------------------------- ! ! PrescribedSoilMoistureInterp ! !----------------------------------------------------------------------- subroutine PrescribedSoilMoistureInterp(bounds, soilstate_inst, waterstatebulk_inst) ! ! Assign data stream information for prescribed soil moisture. ! ! !USES: ! ! !ARGUMENTS: implicit none type(bounds_type) , intent(in) :: bounds type(soilstate_type) , intent(in) :: soilstate_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst ! ! !LOCAL VARIABLES: integer :: rc integer :: c, g, j, ig, n real(r8) :: soilm_liq_frac ! liquid fraction of soil moisture real(r8) :: soilm_ice_frac ! ice fraction of soil moisture real(r8) :: moisture_increment ! soil moisture adjustment increment real(r8) :: h2osoi_vol_initial ! initial vwc value real(r8), pointer :: dataptr2d(:,:) ! first dimension is level, second is data on that level character(*), parameter :: subName = "('PrescribedSoilMoistureInterp')" !----------------------------------------------------------------------- SHR_ASSERT_FL( (lbound(g_to_ig,1) <= bounds%begg ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(g_to_ig,1) >= bounds%endg ), sourcefile, __LINE__) associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input/Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input/Output: [real(r8) (:,:) ] ice water (kg/m2) h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Output: volumetric soil water (m3/m3) h2osoi_vol_prs => waterstatebulk_inst%h2osoi_vol_prs_grc & ! Output: prescribed volumetric soil water (m3/m3) ) SHR_ASSERT_FL( (lbound(h2osoi_vol,1) <= bounds%begc ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_vol,1) >= bounds%endc ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(h2osoi_vol,2) == 1 ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_vol,2) >= nlevsoi ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(dz,1) <= bounds%begc ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(dz,1) >= bounds%endc ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(dz,2) <= 1 ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(dz,2) >= nlevsoi ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(watsat,1) <= bounds%begc ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(watsat,1) >= bounds%endc ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(watsat,2) <= 1 ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(watsat,2) >= nlevsoi ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(h2osoi_liq,1) <= bounds%begc ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_liq,1) >= bounds%endc ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(h2osoi_liq,2) <= 1 ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_liq,2) >= nlevsoi ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(h2osoi_ice,1) <= bounds%begc ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_ice,1) >= bounds%endc ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(h2osoi_ice,2) <= 1 ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_ice,2) >= nlevsoi ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(h2osoi_vol_prs,1) <= bounds%begg ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_vol_prs,1) >= bounds%endg ), sourcefile, __LINE__) SHR_ASSERT_FL( (lbound(h2osoi_vol_prs,2) == 1 ), sourcefile, __LINE__) SHR_ASSERT_FL( (ubound(h2osoi_vol_prs,2) >= nlevsoi ), sourcefile, __LINE__) ! Get pointer for stream data that is time and spatially interpolate to model time and grid call dshr_fldbun_getFldPtr(sdat_soilm%pstrm(1)%fldbun_model, trim(stream_var_name), fldptr2=dataptr2d, rc=rc) if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) then call ESMF_Finalize(endflag=ESMF_END_ABORT) end if ! Check that inner most dimensino of dataptr2d is equal to nlevsoi if (size(dataptr2d,dim=1) /= nlevsoi) then if (masterproc) then write(iulog,*) 'ERROR: dataptr2d(dim=1) = ',size(dataptr2d,dim=1),& ' and nlevsoi = ',nlevsoi,' must match ' end if call endrun(trim(subname) // & ' ERROR:: The input soil moisture stream does not have levels equal to nlevsoi') end if ! Set the prescribed soil moisture read from the file everywhere do g = bounds%begg, bounds%endg ig = g_to_ig(g) do j = 1, nlevsoi h2osoi_vol_prs(g,j) = dataptr2d(j,ig) ! If soil moiture is being interpolated in time and the result is ! large that probably means one of the two data points is missing (set to spval) if ( h2osoi_vol_prs(g,j) > 10.0_r8 .and. (h2osoi_vol_prs(g,j) /= spval) )then h2osoi_vol_prs(g,j) = spval end if end do end do do c = bounds%begc, bounds%endc ! Set variable for each gridcell/column combination g = col%gridcell(c) ig = g_to_ig(g) ! EBK Jan/2020, also check weights on gridcell (See https://github.com/ESCOMP/CTSM/issues/847) if ( (lun%itype(col%landunit(c)) == istsoil) .or. & (lun%itype(col%landunit(c)) == istcrop) .and. (col%wtgcell(c) /= 0._r8) ) then ! this is a 2d field (gridcell/nlevsoi) ! do j = 1, nlevsoi ! if soil water is zero, liq/ice fractions cannot be calculated if((h2osoi_liq(c, j) + h2osoi_ice(c, j)) > 0._r8) then ! save original soil moisture value h2osoi_vol_initial = h2osoi_vol(c,j) ! Check if the vegetated land mask from the dataset on the ! file is different if ( (h2osoi_vol_prs(g,j) == spval) .and. (h2osoi_vol_initial /= spval) )then if ( soilm_ignore_data_if_missing )then cycle else write(iulog,*) 'Input soil moisture dataset is not vegetated as expected: gridcell=', & g, ' active = ', col%active(c) call endrun(subname //& ' ERROR:: The input soil moisture stream is NOT vegetated for one of the land points' ) end if end if ! update volumetric soil moisture from data prescribed from the file h2osoi_vol(c,j) = h2osoi_vol_prs(g,j) ! calculate liq/ice mass fractions soilm_liq_frac = h2osoi_liq(c, j) /(h2osoi_liq(c, j) + h2osoi_ice(c, j)) soilm_ice_frac = h2osoi_ice(c, j) /(h2osoi_liq(c, j) + h2osoi_ice(c, j)) ! calculate moisture increment moisture_increment = h2osoi_vol(c,j) - h2osoi_vol_initial ! add limitation check moisture_increment = min((watsat(c,j) - h2osoi_vol_initial),max(-(h2osoi_vol_initial-watmin),moisture_increment)) ! update liq/ice water mass due to (volumetric) moisture increment h2osoi_liq(c,j) = h2osoi_liq(c,j) + (soilm_liq_frac * moisture_increment * dz(c, j) * denh2o) h2osoi_ice(c,j) = h2osoi_ice(c,j) + (soilm_ice_frac * moisture_increment * dz(c, j) * denice) else call endrun(subname // ':: ERROR h2osoil liquid plus ice is zero') endif enddo endif end do end associate end subroutine PrescribedSoilMoistureInterp end module SoilMoistureStreamMod