diff --git a/ldt/DAobs/LISlsmPrecipobs/LISlsmPrecip_obsMod.F90 b/ldt/DAobs/LISlsmPrecipobs/LISlsmPrecip_obsMod.F90 new file mode 100644 index 000000000..649bb2443 --- /dev/null +++ b/ldt/DAobs/LISlsmPrecipobs/LISlsmPrecip_obsMod.F90 @@ -0,0 +1,340 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.3 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +module LISlsmPrecip_obsMod +!BOP +! +! !MODULE: LISlsmPrecip_obsMod +! +! !DESCRIPTION: +! This module handles the use of a LIS model simulation +! output as "observations" for data assimilation. This +! plugin is typically used to handle the computations of +! scaling factors such as cumulative distribution function +! (CDF) for use in DA +! +! !REVISION HISTORY: +! 16 Nov 2021 Mahdi Navari Initial Specification (based on LISlsmSM_obsMod) +! + + PRIVATE +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + PUBLIC :: LISlsmPrecip_obsInit +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + PUBLIC :: lsmprecipobs +! +!EOP + + type, public :: lsmprecipobsdec + integer :: nvars + integer :: nest + integer :: nc,nr + real :: datares + real :: run_dd(8) + character*50 :: map_proj + character*50 :: format + character*50 :: wstyle + character*50 :: wopt + character*100 :: odir + character*20 :: security_class + character*20 :: distribution_class + character*20 :: data_category + character*20 :: area_of_data + character*20 :: write_interval +!-------------------------------------------------------- +! interpolation/upscaling weights +!-------------------------------------------------------- + integer, allocatable :: n11(:) + integer, allocatable :: n12(:) + integer, allocatable :: n21(:) + integer, allocatable :: n22(:) + real, allocatable :: w11(:) + real, allocatable :: w12(:) + real, allocatable :: w21(:) + real, allocatable :: w22(:) + + end type lsmprecipobsdec + + type(lsmprecipobsdec) :: lsmprecipobs + +contains + +!BOP +! !ROUTINE: LISlsmPrecip_obsInit +! \label{LISlsmPrecip_obsInit} +! +! !INTERFACE: + subroutine LISlsmPrecip_obsInit() +! !USES: + use ESMF + use LDT_coreMod + use LDT_DAobsDataMod + use LDT_logMod + + implicit none +! +! !DESCRIPTION: +! This routine initializes the structures required for the handling of a +! land surface model output (from a LIS simulation) as observations. +! +!EOP + integer :: n + integer :: rc + real :: gridDesci(20) + + n = 1 + + lsmprecipobs%run_dd = LDT_rc%udef + + lsmprecipobs%security_class = '' + lsmprecipobs%distribution_class = '' + lsmprecipobs%data_category = '' + lsmprecipobs%area_of_data = '' + lsmprecipobs%write_interval = '' + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%format, & + label="LIS precipitation output format:",rc=rc) + call LDT_verify(rc,'LIS precipitation output format: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%wopt, & + label="LIS precipitation output methodology:",rc=rc) + call LDT_verify(rc,'LIS precipitation output methodology: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%wstyle, & + label="LIS precipitation output naming style:",rc=rc) + call LDT_verify(rc,'LIS precipitation output naming style: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%map_proj, & + label="LIS precipitation output map projection:",rc=rc) + call LDT_verify(rc,'LIS precipitation output map projection: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%nest, & + label="LIS precipitation output nest index:",rc=rc) + call LDT_verify(rc,'LIS precipitation output nest index: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%odir, & + label="LIS precipitation output directory:",rc=rc) + call LDT_verify(rc,'LIS precipitation output directory: not defined') + + ! WMO-convention specific identifiers + if ( lsmprecipobs%wstyle == "WMO convention") then + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%security_class, & + label="LIS precipitation security class:",rc=rc) + call LDT_verify(rc,'LIS precipitation security class: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%distribution_class, & + label="LIS precipitation distribution class:",rc=rc) + call LDT_verify(rc,'LIS precipitation distribution class: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%data_category, & + label="LIS precipitation data category:",rc=rc) + call LDT_verify(rc,'LIS precipitation data category: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%area_of_data, & + label="LIS precipitation area of data:",rc=rc) + call LDT_verify(rc,'LIS precipitation area of data: not defined') + + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%write_interval, & + label="LIS precipitation write interval:",rc=rc) + call LDT_verify(rc,'LIS precipitation write interval: not defined') + endif + + if(lsmprecipobs%map_proj.eq."latlon") then + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain lower left lat:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(1),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain lower left lon:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(2),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain upper right lat:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(3),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain upper right lon:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(4),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain resolution (dx):",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(5),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain resolution (dy):",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(6),rc=rc) + + lsmprecipobs%datares = min(lsmprecipobs%run_dd(5),lsmprecipobs%run_dd(6)) + + lsmprecipobs%nc = (nint((lsmprecipobs%run_dd(4)-lsmprecipobs%run_dd(2))/& + lsmprecipobs%run_dd(5))) + 1 + lsmprecipobs%nr = (nint((lsmprecipobs%run_dd(3)-lsmprecipobs%run_dd(1))/& + lsmprecipobs%run_dd(6))) + 1 + + gridDesci = 0 + gridDesci(1) = 0 + gridDesci(2) = lsmprecipobs%nc + gridDesci(3) = lsmprecipobs%nr + gridDesci(4) = lsmprecipobs%run_dd(1) + gridDesci(5) = lsmprecipobs%run_dd(2) + gridDesci(6) = 128 + gridDesci(7) = lsmprecipobs%run_dd(3) + gridDesci(8) = lsmprecipobs%run_dd(4) + gridDesci(9) = lsmprecipobs%run_dd(5) + gridDesci(10) = lsmprecipobs%run_dd(6) + gridDesci(20) = 64 + + elseif(lsmprecipobs%map_proj.eq."lambert") then + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain lower left lat:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(1),rc=rc) + call LDT_verify(rc,'LIS precipitation domain lower left lat: not defined') + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain lower left lon:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(2),rc=rc) + call LDT_verify(rc,'LIS precipitation domain lower left lon: not defined') + + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain true lat1:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(3),rc=rc) + call LDT_verify(rc,'LIS precipitation domain true lat1: not defined') + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain true lat2:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(4),rc=rc) + call LDT_verify(rc,'LIS precipitation domain true lat2: not defined') + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain standard lon:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(5),rc=rc) + call LDT_verify(rc,'LIS precipitation domain standard lon: not defined') + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain resolution:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(6),rc=rc) + call LDT_verify(rc,'LIS precipitation domain resolution: not defined') + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain x-dimension size:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(7),rc=rc) + call LDT_verify(rc,'LIS precipitation domain x-dimension size: not defined') + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain y-dimension size:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(8),rc=rc) + call LDT_verify(rc,'LIS precipitation domain y-dimension size: not defined') + + lsmprecipobs%datares = lsmprecipobs%run_dd(6)/100.0 + + lsmprecipobs%nc = lsmprecipobs%run_dd(7) + lsmprecipobs%nr = lsmprecipobs%run_dd(8) + + gridDesci = 0 + gridDesci(1) = 3 + gridDesci(2) = lsmprecipobs%nc + gridDesci(3) = lsmprecipobs%nr + gridDesci(4) = lsmprecipobs%run_dd(1) + gridDesci(5) = lsmprecipobs%run_dd(2) + gridDesci(6) = 8 + gridDesci(7) = lsmprecipobs%run_dd(4) + gridDesci(8) = lsmprecipobs%run_dd(6) + gridDesci(9) = lsmprecipobs%run_dd(6) + gridDesci(10) = lsmprecipobs%run_dd(3) + gridDesci(11) = lsmprecipobs%run_dd(5) + gridDesci(20) = 64 + + elseif(lsmprecipobs%map_proj.eq."polar") then + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain lower left lat:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(1),rc=rc) + + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain lower left lon:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(2),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain true lat1:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(3),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain true lat2:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(4),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain standard lon:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(5),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain resolution:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(6),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain x-dimension size:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(7),rc=rc) + + call ESMF_ConfigFindLabel(LDT_config,& + "LIS precipitation domain y-dimension size:",rc=rc) + call ESMF_ConfigGetAttribute(LDT_config,lsmprecipobs%run_dd(8),rc=rc) + + endif + +!------------------------------------------------------------------- +! if the LIS output (obs) is at a coarser resolution than the +! LDT grid, then setup the weights for interpolation. Else +! setup the weights for upscaling. +!------------------------------------------------------------------- + if(LDT_isLDTatAfinerResolution(n,lsmprecipobs%datares)) then + + allocate(lsmprecipobs%n11(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(lsmprecipobs%n12(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(lsmprecipobs%n21(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(lsmprecipobs%n22(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + + allocate(lsmprecipobs%w11(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(lsmprecipobs%w12(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(lsmprecipobs%w21(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + allocate(lsmprecipobs%w22(LDT_rc%lnc(n)*LDT_rc%lnr(n))) + + call bilinear_interp_input(n, gridDesci, & + lsmprecipobs%n11, & + lsmprecipobs%n12, lsmprecipobs%n21, & + lsmprecipobs%n22, lsmprecipobs%w11, & + lsmprecipobs%w12, lsmprecipobs%w21, & + lsmprecipobs%w22) + + else + + allocate(lsmprecipobs%n11(& + lsmprecipobs%nc*& + lsmprecipobs%nr)) + + call upscaleByAveraging_input(& + gridDesci,& + LDT_rc%gridDesc(n,:),& + lsmprecipobs%nc*lsmprecipobs%nr,& + LDT_rc%lnc(n)*LDT_rc%lnr(n),& + lsmprecipobs%n11) + + endif + +! which variable we want in the DA obs computations. + call LDT_initializeDAobsEntry(LDT_DAobsData(1)%totalprecip_obs, "kg m-2",1,1) + LDT_DAobsData(1)%totalprecip_obs%selectStats = 1 + + end subroutine LISlsmPrecip_obsInit + +end module LISlsmPrecip_obsMod diff --git a/ldt/DAobs/LISlsmPrecipobs/readLISlsmPrecipobs.F90 b/ldt/DAobs/LISlsmPrecipobs/readLISlsmPrecipobs.F90 new file mode 100644 index 000000000..109685ceb --- /dev/null +++ b/ldt/DAobs/LISlsmPrecipobs/readLISlsmPrecipobs.F90 @@ -0,0 +1,579 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.3 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LDT_misc.h" +!BOP +! +! !ROUTINE: readLISlsmPrecipobs +! \label{readLISlsmPrecipobs} +! +! !INTERFACE: +! !REVISION HISTORY: +! 2Dec2021: Mahdi Navari ; Initial Specification (based on readLISlsmSM) +subroutine readLISlsmPrecipobs(n) +! !USES: +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif +#if (defined USE_GRIBAPI) + use grib_api +#endif + use LDT_coreMod + use LDT_DAobsDataMod + use LDT_historyMod + use LDT_logMod + use LISlsmPrecip_obsMod, only : lsmprecipobs +! +! !DESCRIPTION: +! This routine reads the total precipitation fields from a LIS model +! simulation. +! +!EOP + implicit none + + integer, intent(in) :: n + + character*200 :: fname + logical :: file_exists + real :: precip_data(LDT_rc%lnc(n),LDT_rc%lnr(n)) + + integer :: t, index + integer :: ftn + integer :: iret + real :: topLev,botLev,param_num,trange + integer :: igrib,nvars + real :: precipvalue1d(lsmprecipobs%nc*lsmprecipobs%nr) + real :: precipvalue2d(lsmprecipobs%nc, lsmprecipobs%nr) + integer :: c,r + character*20 :: vname + integer :: varid + + interface + subroutine create_lsm_output_fname(n, form, fname, odir, wstyle, wopt, & + run_dd, map_proj, security_class, & + distribution_class, data_category, & + area_of_data, write_interval) + integer, intent(IN) :: n + character(len=*) :: fname + character(len=*) :: form + character(len=*) :: odir + character(len=*) :: wstyle + character(len=*) :: wopt + real, dimension(8), optional :: run_dd + character(len=*), optional :: map_proj + character(len=*), optional :: security_class + character(len=*), optional :: distribution_class + character(len=*), optional :: data_category + character(len=*), optional :: area_of_data + character(len=*), optional :: write_interval + end subroutine create_lsm_output_fname + end interface + + +#if (defined USE_GRIBAPI) + precip_data = LDT_rc%udef + + call create_lsm_output_fname(lsmprecipobs%nest, & + lsmprecipobs%format, & + fname, & + lsmprecipobs%odir, & + lsmprecipobs%wstyle, & + lsmprecipobs%wopt, & + lsmprecipobs%run_dd, & + lsmprecipobs%map_proj, & + lsmprecipobs%security_class, & + lsmprecipobs%distribution_class, & + lsmprecipobs%data_category, & + lsmprecipobs%area_of_data, & + lsmprecipobs%write_interval) + + inquire(file=trim(fname),exist=file_exists) + + if(file_exists) then + write(LDT_logunit,*) '[INFO] reading LSM output ',trim(fname) + if(lsmprecipobs%format.eq."binary") then + write(LDT_logunit,*) '[ERR] DA preprocessing on the binary format is not ' + write(LDT_logunit,*) '[ERR] currently supported. Program stopping....' + call LDT_endrun() + + elseif(lsmprecipobs%format.eq."grib1") then + if(lsmprecipobs%wstyle.ne."WMO convention") then + write(LDT_logunit,*) '[ERR] LDT currently does not support this style of grib output' + call LDT_endrun + endif + + call grib_open_file(ftn,trim(fname),'r',iret) + if(iret.ne.0) then + write(LDT_logunit,*) '[ERR] Could not open file: ',trim(fname) + call LDT_endrun() + endif + call grib_multi_support_on + + do + call grib_new_from_file(ftn,igrib,iret) + if(iret==GRIB_END_OF_FILE) then + exit + endif + + call grib_get(igrib, 'indicatorOfParameter',param_num, iret) + call LDT_verify(iret, & + 'grib_get: indicatorOfParameter failed in readLISlsmPrecipObs') + + call grib_get(igrib, 'topLevel',topLev, iret) + call LDT_verify(iret, & + 'grib_get: topLevel failed in readLISlsmPrecipObs') + + call grib_get(igrib, 'bottomLevel',botLev, iret) + call LDT_verify(iret, & + 'grib_get: bottomLevel failed in readLISlsmPrecipObs') + + call grib_get(igrib, 'timeRangeIndicator',trange, iret) + call LDT_verify(iret, & + 'grib_get: timeRangeIndicator failed in readLISlsmPrecipObs') + +!right now specifically geared for AFWA outputs. + if(param_num.eq.201.and.topLev.eq.0.and.botlev.eq.10.and.& + trange.eq.1) then + + call grib_get(igrib,'values',precipvalue1d,iret) + call LDT_verify(iret,& + 'grib_get: values failed in readLISlsmSMobs') + + call transformPrecipDataToLDTgrid(n,precipvalue1d,precip_data) + + endif + + call grib_release(igrib,iret) + call LDT_verify(iret, 'error in grib_release in readLISlsmSMObs') + enddo + + call grib_close_file(ftn) + + elseif(lsmprecipobs%format.eq."netcdf") then +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + + iret = nf90_open(path=trim(fname),mode=nf90_nowrite, ncid=ftn) + call LDT_verify(iret, 'Error opening file '//trim(fname)) + +! The code looks for instantaneous variables. If it doesn't exist, +! the time averaged data fields will be read in. + + iret = nf90_inq_varid(ftn, 'TotalPrecip_inst', varid) + vname = 'TotalPrecip_inst' + if(iret.ne.0) then + vname = 'TotalPrecip_tavg' + endif + + if(lsmprecipobs%datares.eq.LDT_rc%gridDesc(n,10)) then + call LDT_readLISSingleNetcdfVar(n,ftn, vname,& + 1,lsmprecipobs%nc, lsmprecipobs%nr, precipvalue2d) + else + call LDT_readLISSingleNetcdfVar(n,ftn, vname,& + 1,lsmprecipobs%nc, lsmprecipobs%nr, precipvalue2d) + endif + + iret = nf90_close(ftn) + call LDT_verify(iret,'Error in nf90_close') + + do r=1,lsmprecipobs%nr + do c=1, lsmprecipobs%nc + precipvalue1d(c+(r-1)*lsmprecipobs%nc) = precipvalue2d(c,r) + enddo + enddo + + call transformPrecipDataToLDTgrid(n,precipvalue1d,precip_data) + +#endif + endif + else + write(LDT_logunit,*) '[WARN] LIS file '//trim(fname) + write(LDT_logunit,*) '[WARN] not found ...' + precip_data = LDT_rc%udef + endif + + call LDT_logSingleDAobs(n,LDT_DAobsData(1)%totalprecip_obs,& + precip_data,vlevel=1) + +#endif +end subroutine readLISlsmPrecipobs + + +!BOP +! +! !ROUTINE: create_lsm_output_fname +! \label{create_lsm_output_fname} +! +! !INTERFACE: +subroutine create_lsm_output_fname(n, form, fname, odir, wstyle, wopt, & + run_dd, map_proj, security_class, & + distribution_class, data_category, & + area_of_data, write_interval) +! !USES: + use LDT_coreMod, only : LDT_rc + use LDT_logMod + + implicit none +! !ARGUMENTS: + integer, intent(IN) :: n + character(len=*) :: fname + character(len=*) :: form + character(len=*) :: odir + character(len=*) :: wstyle + character(len=*) :: wopt + real, dimension(8), optional :: run_dd + character(len=*), optional :: map_proj + character(len=*), optional :: security_class + character(len=*), optional :: distribution_class + character(len=*), optional :: data_category + character(len=*), optional :: area_of_data + character(len=*), optional :: write_interval +! +! !DESCRIPTION: +! Create the file name for the output data files. It creates both the GSWP +! style of output filenames and the standard LIS style. The convention used +! in LIS creates a filename in a hierarchical style (output directory, +! model name, date, file extention) +! +! 2 level hierarchy +! \begin{verbatim} +! //LIS_HIST_. +! \end{verbatim} +! 3 level hierarchy +! \begin{verbatim} +! ///LIS_HIST_. +! \end{verbatim} +! 4 level hierarchy +! \begin{verbatim} +! ////LIS_HIST_. +! \end{verbatim} +! WMO convention +! \begin{verbatim} +! / +! \end{verbatim} +! A filename in the convention of weather products (such as): \newline +! {\small +! PS.AFWA\_SC.U\_DI.C\_DC.ANLYS\_GP.LIS\_GR.C0P25DEG\_AR.GLOBAL\_PA.03-HR-SUM\_DD.YYYYMMDD\_DT.HH00\_DF.GR1 \newline +! } +! where \newline +! PS = Product source \newline +! SC = security classification \newline +! DI = distribution classification \newline +! DC = data category \newline +! GP = generating process \newline +! GR = grid \newline +! AR = area of data \newline +! PA = parameter \newline +! DD = date \newline +! DT = data time \newline +! DF = data format \newline +! +! The arguments are: +! \begin{description} +! \item [n] +! index of the domain or nest +! \item [fname] +! the created file name. +! \item [model\_name] +! string describing the name of the model +! \item [writeint] +! output writing interval of the model +! \item [style] +! style option as described above +! \end{description} +!EOP + character(len=8) :: date + character(len=10) :: time + character(len=5) :: zone + integer, dimension(8) :: values + character(len=20) :: mname + character(len=10) :: cdate + character(len=14) :: cdate1 + character(len=2) :: fint + character(len=10) :: fres + character(len=10) :: fres2 + character(len=10) :: fres3 + character*1 :: fres1(10) + character(len=1) :: fproj + integer :: curr_mo = 0 + character(len=200) :: dname + character(len=200), save :: out_fname + integer :: i, c + + mname = 'SURFACEMODEL' + if(wstyle.eq."4 level hierarchy") then + write(unit=cdate1, fmt='(i4.4, i2.2, i2.2, i2.2, i2.2)') & + LDT_rc%yr, LDT_rc%mo, & + LDT_rc%da, LDT_rc%hr,LDT_rc%mn + + dname = trim(odir)//'/' + dname = trim(dname)//trim(mname)//'/' + + write(unit=cdate, fmt='(i4.4)') LDT_rc%yr + dname = trim(dname)//trim(cdate)//'/' + + write(unit=cdate, fmt='(i4.4, i2.2, i2.2)') LDT_rc%yr, LDT_rc%mo, LDT_rc%da + dname = trim(dname)//trim(cdate) + + out_fname = trim(dname)//'/LIS_HIST_'//trim(cdate1) + + write(unit=cdate, fmt='(a2,i2.2)') '.d',n + out_fname = trim(out_fname)//trim(cdate) + + select case ( form ) + case ( "binary" ) + if(wopt.eq."1d tilespace") then + out_fname = trim(out_fname)//'.ts4r' + elseif(wopt.eq."2d gridspace") then + out_fname = trim(out_fname)//'.gs4r' + elseif(wopt.eq."1d gridspace") then + out_fname = trim(out_fname)//'.gs4r' + endif + case ("grib1") + out_fname = trim(out_fname)//'.grb' + case ("netcdf") + out_fname = trim(out_fname)//'.nc' + case ("grib2") + out_fname = trim(out_fname)//'.gr2' + case default + call ldt_log_msg('ERR: create_lsm_output_fname -- '// & + 'Unrecognized output format') + call LDT_endrun + endselect + elseif(wstyle.eq."3 level hierarchy") then + write(unit=cdate1, fmt='(i4.4, i2.2, i2.2, i2.2, i2.2)') & + LDT_rc%yr, LDT_rc%mo, & + LDT_rc%da, LDT_rc%hr,LDT_rc%mn + + dname = trim(odir)//'/' + dname = trim(dname)//trim(mname)//'/' + + write(unit=cdate, fmt='(i4.4, i2.2)') LDT_rc%yr, LDT_rc%mo + dname = trim(dname)//trim(cdate)//'/' + + out_fname = trim(dname)//'LIS_HIST_'//trim(cdate1) + + write(unit=cdate, fmt='(a2,i2.2)') '.d',n + out_fname = trim(out_fname)//trim(cdate) + + select case ( form ) + case ("binary") + if(wopt.eq."1d tilespace") then + out_fname = trim(out_fname)//'.ts4r' + elseif(wopt.eq."2d gridspace") then + out_fname = trim(out_fname)//'.gs4r' + elseif(wopt.eq."1d gridspace") then + out_fname = trim(out_fname)//'.gs4r' + endif + case ("grib1") + out_fname = trim(out_fname)//'.grb' + case ("netcdf") + out_fname = trim(out_fname)//'.nc' + case ("grib2") + out_fname = trim(out_fname)//'.gr2' + case default + call ldt_log_msg('ERR: create_lsm_output_fname -- '// & + 'Unrecognized form value') + call LDT_endrun + endselect + elseif(wstyle.eq."2 level hierarchy") then + write(unit=cdate1, fmt='(i4.4, i2.2, i2.2, i2.2, i2.2)') & + LDT_rc%yr, LDT_rc%mo, & + LDT_rc%da, LDT_rc%hr,LDT_rc%mn + + dname = trim(odir)//'/' + dname = trim(dname)//trim(mname)//'/' + + out_fname = trim(dname)//'LIS_HIST_'//trim(cdate1) + + write(unit=cdate, fmt='(a2,i2.2)') '.d',n + out_fname = trim(out_fname)//trim(cdate) + + select case ( form ) + case ("binary") + if(wopt.eq."1d tilespace") then + out_fname = trim(out_fname)//'.ts4r' + elseif(wopt.eq."2d gridspace") then + out_fname = trim(out_fname)//'.gs4r' + elseif(wopt.eq."1d gridspace") then + out_fname = trim(out_fname)//'.gs4r' + endif + case ("grib1") + out_fname = trim(out_fname)//'.grb' + case ("netcdf") + out_fname = trim(out_fname)//'.nc' + case ("grib2") + out_fname = trim(out_fname)//'.gr2' + case default + call ldt_log_msg('ERR: create_lsm_output_fname -- '// & + 'Unrecognized form value') + call LDT_endrun + endselect + elseif(wstyle.eq."WMO convention") then + if ( .not. present(run_dd) .or. & + .not. present(security_class) .or. & + .not. present(distribution_class) .or. & + .not. present(data_category) .or. & + .not. present(area_of_data) .or. & + .not. present(write_interval) ) then + call ldt_log_msg('ERR: create_lsm_output_fname -- '// & + 'missing WMO convention identifiers') + call LDT_endrun + endif + + write(unit=cdate1, fmt='(i4.4, i2.2, i2.2)') & + LDT_rc%yr, LDT_rc%mo, LDT_rc%da + + write(unit=cdate, fmt='(i2.2, i2.2)') LDT_rc%hr, LDT_rc%mn + + if(map_proj.eq."polar") then + fproj = 'P' + print *,"fres ",run_dd(6) + if (run_dd(6) .ge. 10.) then + write(unit=fres, fmt='(i2)') nint(run_dd(6)) + else + write(unit=fres, fmt='(i1)') nint(run_dd(6)) + endif + fres2 = trim(fres)//'KM' + elseif(map_proj.eq."lambert") then + fproj = 'L' + print *,"fres ",run_dd(6) +! write(unit=fres, fmt='(f2.0)') run_dd(6) + write(unit=fres, fmt='(f3.0)') run_dd(6) + if (run_dd(6) .ge. 10.) then + write(unit=fres, fmt='(i2)') nint(run_dd(6)) + else + write(unit=fres, fmt='(i1)') nint(run_dd(6)) + endif + fres2 = trim(fres)//'KM' + elseif(map_proj.eq."mercator") then + fproj = 'M' + write(unit=fres, fmt='(i2.2)') run_dd(6) + fres = trim(fres)//'KM' + elseif(map_proj.eq."gaussian") then + fproj = 'G' + write(unit=fres, fmt='(i2.2)') run_dd(5)*100 + fres2 = '0P'//trim(fres)//'DEG' + else + fproj = 'C' + write(unit=fres, fmt='(i10)') nint(run_dd(6)*100) + read(unit=fres,fmt='(10a1)') (fres1(i),i=1,10) + c = 0 + do i=1,10 + if(fres1(i).ne.' '.and.c==0) c = i + enddo + if (run_dd(6) .lt. 0.1) then + fres3 = '0P0' + else + fres3 = '0P' + end if + fres2 = fres3 + do i=c,10 + fres2 = trim(fres2)//trim(fres1(i)) + enddo + fres2 = trim(fres2)//'DEG' + endif + + out_fname = trim(odir)//'/'//& + '/PS.AFWA_SC.'//trim(security_class)//& + '_DI.'//trim(distribution_class)//& + '_DC.'//trim(data_category)//& + '_GP.LIS_GR.'//& + trim(fproj)//trim(fres2)//& + '_AR.'//trim(area_of_data)//& + '_PA.'//trim(write_interval)//'-HR-SUM_DD.'//& + trim(cdate1)//'_DT.'//trim(cdate)//'_DF.GR1' + endif + fname = out_fname + end subroutine create_lsm_output_fname + +!BOP +! +! !ROUTINE: transformPrecipDataToLDTgrid +! \label{trasnformDataToLDTgrid} +! +! !INTERFACE: + subroutine transformPrecipDataToLDTgrid(n, precip_inp, precip_out) +! !USES: + use LDT_coreMod + use LISlsmPrecip_obsMod + + implicit none +! !ARGUMENTS: + integer :: n + real :: precip_inp(lsmprecipobs%nc*lsmprecipobs%nr) + real :: precip_out(LDT_rc%lnc(n),LDT_rc%lnr(n)) +! +! !DESCRIPTION: +! This routine interpolates or upscales the input data to +! the LDT grid. If the input data is finer than the LDT +! grid, the input data is upscaled. If the input data is +! coarser, then it is interpolated to the LDT grid. +! +!EOP + integer :: ios + integer :: c,r + logical*1 :: precip_data_b(lsmprecipobs%nc*lsmprecipobs%nr) + real :: precipobs_ip(LDT_rc%lnc(n)*LDT_rc%lnr(n)) + logical*1 :: precipobs_b_ip(lsmprecipobs%nc*lsmprecipobs%nr) + + do r=1,lsmprecipobs%nr + do c=1, lsmprecipobs%nc + if(precip_inp(c+(r-1)*lsmprecipobs%nc).ne.LDT_rc%udef) then + precip_data_b(c+(r-1)*lsmprecipobs%nc) = .true. + else + precip_data_b(c+(r-1)*lsmprecipobs%nc) = .false. + endif + if(precip_inp(c+(r-1)*lsmprecipobs%nc).gt.1) then + precip_inp(c+(r-1)*lsmprecipobs%nc) = LDT_rc%udef + precip_data_b(c+(r-1)*lsmprecipobs%nc) = .false. + endif + enddo + enddo + + if(LDT_isLDTatAfinerResolution(n,lsmprecipobs%datares)) then + +!-------------------------------------------------------------------------- +! Interpolate to the LDT running domain +!-------------------------------------------------------------------------- + call bilinear_interp(LDT_rc%gridDesc(n,:),& + precip_data_b, precip_inp, precipobs_b_ip, precipobs_ip, & + lsmprecipobs%nc*lsmprecipobs%nr, & + LDT_rc%lnc(n)*LDT_rc%lnr(n), & + LDT_domain(n)%lat, LDT_domain(n)%lon,& + lsmprecipobs%w11, lsmprecipobs%w12, & + lsmprecipobs%w21, lsmprecipobs%w22, & + lsmprecipobs%n11, lsmprecipobs%n12, & + lsmprecipobs%n21, lsmprecipobs%n22, & + LDT_rc%udef, ios) + + call neighbor_interp(LDT_rc%gridDesc(n,:),& + precip_data_b, precip_inp, precipobs_b_ip, precipobs_ip, & + lsmprecipobs%nc*lsmprecipobs%nr, & + LDT_rc%lnc(n)*LDT_rc%lnr(n), & + LDT_domain(n)%lat, LDT_domain(n)%lon,& + lsmprecipobs%n11, LDT_rc%udef, ios) + else + call upscaleByAveraging(& + lsmprecipobs%nc*lsmprecipobs%nr,& + LDT_rc%lnc(n)*LDT_rc%lnr(n),LDT_rc%udef, & + lsmprecipobs%n11,precip_data_b, precip_inp, precipobs_b_ip,precipobs_ip) + + endif + + do r=1,LDT_rc%lnr(n) + do c=1,LDT_rc%lnc(n) + if(precipobs_b_ip(c+(r-1)*LDT_rc%lnc(n))) then + precip_out(c,r) = precipobs_ip(c+(r-1)*LDT_rc%lnc(n)) + else + precip_out(c,r) = LDT_rc%udef + endif + enddo + enddo + + end subroutine transformPrecipDataToLDTgrid diff --git a/ldt/DAobs/LISlsmPrecipobs/read_Drange.F90 b/ldt/DAobs/LISlsmPrecipobs/read_Drange.F90 new file mode 100644 index 000000000..5d4f0b7af --- /dev/null +++ b/ldt/DAobs/LISlsmPrecipobs/read_Drange.F90 @@ -0,0 +1,112 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.3 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LDT_misc.h" +#include "LDT_NetCDF_inc.h" +!BOP +! !ROUTINE: read_Drange +! \label{read_Drange} + +! !REVISION HISTORY: +! 2Dec2021: Mahdi Navari ; Initial Specification +! +! !INTERFACE: +subroutine read_Drange(ngrid, filename, varname, xrange) + + use LDT_coreMod + use LDT_logMod + use LDT_historyMod +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !ARGUMENTS: + integer, intent(in) :: ngrid + character(len=*) :: filename + character(len=*) :: varname + real :: xrange(ngrid,2) + +! +! !DESCRIPTION: +! This routine reads the input CDF file (generated by LDT in NETCDF format) +! The xrange values and the corresponding CDFs are read for each grid point. +! Both these fields are expected to be in the 1-d grid vector dimension. +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[nbins] number of bins used to compute the model and obs CDFs +! \item[filename] name of the CDF file +! \item[varname] name of the variable being extracted. +! \item[xrange] x-axis values corresponding to the CDF +! \item[cdf] y-axis (CDF) values corresponding to the CDF +! \end{description} +!EOP + integer :: j,kk + integer :: ngridId, nbinsId, nlevsId,ntimesId + integer :: ngrid_file, nbins_file, nlevs_file, ntimes_file + integer :: xid, cdfid + real, allocatable :: xrange_file(:,:,:) + integer :: nid + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + write(LDT_logunit,*) '[INFO] Reading Drange form CDF file ',trim(filename) + call LDT_verify(nf90_open(path=trim(filename),mode=NF90_NOWRITE,& + ncid=nid),'failed to open file '//trim(filename)) + + call LDT_verify(nf90_inq_dimid(nid,"ngrid",ngridId), & + 'nf90_inq_dimid failed for ngrid') + call LDT_verify(nf90_inq_dimid(nid,"nbins",nbinsId), & + 'nf90_inq_dimid failed for nbins') + call LDT_verify(nf90_inq_dimid(nid,trim(varname)//"_levels",nlevsId), & + 'nf90_inq_dimid failed for '//trim(varname)//"_levels") + call LDT_verify(nf90_inq_dimid(nid,"ntimes",ntimesId), & + 'nf90_inq_dimid failed for ntimes') + + call LDT_verify(nf90_inquire_dimension(nid,ngridId, len=ngrid_file),& + 'nf90_inquire_dimension failed for ngrid') + call LDT_verify(nf90_inquire_dimension(nid,nbinsId, len=nbins_file),& + 'nf90_inquire_dimension failed for nbins') + call LDT_verify(nf90_inquire_dimension(nid,nlevsId, len=nlevs_file),& + 'nf90_inquire_dimension failed for nbins') + call LDT_verify(nf90_inquire_dimension(nid,ntimesId, len=ntimes_file),& + 'nf90_inquire_dimension failed for ntimes') + + if (ntimes_file .gt. 1) then + write(LDT_logunit,*) '[ERR] The number of times specified in the file '//& + trim(filename) + write(LDT_logunit,*) '[ERR] (',ntimes_file, & + ') should be 1 set the Temporal resolution of precipitation CDFs to "yearly" ' + call LDT_endrun() + endif + + allocate(xrange_file(ngrid_file,nlevs_file, nbins_file)) + + do j=1,ntimes_file + call LDT_verify(nf90_inq_varid(nid,trim(varname)//'_xrange',xid),& + 'nf90_inq_varid failed for for '//trim(varname)//'_xrange') + + call LDT_verify(nf90_get_var(nid,xid,xrange_file, & + start=(/1,j,1,1/), count=(/ngrid_file,1,nlevs_file,nbins_file/)),& + 'nf90_get_var failed for '//trim(varname)//'_xrange') + enddo + + xrange(:,1) = xrange_file(:,1,1) ! min value for each gridcell + xrange(:,2) = xrange_file(:,1,nbins_file) ! max value for each gridcell + + deallocate(xrange_file) + + call LDT_verify(nf90_close(nid),& + 'failed to close file '//trim(filename)) + write(LDT_logunit,*) '[INFO] Successfully read CDF file ',trim(filename) +#endif + +end subroutine read_Drange + diff --git a/ldt/DAobs/LISlsmPrecipobs/read_Precip_climo.F90 b/ldt/DAobs/LISlsmPrecipobs/read_Precip_climo.F90 new file mode 100644 index 000000000..34fb4f906 --- /dev/null +++ b/ldt/DAobs/LISlsmPrecipobs/read_Precip_climo.F90 @@ -0,0 +1,116 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.3 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LDT_misc.h" +#include "LDT_NetCDF_inc.h" +!BOP +! !ROUTINE: read_Precip_climo +! \label{read_Precip_climo} + +! !REVISION HISTORY: +! 19Jan2022: Mahdi Navari ; Initial Specification +! +! !INTERFACE: +subroutine read_Precip_climo(ncol, nrow, filename, precip) + + use LDT_coreMod + use LDT_logMod + use LDT_historyMod +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none + ! !ARGUMENTS: + integer, intent(in) :: ncol, nrow + character(len=*), intent(in) :: filename + real :: precip(ncol,nrow,12) + logical :: file_exists +! +! !DESCRIPTION: +! This routine reads the input CDF file (generated by LDT in NETCDF format) +! The xrange values and the corresponding CDFs are read for each grid point. +! Both these fields are expected to be in the 1-d grid vector dimension. +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[nbins] number of bins used to compute the model and obs CDFs +! \item[filename] name of the CDF file +! \item[varname] name of the variable being extracted. +! \item[xrange] x-axis values corresponding to the CDF +! \item[cdf] y-axis (CDF) values corresponding to the CDF +! \end{description} +!EOP + real, allocatable :: precip_climo(:,:,:) + integer :: ios,nid,ncId,nrId,varId + integer :: nc,nr,i,c,r + character*3 :: month_name(12) + + month_name = (/"JAN","FEB","MAR","APR","MAY","JUN",& + "JUL","AUG","SEP","OCT","NOV","DEC"/) + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + inquire(file=trim(filename), exist=file_exists) + + if (.not. file_exists) then + write(LDT_logunit,*) '[ERR] (',filename, & + ') does not exist. Please provide a monthly precipitation climatology in CDF format ' + call LDT_endrun() + end if + + + write(LDT_logunit,*) '[INFO] Reading Precipitation climatology from CDF file ',trim(filename) + call LDT_verify(nf90_open(path=trim(filename),mode=NF90_NOWRITE,& + ncid=nid),'failed to open file '//trim(filename)) + + ios = nf90_inq_dimid(nid,"east_west",ncId) + call LDT_verify(ios,'Error in nf90_inq_dimid in read_precip_climo') + + ios = nf90_inq_dimid(nid,"north_south",nrId) + call LDT_verify(ios,'Error in nf90_inq_dimid in read_precip_climo') + + ios = nf90_inquire_dimension(nid,ncId, len=nc) + call LDT_verify(ios,'Error in nf90_inquire_dimension in read_precip_climo') + + ios = nf90_inquire_dimension(nid,nrId, len=nr) + call LDT_verify(ios,'Error in nf90_inquire_dimension in read_precip_climo') + + if (ncol .ne. nc .or. nrow .ne.nr) then + write(LDT_logunit,*) & + '[ERR] The number of columns or rows specified in the file '//& + trim(filename) + write(LDT_logunit,*) '[ERR] (', nc, nr, & + ') is different from the number of columns and rows specified' + write(LDT_logunit,*) '[ERR] in the ldt.config file (', ncol, nrow, ')' + call LDT_endrun() + endif + + allocate(precip_climo(nc,nr,12)) + + do i=1,12 + ios = nf90_inq_varid(nid,'TotalPrecip_'//trim(month_name(i)),varId) + call LDT_verify(ios,'Precipitation climo field not found in the file') + + ios = nf90_get_var(nid,varId,precip_climo(:,:,i)) + call LDT_verify(ios,'Error in nf90_get_var in LDT_procDataForREAL') + enddo + + ios = nf90_close(nid) + call LDT_verify(ios,'Error in nf90_close in readldtparam_real_2d') + write(LDT_logunit,*) '[INFO] Successfully read Precipitation climo file ',trim(filename) + + precip = precip_climo + deallocate(precip_climo) + +! end USE_NETCDF4 +#endif + +end subroutine read_Precip_climo + diff --git a/ldt/configs/ldt.config.adoc b/ldt/configs/ldt.config.adoc index aee7c0648..0ca2343ba 100644 --- a/ldt/configs/ldt.config.adoc +++ b/ldt/configs/ldt.config.adoc @@ -3703,6 +3703,54 @@ landcover.1gd4r + CDF grouping attributes file: cdf_grouping.txt .... +`Stratify CDFs by external data:` specifies whether to stratify CDFs for each pixel by an externally specified categorical map; for example, by precipitation climatology. +This option generates a stratified CDF for each grid cell. That means if the domain contains 10000 grid cells, the CDF file also has 10000 CDFs. But CDFs in grid cells with a similar precipitation climatology are the same. +Acceptable values are: + +[cols="<,<",] +|=== +|Value |Description + +|0 |do not stratify by external data +|1 |stratify by external data +|=== + +.Example _ldt.config_ entry +.... +Stratify CDFs by external data: 1 +.... + +`Number of bins to use for stratification:` specifies the number of stratification bins to use for computing the CDF. + +.Example _ldt.config_ entry +.... +Number of bins to use for stratification: 15 +.... + +`Stratification data source:` specifies the name of the stratification data source + +.Example _ldt.config_ entry +.... +Stratification data source: "LIS LSM total precipitation" +.... + +`External stratification file:` Specifies the name of a NetCDF file that contains the monthly precipitation climatology. The user needs to run an LSM and generate daily total precipitation output for a long period. Then the user needs to run LVT to generate monthly precipitation climatology. + + +.Example _ldt.config_ entry +.... +External stratification file: "LVT_MEAN_FINAL.202201010000.d01.nc" +.... + + +`Write stratified geolocation independent CDFs:` specifies whether to write the stratified CDF or not. +This option eliminates the similar CDFs and collapses the stratified CDFs into a number of stratified bins. For example, if the domain contains 10000 grid cells and you chose to stratify the CDFs to 15 precipitation bins. This option collapses the 10000 CDFs to 15 CDFs. This CDF is a geolocation independent CDFs. + +.Example _ldt.config_ entry +.... +Write stratified geolocation independent CDFs: 1 +.... + `Temporal averaging interval:` specifies temporal averaging interval to be used while computing the CDF. .Example _ldt.config_ entry diff --git a/ldt/core/LDT_CDFMod.F90 b/ldt/core/LDT_CDFMod.F90 index 8ffa89441..a3d68e13c 100644 --- a/ldt/core/LDT_CDFMod.F90 +++ b/ldt/core/LDT_CDFMod.F90 @@ -17,6 +17,8 @@ module LDT_CDFMod ! ! !REVISION HISTORY: ! 2 Oct 2008 Sujay Kumar Initial Specification +! 2 Dec 2021: Mahdi Navari; modified to stratify CDF based on precipitation +! and save the stratified CDF ! !EOP use LDT_DAobsDataMod @@ -214,14 +216,14 @@ subroutine computeSingleCDF(n,obs, metrics) integer :: t,i,j,k integer :: c,r - integer :: sindex + integer :: sindex,sindex0,sindex1 integer, allocatable :: strat_bincounts(:,:,:,:) real, allocatable :: strat_cdf(:,:,:,:) ! real :: datamask(LDT_rc%ngrid(n)) if(LDT_rc%endtime.eq.1) then if(obs%selectOpt.eq.1.and.metrics%selectOpt.eq.1) then - if(LDT_rc%group_cdfs.eq.0) then + if(LDT_rc%group_cdfs.eq.0 .and. LDT_rc%strat_cdfs.eq.0) then do t=1,LDT_rc%ngrid(n) do j=1,LDT_rc%cdf_ntimes do k=1,obs%vlevels @@ -237,7 +239,7 @@ subroutine computeSingleCDF(n,obs, metrics) enddo enddo enddo - elseif(LDT_rc%group_cdfs.eq.1) then + elseif(LDT_rc%group_cdfs.eq.1 .and. LDT_rc%strat_cdfs.eq.0) then allocate(strat_bincounts(LDT_rc%group_cdfs_nbins, & LDT_rc%cdf_ntimes, & @@ -294,7 +296,136 @@ subroutine computeSingleCDF(n,obs, metrics) enddo ! endif enddo - + ! MN: save stratified CDF + metrics%strat_cdf = strat_cdf + deallocate(strat_bincounts) + deallocate(strat_cdf) + +!MN: Startification based on monthly precipitation climatology +! monthly total precipitation climatology for each pixel are stored in LDT_rc%stratification_data + + elseif(LDT_rc%group_cdfs.eq.0 .and. LDT_rc%strat_cdfs.eq.1 ) then + + allocate(strat_bincounts(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels,& + LDT_rc%cdf_nbins)) + allocate(strat_cdf(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels,& + LDT_rc%cdf_nbins)) + + strat_bincounts = 0 + strat_cdf = 0 + + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + do i=1,LDT_rc%cdf_nbins + sindex = LDT_rc%stratification_data(t,j) + strat_bincounts(sindex,j,k,i) = & + strat_bincounts(sindex,j,k,i) + & + metrics%cdf_bincounts(t,j,k,i) + + enddo + enddo + enddo + enddo + + do t=1,LDT_rc%strat_cdfs_nbins + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + if(sum(strat_bincounts(t,j,k,:)).ne.0) then + strat_cdf(t,j,k,1) = & + strat_bincounts(t,j,k,1)/& ! strat_cdf(t,j,k,1)/& + float(sum(strat_bincounts(t,j,k,:))) + do i=2,LDT_rc%cdf_nbins + strat_cdf(t,j,k,i) = strat_cdf(t,j,k,i-1) + & + strat_bincounts(t,j,k,i)/& + float(sum(strat_bincounts(t,j,k,:))) + enddo + endif + enddo + enddo + enddo + + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + do i=1,LDT_rc%cdf_nbins + sindex = LDT_rc%stratification_data(t,j) + metrics%cdf(t,j,k,i) = strat_cdf(sindex,j,k,i) + enddo + enddo + enddo + enddo + ! MN: save stratified CDF + metrics%strat_cdf = strat_cdf + deallocate(strat_bincounts) + deallocate(strat_cdf) + + elseif(LDT_rc%group_cdfs.eq.1 .and. LDT_rc%strat_cdfs.eq.1 ) then + + allocate(strat_bincounts(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels,& + LDT_rc%cdf_nbins)) + allocate(strat_cdf(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels,& + LDT_rc%cdf_nbins)) + + strat_bincounts = 0 + strat_cdf = 0 + + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + do i=1,LDT_rc%cdf_nbins + sindex0 = LDT_rc%stratification_data(t,j) + sindex1 = LDT_rc%cdf_strat_data(t) + sindex = sindex0 + (sindex1 - 1)*LDT_rc%strat_cdfs_nbins + strat_bincounts(sindex,j,k,i) = & + strat_bincounts(sindex,j,k,i) + & + metrics%cdf_bincounts(t,j,k,i) + + enddo + enddo + enddo + enddo + + do t=1,LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + if(sum(strat_bincounts(t,j,k,:)).ne.0) then + strat_cdf(t,j,k,1) = & + strat_bincounts(t,j,k,1)/& ! strat_cdf(t,j,k,1)/& + float(sum(strat_bincounts(t,j,k,:))) + do i=2,LDT_rc%cdf_nbins + strat_cdf(t,j,k,i) = strat_cdf(t,j,k,i-1) + & + strat_bincounts(t,j,k,i)/& + float(sum(strat_bincounts(t,j,k,:))) + enddo + endif + enddo + enddo + enddo + + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + do i=1,LDT_rc%cdf_nbins + sindex0 = LDT_rc%stratification_data(t,j) + sindex1 = LDT_rc%cdf_strat_data(t) + sindex = sindex0 + (sindex1 - 1)*LDT_rc%strat_cdfs_nbins + metrics%cdf(t,j,k,i) = strat_cdf(sindex,j,k,i) + enddo + enddo + enddo + enddo + + ! MN: save stratified CDF + metrics%strat_cdf = strat_cdf deallocate(strat_bincounts) deallocate(strat_cdf) diff --git a/ldt/core/LDT_DAmetricsDataMod.F90 b/ldt/core/LDT_DAmetricsDataMod.F90 index dadb62de4..23fc1e607 100644 --- a/ldt/core/LDT_DAmetricsDataMod.F90 +++ b/ldt/core/LDT_DAmetricsDataMod.F90 @@ -18,6 +18,7 @@ module LDT_DAmetricsDataMod ! ! !REVISION HISTORY: ! 02 Oct 2008: Sujay Kumar; Initial version +! 28 Feb 2022: Mahdi navari modified to save stratified CDFs ! ! !USES: use ESMF @@ -36,6 +37,9 @@ module LDT_DAmetricsDataMod real, allocatable :: delta(:,:,:) real, allocatable :: xrange(:,:,:,:) real, allocatable :: cdf(:,:,:,:) + real, allocatable :: strat_xrange(:,:,:,:) + real, allocatable :: strat_cdf(:,:,:,:) + real, allocatable :: sx_mu(:,:,:) real, allocatable :: mu(:,:,:) diff --git a/ldt/core/LDT_DAmetricsMod.F90 b/ldt/core/LDT_DAmetricsMod.F90 index 78a76dc27..ce02ea2d9 100644 --- a/ldt/core/LDT_DAmetricsMod.F90 +++ b/ldt/core/LDT_DAmetricsMod.F90 @@ -18,6 +18,7 @@ module LDT_DAmetricsMod ! ! !REVISION HISTORY: ! 02 Oct 2008: Sujay Kumar; Initial version +! 2 Dec 2021: Mahdi Navari; modified to compute CDF for precipitation ! use ESMF use LDT_DAmetricsDataMod @@ -99,6 +100,8 @@ subroutine LDT_DAmetricsInit LDT_DAobsData(n)%vod_obs,LDT_DAmetrics%vod) call registerMetricsEntry(LDT_DA_MOC_LAI,nsize,& LDT_DAobsData(n)%lai_obs,LDT_DAmetrics%lai) + call registerMetricsEntry(LDT_DA_MOC_TOTALPRECIP,nsize,& + LDT_DAobsData(n)%totalprecip_obs,LDT_DAmetrics%totalprecip) call registerMetricsEntry(LDT_DA_MOC_GVF,nsize,& LDT_DAobsData(n)%gvf_obs,LDT_DAmetrics%gvf) !Y.Kwon !------------------------------------------------------------------------ @@ -597,6 +600,33 @@ subroutine initMetricsEntry(nsize, obs, metrics) LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) allocate(metrics%cdf(nsize,& LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) + + if(LDT_rc%write_strat_cdfs.eq.1) then + if(LDT_rc%group_cdfs.eq.1 .and. LDT_rc%strat_cdfs.eq.0) then + allocate(metrics%strat_xrange(LDT_rc%group_cdfs_nbins,& + LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) + allocate(metrics%strat_cdf(LDT_rc%group_cdfs_nbins,& + LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) + LDT_rc%stratified_cdfs_nbins = LDT_rc%group_cdfs_nbins + elseif(LDT_rc%group_cdfs.eq.0 .and. LDT_rc%strat_cdfs.eq.1 ) then + allocate(metrics%strat_xrange(LDT_rc%strat_cdfs_nbins,& + LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) + allocate(metrics%strat_cdf(LDT_rc%strat_cdfs_nbins,& + LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) + LDT_rc%stratified_cdfs_nbins = LDT_rc%strat_cdfs_nbins + elseif(LDT_rc%group_cdfs.eq.1 .and. LDT_rc%strat_cdfs.eq.1 ) then + allocate(metrics%strat_xrange(LDT_rc%group_cdfs_nbins*& + LDT_rc%strat_cdfs_nbins,& + LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) + allocate(metrics%strat_cdf(LDT_rc%group_cdfs_nbins*& + LDT_rc%strat_cdfs_nbins,& + LDT_rc%cdf_ntimes, obs%vlevels,LDT_rc%cdf_nbins)) + LDT_rc%stratified_cdfs_nbins = LDT_rc%group_cdfs_nbins*& + LDT_rc%strat_cdfs_nbins + endif + metrics%strat_xrange = 0 + metrics%strat_cdf = 0 + endif metrics%cdf_bincounts = 0 metrics%delta = 0 @@ -735,7 +765,7 @@ subroutine LDT_computeDAobsMetrics(n, pass) implicit none - character(len=LDT_CONST_PATH_LEN) :: fname_cdf + character(len=LDT_CONST_PATH_LEN) :: fname_cdf, fname_strat_cdf character(len=LDT_CONST_PATH_LEN) :: fname_domain integer :: pass integer :: rc @@ -782,6 +812,14 @@ subroutine LDT_computeDAobsMetrics(n, pass) write(LDT_logunit,*) 'Writing CDF domain file ',trim(fname_domain) iret=nf90_create(path=trim(fname_domain),cmode=nf90_clobber,& ncid=LDT_rc%ftn_DAobs_domain) + + if(LDT_rc%write_strat_cdfs.eq.1) then + fname_strat_cdf = trim(LDT_rc%odir)//'/stratified_'//& + trim(LDT_rc%dapreprocfile)//'.nc' + write(LDT_logunit,*) 'Writing stratified CDF file ',trim(fname_strat_cdf) + iret=nf90_create(path=trim(fname_strat_cdf),cmode=nf90_clobber,& + ncid=LDT_rc%ftn_strat_cdf) + endif endif #endif #if (defined USE_NETCDF4) @@ -796,7 +834,16 @@ subroutine LDT_computeDAobsMetrics(n, pass) write(LDT_logunit,*) 'Writing CDF domain file ',trim(fname_domain) iret=nf90_create(path=trim(fname_domain),cmode=nf90_netcdf4,& ncid=LDT_rc%ftn_DAobs_domain) + + if(LDT_rc%write_strat_cdfs.eq.1) then + fname_strat_cdf = trim(LDT_rc%odir)//'/stratified_'//& + trim(LDT_rc%dapreprocfile)//'.nc' + write(LDT_logunit,*) 'Writing stratified CDF file ',trim(fname_strat_cdf) + iret=nf90_create(path=trim(fname_strat_cdf),cmode=nf90_netcdf4,& + ncid=LDT_rc%ftn_strat_cdf) + endif endif + #endif call outputFinalMetrics(n,pass) @@ -806,6 +853,10 @@ subroutine LDT_computeDAobsMetrics(n, pass) write(LDT_logunit,*) 'Successfully wrote CDF file ',trim(fname_cdf) iret=nf90_close(LDT_rc%ftn_DAobs_domain) write(LDT_logunit,*) 'Successfully wrote CDF file ',trim(fname_domain) + if(LDT_rc%write_strat_cdfs.eq.1) then + iret=nf90_close(LDT_rc%ftn_strat_cdf) + write(LDT_logunit,*) 'Successfully wrote geolocation independent stratified CDF file ',trim(fname_strat_cdf) + endif endif #endif @@ -840,7 +891,7 @@ subroutine outputFinalMetrics(n,pass) !EOP integer :: index integer :: c,r - integer :: dimID(4) + integer :: dimID(4), dimID_strat(4) integer :: bdimID(3) character(len=8) :: date character(len=10) :: time @@ -902,6 +953,52 @@ subroutine outputFinalMetrics(n,pass) LDT_DAobsDataPtr(1,index)%dataEntryPtr) enddo + +! geolocation independent stratified CDF + if(LDT_rc%write_strat_cdfs.eq.1) then +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + call LDT_verify(nf90_def_dim(LDT_rc%ftn_strat_cdf,'n_strat_bins',& + LDT_rc%stratified_cdfs_nbins,dimID_strat(1)),'nf90_def_dim failed for n_strat_bins') + call LDT_verify(nf90_def_dim(LDT_rc%ftn_strat_cdf,'ntimes',& + LDT_rc%cdf_ntimes,dimID_strat(2)),'nf90_def_dim failed for ntimes') + call LDT_verify(nf90_def_dim(LDT_rc%ftn_strat_cdf,'nbins',& + LDT_rc%cdf_nbins,dimID_strat(4)),'nf90_def_dim failed for nbins') + + call LDT_verify(nf90_put_att(LDT_rc%ftn_strat_cdf,NF90_GLOBAL,& + "missing_value", -9999.0)) + call LDT_verify(nf90_put_att(LDT_rc%ftn_strat_cdf,NF90_GLOBAL,& + "temporal_resolution_CDF", LDT_rc%cdf_ntimes)) + call LDT_verify(nf90_put_att(LDT_rc%ftn_strat_cdf,NF90_GLOBAL,& + "title", & + "Land Data Toolkit (LDT) output")) + call LDT_verify(nf90_put_att(LDT_rc%ftn_strat_cdf,NF90_GLOBAL,& + "institution", & + "NASA GSFC Hydrological Sciences Laboratory")) + call LDT_verify(nf90_put_att(LDT_rc%ftn_strat_cdf,NF90_GLOBAL,& + "history", & + "created on date: "//date(1:4)//"-"//date(5:6)//"-"//& + date(7:8)//"T"//time(1:2)//":"//time(3:4)//":"//time(5:10))) + !call LDT_verify(nf90_put_att(LDT_rc%ftn_cdf,NF90_GLOBAL,"references", & + ! "Arsenault_etal_GMD_2018, Kumar_etal_EMS_2006")) + call LDT_verify(nf90_put_att(LDT_rc%ftn_strat_cdf,NF90_GLOBAL,"comment", & + "website: http://lis.gsfc.nasa.gov/")) +#endif + do index=1,LDT_DA_MOC_COUNT + call writeFinalSingleStratifiedCDFEntryHeader(LDT_DAmetricsPtr(index)%dataEntryPtr,& + LDT_DAobsDataPtr(1,index)%dataEntryPtr, dimID_strat) + enddo + +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + call LDT_verify(nf90_enddef(LDT_rc%ftn_strat_cdf)) +#endif + + do index=1,LDT_DA_MOC_COUNT + call writeFinalSingleStratifiedCDFEntry(pass,LDT_DAmetricsPtr(index)%dataEntryPtr,& + LDT_DAobsDataPtr(1,index)%dataEntryPtr) + enddo + + endif + !domain file #if (defined USE_NETCDF3 || defined USE_NETCDF4) @@ -1440,6 +1537,73 @@ subroutine writeFinalSingleEntryHeader(metrics, obs, dimID) end subroutine writeFinalSingleEntryHeader +!BOP +! !ROUTINE: writeFinalSingleStratifiedCDFEntryHeader +! \label{writeFinalSingleStratifiedCDFEntryHeader} +! +! !INTERFACE: + subroutine writeFinalSingleStratifiedCDFEntryHeader(metrics, obs, dimID) +! !USES: + use LDT_coreMod, only : LDT_rc + use LDT_historyMod, only : LDT_writevar_gridded + + implicit none +! !ARGUMENTS: + type(DAmetricsEntry) :: metrics + type(LDT_DAmetadataEntry) :: obs + integer :: dimID(4) +! +! !DESCRIPTION: +! This routine writes the specified set of statistics for a +! single variable at the end of the analysis. +!EOP + + integer :: k + integer :: varid1, varid2 + integer :: i,c,r + integer :: n + character*100 :: vname + integer :: shuffle, deflate, deflate_level + +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + n = 1 + shuffle = NETCDF_shuffle + deflate = NETCDF_deflate + deflate_level =NETCDF_deflate_level + + if(obs%selectOpt.eq.1.and.metrics%selectOpt.eq.1) then + vname = trim(obs%standard_name)//'_levels' + call LDT_verify(nf90_def_dim(LDT_rc%ftn_strat_cdf,trim(vname),& + obs%vlevels,dimID(3))) + + if(LDT_rc%write_strat_cdfs.eq.1) then + call LDT_verify(nf90_def_var(LDT_rc%ftn_strat_cdf,& + trim(obs%standard_name)//'_xrange_stratified',& + nf90_float, dimids = dimID, varid=obs%varID(1))) + +#if (defined USE_NETCDF4) + call LDT_verify(nf90_def_var_deflate(LDT_rc%ftn_strat_cdf,& + obs%varID(1), shuffle, deflate, deflate_level),& + 'nf90_def_var_deflate failed in LDT_DAmetricsMod') +#endif + + call LDT_verify(nf90_def_var(LDT_rc%ftn_strat_cdf,& + trim(obs%standard_name)//'_CDF_stratified',& + nf90_float, dimids = dimID, varid=obs%varID(4))) + +#if (defined USE_NETCDF4) + call LDT_verify(nf90_def_var_deflate(LDT_rc%ftn_strat_cdf,& + obs%varID(4), shuffle, deflate, deflate_level),& + 'nf90_def_var_deflate failed in LDT_DAmetricsMod') +#endif + endif + + endif +#endif + + end subroutine writeFinalSingleStratifiedCDFEntryHeader + + !BOP ! !ROUTINE: writeFinalSingleEntry @@ -1519,6 +1683,61 @@ subroutine writeFinalSingleEntry(pass,metrics, obs) end subroutine writeFinalSingleEntry +!BOP +! !ROUTINE: writeFinalSingleStratifiedCDFEntry +! \label{writeFinalSingleStratifiedCDFEntry} +! +! !INTERFACE: + subroutine writeFinalSingleStratifiedCDFEntry(pass,metrics, obs) +! !USES: + use LDT_coreMod, only : LDT_rc + use LDT_historyMod, only : LDT_writevar_gridded + + implicit none +! !ARGUMENTS: + integer :: pass + type(DAmetricsEntry) :: metrics + type(LDT_DAmetadataEntry) :: obs +! +! !DESCRIPTION: +! This routine writes the specified set of statistics for a +! single variable at the end of the analysis. +!EOP + + integer :: k + integer :: varid, varid2 + integer :: i,c,r + integer :: n + integer :: iret + + n = 1 + if(obs%selectOpt.eq.1.and.metrics%selectOpt.eq.1) then + + if(pass.eq.2.and.LDT_rc%comp_cdf.eq.1) then +! MN: The LDT_writevar_gridded is not generic it assumes the +! first dimension of CDF related vaialbes is LDT_rc%glbngrid(n) +! But for geolocation independent CDFs the first dimension is +! LDT_rc%stratified_cdfs_nbins +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + if(LDT_rc%write_strat_cdfs.eq.1) then + iret = nf90_put_var(LDT_rc%ftn_strat_cdf,obs%varID(1),metrics%strat_xrange,(/1,1,1,1/),& + (/LDT_rc%stratified_cdfs_nbins,LDT_rc%cdf_ntimes, & + obs%vlevels, LDT_rc%cdf_nbins/)) + call LDT_verify(iret, 'nf90_put_var failed in writevar_gridded_real_4d') + + iret = nf90_put_var(LDT_rc%ftn_strat_cdf,obs%varID(4),metrics%strat_cdf,(/1,1,1,1/),& + (/LDT_rc%stratified_cdfs_nbins, LDT_rc%cdf_ntimes, & + obs%vlevels, LDT_rc%cdf_nbins/)) + call LDT_verify(iret, 'nf90_put_var failed in writevar_gridded_real_4d') + endif +#endif + + endif + endif + + end subroutine writeFinalSingleStratifiedCDFEntry + + !BOP ! !ROUTINE: writeFinalSingleDomainEntryHeader ! \label{writeFinalSingleDomainEntryHeader} diff --git a/ldt/core/LDT_DAobsDataMod.F90 b/ldt/core/LDT_DAobsDataMod.F90 index d3d0431f2..27392061d 100644 --- a/ldt/core/LDT_DAobsDataMod.F90 +++ b/ldt/core/LDT_DAobsDataMod.F90 @@ -21,6 +21,7 @@ module LDT_DAobsDataMod ! ! !REVISION HISTORY: ! 2 Oct 2008 Sujay Kumar Initial Specification +! 2 Dec 2021: Mahdi Navari; modified to compute CDF for precipitation ! !USES: PRIVATE @@ -49,6 +50,7 @@ module LDT_DAobsDataMod public :: LDT_DA_MOC_LAI public :: LDT_DA_MOC_GVF !Y.Kwon public :: LDT_DA_MOC_COUNT + public :: LDT_DA_MOC_TOTALPRECIP ! public :: LDT_MOC_GRIB_COUNT ! ALMA ENERGY BALANCE COMPONENTS @@ -60,8 +62,9 @@ module LDT_DAobsDataMod integer, parameter :: LDT_DA_MOC_LAI = 6 integer, parameter :: LDT_DA_MOC_GVF = 7 !Y.Kwon integer, parameter :: LDT_DA_MOC_SOILTEFF = 8 !Y.Kwon + integer, parameter :: LDT_DA_MOC_TOTALPRECIP= 9 ! READ ABOVE NOTE ABOUT SPECIAL CASE INDICES - integer, parameter :: LDT_DA_MOC_COUNT = 8 !Y.Kwon + integer, parameter :: LDT_DA_MOC_COUNT = 9 !Y.Kwon ! Add the special cases. LDT_MOC_GRIB_COUNT should be used only in ! LDT_gribMod.F90. ! integer, parameter :: LDT_MOC_GRIB_COUNT = 100 @@ -103,6 +106,7 @@ module LDT_DAobsDataMod type(LDT_DAmetadataEntry) :: teff ! effective soil temperature (Y.Kwon) type(LDT_DAmetadataEntry) :: vod type(LDT_DAmetadataEntry) :: lai + type(LDT_DAmetadataEntry) :: totalprecip type(LDT_DAmetadataEntry) :: gvf !Y.Kwon end type output_meta @@ -119,6 +123,7 @@ module LDT_DAobsDataMod type(LDT_DAmetadataEntry) :: tws_obs type(LDT_DAmetadataEntry) :: vod_obs type(LDT_DAmetadataEntry) :: lai_obs + type(LDT_DAmetadataEntry) :: totalprecip_obs type(LDT_DAmetadataEntry) :: gvf_obs !Y.Kwon end type obs_list_dec @@ -209,6 +214,8 @@ subroutine LDT_DAobsEntryInit(i, nsize) LDT_DAobsData(i)%vod_obs,1,nsize,(/"-"/)) call register_obsDataEntry(i,LDT_DA_MOC_LAI ,& LDT_DAobsData(i)%lai_obs,1,nsize,(/"-"/)) + call register_obsDataEntry(i,LDT_DA_MOC_TOTALPRECIP ,& + LDT_DAobsData(i)%totalprecip_obs,1,nsize,(/"kg/m2"/)) call register_obsDataEntry(i,LDT_DA_MOC_GVF ,& LDT_DAobsData(i)%gvf_obs,1,nsize,(/"-"/)) !Y.Kwon end subroutine LDT_DAobsEntryInit diff --git a/ldt/core/LDT_DAobservationsMod.F90 b/ldt/core/LDT_DAobservationsMod.F90 index 35b1f76e9..9949bb8bc 100644 --- a/ldt/core/LDT_DAobservationsMod.F90 +++ b/ldt/core/LDT_DAobservationsMod.F90 @@ -18,6 +18,7 @@ module LDT_DAobservationsMod ! ! !REVISION HISTORY: ! 02 Oct 2008 Sujay Kumar Initial Specification +! 2 Dec 2021: Mahdi Navari; modified to compute CDF for precipitation ! ! !USES: @@ -71,6 +72,7 @@ subroutine LDT_DAobsDataInit call default_init_obsEntry(LDT_DAobsData(i)%vod_obs, "VOD") call default_init_obsEntry(LDT_DAobsData(i)%lai_obs, "LAI") call default_init_obsEntry(LDT_DAobsData(i)%gvf_obs, "GVF") !Y.Kwon + call default_init_obsEntry(LDT_DAobsData(i)%totalprecip_obs, "TotalPrecip") enddo call daobservationsetup(trim(LDT_rc%obs_src)//char(0)) diff --git a/ldt/core/LDT_DApreprocMod.F90 b/ldt/core/LDT_DApreprocMod.F90 index 9da0da58a..bce21191e 100644 --- a/ldt/core/LDT_DApreprocMod.F90 +++ b/ldt/core/LDT_DApreprocMod.F90 @@ -16,6 +16,7 @@ module LDT_DApreprocMod ! ! !REVISION HISTORY: ! 24 Nov 2008 Sujay Kumar Initial Specification +! 2 Dec 2021: Mahdi Navari; modified to save stratify CDF ! use ESMF @@ -53,13 +54,13 @@ subroutine LDT_DApreprocInit() character*20 :: tres integer :: rc integer :: n - integer :: c,r + integer :: c,r,j integer :: ios1 integer :: ftn character*50 :: preprocMethod real :: delta real, allocatable :: cdf_strat_data(:,:) - + real, allocatable :: stratification_data(:,:,:) n = 1 call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%obs_src,& label="DA observation source:",rc=rc) @@ -86,7 +87,7 @@ subroutine LDT_DApreprocInit() call LDT_verify(rc,'Name of the preprocessed DA file: not defined') if(LDT_rc%comp_cdf.gt.0) then - + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%cdf_nbins,& label="Number of bins to use in the CDF:",rc=rc) call LDT_verify(rc,'Number of bins to use in the CDF: not defined') @@ -166,11 +167,17 @@ subroutine LDT_DApreprocInit() do r=1,LDT_rc%lnr(n) do c=1,LDT_rc%lnc(n) - if(LDT_domain(n)%gindex(c,r).ne.-1) then - - LDT_rc%cdf_strat_data(LDT_domain(n)%gindex(c,r)) = & + if(LDT_domain(n)%gindex(c,r).ne.-1) then + + LDT_rc%cdf_strat_data(LDT_domain(n)%gindex(c,r)) = & nint((cdf_strat_data(c,r) - LDT_rc%group_cdfs_min)/& delta)+1 + if (LDT_rc%cdf_strat_data(LDT_domain(n)%gindex(c,r)) .gt. LDT_rc%group_cdfs_nbins) then + write(LDT_logunit,*) '[INFO] Group bins is larger then Max Group bins',& + LDT_rc%cdf_strat_data(LDT_domain(n)%gindex(c,r)), 'vs.', LDT_rc%group_cdfs_nbins ,& + 'Value adjusted the Max Group bins' + LDT_rc%cdf_strat_data(LDT_domain(n)%gindex(c,r)) = LDT_rc%group_cdfs_nbins + endif endif enddo enddo @@ -179,49 +186,108 @@ subroutine LDT_DApreprocInit() write(LDT_logunit,*) '[INFO] Finished reading ',& trim(LDT_rc%group_cdfs_strat_file) endif + +!This part reads the monthly total precipitation climatology and generates +! stratification input data LDT_rc%stratification_data(LDT_rc%ngrid(n),LDT_rc%cdf_ntimes). + + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%strat_cdfs,& + label="Stratify CDFs by external data:",default=0, rc=rc) + call LDT_verify(rc,"Stratify CDFs by external data: not defined") + + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%write_strat_cdfs,& + label="Write stratified geolocation independent CDFs:",default=0, rc=rc) + call LDT_verify(rc,"Write stratify geolocation independent CDFs:: not defined") + + if(LDT_rc%strat_cdfs.gt.0) then + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%strat_src,& + label="Stratification data source:", rc=rc) + call LDT_verify(rc,"Stratification data source: not defined") + + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%strat_cdfs_nbins,& + label="Number of bins to use for stratification:",rc=rc) + call LDT_verify(rc,"Number of bins to use for stratification: not defined") + + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%strat_file,& + label="External stratification file:",rc=rc) + call LDT_verify(rc,"External stratification file: not defined") + + allocate(LDT_rc%stratification_data(LDT_rc%ngrid(n),LDT_rc%cdf_ntimes)) + allocate(stratification_data(LDT_rc%lnc(n),LDT_rc%lnr(n),LDT_rc%cdf_ntimes)) + + call read_Precip_climo (LDT_rc%lnc(n), LDT_rc%lnr(n), LDT_rc%strat_file, stratification_data) ! + + do j=1,LDT_rc%cdf_ntimes + LDT_rc%strat_cdfs_min = 0. !minval returns -9999. minval(stratification_data(:,:,j)) ! min value over the entire domain + LDT_rc%strat_cdfs_max = maxval(stratification_data(:,:,j)) ! max value over the entire domain + delta = (LDT_rc%strat_cdfs_max-LDT_rc%strat_cdfs_min)/& + LDT_rc%strat_cdfs_nbins + do r=1,LDT_rc%lnr(n) + do c=1,LDT_rc%lnc(n) + if(LDT_domain(n)%gindex(c,r).ne.-1) then + if(stratification_data(c,r,j) .gt. 0) then + LDT_rc%stratification_data(LDT_domain(n)%gindex(c,r), j) = & + nint((stratification_data(c,r,j) - LDT_rc%strat_cdfs_min)/& + delta)+1 + if (LDT_rc%stratification_data(LDT_domain(n)%gindex(c,r), j) .gt. LDT_rc%strat_cdfs_nbins) then + write(LDT_logunit,*) '[INFO] Startification bins is larger then Max Startification bins',& + LDT_rc%stratification_data(LDT_domain(n)%gindex(c,r), j) , 'vs.', LDT_rc%strat_cdfs_nbins,& + 'Value adjusted to the Max Startification bins' + LDT_rc%stratification_data(LDT_domain(n)%gindex(c,r), j) = LDT_rc%strat_cdfs_nbins + endif + else + write(LDT_logunit,*) '[INFO] Total precip is zero or undefined',& + stratification_data(c,r,j) ,c,r,j,& + 'Value adjusted to the Min Startification bins' + LDT_rc%stratification_data(LDT_domain(n)%gindex(c,r), j) = 1 + endif + endif + enddo + enddo + enddo + endif endif - - if(LDT_rc%comp_obsgrid.eq.1) then + + if(LDT_rc%comp_obsgrid.eq.1) then LDT_rc%pass = 0 else - - if(LDT_rc%tavgInterval.lt.LDT_rc%ts) then + + if(LDT_rc%tavgInterval.lt.LDT_rc%ts) then write(LDT_logunit,*) '[ERR] Temporal averaging interval must be greater than' write(LDT_logunit,*) '[ERR] or equal to the LIS output timestep. ' write(LDT_logunit,*) '[ERR] Program stopping....' call LDT_endrun() endif - + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%obsCountThreshold,& label="Observation count threshold:",& rc=rc) - call LDT_verify(rc,'Observation count threshold: not defined') - - if(LDT_rc%comp_cdf.eq.1) then + call LDT_verify(rc,'Observation count threshold: not defined') + + if(LDT_rc%comp_cdf.eq.1) then LDT_rc%pass = 2 - elseif(LDT_rc%anomalyObsProc.eq.1) then + elseif(LDT_rc%anomalyObsProc.eq.1) then LDT_rc%pass = 2 endif - + endif call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%applyMask,& label="Apply external mask:",& rc=rc) call LDT_verify(rc,'Apply external mask: not defined') - + call ESMF_ConfigGetAttribute(LDT_config,LDT_rc%maskdir,& label="External mask directory:",& rc=rc) call LDT_verify(rc,'External mask directory: not defined') - + call LDT_DAobs_plugin - + LDT_rc%nobs = 1 - + ! This flag is set to prompt the writing of the domain file - if(LDT_rc%anomalyObsProc.eq.1) then - LDT_rc%comp_obsgrid = 1 + if(LDT_rc%anomalyObsProc.eq.1) then + LDT_rc%comp_obsgrid = 1 endif end subroutine LDT_DApreprocInit diff --git a/ldt/core/LDT_DrangeMod.F90 b/ldt/core/LDT_DrangeMod.F90 index 53f97360c..b08e3b18a 100644 --- a/ldt/core/LDT_DrangeMod.F90 +++ b/ldt/core/LDT_DrangeMod.F90 @@ -17,6 +17,7 @@ module LDT_DrangeMod ! ! !REVISION HISTORY: ! 2 Oct 2008 Sujay Kumar Initial Specification +! 16 Feb 2022 Mahdi Navari; modified to save stratify CDF ! !EOP use LDT_DAobsDataMod @@ -209,7 +210,7 @@ subroutine computeSingleDrange(n,obs, metrics) type(DAmetricsEntry) :: metrics integer :: t,i,j,k,c,r - integer :: sindex + integer :: sindex,sindex0,sindex1 real, allocatable :: strat_xrange(:,:,:,:) real, allocatable :: strat_delta(:,:,:) real, allocatable :: strat_mask(:,:,:) @@ -221,7 +222,7 @@ subroutine computeSingleDrange(n,obs, metrics) if(LDT_rc%endtime.eq.1) then if(obs%selectOpt.eq.1.and.metrics%selectOpt.eq.1) then - if(LDT_rc%group_cdfs.eq.0) then + if(LDT_rc%group_cdfs.eq.0 .and. LDT_rc%strat_cdfs.eq.0) then do t=1,LDT_rc%ngrid(n) do j=1,LDT_rc%cdf_ntimes do k=1,obs%vlevels @@ -246,7 +247,7 @@ subroutine computeSingleDrange(n,obs, metrics) enddo enddo enddo - elseif(LDT_rc%group_cdfs.eq.1) then + elseif(LDT_rc%group_cdfs.eq.1 .and. LDT_rc%strat_cdfs.eq.0) then allocate(strat_drange_total(LDT_rc%group_cdfs_nbins, & LDT_rc%cdf_ntimes, & obs%vlevels)) @@ -337,12 +338,219 @@ subroutine computeSingleDrange(n,obs, metrics) enddo enddo enddo + metrics%strat_xrange = strat_xrange + deallocate(strat_drange_total) + deallocate(strat_maxval) + deallocate(strat_minval) + deallocate(strat_mask) + deallocate(strat_delta) + +!MN: Startification based on monthly precipitation climatology +! monthly total precipitation climatology for each pixel are stored in LDT_rc%stratification_data + + elseif(LDT_rc%group_cdfs.eq.0 .and. LDT_rc%strat_cdfs.eq.1) then + allocate(strat_drange_total(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_maxval(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_minval(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_mask(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_delta(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_xrange(LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels, & + LDT_rc%cdf_nbins)) + + strat_drange_total = 0 + strat_maxval = -1000000.0 + strat_minval = 1000000.0 + strat_mask = 0 + strat_delta = 0 + strat_xrange = 0 + + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + sindex = LDT_rc%stratification_data(t,j) + strat_drange_total(sindex,j,k) = & + strat_drange_total(sindex,j,k) + & + metrics%count_drange_total(t,j,k) + if(metrics%maxval(t,j,k).gt.& + strat_maxval(sindex,j,k)) then + strat_maxval(sindex,j,k) = metrics%maxval(t,j,k) + endif + if(metrics%minval(t,j,k).lt.& + strat_minval(sindex,j,k)) then + strat_minval(sindex,j,k) = metrics%minval(t,j,k) + endif + enddo + enddo + enddo + + do t=1,LDT_rc%strat_cdfs_nbins + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + if(strat_drange_total(t,j,k).le.& + LDT_rc%obsCountThreshold) then + strat_maxval(t,j,k) = LDT_rc%udef + strat_minval(t,j,k) = LDT_rc%udef + strat_mask(t,j,k) = LDT_rc%udef + else + strat_mask(t,j,k) = strat_drange_total(t,j,k) + strat_delta(t,j,k) = & + (strat_maxval(t,j,k) - & + strat_minval(t,j,k))/& + (LDT_rc%cdf_nbins-1) + strat_xrange(t,j,k,1) = strat_minval(t,j,k) + do i=2, LDT_rc%cdf_nbins + strat_xrange(t,j,k,i) = & + strat_xrange(t,j,k,i-1) + & + strat_delta(t,j,k) + + enddo + endif + enddo + enddo + enddo + + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + sindex = LDT_rc%stratification_data(t,j) + if(strat_mask(sindex,j,k).eq.LDT_rc%udef) then + metrics%maxval(t,j,k) = LDT_rc%udef + metrics%minval(t,j,k) = LDT_rc%udef + metrics%mask(t,j,k) = LDT_rc%udef + else + metrics%mask(t,j,k) = strat_mask(sindex,j,k) + metrics%delta(t,j,k) = strat_delta(sindex,j,k) + metrics%xrange(t,j,k,:) =strat_xrange(sindex,j,k,:) + + metrics%maxval(t,j,k) = strat_maxval(sindex,j,k) + metrics%minval(t,j,k) = strat_minval(sindex,j,k) + endif + enddo + enddo + enddo + metrics%strat_xrange = strat_xrange + deallocate(strat_drange_total) + deallocate(strat_maxval) + deallocate(strat_minval) + deallocate(strat_mask) + deallocate(strat_delta) + deallocate(strat_xrange) + + elseif(LDT_rc%group_cdfs.eq.1 .and. LDT_rc%strat_cdfs.eq.1) then + allocate(strat_drange_total(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_maxval(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_minval(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_mask(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_delta(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels)) + allocate(strat_xrange(LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins, & + LDT_rc%cdf_ntimes, & + obs%vlevels, & + LDT_rc%cdf_nbins)) + + strat_drange_total = 0 + strat_maxval = -1000000.0 + strat_minval = 1000000.0 + strat_mask = 0 + strat_delta = 0 + strat_xrange = 0 + + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + sindex0 = LDT_rc%stratification_data(t,j) + sindex1 = LDT_rc%cdf_strat_data(t) + sindex = sindex0 + (sindex1 - 1)*LDT_rc%strat_cdfs_nbins + strat_drange_total(sindex,j,k) = & + strat_drange_total(sindex,j,k) + & + metrics%count_drange_total(t,j,k) + if(metrics%maxval(t,j,k).gt.& + strat_maxval(sindex,j,k)) then + strat_maxval(sindex,j,k) = metrics%maxval(t,j,k) + endif + if(metrics%minval(t,j,k).lt.& + strat_minval(sindex,j,k)) then + strat_minval(sindex,j,k) = metrics%minval(t,j,k) + endif + enddo + enddo + enddo + + do t=1,LDT_rc%group_cdfs_nbins*LDT_rc%strat_cdfs_nbins + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + if(strat_drange_total(t,j,k).le.& + LDT_rc%obsCountThreshold) then + strat_maxval(t,j,k) = LDT_rc%udef + strat_minval(t,j,k) = LDT_rc%udef + strat_mask(t,j,k) = LDT_rc%udef + else + strat_mask(t,j,k) = strat_drange_total(t,j,k) + strat_delta(t,j,k) = & + (strat_maxval(t,j,k) - & + strat_minval(t,j,k))/& + (LDT_rc%cdf_nbins-1) + strat_xrange(t,j,k,1) = strat_minval(t,j,k) + do i=2, LDT_rc%cdf_nbins + strat_xrange(t,j,k,i) = & + strat_xrange(t,j,k,i-1) + & + strat_delta(t,j,k) + + enddo + endif + enddo + enddo + enddo + do t=1,LDT_rc%ngrid(n) + do j=1,LDT_rc%cdf_ntimes + do k=1,obs%vlevels + sindex0 = LDT_rc%stratification_data(t,j) + sindex1 = LDT_rc%cdf_strat_data(t) + sindex = sindex0 + (sindex1 - 1)*LDT_rc%strat_cdfs_nbins + if(strat_mask(sindex,j,k).eq.LDT_rc%udef) then + metrics%maxval(t,j,k) = LDT_rc%udef + metrics%minval(t,j,k) = LDT_rc%udef + metrics%mask(t,j,k) = LDT_rc%udef + else + metrics%mask(t,j,k) = strat_mask(sindex,j,k) + metrics%delta(t,j,k) = strat_delta(sindex,j,k) + metrics%xrange(t,j,k,:) =strat_xrange(sindex,j,k,:) + + metrics%maxval(t,j,k) = strat_maxval(sindex,j,k) + metrics%minval(t,j,k) = strat_minval(sindex,j,k) + endif + enddo + enddo + enddo + metrics%strat_xrange = strat_xrange deallocate(strat_drange_total) deallocate(strat_maxval) deallocate(strat_minval) deallocate(strat_mask) deallocate(strat_delta) + deallocate(strat_xrange) endif endif diff --git a/ldt/core/LDT_PRIV_rcMod.F90 b/ldt/core/LDT_PRIV_rcMod.F90 index 6c419e1b5..318e20217 100644 --- a/ldt/core/LDT_PRIV_rcMod.F90 +++ b/ldt/core/LDT_PRIV_rcMod.F90 @@ -363,6 +363,18 @@ module LDT_PRIV_rcMod integer :: group_cdfs_nbins integer :: daily_interp_switch !0:on; 1:off (Y.Kwon) + integer :: strat_cdfs + integer :: write_strat_cdfs + character*50 :: strat_src + character*50 :: strat_file + !character*50 :: strat_cdfs_attrib_file + integer :: strat_cdfs_nbins + integer :: stratified_cdfs_nbins + real :: strat_cdfs_min + real :: strat_cdfs_max + integer, allocatable :: stratification_data(:,:) + + integer :: sp_sampl_cdfs integer :: sp_sample_cdf_rad @@ -377,6 +389,7 @@ module LDT_PRIV_rcMod integer :: pass_id integer :: ftn_cdf + integer :: ftn_strat_cdf integer :: ftn_DAobs_domain character*100 :: institution = 'NASA GSFC' diff --git a/ldt/make/Filepath b/ldt/make/Filepath index a5237a862..c8a408916 100644 --- a/ldt/make/Filepath +++ b/ldt/make/Filepath @@ -1 +1 @@ -dirs := . ../main ../core ../interp ../plugins ../lib/bil ../lib/xmrg ../params/mask ../params/landcover ../params/soils ../params/topo ../params/albedo ../params/gfrac ../params/LAISAI ../params/slopetype ../params/tbot ../params/pet ../params/climateparms ../params/crop_parms ../params/irrigation ../params/catchment ../params/CLM2 ../params/GeoWRSI ../params/FLAKE ../params/HYMAP ../params/JULES ../params/Mosaic ../params/Noah ../params/RDHM ../params/RUC ../params/SACHTET ../params/Snow17 ../params/SnowModel ../params/SiB2 ../params/VIC ../params/glacier ../runmodes/ANNproc ../runmodes/DApreproc ../runmodes/LSMparamproc ../runmodes/EnsRstpreproc ../runmodes/NUWRFpreproc ../runmodes/Metforcproc ../runmodes/MetTimeDScale ../runmodes/climoRstproc ../runmodes/rstTransformProc ../runmodes/StatDscaleMetForc ../runmodes/OPTUEproc/ ../domains/latlon ../domains/gaussian ../domains/hrap ../domains/lambert ../domains/polar ../domains/merc ../domains/easev2 ../DAobs/MCD15A2H_LAI ../DAobs/LISlsmSMobs ../DAobs/LISlsmTEFFobs ../DAobs/synthetic_sm ../DAobs/NASA_AMSRE_sm ../DAobs/LPRM_AMSREsm ../DAobs/LPRMvod ../DAobs/ESACCI_sm ../DAobs/WindSat_sm ../DAobs/GRACE_tws ../DAobs/GRACEQL_tws ../DAobs/SMOPS ../DAobs/ASCAT_TUW ../DAobs/SMOS_L2sm ../DAobs/SMOS_NESDIS ../DAobs/GCOMW_AMSR2L3snd ../DAobs/GCOMW_AMSR2L3sm ../DAobs/Aquarius_L2sm ../DAobs/simGRACE_JPL ../DAobs/SMMR_SNWD ../DAobs/SSMI_SNWD ../DAobs/ANSA_SNWD ../DAobs/NASA_SMAPsm ../DAobs/SMOS_NRTNN_L2sm ../DAobs/NASA_SMAPvod ../DAobs/GLASSlai ../DAobs/THySM ../DAobs/GEOS_FP_TEFFobs ../DAobs/SMAP_E_OPLsm ../DAobs/VIIRS_GVF ../DAobs/CDFS_GVF ../RESTRICTED/usaf/usaf ../metforcing/agrradps ../metforcing/cmap ../metforcing/cmorph ../metforcing/ecmwf ../metforcing/gdas ../metforcing/geos5fcst ../metforcing/gfs ../metforcing/gldas ../metforcing/gswp1 ../metforcing/gswp2 ../metforcing/merra2 ../metforcing/era5 ../metforcing/nam242 ../metforcing/narr ../metforcing/princeton ../metforcing/nldas2 ../metforcing/WRFoutv2 ../metforcing/WRFAKdom ../metforcing/stg2 ../metforcing/stg4 ../metforcing/3B42RTV7 ../metforcing/3B42V6 ../metforcing/3B42V7 ../metforcing/RFE2Daily ../metforcing/RFE2gdas ../metforcing/chirps2 ../MetforcScale ../ANNdata/LISlsmSM ../ANNdata/synthetic_sm ../ANNdata/LPRM_AMSREsm ../ANNdata/MODIS_LST ../ANNdata/MOD10A1 ../ANNdata/GHCNsnwd ../ANNdata/GCOMW_AMSR2_TB ../statDscale/BayesianMerging ../statDscale/Climo ../params/CLM45 ../USAFSI ../SMAP_E_OPL ../runmodes/USAFSI ../runmodes/LISHydropreproc ../runmodes/SMAP_E_OPL ../params/Crocus ../runmodes/obsSim ../obsSim/NatureRun/LISout/ ../obsSim/OSSEmask/LISout/ ../obsSim/OSSEmask/AMSR2/ ../obsSim/OSSEmask/TSMM/ ../obsSim/OSSEmask/MODIS/ ../obsSim/OSSEmask/Sentinel1A/ +dirs := . ../main ../core ../interp ../plugins ../lib/bil ../lib/xmrg ../params/mask ../params/landcover ../params/soils ../params/topo ../params/albedo ../params/gfrac ../params/LAISAI ../params/slopetype ../params/tbot ../params/pet ../params/climateparms ../params/crop_parms ../params/irrigation ../params/catchment ../params/CLM2 ../params/GeoWRSI ../params/FLAKE ../params/HYMAP ../params/JULES ../params/Mosaic ../params/Noah ../params/RDHM ../params/RUC ../params/SACHTET ../params/Snow17 ../params/SnowModel ../params/SiB2 ../params/VIC ../params/glacier ../runmodes/ANNproc ../runmodes/DApreproc ../runmodes/LSMparamproc ../runmodes/EnsRstpreproc ../runmodes/NUWRFpreproc ../runmodes/Metforcproc ../runmodes/MetTimeDScale ../runmodes/climoRstproc ../runmodes/rstTransformProc ../runmodes/StatDscaleMetForc ../runmodes/OPTUEproc/ ../domains/latlon ../domains/gaussian ../domains/hrap ../domains/lambert ../domains/polar ../domains/merc ../domains/easev2 ../DAobs/MCD15A2H_LAI ../DAobs/LISlsmSMobs ../DAobs/LISlsmPrecipobs ../DAobs/LISlsmTEFFobs ../DAobs/synthetic_sm ../DAobs/NASA_AMSRE_sm ../DAobs/LPRM_AMSREsm ../DAobs/LPRMvod ../DAobs/ESACCI_sm ../DAobs/WindSat_sm ../DAobs/GRACE_tws ../DAobs/GRACEQL_tws ../DAobs/SMOPS ../DAobs/ASCAT_TUW ../DAobs/SMOS_L2sm ../DAobs/SMOS_NESDIS ../DAobs/GCOMW_AMSR2L3snd ../DAobs/GCOMW_AMSR2L3sm ../DAobs/Aquarius_L2sm ../DAobs/simGRACE_JPL ../DAobs/SMMR_SNWD ../DAobs/SSMI_SNWD ../DAobs/ANSA_SNWD ../DAobs/NASA_SMAPsm ../DAobs/SMOS_NRTNN_L2sm ../DAobs/NASA_SMAPvod ../DAobs/GLASSlai ../DAobs/THySM ../DAobs/GEOS_FP_TEFFobs ../DAobs/SMAP_E_OPLsm ../DAobs/VIIRS_GVF ../DAobs/CDFS_GVF ../RESTRICTED/usaf/usaf ../metforcing/agrradps ../metforcing/cmap ../metforcing/cmorph ../metforcing/ecmwf ../metforcing/gdas ../metforcing/geos5fcst ../metforcing/gfs ../metforcing/gldas ../metforcing/gswp1 ../metforcing/gswp2 ../metforcing/merra2 ../metforcing/era5 ../metforcing/nam242 ../metforcing/narr ../metforcing/princeton ../metforcing/nldas2 ../metforcing/WRFoutv2 ../metforcing/WRFAKdom ../metforcing/stg2 ../metforcing/stg4 ../metforcing/3B42RTV7 ../metforcing/3B42V6 ../metforcing/3B42V7 ../metforcing/RFE2Daily ../metforcing/RFE2gdas ../metforcing/chirps2 ../MetforcScale ../ANNdata/LISlsmSM ../ANNdata/synthetic_sm ../ANNdata/LPRM_AMSREsm ../ANNdata/MODIS_LST ../ANNdata/MOD10A1 ../ANNdata/GHCNsnwd ../ANNdata/GCOMW_AMSR2_TB ../statDscale/BayesianMerging ../statDscale/Climo ../params/CLM45 ../USAFSI ../SMAP_E_OPL ../runmodes/USAFSI ../runmodes/LISHydropreproc ../runmodes/SMAP_E_OPL ../params/Crocus ../runmodes/obsSim ../obsSim/NatureRun/LISout/ ../obsSim/OSSEmask/LISout/ ../obsSim/OSSEmask/AMSR2/ ../obsSim/OSSEmask/TSMM/ ../obsSim/OSSEmask/MODIS/ ../obsSim/OSSEmask/Sentinel1A/ diff --git a/ldt/plugins/LDT_DAobs_pluginMod.F90 b/ldt/plugins/LDT_DAobs_pluginMod.F90 index 93b91724d..105b70d7c 100644 --- a/ldt/plugins/LDT_DAobs_pluginMod.F90 +++ b/ldt/plugins/LDT_DAobs_pluginMod.F90 @@ -72,6 +72,7 @@ subroutine LDT_DAobs_plugin use LPRMvod_obsMod, only : LPRMvod_obsinit use MCD15A2Hlai_obsMod, only : MCD15A2Hlai_obsinit use THySM_obsMod, only : THySM_obsinit + use LISlsmPrecip_obsMod, only : LISlsmPrecip_obsInit use VIIRSGVFobsMod, only : VIIRSGVFobsinit !Y.Kwon use CDFSGVFobsMod, only : CDFSGVFobsinit !Y.Kwon use GEOSTEFF_obsMod, only : GEOSTeffobsinit !Y.Kwon @@ -104,6 +105,7 @@ subroutine LDT_DAobs_plugin external readLPRMvodObs external readMCD15A2HlaiObs external readTHySMobs + external readLISlsmPrecipObs external readVIIRS_GVFObs !Y.Kwon external readCDFS_GVFObs !Y.Kwon external readGEOSTEFFObs !Y.Kwon @@ -112,6 +114,9 @@ subroutine LDT_DAobs_plugin call registerdaobssetup(trim(LDT_LISlsmSMobsId)//char(0), LISlsmSM_obsInit) call registerdaobsread(trim(LDT_LISlsmSMobsId)//char(0), readLISlsmSMObs) + call registerdaobssetup(trim(LDT_LISlsmPrecipobsId)//char(0), LISlsmPrecip_obsInit) + call registerdaobsread(trim(LDT_LISlsmPrecipobsId)//char(0), readLISlsmPrecipObs) + !Y.Kwon call registerdaobssetup(trim(LDT_LISlsmTEFFobsId)//char(0), LISlsmTEFF_obsInit) call registerdaobsread(trim(LDT_LISlsmTEFFobsId)//char(0), readLISlsmTEFFObs) diff --git a/ldt/plugins/LDT_pluginIndices.F90 b/ldt/plugins/LDT_pluginIndices.F90 index eda68dcff..3d2c4f8eb 100644 --- a/ldt/plugins/LDT_pluginIndices.F90 +++ b/ldt/plugins/LDT_pluginIndices.F90 @@ -120,9 +120,10 @@ module LDT_pluginIndices = "GLASS LAI" character*50, public, parameter :: LDT_LPRMvodobsId & = "LPRM vegetation optical depth" - character*50, public, parameter :: LDT_MCD15A2HlaiobsId & = "MCD15A2H LAI" + character*50, public, parameter :: LDT_LISlsmPrecipobsId & + = "LIS LSM total precipitation" character*50, public, parameter :: LDT_VIIRSgvfobsId & = "VIIRS GVF" !Y.Kwon character*50, public, parameter :: LDT_CDFSgvfobsId & diff --git a/ldt/testcases/DAobs_LSM_pclimo_transfer/README b/ldt/testcases/DAobs_LSM_pclimo_transfer/README new file mode 100644 index 000000000..8a261ad3b --- /dev/null +++ b/ldt/testcases/DAobs_LSM_pclimo_transfer/README @@ -0,0 +1,15 @@ +LSM Soil Moisture Test Case for Precip Climo and Relocating CDFs + +This testcase processes LIS LSM data (Noah.4.0.1) to generate a CDF. + +This directory contains: + +* This README file. +* The ldt.config file used for this test case. (This file should be + edited to make sure that the locations of the parameter and forcing + files are specified correctly.) + +To run this test case: +* Generate the LDT executable. +* Run the LDT executable using the ldt.config file and the sample input + data. diff --git a/ldt/testcases/DAobs_LSM_pclimo_transfer/ldt.config b/ldt/testcases/DAobs_LSM_pclimo_transfer/ldt.config new file mode 100755 index 000000000..0ed463e6e --- /dev/null +++ b/ldt/testcases/DAobs_LSM_pclimo_transfer/ldt.config @@ -0,0 +1,308 @@ +# == LDT Main Entry Options == + +LDT running mode: "DA preprocessing" # LDT type of run-mode (top-level option) +Processed LSM parameter filename: ../procLSM/lis_input.d01.nc # Final output file read by LIS-7 + +LIS number of nests: 1 # Total number of nests run by LIS +Number of surface model types: 1 # Total number of desired surface model types +Surface model types: "LSM" # Surface models: LSM | Openwater +Land surface model: "Noah-MP.4.0.1" # Enter LSM(s) of choice +Routing model: "none" # "HYMAP" # "none" +Lake model: "none" # Enter Lake model(s) of choice +Water fraction cutoff value: 0.5 # Fraction at which gridcell is designated as 'water' +#Glacier fraction cutoff value: 0.2 + +Number of met forcing sources: 0 # Enter number of forcing types +Met forcing sources: "NLDAS2" # Enter 'none' if no forcing selected +Met spatial transform methods: neighbor # bilinear | budget-bilinear | neighbor | average +Topographic correction method (met forcing): "none" # none | lapse-rate + +LDT diagnostic file: ldtlog.strat.p.climo.15bins # Log-based diagnostic output file +Undefined value: -9999.0 # Universal undefined value +LDT output directory: LDTOUT.strat.p.climo.15bins # If metrics or stats are written out +Number of ensembles per tile: 1 # The number of ensemble members per tile + +#LIS domain: +Map projection of the LIS domain: latlon +Run domain lower left lat: 25.0625 +Run domain lower left lon: -124.9375 +Run domain upper right lat: 52.9375 +Run domain upper right lon: -67.0625 +Run domain resolution (dx): 0.125 +Run domain resolution (dy): 0.125 + +# == Landcover, Landmask and Soil Texture Parameters == + +#Landcover/Mask Parameter Inputs # +Landcover data source: "AVHRR" #"MODIS_Native" # +Landcover classification: "UMD" # "IGBPNCEP" # # Enter land cover classification type +Landcover file: ../input/LS_PARAMETERS/NLDAS_0.125/umdveg_dom_nldas.1gd4r # ./input/LS_PARAMETERS/noah_2dparms/igbp.bin # # Landcover map path +Landcover spatial transform: none +Landcover fill option: neighbor +Landcover fill radius: 5 +Landcover fill value: 4 +Landcover map projection: latlon +Landcover lower left lat: 25.0625 +Landcover lower left lon: -124.9375 +Landcover upper right lat: 52.9375 +Landcover upper right lon: -67.0625 +Landcover resolution (dx): 0.125 +Landcover resolution (dy): 0.125 + +Create or readin landmask: "readin" +Landmask data source: "AVHRR" # If 'created', recommended to put Landcover source name here +Landmask file: ../input/LS_PARAMETERS/irrigation/conus_modis/IRR_N125_mask.1gd4r +Landmask spatial transform: none +Landmask map projection: latlon +Landmask lower left lat: 25.0625 +Landmask lower left lon: -124.9375 +Landmask upper right lat: 52.9375 +Landmask upper right lon: -67.0625 +Landmask resolution (dx): 0.125 +Landmask resolution (dy): 0.125 + + +#Soil texture map: +Soil texture data source: ISRIC # STATSGOFAO_LIS +Soil texture map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/TEXMHT_M_sl1_250m.tif # ./input/LS_PARAMETERS/NLDAS_0.125/soil_texture_statsgo_nldas.1gd4r # Enter soil texture map +Soil texture spatial transform: mode # none # none | mode | neighbor | tile +Soil texture fill option: neighbor # none | neighbor +Soil texture fill radius: 5 # Number of pixels to search for neighbor +Soil texture fill value: 6 # Static value to fill where missing +Soil texture map projection: latlon +Soil texture lower left lat: 25.0625 +Soil texture lower left lon: -124.9375 +Soil texture upper right lat: 52.9375 +Soil texture upper right lon: -67.0625 +Soil texture resolution (dx): 0.125 +Soil texture resolution (dy): 0.125 + + +# Soil fraction +#Porosity map: # ../LS_PARAMETERS/UMD/25KM/porosity_FAO +Porosity data source: CONSTANT +Porosity fill value: 0.32 # Porosity fill is currently locked with Soils fill options + +Soil fraction data source: none # ISRIC # none +Soil fraction number of bands: 1 +Sand fraction map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/SNDPPT_M_sl1_250m_ll.tif +Clay fraction map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/CLYPPT_M_sl1_250m_ll.tif +Silt fraction map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/SLTPPT_M_sl1_250m_ll.tif +Soils spatial transform: average +Soils fill option: neighbor # none +Soils fill radius: 1 +Soils fill value: 0.333 +Sand fraction fill value: 0.30 +Clay fraction fill value: 0.35 +Silt fraction fill value: 0.35 + +Soils map projection: latlon +Soils lower left lat: -56.008104 +Soils lower left lon: -180.00 +Soils upper right lat: 83.9991672 +Soils upper right lon: 179.9999424 +Soils resolution (dx): 0.0020833 +Soils resolution (dy): 0.0020833 + + +#Topography maps +#Elevation data source: GTOPO30_Native +#Elevation number of bands: 1 +#Elevation map: ./input/LS_PARAMETERS/topo_parms/GTOPO30/raw_updated +#Elevation fill option: neighbor +#Elevation fill radius: 3 +#Elevation fill value: 287 +# SRTM Elevation data entries: +Elevation data source: "SRTM_LIS" +Elevation map: ../input/LS_PARAMETERS/topo_parms/SRTM/SRTM30/srtm_elev1km.1gd4r +Elevation number of bands: 1 +Slope data source: "SRTM_LIS" +Slope map: ../input/LS_PARAMETERS/topo_parms/SRTM/SRTM30/srtm_slope1km.1gd4r +Slope number of bands: 1 +Aspect data source: "SRTM_LIS" +Aspect map: ../input/LS_PARAMETERS/topo_parms/SRTM/SRTM30/srtm_aspect1km.1gd4r +Aspect number of bands: 1 + + +Topography spatial transform: average +Elevation fill option: average +Elevation fill radius: 5 +Elevation fill value: 0 +Slope fill option: average +Slope fill radius: 5 +Slope fill value: 0.1 +Aspect fill option: average +Aspect fill radius: 5 +Aspect fill value: 0 +Topography map projection: latlon +Topography lower left lat: -59.995 +Topography lower left lon: -179.995 +Topography upper right lat: 89.995 +Topography upper right lon: 179.995 +Topography resolution (dx): 0.01 +Topography resolution (dy): 0.01 + + +# Albedo maps: +Albedo data source: NCEP_Native # NCEP_LIS +Albedo map: ../input/LS_PARAMETERS/noah_2dparms/albedo # ./input/LS_PARAMETERS/NLDAS_0.125/albedo_nldas # Albedo files +Albedo climatology interval: monthly # monthly | quarterly +Albedo spatial transform: bilinear # average | neighbor | bilinear | budget-bilinear +Albedo fill option: neighbor # none | neighbor | average +Albedo fill radius: 1 # Number of pixels to search for neighbor +Albedo fill value: 0.14 # Static value to fill where missing +Albedo map projection: latlon +Albedo lower left lat: 25.0625 +Albedo lower left lon: -124.9375 +Albedo upper right lat: 52.9375 +Albedo upper right lon: -67.0625 +Albedo resolution (dx): 0.125 +Albedo resolution (dy): 0.12.5 + +Max snow albedo data source: NCEP_Native # NCEP_LIS +Max snow albedo map: ../input/LS_PARAMETERS/noah_2dparms/maxsnoalb.asc #./input/LS_PARAMETERS/NLDAS_0.125/maxsnalb_nldas.1gd4r # Max. snow albedo map +Max snow albedo spatial transform: budget-bilinear # none # average | neighbor | bilinear | budget-bilinear +Max snow albedo fill option: neighbor # none | neighbor | average +Max snow albedo fill radius: 3 # Number of pixels to search for neighbor +Max snow albedo fill value: 0.65 # Static value to fill where missing +Max snow albedo map projection: latlon +Max snow albedo lower left lat: 25.0625 +Max snow albedo lower left lon: -124.9375 +Max snow albedo upper right lat: 52.9375 +Max snow albedo upper right lon: -67.0625 +Max snow albedo resolution (dx): 0.125 +Max snow albedo resolution (dy): 0.125 + +# Greenness fraction maps: +Greenness data source: NCEP_Native #NCEP_LIS +Greenness fraction map: ../input/LS_PARAMETERS/noah_2dparms/gfrac # ./input/LS_PARAMETERS/NLDAS_0.125/gfrac_nldas # Greenness fraction map +Greenness climatology interval: monthly # monthly +Calculate min-max greenness fraction: .false. +Greenness maximum map: ../input/LS_PARAMETERS/noah_2dparms/gfrac_max.asc # ./input/LS_PARAMETERS/NLDAS_0.125/gfrac_nldas.MAX.1gd4r # Maximum greenness fraction map +Greenness minimum map: ../input/LS_PARAMETERS/noah_2dparms/gfrac_min.asc # ./input/LS_PARAMETERS/NLDAS_0.125/gfrac_nldas.MIN.1gd4r # Minimum greenness fraction map +Greenness spatial transform: bilinear #none # average | neighbor | bilinear | budget-bilinear +Greenness fill option: neighbor # none | neighbor | average +Greenness fill radius: 1 # Number of pixels to search for neighbor +Greenness fill value: 0.30 # Static value to fill where missing +Greenness maximum fill value: 0.40 # Static value to fill where missing +Greenness minimum fill value: 0.20 # Static value to fill where missing +Greenness map projection: latlon +Greenness lower left lat: 25.0625 +Greenness lower left lon: -124.9375 +Greenness upper right lat: 52.9375 +Greenness upper right lon: -67.0625 +Greenness resolution (dx): 0.125 +Greenness resolution (dy): 0.125 + +# Slope type map: +Slope type data source: NCEP_Native # NCEP_LIS +Slope type map: ../input/LS_PARAMETERS/noah_2dparms/islope #./input/LS_PARAMETERS/NLDAS_0.125/noah_slope_nldas.1gd4r # Slope type map +Slope type spatial transform: neighbor # none # none | neighbor | mode +Slope type fill option: neighbor # none | neighbor +Slope type fill radius: 1 # Number of pixels to search for neighbor +Slope type fill value: 3. # Static value to fill where missing +Slope type map projection: latlon +Slope type lower left lat: 25.0625 +Slope type lower left lon: -124.9375 +Slope type upper right lat: 52.9375 +Slope type upper right lon: -67.0625 +Slope type resolution (dx): 0.125 +Slope type resolution (dy): 0.125 + +# Bottom temperature map (lapse-rate correction option): +Bottom temperature data source: NCEP_LIS +Bottom temperature map: ../input/LS_PARAMETERS/NLDAS_0.125/tbot_nldas.1gd4r +Bottom temperature spatial transform: none # none | average | neighbor | bilinear | budget-bilinear +Bottom temperature fill option: none # none | average | neighbor +Bottom temperature fill radius: 1 # Number of pixels to search for neighbor +Bottom temperature fill value: 287. # Static value to fill where missing +Bottom temperature topographic downscaling: "none" # none | lapse-rate +Bottom temperature map projection: latlon # Projection type +Bottom temperature lower left lat: 25.0625 +Bottom temperature lower left lon: -124.9375 +Bottom temperature upper right lat: 52.9375 +Bottom temperature upper right lon: -67.0625 +Bottom temperature resolution (dx): 0.125 +Bottom temperature resolution (dy): 0.125 + +# ======================================================= +#Forcing elevation +NLDAS2 forcing directory: /discover/nobackup/projects/lis/MET_FORCING/NLDAS2.FORCING +NLDAS2 data center source: "GES-DISC" +NLDAS2 use model level data: 0 +NLDAS2 use model based swdown: 0 +NLDAS2 use model based precip: 0 +NLDAS2 use model based pressure: 0 +Forcing variables list file: ./ldt_forcing_vars.txt + +NLDAS2 elevation difference map: ../input/LS_PARAMETERS/NLDAS_0.125/NARR_elev-diff.1gd4r +NARR terrain height map: ../input/LS_PARAMETERS/NLDAS_0.125/NARR_elevation.1gd4r + +#Noah-MP LSM inputs +Noah-MP PBL Height Value: 900. + +# ----------------------------------------------------------------------------- +# Model-based CDF file generation: +#------------------------------------------------------------------------------ +DA preprocessing method: "CDF generation" +DA observation source: "LIS LSM soil moisture" # "LSM soil moisture" +Name of the preprocessed DA file: "cdf_noahmp401" +Apply anomaly correction to obs: 0 +Temporal resolution of CDFs: "monthly" # yearly | monthly +Number of bins to use in the CDF: 100 +Observation count threshold: 30 +Temporal averaging interval: "1da" +Apply external mask: 0 +External mask directory: none + +LIS soil moisture output model name: "Noah.4.0.1" +LIS soil moisture output directory: ../sim/output.ol/ +LIS soil moisture output format: "netcdf" +LIS soil moisture output methodology: "2d gridspace" +LIS soil moisture output naming style: "3 level hierarchy" +LIS soil moisture output nest index: 1 +LIS soil moisture output map projection: "latlon" +LIS soil moisture domain lower left lat: 25.0625 +LIS soil moisture domain upper right lat: 52.9375 +LIS soil moisture domain lower left lon: -124.9375 +LIS soil moisture domain upper right lon: -67.0625 +LIS soil moisture domain resolution (dx): 0.125 +LIS soil moisture domain resolution (dy): 0.125 + + +Starting year: 2015 +Starting month: 03 +Starting day: 01 +Starting hour: 0 +Starting minute: 0 +Starting second: 0 +Ending year: 2022 +Ending month: 01 +Ending day: 01 +Ending hour: 0 +Ending minute: 0 +Ending second: 0 +LIS output timestep: "6hr" + + +Maximum number of surface type tiles per grid: 1 +Minimum cutoff percentage (surface type tiles): 0.05 +Maximum number of soil texture tiles per grid: 1 +Minimum cutoff percentage (soil texture tiles): 0.05 +Maximum number of soil fraction tiles per grid: 1 +Minimum cutoff percentage (soil fraction tiles): 0.05 +Maximum number of elevation bands per grid: 1 +Minimum cutoff percentage (elevation bands): 0.05 +Maximum number of slope bands per grid: 1 +Minimum cutoff percentage (slope bands): 0.05 +Maximum number of aspect bands per grid: 1 +Minimum cutoff percentage (aspect bands): 0.05 + +# ----------------------------------------------------------------------------- +# Stratification based on precipitation: +#------------------------------------------------------------------------------ +Stratify CDFs by external data: 1 +Number of bins to use for stratification: 15 +Stratification data source: "LIS LSM total precipitation" +External stratification file: "LVT_MEAN_FINAL.202201010000.d01.nc" # "../Precip_climatology/Precip.climo.us.nldas2/LVT_MEAN_FINAL.202109010000.d01.nc" +Write stratified geolocation independent CDFs: 1 diff --git a/ldt/testcases/DAobs_SMAP_pclimo_transfer/README b/ldt/testcases/DAobs_SMAP_pclimo_transfer/README new file mode 100644 index 000000000..aca2520d4 --- /dev/null +++ b/ldt/testcases/DAobs_SMAP_pclimo_transfer/README @@ -0,0 +1,15 @@ +SMAP Soil Moisture Test Case for Precip Climo and Relocating CDFs + +This testcase processes SMAP data to generate a CDF. + +This directory contains: + +* This README file. +* The ldt.config file used for this test case. (This file should be + edited to make sure that the locations of the parameter and forcing + files are specified correctly.) + +To run this test case: +* Generate the LDT executable. +* Run the LDT executable using the ldt.config file and the sample input + data. diff --git a/ldt/testcases/DAobs_SMAP_pclimo_transfer/README~ b/ldt/testcases/DAobs_SMAP_pclimo_transfer/README~ new file mode 100644 index 000000000..8a261ad3b --- /dev/null +++ b/ldt/testcases/DAobs_SMAP_pclimo_transfer/README~ @@ -0,0 +1,15 @@ +LSM Soil Moisture Test Case for Precip Climo and Relocating CDFs + +This testcase processes LIS LSM data (Noah.4.0.1) to generate a CDF. + +This directory contains: + +* This README file. +* The ldt.config file used for this test case. (This file should be + edited to make sure that the locations of the parameter and forcing + files are specified correctly.) + +To run this test case: +* Generate the LDT executable. +* Run the LDT executable using the ldt.config file and the sample input + data. diff --git a/ldt/testcases/DAobs_SMAP_pclimo_transfer/ldt.config b/ldt/testcases/DAobs_SMAP_pclimo_transfer/ldt.config new file mode 100755 index 000000000..456ea0382 --- /dev/null +++ b/ldt/testcases/DAobs_SMAP_pclimo_transfer/ldt.config @@ -0,0 +1,316 @@ +# == LDT Main Entry Options == + +LDT running mode: "DA preprocessing" # LDT type of run-mode (top-level option) +Processed LSM parameter filename: ../procLSM/lis_input.d01.nc # Final output file read by LIS-7 + +LIS number of nests: 1 # Total number of nests run by LIS +Number of surface model types: 1 # Total number of desired surface model types +Surface model types: "LSM" # Surface models: LSM | Openwater +Land surface model: "Noah-MP.4.0.1" # Enter LSM(s) of choice +Routing model: "none" # "HYMAP" # "none" +Lake model: "none" # Enter Lake model(s) of choice +Water fraction cutoff value: 0.5 # Fraction at which gridcell is designated as 'water' +#Glacier fraction cutoff value: 0.2 + +Number of met forcing sources: 0 # Enter number of forcing types +Met forcing sources: "NLDAS2" # Enter 'none' if no forcing selected +Met spatial transform methods: neighbor # bilinear | budget-bilinear | neighbor | average +Topographic correction method (met forcing): "none" # none | lapse-rate + +LDT diagnostic file: ldtlog.strat.p.climo.15bins # Log-based diagnostic output file +Undefined value: -9999.0 # Universal undefined value +LDT output directory: LDTOUT.strat.p.climo.15bins # If metrics or stats are written out +Number of ensembles per tile: 1 # The number of ensemble members per tile + +#LIS domain: +Map projection of the LIS domain: latlon +Run domain lower left lat: 25.0625 +Run domain lower left lon: -124.9375 +Run domain upper right lat: 52.9375 +Run domain upper right lon: -67.0625 +Run domain resolution (dx): 0.125 +Run domain resolution (dy): 0.125 + +# == Landcover, Landmask and Soil Texture Parameters == + +#Landcover/Mask Parameter Inputs # +Landcover data source: "AVHRR" #"MODIS_Native" # +Landcover classification: "UMD" # "IGBPNCEP" # # Enter land cover classification type +Landcover file: ../input/LS_PARAMETERS/NLDAS_0.125/umdveg_dom_nldas.1gd4r # ./input/LS_PARAMETERS/noah_2dparms/igbp.bin # # Landcover map path +Landcover spatial transform: none +Landcover fill option: neighbor +Landcover fill radius: 5 +Landcover fill value: 4 +Landcover map projection: latlon +Landcover lower left lat: 25.0625 +Landcover lower left lon: -124.9375 +Landcover upper right lat: 52.9375 +Landcover upper right lon: -67.0625 +Landcover resolution (dx): 0.125 +Landcover resolution (dy): 0.125 + +Create or readin landmask: "readin" +Landmask data source: "AVHRR" # If 'created', recommended to put Landcover source name here +Landmask file: ../input/LS_PARAMETERS/irrigation/conus_modis/IRR_N125_mask.1gd4r +Landmask spatial transform: none +Landmask map projection: latlon +Landmask lower left lat: 25.0625 +Landmask lower left lon: -124.9375 +Landmask upper right lat: 52.9375 +Landmask upper right lon: -67.0625 +Landmask resolution (dx): 0.125 +Landmask resolution (dy): 0.125 + + +#Soil texture map: +Soil texture data source: ISRIC # STATSGOFAO_LIS +Soil texture map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/TEXMHT_M_sl1_250m.tif # ./input/LS_PARAMETERS/NLDAS_0.125/soil_texture_statsgo_nldas.1gd4r # Enter soil texture map +Soil texture spatial transform: mode # none # none | mode | neighbor | tile +Soil texture fill option: neighbor # none | neighbor +Soil texture fill radius: 5 # Number of pixels to search for neighbor +Soil texture fill value: 6 # Static value to fill where missing +Soil texture map projection: latlon +Soil texture lower left lat: 25.0625 +Soil texture lower left lon: -124.9375 +Soil texture upper right lat: 52.9375 +Soil texture upper right lon: -67.0625 +Soil texture resolution (dx): 0.125 +Soil texture resolution (dy): 0.125 + + +# Soil fraction +#Porosity map: # ../LS_PARAMETERS/UMD/25KM/porosity_FAO +Porosity data source: CONSTANT +Porosity fill value: 0.32 # Porosity fill is currently locked with Soils fill options + +Soil fraction data source: none # ISRIC # none +Soil fraction number of bands: 1 +Sand fraction map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/SNDPPT_M_sl1_250m_ll.tif +Clay fraction map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/CLYPPT_M_sl1_250m_ll.tif +Silt fraction map: ../input/LS_PARAMETERS/soil_parms/ISRIC/v2017/SLTPPT_M_sl1_250m_ll.tif +Soils spatial transform: average +Soils fill option: neighbor # none +Soils fill radius: 1 +Soils fill value: 0.333 +Sand fraction fill value: 0.30 +Clay fraction fill value: 0.35 +Silt fraction fill value: 0.35 + +Soils map projection: latlon +Soils lower left lat: -56.008104 +Soils lower left lon: -180.00 +Soils upper right lat: 83.9991672 +Soils upper right lon: 179.9999424 +Soils resolution (dx): 0.0020833 +Soils resolution (dy): 0.0020833 + + +#Topography maps +#Elevation data source: GTOPO30_Native +#Elevation number of bands: 1 +#Elevation map: ./input/LS_PARAMETERS/topo_parms/GTOPO30/raw_updated +#Elevation fill option: neighbor +#Elevation fill radius: 3 +#Elevation fill value: 287 +# SRTM Elevation data entries: +Elevation data source: "SRTM_LIS" +Elevation map: ../input/LS_PARAMETERS/topo_parms/SRTM/SRTM30/srtm_elev1km.1gd4r +Elevation number of bands: 1 +Slope data source: "SRTM_LIS" +Slope map: ../input/LS_PARAMETERS/topo_parms/SRTM/SRTM30/srtm_slope1km.1gd4r +Slope number of bands: 1 +Aspect data source: "SRTM_LIS" +Aspect map: ../input/LS_PARAMETERS/topo_parms/SRTM/SRTM30/srtm_aspect1km.1gd4r +Aspect number of bands: 1 + + +Topography spatial transform: average +Elevation fill option: average +Elevation fill radius: 5 +Elevation fill value: 0 +Slope fill option: average +Slope fill radius: 5 +Slope fill value: 0.1 +Aspect fill option: average +Aspect fill radius: 5 +Aspect fill value: 0 +Topography map projection: latlon +Topography lower left lat: -59.995 +Topography lower left lon: -179.995 +Topography upper right lat: 89.995 +Topography upper right lon: 179.995 +Topography resolution (dx): 0.01 +Topography resolution (dy): 0.01 + + +# Albedo maps: +Albedo data source: NCEP_Native # NCEP_LIS +Albedo map: ../input/LS_PARAMETERS/noah_2dparms/albedo # ./input/LS_PARAMETERS/NLDAS_0.125/albedo_nldas # Albedo files +Albedo climatology interval: monthly # monthly | quarterly +Albedo spatial transform: bilinear # average | neighbor | bilinear | budget-bilinear +Albedo fill option: neighbor # none | neighbor | average +Albedo fill radius: 1 # Number of pixels to search for neighbor +Albedo fill value: 0.14 # Static value to fill where missing +Albedo map projection: latlon +Albedo lower left lat: 25.0625 +Albedo lower left lon: -124.9375 +Albedo upper right lat: 52.9375 +Albedo upper right lon: -67.0625 +Albedo resolution (dx): 0.125 +Albedo resolution (dy): 0.12.5 + +Max snow albedo data source: NCEP_Native # NCEP_LIS +Max snow albedo map: ../input/LS_PARAMETERS/noah_2dparms/maxsnoalb.asc #./input/LS_PARAMETERS/NLDAS_0.125/maxsnalb_nldas.1gd4r # Max. snow albedo map +Max snow albedo spatial transform: budget-bilinear # none # average | neighbor | bilinear | budget-bilinear +Max snow albedo fill option: neighbor # none | neighbor | average +Max snow albedo fill radius: 3 # Number of pixels to search for neighbor +Max snow albedo fill value: 0.65 # Static value to fill where missing +Max snow albedo map projection: latlon +Max snow albedo lower left lat: 25.0625 +Max snow albedo lower left lon: -124.9375 +Max snow albedo upper right lat: 52.9375 +Max snow albedo upper right lon: -67.0625 +Max snow albedo resolution (dx): 0.125 +Max snow albedo resolution (dy): 0.125 + +# Greenness fraction maps: +Greenness data source: NCEP_Native #NCEP_LIS +Greenness fraction map: ../input/LS_PARAMETERS/noah_2dparms/gfrac # ./input/LS_PARAMETERS/NLDAS_0.125/gfrac_nldas # Greenness fraction map +Greenness climatology interval: monthly # monthly +Calculate min-max greenness fraction: .false. +Greenness maximum map: ../input/LS_PARAMETERS/noah_2dparms/gfrac_max.asc # ./input/LS_PARAMETERS/NLDAS_0.125/gfrac_nldas.MAX.1gd4r # Maximum greenness fraction map +Greenness minimum map: ../input/LS_PARAMETERS/noah_2dparms/gfrac_min.asc # ./input/LS_PARAMETERS/NLDAS_0.125/gfrac_nldas.MIN.1gd4r # Minimum greenness fraction map +Greenness spatial transform: bilinear #none # average | neighbor | bilinear | budget-bilinear +Greenness fill option: neighbor # none | neighbor | average +Greenness fill radius: 1 # Number of pixels to search for neighbor +Greenness fill value: 0.30 # Static value to fill where missing +Greenness maximum fill value: 0.40 # Static value to fill where missing +Greenness minimum fill value: 0.20 # Static value to fill where missing +Greenness map projection: latlon +Greenness lower left lat: 25.0625 +Greenness lower left lon: -124.9375 +Greenness upper right lat: 52.9375 +Greenness upper right lon: -67.0625 +Greenness resolution (dx): 0.125 +Greenness resolution (dy): 0.125 + +# Slope type map: +Slope type data source: NCEP_Native # NCEP_LIS +Slope type map: ../input/LS_PARAMETERS/noah_2dparms/islope #./input/LS_PARAMETERS/NLDAS_0.125/noah_slope_nldas.1gd4r # Slope type map +Slope type spatial transform: neighbor # none # none | neighbor | mode +Slope type fill option: neighbor # none | neighbor +Slope type fill radius: 1 # Number of pixels to search for neighbor +Slope type fill value: 3. # Static value to fill where missing +Slope type map projection: latlon +Slope type lower left lat: 25.0625 +Slope type lower left lon: -124.9375 +Slope type upper right lat: 52.9375 +Slope type upper right lon: -67.0625 +Slope type resolution (dx): 0.125 +Slope type resolution (dy): 0.125 + +# Bottom temperature map (lapse-rate correction option): +Bottom temperature data source: NCEP_LIS +Bottom temperature map: ../input/LS_PARAMETERS/NLDAS_0.125/tbot_nldas.1gd4r +Bottom temperature spatial transform: none # none | average | neighbor | bilinear | budget-bilinear +Bottom temperature fill option: none # none | average | neighbor +Bottom temperature fill radius: 1 # Number of pixels to search for neighbor +Bottom temperature fill value: 287. # Static value to fill where missing +Bottom temperature topographic downscaling: "none" # none | lapse-rate +Bottom temperature map projection: latlon # Projection type +Bottom temperature lower left lat: 25.0625 +Bottom temperature lower left lon: -124.9375 +Bottom temperature upper right lat: 52.9375 +Bottom temperature upper right lon: -67.0625 +Bottom temperature resolution (dx): 0.125 +Bottom temperature resolution (dy): 0.125 + +# ======================================================= +#Forcing elevation +NLDAS2 forcing directory: /discover/nobackup/projects/lis/MET_FORCING/NLDAS2.FORCING +NLDAS2 data center source: "GES-DISC" +NLDAS2 use model level data: 0 +NLDAS2 use model based swdown: 0 +NLDAS2 use model based precip: 0 +NLDAS2 use model based pressure: 0 +Forcing variables list file: ./ldt_forcing_vars.txt + +NLDAS2 elevation difference map: ../input/LS_PARAMETERS/NLDAS_0.125/NARR_elev-diff.1gd4r +NARR terrain height map: ../input/LS_PARAMETERS/NLDAS_0.125/NARR_elevation.1gd4r + +#Noah-MP LSM inputs +Noah-MP PBL Height Value: 900. + +# ----------------------------------------------------------------------------- +# SMAP SM L3 CDF Generation: +#------------------------------------------------------------------------------ +DA preprocessing method: "CDF generation" +DA observation source: "NASA SMAP soil moisture" +Name of the preprocessed DA file: "cdf_smapobs" + +Apply anomaly correction to obs: 0 +Temporal resolution of CDFs: "monthly" # monthly | yearly +Number of bins to use in the CDF: 100 +Observation count threshold: 30 +Temporal averaging interval: "1da" +Apply external mask: 0 +External mask directory: none + +NASA SMAP soil moisture observation directory: ../input/RS_DATA/SMAP/SPL3SMP.007 +NASA SMAP soil moisture data designation: SPL3SMP +NASA SMAP search radius for openwater proximity detection: 3 +SMAP(NASA) soil moisture Composite Release ID (e.g., R16): "R17" + + +Starting year: 2015 +Starting month: 03 +Starting day: 01 +Starting hour: 0 +Starting minute: 0 +Starting second: 0 +Ending year: 2022 +Ending month: 01 +Ending day: 01 +Ending hour: 0 +Ending minute: 0 +Ending second: 0 +LIS output timestep: "6hr" + + +Maximum number of surface type tiles per grid: 1 +Minimum cutoff percentage (surface type tiles): 0.05 +Maximum number of soil texture tiles per grid: 1 +Minimum cutoff percentage (soil texture tiles): 0.05 +Maximum number of soil fraction tiles per grid: 1 +Minimum cutoff percentage (soil fraction tiles): 0.05 +Maximum number of elevation bands per grid: 1 +Minimum cutoff percentage (elevation bands): 0.05 +Maximum number of slope bands per grid: 1 +Minimum cutoff percentage (slope bands): 0.05 +Maximum number of aspect bands per grid: 1 +Minimum cutoff percentage (aspect bands): 0.05 + +# ----------------------------------------------------------------------------- +# Stratification based on precipitation: +#------------------------------------------------------------------------------ +Stratify CDFs by external data: 1 +Number of bins to use for stratification: 15 +Stratification data source: "LIS LSM total precipitation" +External stratification file: "LVT_MEAN_FINAL.202201010000.d01.nc" # "../Precip_climatology/Precip.climo.eu.merra/LVT_MEAN_FINAL.202201010000.d01.nc" +Write stratified geolocation independent CDFs: 1 + +#LIS precipitation output model name: "Noah.4.0.1" +#LIS precipitation output directory: output.precip +#LIS precipitation output format: "netcdf" +#LIS precipitation output methodology: "2d gridspace" +#LIS precipitation output naming style: "3 level hierarchy" +#LIS precipitation output nest index: 1 +#LIS precipitation output map projection: "lambert" +#LIS precipitation domain lower left lat: 24.01347 #25.20226 +#LIS precipitation domain lower left lon: -13.1299 #-12.5475 +#LIS precipitation domain true lat1: 66.0 +#LIS precipitation domain true lat2: 29.0 +#LIS precipitation domain standard lon: 15.108 +#LIS precipitation domain resolution: 15.0 #5.0 # 1.0 +#LIS precipitation domain x-dimension size: 385 #1114 +#LIS precipitation domain y-dimension size: 289 #826 # 780 diff --git a/lis/configs/lis.config.adoc b/lis/configs/lis.config.adoc index 82382f346..271accd18 100644 --- a/lis/configs/lis.config.adoc +++ b/lis/configs/lis.config.adoc @@ -3062,10 +3062,12 @@ whether the observation error standard deviation is to be scaled using model and observation standard deviation. `SMAP(NASA) model CDF file:` specifies the name of the model CDF file -(observations will be scaled into this climatology). +(observations will be scaled into this climatology). +Note: Soil moisture CDF grouped (stratified) by land cover or precipitation climatology or both simultaneously also can be used here. `SMAP(NASA) observation CDF file:` specifies the name of the observation CDF file. +Note: Soil moisture CDF grouped (stratified) by land cover or precipitation climatology or both simultaneously also can be used here. `SMAP(NASA) soil moisture number of bins in the CDF:` specifies the number of bins in the CDF. @@ -3141,6 +3143,33 @@ SMAP(NRT) CDF read option: .... +[[sssec_cdftransfersmda,Transfering stratified CDFs from one domain to another]] +==== Transfering stratified CDFs from one domain to another + +`Use CDF transfer for soil moisture data assimilation:` specifies +whether to use CDF transfer method. + +`Reference domain model CDF file:` specifies the reference domain model CDF name and data directory. + +`Reference domain obs CDF file:` specifies the reference domain obs CDF name and data directory. + +`Number of bins in the soil moisture CDF:` specifies the number of bins in the CDF. + +`Reference domain precipitation climatology data source:` specifies the reference domain precipitation climatology generated by LVT. + +`Target domain precipitation climatology data source:` specifies the target domain precipitation climatology generated by LVT. + + +.Example _lis.config_ entry +.... +Use CDF transfer for soil moisture data assimilation: +Reference domain model CDF file: ref_model_cdf_forcing_diff/stratified_cdf_noahmp401.nc +Reference domain obs CDF file: ref_obs_cdf_forcing_diff/stratified_cdf_smapobs.nc +Number of bins in the soil moisture CDF: 100 +Reference domain precipitation climatology data source: Precip.climo.us.nldas2/LVT_MEAN_FINAL.202201010000.d01.nc +Target domain precipitation climatology data source: Precip.climo.eu.merra/LVT_MEAN_FINAL.202201010000.d01.nc +.... + [[sssec_smapnasavodda,SMAP (NASA) vegetation optical depth assimilation]] ==== SMAP (NASA) vegetation optical depth assimilation diff --git a/lis/core/LIS_dataAssimMod.F90 b/lis/core/LIS_dataAssimMod.F90 index 4a3a6bbb3..1b70c12c3 100644 --- a/lis/core/LIS_dataAssimMod.F90 +++ b/lis/core/LIS_dataAssimMod.F90 @@ -52,6 +52,11 @@ module LIS_dataAssimMod ! 7 Sep 2017 Mahdi Navari; set a condition to deal the skewness of the both model and obs CDF ! 1 Apr 2019 Yonghwan Kwon; include an option to read soil moisture CDF information for each month separately ! 5 Mar 2021 Eric Kemp; reduced memory usage in read_CDFdata_month +! 25 Feb 2022 Mahdi Navari; added LIS_getCDFtransferattributes +! LIS_rescale_with_stratified_CDF +! read_CDFtransferdata_all +! read_Precip_climo_maxval +! read_Precip_climo ! 1 Oct 2022 Yeosang Yoon; excluded RAPID from 'Routing_DAinst' ! ! !USES: @@ -60,6 +65,7 @@ module LIS_dataAssimMod use LIS_logMod use LIS_historyMod use LIS_DAobservationsMod + use LIS_mpiMod #if(defined USE_NETCDF3 || defined USE_NETCDF4) use netcdf @@ -83,6 +89,11 @@ module LIS_dataAssimMod public :: LIS_readCDFdata !read the CDF files generated by LDT public :: LIS_readMeanSigmaData public :: LIS_forwardEstimate_with_ANN + public :: LIS_getCDFtransferattributes + public :: LIS_rescale_with_stratified_CDF + public :: read_CDFtransferdata_all + public :: read_Precip_climo + public :: read_Precip_climo_maxval !EOP public :: LIS_DA_struc @@ -987,6 +998,174 @@ subroutine LIS_rescale_with_anomaly(& TRACE_EXIT("DA_rescaleAno") end subroutine LIS_rescale_with_anomaly !------------------------------------------------------------------kyh20210422 +! +!MN 2022.02.24 +!BOP +! +! !ROUTINE: LIS_rescale_with_stratified_CDF +! \label{LIS_rescale_with_stratified_CDF} + +! !REVISION HISTORY: +! 3 March 2022: Mahdi Navari ; Initial Specification +! +! !INTERFACE: + subroutine LIS_rescale_with_stratified_CDF(& + n, & + k, & + nbins, & + ntimes, & + max_obs_value, & + min_obs_value, & + model_xrange, & + obs_xrange, & + model_cdf, & + obs_cdf, & + ref_precip_climo_maxval, & + target_precip_climo,& + n_strat_bins,& + obs_value) + + + implicit none +! +! !ARGUMENTS: + integer :: n + integer :: k + integer :: nbins + integer :: n_strat_bins + integer :: ntimes + real :: max_obs_value + real :: min_obs_value + real :: model_xrange(n_strat_bins,ntimes, nbins) + real :: obs_xrange(n_strat_bins,ntimes, nbins) + real :: model_cdf(n_strat_bins,ntimes, nbins) + real :: obs_cdf(n_strat_bins,ntimes, nbins) + real :: obs_value(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + real :: ref_precip_climo_maxval(ntimes) + real :: target_precip_climo(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k),ntimes) +! +! !DESCRIPTION: +! +! This routine rescales the input observation data to the model's +! climatology so that the cumulative distribution functions (CDFs) +! of the observations and the model match (But we use geolocation +! independent CDF, therefore we match the input observation data to +! model's climatology for each stratified bins not for each grid point). +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[nbins] number of bins used to compute the model and obs CDFs +! \item[max\_obs\_value] maximum allowable value of observation +! \item[min\_obs\_value] minimum allowable value of observation +! \item[model\_xrange] x-axis values corresponding to the model CDF +! \item[obs\_xrange] x-axis values corresponding to the obs CDF +! \item[model\_cdf] y-axis (CDF) values corresponding to the model CDF +! \item[obs\_cdf] y-axis (CDF) values corresponding to the obs CDF +! \item[obs\_value] observation value to be rescaled. +! \end{description} +!EOP + + real :: model_delta(n_strat_bins) + real :: obs_delta(n_strat_bins) + real :: ref_p_min, ref_p_max,delta_p_ref + integer :: t,i,kk + integer :: binval, strat_binval + integer :: col,row + real :: cdf_obsval + real :: obs_tmp, obs_in + integer, dimension (1) :: index_25 , index_75 + real :: Lb_xrange, Ub_xrange, iqr_obs, iqr_model + + if(ntimes.gt.1) then + kk = LIS_rc%mo + else + kk = 1 + endif + do t=1,n_strat_bins + model_delta(t) = model_xrange(t,kk,2)-model_xrange(t,kk,1) + obs_delta(t) = obs_xrange(t,kk,2)-obs_xrange(t,kk,1) + enddo + do t=1,LIS_rc%obs_ngrid(k) + + col = LIS_obs_domain(n,k)%col(t) + row = LIS_obs_domain(n,k)%row(t) + + ! find the reference CDF stratification bin value for each grid cell + if(target_precip_climo(col,row,kk).ne.-9999.0) then + ref_p_min = 0. !minval returns -9999. minval(ref_precip_climo(:,:,j)) ! min value over the entire domain + ref_p_max = ref_precip_climo_maxval(kk) ! max value over the entire domain + delta_p_ref = (ref_p_max-ref_p_min)/& + n_strat_bins + strat_binval = nint( (target_precip_climo(col,row,kk)-ref_p_min)/& + delta_p_ref)+1 + if(strat_binval.gt.n_strat_bins) strat_binval = n_strat_bins + if(strat_binval.le.0) strat_binval = 1 + endif + + index_25 = minloc(abs(obs_cdf(strat_binval,kk,:) - 0.25)) + index_75 = minloc(abs(obs_cdf(strat_binval,kk,:) - 0.75)) + Lb_xrange = obs_xrange(strat_binval,kk,index_25(1)) + Ub_xrange = obs_xrange(strat_binval,kk,index_75(1)) + iqr_obs = Ub_xrange - Lb_xrange + + index_25 = minloc(abs(model_cdf(strat_binval,kk,:) - 0.25)) + index_75 = minloc(abs(model_cdf(strat_binval,kk,:) - 0.75)) + Lb_xrange = model_xrange(strat_binval,kk,index_25(1)) + Ub_xrange = model_xrange(strat_binval,kk,index_75(1)) + iqr_model = Ub_xrange - Lb_xrange + + if( iqr_obs .lt. 0.05 * (obs_xrange(strat_binval,kk,nbins)-obs_xrange(strat_binval,kk,1)) .or. & + iqr_model .lt. 0.05 * (model_xrange(strat_binval,kk,nbins)-model_xrange(strat_binval,kk,1)) ) then + obs_delta(strat_binval) = 0 + endif + + if(obs_value(col,row).ne.-9999.0) then + obs_in = obs_value(col,row) + if(obs_delta(strat_binval).gt.0) then + binval = nint((obs_value(col,row)-obs_xrange(strat_binval,kk,1))/& + obs_delta(strat_binval))+1 + if(binval.gt.nbins) binval = nbins + if(binval.le.0) binval = 1 + cdf_obsval = obs_cdf(strat_binval,kk,binval) + if(cdf_obsval.gt.1.0) cdf_obsval = 1.0 + if(cdf_obsval .gt.0)then + i=1 + do while((model_cdf(strat_binval,kk,i).lt.cdf_obsval).and.& + (i.le.nbins)) + i = i+1 + if(i.gt.nbins) exit + enddo + if(i.gt.nbins) i = i-1 + obs_tmp = model_xrange(strat_binval,kk,i) + + if(obs_tmp.gt.max_obs_value) then + obs_tmp = LIS_rc%udef + endif + + if(obs_tmp.le.min_obs_value) then + obs_tmp = LIS_rc%udef + endif + obs_value(col,row) = obs_tmp + else + obs_value(col,row) = LIS_rc%udef + endif + else + obs_value(col,row) = LIS_rc%udef + endif + if(obs_value(col,row).le.min_obs_value.and.& + obs_value(col,row).ne.-9999.0) then + write(LIS_logunit,*) '[ERR] Problem in CDF scaling of observations in the DA instance ',k + write(LIS_logunit,*) '[ERR] ',col,row,obs_value(col,row), obs_in + call LIS_endrun() + endif + else + obs_value(col,row) = LIS_rc%udef + endif + enddo + !TRACE_EXIT("DA_rescaleCDF") + end subroutine LIS_rescale_with_stratified_CDF +!MN 2022.02.24 !BOP @@ -1046,6 +1225,58 @@ subroutine LIS_getCDFattributes(k,filename, ntimes, ngrid) end subroutine LIS_getCDFattributes +!BOP +! !ROUTINE: LIS_getCDFtransferattributes +! \label{LIS_getCDFtransferattributes} + +! !REVISION HISTORY: +! 3 March 2022: Mahdi Navari ; Initial Specification +! +! !INTERFACE: + subroutine LIS_getCDFtransferattributes(k,filename, n_strat_bins, ntimes) + + implicit none +! !ARGUMENTS: + integer :: k + character(len=*), intent(in) :: filename + integer :: ntimes + integer :: n_strat_bins +! +! !DESCRIPTION: +! This routine reads the number of stratification bins abd other CDF attribute from the +! CDF file +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[ntimes] temporal resolution of the CDFs +! \end{description} +!EOP + + integer :: nid, binid + + !TRACE_ENTER("DA_getCDFatt") +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + call LIS_verify(nf90_open(path=trim(filename),mode=NF90_NOWRITE,& + ncid=nid),'failed to open file '//trim(filename)) + + call LIS_verify(nf90_inq_dimid(nid, 'n_strat_bins',binId), & + 'Error nf90_inq_dimid: n_strat_bins') + + call LIS_verify(nf90_inquire_dimension(nid, binId, len=n_strat_bins), & + 'Error nf90_inquire_dimension:n_strat_bins') + + call LIS_verify(nf90_get_att(nid, NF90_GLOBAL, & + 'temporal_resolution_CDF', & + ntimes), & + 'Error in nf90_get_att: temporal_resolution_CDF') + + call LIS_verify(nf90_close(nid)) +#endif + !TRACE_EXIT("DA_getCDFatt") + + end subroutine LIS_getCDFtransferattributes + !BOP ! !ROUTINE: read_CDFdata_all @@ -1472,6 +1703,286 @@ subroutine read_MeanSigmaData_month(n, k, ntimes, ngrid, filename, varname, mu, TRACE_EXIT("DA_readSigma") end subroutine read_MeanSigmaData_month + +!BOP +! !ROUTINE: read_CDFtransferdata_all +! \label{read_CDFtransferdata_all} +! !REVISION HISTORY: +! 1 March 2022: Mahdi Navari ; Initial Specification +! +! !INTERFACE: + subroutine read_CDFtransferdata_all(n, k, nbins, ntimes, n_strat_bins, & + filename, varname, xrange, cdf) + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + integer, intent(in) :: nbins + integer, intent(in) :: ntimes + integer, intent(in) :: n_strat_bins + character(len=*) :: filename + character(len=*) :: varname + real :: xrange(n_strat_bins,ntimes, nbins) + real :: cdf(n_strat_bins,ntimes, nbins) + +! +! !DESCRIPTION: +! This routine reads the input CDF file (generated by LDT in NETCDF format) +! The xrange values and the corresponding CDFs are read for each stratification bin. +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[nbins] number of bins used to compute the model and obs CDFs +! \item[filename] name of the CDF file +! \item[varname] name of the variable being extracted. +! \item[xrange] x-axis values corresponding to the CDF +! \item[cdf] y-axis (CDF) values corresponding to the CDF +! \end{description} +!EOP + integer :: j,kk + integer :: nstratbinId, nbinsId, nlevsId + integer :: nstratbin_file, nbins_file, nlevs_file + integer :: xid, cdfid + real, allocatable :: xrange_file(:,:,:,:) + real, allocatable :: cdf_file(:,:,:,:) + integer :: nid + + !TRACE_ENTER("DA_readCDF") +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + write(LIS_logunit,*) '[INFO] Reading stratified geolocation independent reference CDF file ',trim(filename) + call LIS_verify(nf90_open(path=trim(filename),mode=NF90_NOWRITE,& + ncid=nid),'failed to open file '//trim(filename)) + + call LIS_verify(nf90_inq_dimid(nid,"n_strat_bins",nstratbinId), & + 'nf90_inq_dimid failed for n_strat_bins') + call LIS_verify(nf90_inq_dimid(nid,"nbins",nbinsId), & + 'nf90_inq_dimid failed for nbins') + call LIS_verify(nf90_inq_dimid(nid,trim(varname)//"_levels",nlevsId), & + 'nf90_inq_dimid failed for '//trim(varname)//"_levels") + + call LIS_verify(nf90_inquire_dimension(nid,nstratbinId, len=nstratbin_file),& + 'nf90_inquire_dimension failed for n_strat_bins') + call LIS_verify(nf90_inquire_dimension(nid,nbinsId, len=nbins_file),& + 'nf90_inquire_dimension failed for nbins') + call LIS_verify(nf90_inquire_dimension(nid,nlevsId, len=nlevs_file),& + 'nf90_inquire_dimension failed for nlevels') + + if(nbins.ne.nbins_file) then + write(LIS_logunit,*) '[ERR] The number of bins specified in the file '//& + trim(filename) + write(LIS_logunit,*) '[ERR] (',nbins_file, & + ') is different from the number of bins specified' + write(LIS_logunit,*) '[ERR] in the lis.config file (',nbins,')' + call LIS_endrun() + endif + + allocate(xrange_file(nstratbin_file,ntimes,nlevs_file,nbins)) + allocate(cdf_file(nstratbin_file,ntimes,nlevs_file,nbins)) + + call LIS_verify(nf90_inq_varid(nid,trim(varname)//'_xrange_stratified',xid),& + 'nf90_inq_varid failed for for '//trim(varname)//'_xrange_stratified') + call LIS_verify(nf90_inq_varid(nid,trim(varname)//'_CDF_stratified',cdfid),& + 'nf90_inq_varid failed for '//trim(varname)//'_CDF_stratified') + + call LIS_verify(nf90_get_var(nid,xid,xrange_file, & + start=(/1,1,1,1/), count=(/nstratbin_file,ntimes,nlevs_file,nbins/)),& + 'nf90_get_var failed for '//trim(varname)//'_xrange_stratified') + call LIS_verify(nf90_get_var(nid,cdfid,cdf_file,& + start=(/1,1,1,1/), count=(/nstratbin_file,ntimes,nlevs_file,nbins/)),& + 'nf90_get_var failed for '//trim(varname)//'_CDF_stratified') + + xrange = xrange_file(:,:,1,:) + cdf = cdf_file(:,:,1,:) + + deallocate(xrange_file) + deallocate(cdf_file) + + call LIS_verify(nf90_close(nid),& + 'failed to close file '//trim(filename)) + write(LIS_logunit,*) '[INFO] Successfully read CDF file ',trim(filename) +#endif + ! TRACE_EXIT("DA_readCDF") + end subroutine read_CDFtransferdata_all + +!BOP +! !ROUTINE: read_Precip_climo +! \label{read_Precip_climo} + +! !REVISION HISTORY: +! 28 Feb 2022: Mahdi Navari ; Initial Specification +! +! !INTERFACE: +subroutine read_Precip_climo(n,k, filename, precip) + use LIS_coreMod + use LIS_logMod + use LIS_historyMod +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + character(len=*), intent(in) :: filename + real, allocatable, intent(inout) :: precip(:,:,:) + logical :: file_exists +! +! !DESCRIPTION: +! This routine reads the input CDF file (generated by LVT in NETCDF format) +! +! The arguments are: +! \begin{description} +! \item[filename] name of the CDF file +! \item[varname] name of the variable being extracted. +! \end{description} +!EOP + integer :: ios,nid,ncId,nrId,varId + integer :: gnc,gnr,gr,gc,nc,nr,i,c,r + character*3 :: month_name(12) + integer :: gindex,gid + integer :: ierr + + month_name = (/"JAN","FEB","MAR","APR","MAY","JUN",& + "JUL","AUG","SEP","OCT","NOV","DEC"/) + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + inquire(file=trim(filename), exist=file_exists) + if(file_exists) then + write(LIS_logunit,*) '[INFO] Reading Precipitation climatology form CDF file ',trim(filename) + + call LIS_verify(nf90_open(path=trim(filename),mode=NF90_NOWRITE,& + ncid=nid),'failed to open file '//trim(filename)) + + ios = nf90_inq_dimid(nid,"east_west",ncId) + call LIS_verify(ios,'Error in nf90_inq_dimid in read_precip_climo') + + ios = nf90_inq_dimid(nid,"north_south",nrId) + call LIS_verify(ios,'Error in nf90_inq_dimid in read_precip_climo') + + ios = nf90_inquire_dimension(nid,ncId, len=gnc) + call LIS_verify(ios,'Error in nf90_inquire_dimension in read_precip_climo') + + ios = nf90_inquire_dimension(nid,nrId, len=gnr) + call LIS_verify(ios,'Error in nf90_inquire_dimension in read_precip_climo') + + if (gnc .ne. LIS_rc%gnc(n) .or. gnr .ne. LIS_rc%gnr(n)) then + write(LIS_logunit,*) '[ERR] domain dimension mismatch. Number of rows and columns in the ' + write(LIS_logunit,*) 'precipitation climatology file differs from that in the LIS domain ' + call LIS_endrun() + endif + + allocate(precip(gnc,gnr,12)) + do i=1,12 + ios = nf90_inq_varid(nid,'TotalPrecip_'//trim(month_name(i)),varId) + call LIS_verify(ios,'Precipitation climo field not found in the file') + + ios = nf90_get_var(nid,varId,precip(:,:,i)) ! precip_climo(:,:,i)) + call LIS_verify(ios,'Error in nf90_get_var in read_Precip_climo') + enddo + + ios = nf90_close(nid) + call LIS_verify(ios,'Error in nf90_close in read_Precip_climo') + write(LIS_logunit,*) '[INFO] Successfully read Precipitation climo file ',trim(filename) + endif + +! end USE_NETCDF4 +#endif + +end subroutine read_Precip_climo + +!BOP +! !ROUTINE: read_Precip_climo +! \label{read_Precip_climo} + +! !REVISION HISTORY: +! 28 Feb 2022: Mahdi Navari ; Initial Specification +! +! !INTERFACE: +subroutine read_Precip_climo_maxval(n,filename, strat_cdfs_max) + + use LIS_coreMod + use LIS_logMod + use LIS_historyMod +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + character(len=*), intent(in) :: filename + real :: strat_cdfs_max(12) + logical :: file_exists +! +! !DESCRIPTION: +! This routine reads the input CDF file (generated by LVT in NETCDF format) +! +! The arguments are: +! \begin{description} +! \item[filename] name of the CDF file +! \item[varname] name of the variable being extracted. +! \end{description} +!EOP + real, allocatable :: precip_climo(:,:,:) + integer :: ios,nid,ncId,nrId,varId + integer :: gnc,gnr,gr,gc,nc,nr,i,c,r + character*3 :: month_name(12) + integer :: gindex + + month_name = (/"JAN","FEB","MAR","APR","MAY","JUN",& + "JUL","AUG","SEP","OCT","NOV","DEC"/) + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + inquire(file=trim(filename), exist=file_exists) + if(file_exists) then + write(LIS_logunit,*) '[INFO] Reading Precipitation climatology form CDF file ',trim(filename) + + call LIS_verify(nf90_open(path=trim(filename),mode=NF90_NOWRITE,& + ncid=nid),'failed to open file '//trim(filename)) + + ios = nf90_inq_dimid(nid,"east_west",ncId) + call LIS_verify(ios,'Error in nf90_inq_dimid in read_precip_climo_maxval') + + ios = nf90_inq_dimid(nid,"north_south",nrId) + call LIS_verify(ios,'Error in nf90_inq_dimid in read_precip_climo_maxval') + + ios = nf90_inquire_dimension(nid,ncId, len=gnc) + call LIS_verify(ios,'Error in nf90_inquire_dimension in read_precip_climo_maxval') + + ios = nf90_inquire_dimension(nid,nrId, len=gnr) + call LIS_verify(ios,'Error in nf90_inquire_dimension in read_precip_climo_maxval') + + allocate(precip_climo(gnc,gnr,12)) + + do i=1,12 + ios = nf90_inq_varid(nid,'TotalPrecip_'//trim(month_name(i)),varId) + call LIS_verify(ios,'Precipitation climo field not found in the file') + + ios = nf90_get_var(nid,varId,precip_climo(:,:,i)) + call LIS_verify(ios,'Error in nf90_get_var in read_Precip_climo_maxval') + enddo + + ios = nf90_close(nid) + call LIS_verify(ios,'Error in nf90_close in read_Precip_climo_maxval') + write(LIS_logunit,*) '[INFO] Successfully read Precipitation climo file ',trim(filename) + endif + + ! Note: we just need Max and min values for each month (Min is zero) + do i=1,12 ! number of month + strat_cdfs_max(i) = maxval(precip_climo(:,:,i)) ! max value over the entire domain + enddo + + deallocate(precip_climo) + +! end USE_NETCDF4 +#endif + +end subroutine read_Precip_climo_maxval + + !BOP ! ! !ROUTINE: LIS_forwardEstimate_with_ANN diff --git a/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/cdfTransfer_NASASMAPsm_Mod.F90 b/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/cdfTransfer_NASASMAPsm_Mod.F90 new file mode 100644 index 000000000..875dd2207 --- /dev/null +++ b/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/cdfTransfer_NASASMAPsm_Mod.F90 @@ -0,0 +1,537 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.3 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! +! !MODULE: cdfTransfer_NASASMAPsm_Mod +! +! !DESCRIPTION: +! This module contains interfaces and subroutines to +! handle +! +! !REVISION HISTORY: +! 2 Mar 2022 Mahdi Navari; initial specification +! +module cdfTransfer_NASASMAPsm_Mod +! !USES: + use ESMF + use map_utils + + implicit none + + PRIVATE + +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: cdfTransfer_NASASMAPsm_setup +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + public :: cdfT_SMAPsm_struc +!EOP + type, public:: cdfT_SMAPsm_dec + + integer :: useSsdevScal + logical :: startMode + integer :: nc + integer :: nr + real, allocatable :: smobs(:,:) + real*8, allocatable :: smtime(:,:) + + real :: ssdev_inp + integer, allocatable :: n11(:) + real, allocatable :: rlat(:) + real, allocatable :: rlon(:) + + real, allocatable :: model_xrange(:,:,:) + real, allocatable :: obs_xrange(:,:,:) + real, allocatable :: model_cdf(:,:,:) + real, allocatable :: obs_cdf(:,:,:) + real, allocatable :: model_mu(:,:) + real, allocatable :: obs_mu(:,:) + real, allocatable :: model_sigma(:,:) + real, allocatable :: obs_sigma(:,:) + character*20 :: data_designation + character*3 :: release_number + integer :: nbins + integer :: ntimes + + logical :: cdf_read_mon !(for reading monthly CDF when + !LIS_rc%da > 1 but the first model time step, + !e.g., 4/29 13:00:00) + integer :: cdf_read_opt ! 0: read all months at one time + ! 1: read only the current month + character*100 :: modelcdffile + character*100 :: obscdffile + integer :: n_strat_bins + integer :: useCDFtransfer + character*100 :: ref_p_climo_file + character*100 :: target_p_climo_file + real, allocatable :: ref_p_climo_maxval(:) + real, allocatable :: target_p_climo(:,:,:) + + end type cdfT_SMAPsm_dec + + type(cdfT_SMAPsm_dec),allocatable :: cdfT_SMAPsm_struc(:) + +contains + +!BOP +! +! !ROUTINE: cdfTransfer_NASASMAPsm_setup +! \label{cdfTransfer_NASASMAPsm_setup} +! +! !INTERFACE: + subroutine cdfTransfer_NASASMAPsm_setup(k, OBS_State, OBS_Pert_State) +! !USES: + use LIS_coreMod + use LIS_timeMgrMod + use LIS_historyMod + use LIS_dataAssimMod + use LIS_perturbMod + use LIS_DAobservationsMod + use LIS_logmod + + implicit none + +! !ARGUMENTS: + integer :: k + type(ESMF_State) :: OBS_State(LIS_rc%nnest) + type(ESMF_State) :: OBS_Pert_State(LIS_rc%nnest) +! +! !DESCRIPTION: +! +! This routine completes the runtime initializations and +! creation of data strctures required for handling NASASMAP +! soil moisture data. +! +! The arguments are: +! \begin{description} +! \item[OBS\_State] observation state +! \item[OBS\_Pert\_State] observation perturbations state +! \end{description} +!EOP + real, parameter :: minssdev =0.001 + integer :: n,i,t,kk,jj + integer :: ftn + integer :: status + type(ESMF_Field) :: obsField(LIS_rc%nnest) + type(ESMF_ArraySpec) :: intarrspec, realarrspec + type(ESMF_Field) :: pertField(LIS_rc%nnest) + type(ESMF_ArraySpec) :: pertArrSpec + character*100 :: rtsmopssmobsdir + character*100 :: temp + real, allocatable :: obsstd(:) + character*1 :: vid(2) + character*40, allocatable :: vname(:) + real , allocatable :: varmin(:) + real , allocatable :: varmax(:) + type(pert_dec_type) :: obs_pert + real, pointer :: obs_temp(:,:) + real :: gridDesci(50) + real, allocatable :: ssdev(:) + + real, allocatable :: obserr(:,:) + real, allocatable :: lobserr(:,:) + integer :: c,r + real, allocatable :: ssdev_grid(:,:) + integer :: ngrid + real, allocatable :: target_precip_climo(:,:,:) + + allocate(cdfT_SMAPsm_struc(LIS_rc%nnest)) + + call ESMF_ArraySpecSet(intarrspec,rank=1,typekind=ESMF_TYPEKIND_I4,& + rc=status) + call LIS_verify(status) + + call ESMF_ArraySpecSet(realarrspec,rank=1,typekind=ESMF_TYPEKIND_R4,& + rc=status) + call LIS_verify(status) + + call ESMF_ArraySpecSet(pertArrSpec,rank=2,typekind=ESMF_TYPEKIND_R4,& + rc=status) + call LIS_verify(status) + + call ESMF_ConfigFindLabel(LIS_config,"SMAP(NASA) soil moisture data directory:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,rtsmopssmobsdir,& + rc=status) + call LIS_verify(status, 'SMAP(NASA) soil moisture data directory: is missing') + + call ESMF_AttributeSet(OBS_State(n),"Data Directory",& + rtsmopssmobsdir, rc=status) + call LIS_verify(status) + enddo + + call ESMF_ConfigFindLabel(LIS_config,"SMAP(NASA) soil moisture data designation:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,& + cdfT_SMAPsm_struc(n)%data_designation,& + rc=status) + call LIS_verify(status, 'SMAP(NASA) soil moisture data designation: is missing') + enddo + + call ESMF_ConfigFindLabel(LIS_config,"SMAP(NASA) soil moisture Composite Release ID:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,& + cdfT_SMAPsm_struc(n)%release_number,& + rc=status) + call LIS_verify(status, 'SMAP(NASA) soil moisture Composite Release ID: is missing') + enddo + + call ESMF_ConfigFindLabel(LIS_config,"SMAP(NASA) soil moisture use scaled standard deviation model:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,cdfT_SMAPsm_struc(n)%useSsdevScal,& + rc=status) + call LIS_verify(status, 'SMAP(NASA) soil moisture use scaled standard deviation model: is missing') + + enddo + + call ESMF_ConfigFindLabel(LIS_config, "SMAP(NASA) soil moisture number of bins in the CDF:", rc=status) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,cdfT_SMAPsm_struc(n)%nbins, rc=status) + call LIS_verify(status, "SMAP(NASA) soil moisture number of bins in the CDF: not defined") + enddo + + call ESMF_ConfigFindLabel(LIS_config,"Use CDF transfer for soil moisture data assimilation:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,& + cdfT_SMAPsm_struc(n)%useCDFtransfer,& + default=0,rc=status) + call LIS_verify(status, 'Use CDF transfer for soil moisture data assimilation: is missing') + enddo + + call ESMF_ConfigFindLabel(LIS_config,"Reference domain model CDF file:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,& + cdfT_SMAPsm_struc(n)%modelcdffile,& + rc=status) + call LIS_verify(status, 'Reference domain model CDF file: is missing') + enddo + + call ESMF_ConfigFindLabel(LIS_config,"Reference domain obs CDF file:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,& + cdfT_SMAPsm_struc(n)%obscdffile,& + rc=status) + call LIS_verify(status, 'Reference domain obs CDF file: is missing') + enddo + + call ESMF_ConfigFindLabel(LIS_config,"Reference domain precipitation climatology data source:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,& + cdfT_SMAPsm_struc(n)%ref_p_climo_file,& + rc=status) + call LIS_verify(status, 'Reference domain precipitation climatology data source: is missing') + enddo + + call ESMF_ConfigFindLabel(LIS_config,"Target domain precipitation climatology data source:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,& + cdfT_SMAPsm_struc(n)%target_p_climo_file,& + rc=status) + call LIS_verify(status, 'Target domain precipitation climatology data source: is missing') + enddo + + do n=1,LIS_rc%nnest + call ESMF_AttributeSet(OBS_State(n),"Data Update Status",& + .false., rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(OBS_State(n),"Data Update Time",& + -99.0, rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(OBS_State(n),"Data Assimilate Status",& + .false., rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(OBS_State(n),"Number Of Observations",& + LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status) + + enddo + + write(LIS_logunit,*)& + '[INFO] read reference and target domains precipitation climatology data' +!-------------------------------------------------------------------------------- +! Data must be in the local domain +! The stratified CDF is geolocation independent +! Subset the precipition climotology to local space +!------------------------------------------------------------------------------- + do n=1,LIS_rc%nnest + allocate(cdfT_SMAPsm_struc(n)%ref_p_climo_maxval(12)) !cdfT_SMAPsm_struc(n)%ntimes)) + allocate(cdfT_SMAPsm_struc(n)%target_p_climo(LIS_rc%lnc(n), LIS_rc%lnr(n),12)) !cdfT_SMAPsm_struc(n)%ntimes)) + ! Note: we do not need to know the dimenesions of the reference domain precip climatology + ! for allocation so we just extract max value for each month + + call read_Precip_climo_maxval (n,cdfT_SMAPsm_struc(n)%ref_p_climo_file,& + cdfT_SMAPsm_struc(n)%ref_p_climo_maxval) + + call read_Precip_climo (n,k,& + cdfT_SMAPsm_struc(n)%target_p_climo_file,& + target_precip_climo) + + ! subset from the global 2-d space to the local 2-d space + do i=1,12 + cdfT_SMAPsm_struc(n)%target_p_climo(:,:,i) = target_precip_climo(& + LIS_ews_halo_ind(n,LIS_localPet+1):& + LIS_ewe_halo_ind(n,LIS_localPet+1), & + LIS_nss_halo_ind(n,LIS_localPet+1): & + LIS_nse_halo_ind(n,LIS_localPet+1),i) + enddo + enddo + + write(LIS_logunit,*)& + '[INFO] read SMAP(NASA) soil moisture data specifications' + +!---------------------------------------------------------------------------- +! Create the array containers that will contain the observations and +! the perturbations. NASASMAP +! observations are in the grid space. Since there is only one layer +! being assimilated, the array size is LIS_rc%obs_ngrid(k). +! +!---------------------------------------------------------------------------- + + do n=1,LIS_rc%nnest + + write(unit=temp,fmt='(i2.2)') 1 + read(unit=temp,fmt='(2a1)') vid + + obsField(n) = ESMF_FieldCreate(arrayspec=realarrspec,& + grid=LIS_obsvecGrid(n,k),& + name="Observation"//vid(1)//vid(2),rc=status) + call LIS_verify(status) + +!Perturbations State + write(LIS_logunit,*) '[INFO] Opening attributes for observations ',& + trim(LIS_rc%obsattribfile(k)) + ftn = LIS_getNextUnitNumber() + open(ftn,file=trim(LIS_rc%obsattribfile(k)),status='old') + read(ftn,*) + read(ftn,*) LIS_rc%nobtypes(k) + read(ftn,*) + + allocate(vname(LIS_rc%nobtypes(k))) + allocate(varmax(LIS_rc%nobtypes(k))) + allocate(varmin(LIS_rc%nobtypes(k))) + + do i=1,LIS_rc%nobtypes(k) + read(ftn,fmt='(a40)') vname(i) + read(ftn,*) varmin(i),varmax(i) + write(LIS_logunit,*) '[INFO] ',vname(i),varmin(i),varmax(i) + enddo + call LIS_releaseUnitNumber(ftn) + + allocate(ssdev(LIS_rc%obs_ngrid(k))) + + if(trim(LIS_rc%perturb_obs(k)).ne."none") then + allocate(obs_pert%vname(1)) + allocate(obs_pert%perttype(1)) + allocate(obs_pert%ssdev(1)) + allocate(obs_pert%stdmax(1)) + allocate(obs_pert%zeromean(1)) + allocate(obs_pert%tcorr(1)) + allocate(obs_pert%xcorr(1)) + allocate(obs_pert%ycorr(1)) + allocate(obs_pert%ccorr(1,1)) + + call LIS_readPertAttributes(1,LIS_rc%obspertAttribfile(k),& + obs_pert) + +! Set obs err to be uniform (will be rescaled later for each grid point). + ssdev = obs_pert%ssdev(1) + cdfT_SMAPsm_struc(n)%ssdev_inp = obs_pert%ssdev(1) + + pertField(n) = ESMF_FieldCreate(arrayspec=pertArrSpec,& + grid=LIS_obsensOnGrid(n,k),name="Observation"//vid(1)//vid(2),& + rc=status) + call LIS_verify(status) + +! initializing the perturbations to be zero + call ESMF_FieldGet(pertField(n),localDE=0,farrayPtr=obs_temp,rc=status) + call LIS_verify(status) + obs_temp(:,:) = 0 + + call ESMF_AttributeSet(pertField(n),"Perturbation Type",& + obs_pert%perttype(1), rc=status) + call LIS_verify(status) + + if(LIS_rc%obs_ngrid(k).gt.0) then + call ESMF_AttributeSet(pertField(n),"Standard Deviation",& + ssdev,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status) + endif + + call ESMF_AttributeSet(pertField(n),"Std Normal Max",& + obs_pert%stdmax(1), rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(pertField(n),"Ensure Zero Mean",& + obs_pert%zeromean(1),rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(pertField(n),"Temporal Correlation Scale",& + obs_pert%tcorr(1),rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(pertField(n),"X Correlation Scale",& + obs_pert%xcorr(1),rc=status) + + call ESMF_AttributeSet(pertField(n),"Y Correlation Scale",& + obs_pert%ycorr(1),rc=status) + + call ESMF_AttributeSet(pertField(n),"Cross Correlation Strength",& + obs_pert%ccorr(1,:),itemCount=1,rc=status) + + endif + + deallocate(vname) + deallocate(varmax) + deallocate(varmin) + deallocate(ssdev) + + enddo + + write(LIS_logunit,*) & + '[INFO] Created the States to hold the SMAP NASA observations data' + do n=1,LIS_rc%nnest + if(cdfT_SMAPsm_struc(n)%data_designation.eq."SPL3SMP") then + cdfT_SMAPsm_struc(n)%nc = 964 + cdfT_SMAPsm_struc(n)%nr = 406 + elseif(cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP") then + cdfT_SMAPsm_struc(n)%nc = 964 + cdfT_SMAPsm_struc(n)%nr = 406 + elseif(cdfT_SMAPsm_struc(n)%data_designation.eq."SPL3SMP_E") then + cdfT_SMAPsm_struc(n)%nc = 3856 + cdfT_SMAPsm_struc(n)%nr = 1624 + elseif(cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP_E") then + cdfT_SMAPsm_struc(n)%nc = 3856 + cdfT_SMAPsm_struc(n)%nr = 1624 + endif + + allocate(cdfT_SMAPsm_struc(n)%smobs(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k))) + allocate(cdfT_SMAPsm_struc(n)%smtime(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k))) + + cdfT_SMAPsm_struc(n)%smtime = -1 + + enddo + + write(LIS_logunit,*)& + '[INFO] read cdf transfer data specifications' + + do n=1,LIS_rc%nnest + allocate(ssdev(LIS_rc%obs_ngrid(k))) + ssdev = obs_pert%ssdev(1) + + if(cdfT_SMAPsm_struc(n)%useCDFtransfer.gt.0) then ! LIS_rc%dascaloption(k).eq."CDF matching" .and. & + + call LIS_getCDFtransferattributes(k,cdfT_SMAPsm_struc(n)%modelcdffile,& + cdfT_SMAPsm_struc(n)%n_strat_bins , cdfT_SMAPsm_struc(n)%ntimes) + + allocate(cdfT_SMAPsm_struc(n)%model_xrange(& + cdfT_SMAPsm_struc(n)%n_strat_bins, cdfT_SMAPsm_struc(n)%ntimes, & + cdfT_SMAPsm_struc(n)%nbins)) + allocate(cdfT_SMAPsm_struc(n)%obs_xrange(& + cdfT_SMAPsm_struc(n)%n_strat_bins, cdfT_SMAPsm_struc(n)%ntimes, & + cdfT_SMAPsm_struc(n)%nbins)) + allocate(cdfT_SMAPsm_struc(n)%model_cdf(& + cdfT_SMAPsm_struc(n)%n_strat_bins, cdfT_SMAPsm_struc(n)%ntimes, & + cdfT_SMAPsm_struc(n)%nbins)) + allocate(cdfT_SMAPsm_struc(n)%obs_cdf(& + cdfT_SMAPsm_struc(n)%n_strat_bins, cdfT_SMAPsm_struc(n)%ntimes, & + cdfT_SMAPsm_struc(n)%nbins)) + write(LIS_logunit,*)& + '[INFO] Successfully read cdf transfer data specifications' + endif + + if(cdfT_SMAPsm_struc(n)%useSsdevScal.eq.1) then + write(LIS_logunit,*) '[ERR] "use scaled standard deviation model" does not work ' + write(LIS_logunit,*) '[ERR] for CDF transfer method' + call LIS_endrun() + endif + + if(LIS_rc%obs_ngrid(k).gt.0) then + call ESMF_AttributeSet(pertField(n),"Standard Deviation",& + ssdev,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status) + endif + + deallocate(ssdev) + enddo + + + do n=1,LIS_rc%nnest + if(cdfT_SMAPsm_struc(n)%data_designation.eq."SPL3SMP".or.& + cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP") then + gridDesci = 0 + gridDesci(1) = 9 + gridDesci(2) = 964 + gridDesci(3) = 406 + gridDesci(9) = 4 !M36 grid + gridDesci(20) = 64 + gridDesci(10) = 0.36 + gridDesci(11) = 1 !for the global switch + elseif(cdfT_SMAPsm_struc(n)%data_designation.eq."SPL3SMP_E".or.& + cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP_E") then + gridDesci = 0 + gridDesci(1) = 9 + gridDesci(2) = 3856 + gridDesci(3) = 1624 + gridDesci(9) = 5 !M09 grid + gridDesci(20) = 64 + gridDesci(10) = 0.09 + gridDesci(11) = 1 !for the global switch + endif + + allocate(cdfT_SMAPsm_struc(n)%n11(LIS_rc%obs_lnc(k)*& + LIS_rc%obs_lnr(k))) + allocate(cdfT_SMAPsm_struc(n)%rlat(& + LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k))) + allocate(cdfT_SMAPsm_struc(n)%rlon(& + LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k))) + + call neighbor_interp_input_withgrid(gridDesci,& + LIS_rc%obs_gridDesc(k,:),& + LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k),& + cdfT_SMAPsm_struc(n)%rlat, & + cdfT_SMAPsm_struc(n)%rlon, & + cdfT_SMAPsm_struc(n)%n11) + + if(cdfT_SMAPsm_struc(n)%data_designation.eq."SPL3SMP".or.& + cdfT_SMAPsm_struc(n)%data_designation.eq."SPL3SMP_E") then + call LIS_registerAlarm("NASASMAP read alarm",& + 86400.0, 86400.0) + else + call LIS_registerAlarm("NASASMAP read alarm",& + 3600.0, 3600.0) + endif + + cdfT_SMAPsm_struc(n)%startMode = .true. + + call ESMF_StateAdd(OBS_State(n),(/obsField(n)/),rc=status) + call LIS_verify(status) + + call ESMF_StateAdd(OBS_Pert_State(n),(/pertField(n)/),rc=status) + call LIS_verify(status) + + enddo + end subroutine cdfTransfer_NASASMAPsm_setup + +end module cdfTransfer_NASASMAPsm_Mod diff --git a/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/read_cdfTransfer_NASASMAPsm.F90 b/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/read_cdfTransfer_NASASMAPsm.F90 new file mode 100644 index 000000000..d07c4df8c --- /dev/null +++ b/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/read_cdfTransfer_NASASMAPsm.F90 @@ -0,0 +1,1325 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.3 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LIS_misc.h" +!BOP +! !ROUTINE: read_cdfTransfer_NASASMAPsm +! \label{read_cdfTransfer_NASASMAPsm} +! +! !REVISION HISTORY: +! 2 MAR 2022: Mahdi Navari: initial specification based on read_NASASMAPsm.F90 +! +! !INTERFACE: +subroutine read_cdfTransfer_NASASMAPsm(n, k, OBS_State, OBS_Pert_State) +! !USES: + use ESMF + use LIS_mpiMod + use LIS_coreMod + use LIS_logMod + use LIS_timeMgrMod + use LIS_dataAssimMod + use LIS_DAobservationsMod + use map_utils + use LIS_pluginIndices + use cdfTransfer_NASASMAPsm_Mod, only: cdfT_SMAPsm_struc + !use cdfTransfer_NASASMAPsm_Mod, only: cdfT_SMAPsm_struc + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + type(ESMF_State) :: OBS_State + type(ESMF_State) :: OBS_Pert_State +! +! !DESCRIPTION: +! +! reads the AMSRE soil moisture observations +! from NETCDF files and applies the spatial masking for dense +! vegetation, rain and RFI. The data is then rescaled +! to the land surface model's climatology using rescaling +! algorithms. +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[OBS\_State] observations state +! \end{description} +! +!EOP + real, parameter :: minssdev = 0.05 + real, parameter :: maxssdev = 0.11 + real, parameter :: MAX_SM_VALUE = 0.45, MIN_SM_VALUE = 0.0001 + integer :: status + integer :: grid_index + character*100 :: smobsdir + character*100 :: fname + logical :: alarmCheck, file_exists + integer :: t, c, r, jj + real, pointer :: obsl(:) + type(ESMF_Field) :: smfield, pertField + integer :: gid(LIS_rc%obs_ngrid(k)) + integer :: assimflag(LIS_rc%obs_ngrid(k)) + real :: obs_unsc(LIS_rc%obs_ngrid(k)) + logical :: data_update + logical :: data_upd_flag(LIS_npes) + logical :: data_upd_flag_local + logical :: data_upd + real :: smobs(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + real :: smobs_D(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + real :: smobs_A(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + real :: sm_current(LIS_rc%obs_lnc(k), LIS_rc%obs_lnr(k)) + real :: dt + real :: lon + real :: lhour + real :: gmt + integer :: zone + integer :: fnd + real, allocatable :: ssdev(:) + character*4 :: yyyy + character*8 :: yyyymmdd + character*2 :: mm, dd, hh + integer :: yr, mo, da, hr, mn, ss + integer :: cyr, cmo, cda, chr, cmn, css + integer :: nyr, nmo, nda, nhr, nmn, nss + real*8 :: timenow, time1,time2,time3 + integer :: doy + character*200 :: list_files + integer :: mn_ind + integer :: ftn, ierr + integer :: rc + character(len=3) :: CRID + integer, external :: create_filelist ! C function + + + call ESMF_AttributeGet(OBS_State, "Data Directory", & + smobsdir, rc=status) + call LIS_verify(status) + call ESMF_AttributeGet(OBS_State, "Data Update Status", & + data_update, rc=status) + call LIS_verify(status) + + data_upd = .false. + obs_unsc = LIS_rc%udef + + alarmCheck = LIS_isAlarmRinging(LIS_rc, "NASASMAP read alarm") + + smobs_A = LIS_rc%udef + smobs_D = LIS_rc%udef + + cyr = LIS_rc%yr + cmo = LIS_rc%mo + cda = LIS_rc%da + chr = LIS_rc%hr + cmn = LIS_rc%mn + css = LIS_rc%ss + + call LIS_tick(time1,doy,gmt,cyr,cmo,cda,chr,cmn,css,0.0) + nyr = LIS_rc%yr + nmo = LIS_rc%mo + nda = LIS_rc%da + nhr = LIS_rc%hr + nmn = LIS_rc%mn + nss = LIS_rc%ss + + call LIS_tick(time2,doy,gmt,nyr,nmo,nda,nhr,nmn,nss,3600.0) + nyr = LIS_rc%yr + nmo = LIS_rc%mo + nda = LIS_rc%da + nhr = LIS_rc%hr + nmn = LIS_rc%mn + nss = LIS_rc%ss + + call LIS_tick(time3,doy,gmt,nyr,nmo,nda,nhr,nmn,nss,LIS_rc%ts) + + if (alarmCheck .or. cdfT_SMAPsm_struc(n)%startMode) then + cdfT_SMAPsm_struc(n)%startMode = .false. + if ( (cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP_E") .or. & + (cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP") ) then + + cdfT_SMAPsm_struc(n)%smobs = LIS_rc%udef + cdfT_SMAPsm_struc(n)%smtime = -1.0 + + write(yyyymmdd,'(i4.4,2i2.2)') LIS_rc%yr, LIS_rc%mo, LIS_rc%da + write(yyyy,'(i4.4)') LIS_rc%yr + write(mm,'(i2.2)') LIS_rc%mo + write(dd,'(i2.2)') LIS_rc%da + write(hh,'(i2.2)') LIS_rc%hr + + if(LIS_masterproc) then + list_files = trim(smobsdir)//'/'//trim(yyyy)//'.'//trim(mm)//'.'//& + trim(dd)//'/SMAP_L2_*' & + //trim(yyyy)//trim(mm)//trim(dd)//'T'//trim(hh)//'*.h5' + write(LIS_logunit,*) & + '[INFO] Searching for ',trim(list_files) + rc = create_filelist(trim(list_files)//char(0), & + "SMAP_filelist.sm.dat"//char(0)) + if (rc .ne. 0) then + write(LIS_logunit,*) & + '[WARN] Problem encountered when searching for SMAP files' + write(LIS_logunit,*) & + 'Was searching for ',trim(list_files) + write(LIS_logunit,*) & + 'LIS will continue...' + endif + end if +#if (defined SPMD) + call mpi_barrier(lis_mpi_comm,ierr) +#endif + + ftn = LIS_getNextUnitNumber() + open(ftn,file="./SMAP_filelist.sm.dat",status='old',iostat=ierr) + + do while(ierr.eq.0) + read(ftn,'(a)',iostat=ierr) fname + if(ierr.ne.0) then + exit + endif + + mn_ind = index(fname,trim(yyyymmdd)//'T'//trim(hh))+11 + read(fname(mn_ind:mn_ind+1),'(i2.2)') mn + ss=0 + call LIS_tick(timenow,doy,gmt,LIS_rc%yr, LIS_rc%mo, LIS_rc%da, & + LIS_rc%hr, mn, ss, 0.0) + + write(LIS_logunit,*) '[INFO] reading ',trim(fname) + + call read_SMAPL2sm_data_cdfTransfer(n,k,fname,& + cdfT_SMAPsm_struc(n)%smobs,timenow) + enddo + call LIS_releaseUnitNumber(ftn) + + elseif (cdfT_SMAPsm_struc(n)%data_designation .eq. "SPL3SMP_E") then +!--------------------------------------------------------------------------- +! MN: create filename for 9 km product +!--------------------------------------------------------------------------- + + write (yyyy, '(i4.4)') LIS_rc%yr + write (mm, '(i2.2)') LIS_rc%mo + write (dd, '(i2.2)') LIS_rc%da + write (CRID, '(a)') cdfT_SMAPsm_struc(n)%release_number + + ! EMK...Make sure only one PET calls the file system to determine what + ! SMAP files are available. + if (LIS_masterproc) then + + list_files = trim(smobsdir)//'/'//trim(yyyy)//'.'//trim(mm)//'.'// & + trim(dd)//'/SMAP_L3_SM_P_E_' & + //trim(yyyy)//trim(mm)//trim(dd)//'_'// & + trim(CRID)//'*.h5' + write(LIS_logunit,*) & + '[INFO] Searching for ',trim(list_files) + rc = create_filelist(trim(list_files)//char(0), & + "SMAP_filelist.dat"//char(0)) + if (rc .ne. 0) then + write(LIS_logunit,*) & + '[WARN] Problem encountered when searching for SMAP files' + write(LIS_logunit,*) & + 'Was searching for ',trim(list_files) + write(LIS_logunit,*) & + 'LIS will continue...' + endif + end if +#if (defined SPMD) + call mpi_barrier(lis_mpi_comm, ierr) +#endif + + ftn = LIS_getNextUnitNumber() + open (ftn, file="./SMAP_filelist.dat", & + action='read', status='old', iostat=ierr) + +! if multiple files for the same time and orbits are present, the latest +! one will overwrite older ones, though multiple (redundant) reads occur. +! This assumes that the 'ls command' will list the files in that order. + + do while (ierr .eq. 0) + read (ftn, '(a)', iostat=ierr) fname + if (ierr .ne. 0) then + exit + endif + + write (LIS_logunit, *) '[INFO] Reading descending pass ', trim(fname) + call read_NASASMAP_E_data_cdfTransfer(n, k, 'D', fname, smobs_D) + + write (LIS_logunit, *) '[INFO] Reading ascending pass ', trim(fname) + call read_NASASMAP_E_data_cdfTransfer(n, k, 'A', fname, smobs_A) + enddo + + cdfT_SMAPsm_struc(n)%smobs = LIS_rc%udef + cdfT_SMAPsm_struc(n)%smtime = -1 + call LIS_releaseUnitNumber(ftn) +!------------------------------------------------------------------------- +! Ascending pass assumed to be at 6pm localtime and the descending +! pass is assumed to be at 6am local time +!------------------------------------------------------------------------- + do r = 1, LIS_rc%obs_lnr(k) + do c = 1, LIS_rc%obs_lnc(k) + grid_index = LIS_obs_domain(n, k)%gindex(c, r) + if (grid_index .ne. -1) then + + if (smobs_D(c + (r - 1)*LIS_rc%obs_lnc(k)) .ne. -9999.0) then + cdfT_SMAPsm_struc(n)%smobs(c, r) = & + smobs_D(c + (r - 1)*LIS_rc%obs_lnc(k)) + lon = LIS_obs_domain(n, k)%lon(c + (r - 1)*LIS_rc%obs_lnc(k)) + lhour = 6.0 + call LIS_localtime2gmt(gmt, lon, lhour, zone) + cdfT_SMAPsm_struc(n)%smtime(c, r) = gmt + + endif +!------------------------------------------------------------------------- +! The ascending data is used only over locations where descending data +! doesn't exist. +!------------------------------------------------------------------------- + if (smobs_A(c + (r - 1)*LIS_rc%obs_lnc(k)) .ne. -9999.0 .and. & + cdfT_SMAPsm_struc(n)%smobs(c, r) .eq. -9999.0) then + cdfT_SMAPsm_struc(n)%smobs(c, r) = & + smobs_A(c + (r - 1)*LIS_rc%obs_lnc(k)) + lon = LIS_obs_domain(n, k)%lon(c + (r - 1)*LIS_rc%obs_lnc(k)) + lhour = 18.0 + call LIS_localtime2gmt(gmt, lon, lhour, zone) + cdfT_SMAPsm_struc(n)%smtime(c, r) = gmt + endif + endif + enddo + enddo + + elseif (cdfT_SMAPsm_struc(n)%data_designation .eq. "SPL3SMP") then +!------------------------------------------------------------------------- +! MN: create filename for 36km product (SMAP_L3_SM_P_) +!------------------------------------------------------------------------- + + write (yyyy, '(i4.4)') LIS_rc%yr + write (mm, '(i2.2)') LIS_rc%mo + write (dd, '(i2.2)') LIS_rc%da + write (CRID, '(a)') cdfT_SMAPsm_struc(n)%release_number + + ! EMK...Make sure only one PET calls the file system to determine what + ! SMAP files are available. + if (LIS_masterproc) then + list_files = trim(smobsdir)//'/'//trim(yyyy)//'.'//trim(mm)//'.'// & + trim(dd)//'/SMAP_L3_SM_P_' & + //trim(yyyy)//trim(mm)//trim(dd)//'_'// & + trim(CRID)//'*.h5' + write(LIS_logunit,*) & + '[INFO] Searching for ',trim(list_files) + rc = create_filelist(trim(list_files)//char(0), & + "SMAP_filelist.dat"//char(0)) + if (rc .ne. 0) then + write(LIS_logunit,*) & + '[WARN] Problem encountered when searching for SMAP files' + write(LIS_logunit,*) & + 'Was searching for ',trim(list_files) + write(LIS_logunit,*) & + 'LIS will continue...' + endif + end if +#if (defined SPMD) + call mpi_barrier(lis_mpi_comm, ierr) +#endif + + ftn = LIS_getNextUnitNumber() + open (ftn, file="./SMAP_filelist.dat", & + action='read', status='old', iostat=ierr) + +! if multiple files for the same time and orbits are present, the latest +! one will overwrite older ones, though multiple (redundant) reads occur. +! This assumes that the 'ls command' will list the files in that order. + + do while (ierr .eq. 0) + read (ftn, '(a)', iostat=ierr) fname + if (ierr .ne. 0) then + exit + endif + + write (LIS_logunit, *) '[INFO] Reading descending pass ', trim(fname) + call read_NASASMAP_data_cdfTransfer(n, k, 'D', fname, smobs_D) + + write (LIS_logunit, *) '[INFO] Reading ascending pass ', trim(fname) + call read_NASASMAP_data_cdfTransfer(n, k, 'A', fname, smobs_A) + enddo + + cdfT_SMAPsm_struc(n)%smobs = LIS_rc%udef + cdfT_SMAPsm_struc(n)%smtime = -1 + call LIS_releaseUnitNumber(ftn) +!------------------------------------------------------------------------- +! Ascending pass assumed to be at 6pm localtime and the descending +! pass is assumed to be at 6am local time +!------------------------------------------------------------------------- + do r = 1, LIS_rc%obs_lnr(k) + do c = 1, LIS_rc%obs_lnc(k) + grid_index = LIS_obs_domain(n, k)%gindex(c, r) + if (grid_index .ne. -1) then + + if (smobs_D(c + (r - 1)*LIS_rc%obs_lnc(k)) .ne. -9999.0) then + cdfT_SMAPsm_struc(n)%smobs(c, r) = & + smobs_D(c + (r - 1)*LIS_rc%obs_lnc(k)) + lon = LIS_obs_domain(n, k)%lon(c + (r - 1)*LIS_rc%obs_lnc(k)) + lhour = 6.0 + call LIS_localtime2gmt(gmt, lon, lhour, zone) + cdfT_SMAPsm_struc(n)%smtime(c, r) = gmt + + endif +!------------------------------------------------------------------------- +! The ascending data is used only over locations where descending data +! doesn't exist. +!------------------------------------------------------------------------- + if (smobs_A(c + (r - 1)*LIS_rc%obs_lnc(k)) .ne. -9999.0 .and. & + cdfT_SMAPsm_struc(n)%smobs(c, r) .eq. -9999.0) then + cdfT_SMAPsm_struc(n)%smobs(c, r) = & + smobs_A(c + (r - 1)*LIS_rc%obs_lnc(k)) + lon = LIS_obs_domain(n, k)%lon(c + (r - 1)*LIS_rc%obs_lnc(k)) + lhour = 18.0 + call LIS_localtime2gmt(gmt, lon, lhour, zone) + cdfT_SMAPsm_struc(n)%smtime(c, r) = gmt + endif + endif + enddo + enddo + + endif ! sensor + endif ! alram + + call ESMF_StateGet(OBS_State, "Observation01", smfield, & + rc=status) + call LIS_verify(status, 'Error: StateGet Observation01') + + call ESMF_FieldGet(smfield, localDE=0, farrayPtr=obsl, rc=status) + call LIS_verify(status, 'Error: FieldGet') + + fnd = 0 + sm_current = LIS_rc%udef + + ! dt is not defined as absolute value of the time difference to avoid + ! double counting of the data in assimilation. + + if ( (cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP_E") .or. & + (cdfT_SMAPsm_struc(n)%data_designation.eq."SPL2SMP") ) then + + do r=1,LIS_rc%obs_lnr(k) + do c=1,LIS_rc%obs_lnc(k) + if(LIS_obs_domain(n,k)%gindex(c,r).ne.-1) then + grid_index = c+(r-1)*LIS_rc%obs_lnc(k) + dt = (cdfT_SMAPsm_struc(n)%smtime(c,r)-time1) + if(dt.ge.0.and.dt.lt.(time3-time1)) then + sm_current(c,r) = & + cdfT_SMAPsm_struc(n)%smobs(c,r) + if(LIS_obs_domain(n,k)%gindex(c,r).ne.-1) then + obs_unsc(LIS_obs_domain(n,k)%gindex(c,r)) = & + sm_current(c,r) + endif + if(sm_current(c,r).ne.LIS_rc%udef) then + fnd = 1 + endif + endif + endif + enddo + enddo + + else + + do r = 1, LIS_rc%obs_lnr(k) + do c = 1, LIS_rc%obs_lnc(k) + if (LIS_obs_domain(n, k)%gindex(c, r) .ne. -1) then + grid_index = c + (r - 1)*LIS_rc%obs_lnc(k) + + dt = (LIS_rc%gmt - cdfT_SMAPsm_struc(n)%smtime(c, r))*3600.0 + if (dt .ge. 0 .and. dt .lt. LIS_rc%ts) then + sm_current(c, r) = & + cdfT_SMAPsm_struc(n)%smobs(c, r) + if (LIS_obs_domain(n, k)%gindex(c, r) .ne. -1) then + obs_unsc(LIS_obs_domain(n, k)%gindex(c, r)) = & + sm_current(c, r) + endif + if (sm_current(c, r) .ne. LIS_rc%udef) then + fnd = 1 + endif + endif + endif + enddo + enddo + endif + + +!------------------------------------------------------------------------- +! Transform data to the LSM climatology using a geolocation independent +! CDF-scaling approach +!------------------------------------------------------------------------- + if (cdfT_SMAPsm_struc(n)%useCDFtransfer .gt. 0) then + if (LIS_rc%da .eq. 1 .and. LIS_rc%hr .eq. 0 .and. & + LIS_rc%mn .eq. 0 .and. LIS_rc%ss .eq. 0) then + call read_CDFtransferdata_all(n,k,& + cdfT_SMAPsm_struc(n)%nbins,& + cdfT_SMAPsm_struc(n)%ntimes,& + cdfT_SMAPsm_struc(n)%n_strat_bins, & + cdfT_SMAPsm_struc(n)%modelcdffile, & + "SoilMoist",& + cdfT_SMAPsm_struc(n)%model_xrange,& + cdfT_SMAPsm_struc(n)%model_cdf) + + call read_CDFtransferdata_all(n,k,& + cdfT_SMAPsm_struc(n)%nbins,& + cdfT_SMAPsm_struc(n)%ntimes,& + cdfT_SMAPsm_struc(n)%n_strat_bins, & + cdfT_SMAPsm_struc(n)%obscdffile, & + "SoilMoist",& + cdfT_SMAPsm_struc(n)%obs_xrange,& + cdfT_SMAPsm_struc(n)%obs_cdf) + endif + if (fnd .ne. 0) then + + call LIS_rescale_with_stratified_CDF(& + n, & + k, & + cdfT_SMAPsm_struc(n)%nbins, & + cdfT_SMAPsm_struc(n)%ntimes, & + MAX_SM_VALUE, & + MIN_SM_VALUE, & + cdfT_SMAPsm_struc(n)%model_xrange, & + cdfT_SMAPsm_struc(n)%obs_xrange, & + cdfT_SMAPsm_struc(n)%model_cdf, & + cdfT_SMAPsm_struc(n)%obs_cdf, & + cdfT_SMAPsm_struc(n)%ref_p_climo_maxval, & + cdfT_SMAPsm_struc(n)%target_p_climo,& + cdfT_SMAPsm_struc(n)%n_strat_bins,& + sm_current) + endif + endif + + obsl = LIS_rc%udef + do r = 1, LIS_rc%obs_lnr(k) + do c = 1, LIS_rc%obs_lnc(k) + if (LIS_obs_domain(n, k)%gindex(c, r) .ne. -1) then + obsl(LIS_obs_domain(n, k)%gindex(c, r)) = sm_current(c, r) + endif + enddo + enddo + !------------------------------------------------------------------------- + ! Apply LSM based QC and screening of observations + !------------------------------------------------------------------------- + call lsmdaqcobsstate(trim(LIS_rc%lsm)//"+" & + //trim(LIS_CDFTRANSFERNASASMAPsmobsId)//char(0), n, k, OBS_state) + + call LIS_checkForValidObs(n, k, obsl, fnd, sm_current) + + if (fnd .eq. 0) then + data_upd_flag_local = .false. + else + data_upd_flag_local = .true. + endif + +#if (defined SPMD) + call MPI_ALLGATHER(data_upd_flag_local, 1, & + MPI_LOGICAL, data_upd_flag(:), & + 1, MPI_LOGICAL, LIS_mpi_comm, status) + data_upd = any(data_upd_flag) +#else + data_upd = data_upd_flag_local +#endif + + if (data_upd) then + + do t = 1, LIS_rc%obs_ngrid(k) + gid(t) = t + if (obsl(t) .ne. -9999.0) then + assimflag(t) = 1 + else + assimflag(t) = 0 + endif + enddo + + call ESMF_AttributeSet(OBS_State, "Data Update Status", & + .true., rc=status) + call LIS_verify(status) + + if (LIS_rc%obs_ngrid(k) .gt. 0) then + call ESMF_AttributeSet(smField, "Grid Number", & + gid, itemCount=LIS_rc%obs_ngrid(k), rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(smField, "Assimilation Flag", & + assimflag, itemCount=LIS_rc%obs_ngrid(k), rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(smfield, "Unscaled Obs", & + obs_unsc, itemCount=LIS_rc%obs_ngrid(k), rc=status) + call LIS_verify(status, 'Error in setting Unscaled Obs attribute') + + endif + + if(cdfT_SMAPsm_struc(n)%useSsdevScal.eq.1) then + write(LIS_logunit,*) '[ERR] "use scaled standard deviation model" does not work ' + write(LIS_logunit,*) '[ERR] for CDF transfer method' + call LIS_endrun() + endif + else + call ESMF_AttributeSet(OBS_State, "Data Update Status", & + .false., rc=status) + call LIS_verify(status) + endif + +end subroutine read_cdfTransfer_NASASMAPsm + +!BOP +! +! !ROUTINE: read_SMAPL2sm_data_cdfTransfer +! \label{read_SMAPL2sm_data_cdfTransfer} +! +! !INTERFACE: +subroutine read_SMAPL2sm_data_cdfTransfer(n, k,fname, smobs_inp, time) +! +! !USES: + + use LIS_coreMod + use LIS_logMod + use LIS_timeMgrMod + use cdfTransfer_NASASMAPsm_Mod, only : cdfT_SMAPsm_struc + +#if (defined USE_HDF5) + use hdf5 +#endif + + implicit none +! +! !INPUT PARAMETERS: +! + integer :: n + integer :: k + character (len=*) :: fname + real :: smobs_inp(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + real*8 :: time + +! !OUTPUT PARAMETERS: +! +! +! !DESCRIPTION: +! +! +!EOP + +#if (defined USE_HDF5) + + character*100, parameter :: sm_gr_name = "Soil_Moisture_Retrieval_Data" + character*100, parameter :: sm_field_name = "soil_moisture" + character*100, parameter :: sm_qa_name = "retrieval_qual_flag" +!YK + character*100, parameter :: vwc_field_name = "vegetation_water_content" + + integer(hsize_t), dimension(1) :: dims + integer(hsize_t), dimension(1) :: maxdims + integer(hid_t) :: file_id + integer(hid_t) :: dspace_id + integer(hid_t) :: row_id, col_id + integer(hid_t) :: sm_gr_id,sm_field_id, sm_qa_id + integer(hid_t) :: sm_gr_id_A,sm_field_id_A + integer(hid_t) :: vwc_field_id ! YK + real, allocatable :: sm_field(:) + real, allocatable :: vwc_field(:)! YK + integer, allocatable :: sm_qa(:) + integer, allocatable :: ease_row(:) + integer, allocatable :: ease_col(:) + integer :: c,r,t + logical*1 :: sm_data_b(cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr) + logical*1 :: smobs_b_ip(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + real :: sm_data(cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr) + real :: smobs_ip(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + + integer :: status,ios,iret + + call h5open_f(status) + call LIS_verify(status, 'Error opening HDF fortran interface') + + call h5fopen_f(trim(fname),H5F_ACC_RDONLY_F, file_id, status) + call LIS_verify(status, 'Error opening SMAP L2 file ') + + call h5gopen_f(file_id,sm_gr_name,sm_gr_id, status) + call LIS_verify(status, 'Error opening SM group in SMAP L2 file') + + call h5dopen_f(sm_gr_id,sm_field_name,sm_field_id, status) + call LIS_verify(status, 'Error opening SM field in SMAP L2 file') + + call h5dopen_f(sm_gr_id,"EASE_row_index",row_id, status) + call LIS_verify(status, 'Error opening row index field in SMAP L2 file') + + call h5dopen_f(sm_gr_id,"EASE_column_index",col_id, status) + call LIS_verify(status, 'Error opening column index field in SMAP L2 file') + +!YK + call h5dopen_f(sm_gr_id, sm_qa_name,sm_qa_id, status) + call LIS_verify(status, 'Error opening QA field in SMAP L2 file') + +!YK + call h5dopen_f(sm_gr_id, vwc_field_name,vwc_field_id, status) + call LIS_verify(status, 'Error opening Veg water content field in SMAP L2 file') + + call h5dget_space_f(sm_field_id, dspace_id, status) + call LIS_verify(status, 'Error in h5dget_space_f: reaSMAP L2Obs') + +! Size of the arrays +! This routine returns -1 on failure, rank on success. + call h5sget_simple_extent_dims_f(dspace_id, dims, maxdims, status) + if(status.eq.-1) then + call LIS_verify(status, 'Error in h5sget_simple_extent_dims_f: readSMAP L2Obs') + endif + + allocate(sm_field(maxdims(1))) + allocate(sm_qa(maxdims(1))) !YK + allocate(vwc_field(maxdims(1))) !YK + allocate(ease_row(maxdims(1))) + allocate(ease_col(maxdims(1))) + + call h5dread_f(row_id, H5T_NATIVE_INTEGER,ease_row,dims,status) + call LIS_verify(status, 'Error extracting row index from SMAP L2 file') + + call h5dread_f(col_id, H5T_NATIVE_INTEGER,ease_col,dims,status) + call LIS_verify(status, 'Error extracting col index from SMAP L2 file') + + call h5dread_f(sm_field_id, H5T_NATIVE_REAL,sm_field,dims,status) + call LIS_verify(status, 'Error extracting SM field from SMAP L2 file') + +!YK + call h5dread_f(sm_qa_id, H5T_NATIVE_INTEGER,sm_qa,dims,status) + call LIS_verify(status, 'Error extracting SM field from SMAP L2 file') + +!YK get the vegetation water content + call h5dread_f(vwc_field_id, H5T_NATIVE_REAL,vwc_field,dims,status) + call LIS_verify(status, 'Error extracting Veg water content (AM) field from SMAP L2 file') + +!YK + call h5dclose_f(sm_qa_id,status) + call LIS_verify(status,'Error in H5DCLOSE call') + +!YK + call h5dclose_f(vwc_field_id,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5dclose_f(row_id,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5dclose_f(col_id,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5dclose_f(sm_field_id,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5gclose_f(sm_gr_id,status) + call LIS_verify(status,'Error in H5GCLOSE call') + + call h5fclose_f(file_id,status) + call LIS_verify(status,'Error in H5FCLOSE call') + + call h5close_f(status) + call LIS_verify(status,'Error in H5CLOSE call') + + sm_data = LIS_rc%udef + sm_data_b = .false. + +!--------------------------------------------------------------------YK +!grid the data in EASE projection +! The retrieval_quality_field variable's binary representation consists of bits +! that indicate whether retrieval is performed or not at a given grid cell. +! When retrieval is performed, it contains additional bits to further +! indicate the exit status and quality of the retrieval. The first bit +! indicates the recommended quality (0-means retrieval has recommended quality). + do t=1,maxdims(1) + if(ease_col(t).gt.0.and.ease_row(t).gt.0) then + sm_data(ease_col(t) + & + (ease_row(t)-1)*cdfT_SMAPsm_struc(n)%nc) = sm_field(t) + + if(vwc_field(t).gt.5) then !YK Aply QC : if VWC > 5 kg/m2 + sm_data(ease_col(t) + & + (ease_row(t)-1)*cdfT_SMAPsm_struc(n)%nc) = LIS_rc%udef + else + if(sm_data(ease_col(t) + & + (ease_row(t)-1)*cdfT_SMAPsm_struc(n)%nc).ne.-9999.0) then + if(ibits(sm_qa(t),0,1).eq.0) then + sm_data_b(ease_col(t) + & + (ease_row(t)-1)*cdfT_SMAPsm_struc(n)%nc) = .true. + else + sm_data(ease_col(t) + & + (ease_row(t)-1)*cdfT_SMAPsm_struc(n)%nc) = LIS_rc%udef + endif + endif + endif + endif + enddo +!----------------------------------------------------------------------- + +!-------------------------------------------------------------------------- +! Interpolate to the LIS running domain +!-------------------------------------------------------------------------- + call neighbor_interp(LIS_rc%obs_gridDesc(k,:), sm_data_b, sm_data, & + smobs_b_ip, smobs_ip, & + cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr, & + LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k), & + cdfT_SMAPsm_struc(n)%rlat, cdfT_SMAPsm_struc(n)%rlon,& + cdfT_SMAPsm_struc(n)%n11, LIS_rc%udef, ios) + + + deallocate(sm_field) + deallocate(ease_row) + deallocate(ease_col) + +!overwrite the input data + do r=1,LIS_rc%obs_lnr(k) + do c=1,LIS_rc%obs_lnc(k) + if(smobs_ip(c+(r-1)*LIS_rc%obs_lnc(k)).ne.-9999.0) then + smobs_inp(c,r) = & + smobs_ip(c+(r-1)*LIS_rc%obs_lnc(k)) + + cdfT_SMAPsm_struc(n)%smtime(c,r) = & + time + endif + enddo + enddo +#endif + +end subroutine read_SMAPL2sm_data_cdfTransfer + +!BOP +! +! !ROUTINE: read_NASASMAP_E_data_cdfTransfer +! \label{read_NASASMAP_E_data_cdfTransfer} +! +! !INTERFACE: +subroutine read_NASASMAP_E_data_cdfTransfer(n, k, pass, fname, smobs_ip) +! +! !USES: + + use LIS_coreMod, only : LIS_rc, LIS_domain + use LIS_logMod + use LIS_timeMgrMod + use cdfTransfer_NASASMAPsm_Mod, only : cdfT_SMAPsm_struc +#if (defined USE_HDF5) + use hdf5 +#endif + + implicit none +! +! !INPUT PARAMETERS: +! + integer :: n + integer :: k + character (len=*) :: pass + character (len=*) :: fname + real :: smobs_ip(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + + +! !OUTPUT PARAMETERS: +! +! +! !DESCRIPTION: +! This subroutine reads the SMOS NESDIS binary file and applies the data +! quality flags to filter the data. !\normalsize +! +! tb_time_seconds +! Arithmetic average of the same parameters found in the +! fore- and aft-looking groups in the input SPL1CTB granule. +! The resulting parameter thus describes the average of UTC +! acquisition times of SPL1BTB observations whose boresights +! fall within a 36 km EASE-Grid 2.0 cell. The result is then +! expressed in J2000 seconds (the number of seconds since +! 11:58:55.816 on January 1, 2000 UT). +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[fname] name of the RTNASASMAP AMSR-E file +! \item[smobs\_ip] soil moisture data processed to the LIS domain +! \end{description} +! +! +!EOP + +#if (defined USE_HDF5) + + character*100, parameter :: sm_gr_name_D = "Soil_Moisture_Retrieval_Data_AM" + character*100, parameter :: sm_field_name_D = "soil_moisture" + character*100, parameter :: sm_qa_name_D = "retrieval_qual_flag" + character*100, parameter :: sm_gr_name_A = "Soil_Moisture_Retrieval_Data_PM" + character*100, parameter :: sm_field_name_A = "soil_moisture_pm" + character*100, parameter :: sm_qa_name_A = "retrieval_qual_flag_pm" +! MN + character*100, parameter :: vwc_field_name_D = "vegetation_water_content" + character*100, parameter :: vwc_field_name_A = "vegetation_water_content_pm" + + integer(hsize_t), allocatable :: dims(:) + integer(hsize_t), dimension(2) :: dimsm + integer(hsize_t), dimension(2) :: count_file + integer(hsize_t), dimension(2) :: count_mem + integer(hid_t) :: memspace + integer(hid_t) :: dataspace + integer :: memrank = 2 ! scaler--> rank = 0 ; 2D array--> rank = 2 + integer(hsize_t), dimension(2) :: offset_mem = (/0,0/) + integer(hsize_t), dimension(2) :: offset_file = (/0,0/) + integer(hid_t) :: file_id + integer(hid_t) :: sm_gr_id_D,sm_field_id_D,sm_qa_id_D + integer(hid_t) :: sm_gr_id_A,sm_field_id_A,sm_qa_id_A + integer(hid_t) :: vwc_field_id_D ! MN + integer(hid_t) :: vwc_field_id_A ! MN + real, allocatable :: sm_field(:,:) + real, allocatable :: vwc_field(:,:)! MN + integer, allocatable :: sm_qa(:,:) + integer :: c,r,t + logical*1 :: sm_data_b(cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr) + logical*1 :: smobs_b_ip(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + real :: sm_data(cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr) + integer :: status,ios + + dimsm = (/cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr/) + count_file = (/cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr/) + count_mem = (/cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr/) + + allocate(sm_field(cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr)) + allocate(sm_qa(cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr)) + allocate(vwc_field(cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr)) + allocate(dims(2)) + + dims(1) = cdfT_SMAPsm_struc(n)%nc + dims(2) = cdfT_SMAPsm_struc(n)%nr + + call h5open_f(status) + call LIS_verify(status, 'Error opening HDF fortran interface') + + call h5fopen_f(trim(fname),H5F_ACC_RDONLY_F, file_id, status) + call LIS_verify(status, 'Error opening NASASMAP file ') + + if(pass.eq.'D') then + call h5gopen_f(file_id,sm_gr_name_D,sm_gr_id_D, status) + call LIS_verify(status, 'Error opening SM group in NASASMAP file') + + call h5dopen_f(sm_gr_id_D,sm_field_name_D,sm_field_id_D, status) + call LIS_verify(status, 'Error opening SM field in NASASMAP file') + + call h5dget_space_f(sm_field_id_D, dataspace, status) + call LIS_verify(status, 'Error in h5dget_space_f: readNASASMAPObs') + + call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, & + start=offset_file, count=count_file, hdferr=status) + call LIS_verify(status, 'Error setting hyperslab dataspace in readNASASMAPObs') + + call h5screate_simple_f(memrank,dimsm, memspace, status) + call LIS_verify(status, 'Error in h5create_simple_f; read_cdfTransfer_NASASMAPsm') + + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, & + start=offset_mem, count=count_mem, hdferr=status) + call LIS_verify(status, 'Error in h5sselect_hyperslab_f: read_cdfTransfer_NASASMAPsm') + + call h5dread_f(sm_field_id_D, H5T_NATIVE_REAL,sm_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM field from NASASMAPfile') + + call h5dclose_f(sm_field_id_D,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5dopen_f(sm_gr_id_D,sm_qa_name_D,sm_qa_id_D, status) + call LIS_verify(status, 'Error opening SM QA field in NASASMAP file') + + call h5dread_f(sm_qa_id_D, H5T_NATIVE_INTEGER,sm_qa,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM QA field from NASASMAPfile') + + call h5dclose_f(sm_qa_id_D,status) + call LIS_verify(status,'Error in H5DCLOSE call') + +! MN get the vegetation water contnent + call h5dopen_f(sm_gr_id_D,vwc_field_name_D,vwc_field_id_D, status) + call LIS_verify(status, 'Error opening Veg water content field in NASASMAP file') + + call h5dread_f(vwc_field_id_D, H5T_NATIVE_REAL,vwc_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting Veg water content (AM) field from NASASMAPfile') + + call h5dclose_f(vwc_field_id_D,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5gclose_f(sm_gr_id_D,status) + call LIS_verify(status,'Error in H5GCLOSE call') + + else + call h5gopen_f(file_id,sm_gr_name_A,sm_gr_id_A, status) + call LIS_verify(status, 'Error opening SM group in NASASMAP file') + + call h5dopen_f(sm_gr_id_A,sm_field_name_A,sm_field_id_A, status) + call LIS_verify(status, 'Error opening SM field in NASASMAP file') + + call h5dget_space_f(sm_field_id_A, dataspace, status) + call LIS_verify(status, 'Error in h5dget_space_f: readNASASMAPObs') + + call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, & + start=offset_file, count=count_file, hdferr=status) + call LIS_verify(status, 'Error setting hyperslab dataspace in readNASASMAPObs') + + call h5screate_simple_f(memrank,dimsm, memspace, status) + call LIS_verify(status, 'Error in h5create_simple_f; read_cdfTransfer_NASASMAPsm') + + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, & + start=offset_mem, count=count_mem, hdferr=status) + call LIS_verify(status, 'Error in h5sselect_hyperslab_f: read_cdfTransfer_NASASMAPsm') + + call h5dread_f(sm_field_id_A, H5T_NATIVE_REAL,sm_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM field from NASASMAPfile') + + call h5dclose_f(sm_field_id_A,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5dopen_f(sm_gr_id_A,sm_qa_name_A,sm_qa_id_A, status) + call LIS_verify(status, 'Error opening SM QA field in NASASMAP file') + + call h5dread_f(sm_qa_id_A, H5T_NATIVE_INTEGER,sm_qa,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM QA field from NASASMAPfile') + + call h5dclose_f(sm_qa_id_A,status) + call LIS_verify(status,'Error in H5DCLOSE call') + +! MN get the vegetation water contnent + call h5dopen_f(sm_gr_id_A,vwc_field_name_A,vwc_field_id_A, status) + call LIS_verify(status, 'Error opening Veg water content field in NASASMAP file') + + call h5dread_f(vwc_field_id_A, H5T_NATIVE_REAL,vwc_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting Veg water content (AM) field from NASASMAPfile') + + call h5dclose_f(vwc_field_id_A,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5gclose_f(sm_gr_id_A,status) + call LIS_verify(status,'Error in H5GCLOSE call') + + endif + + call h5fclose_f(file_id,status) + call LIS_verify(status,'Error in H5FCLOSE call') + + call h5close_f(status) + call LIS_verify(status,'Error in H5CLOSE call') + + sm_data_b = .false. + t = 1 + +! The retrieval_quality_field variable's binary representation consists of bits +! that indicate whether retrieval is performed or not at a given grid cell. +! When retrieval is performed, it contains additional bits to further +! indicate the exit status and quality of the retrieval. The first bit +! indicates the recommended quality (0-means retrieval has recommended quality). + + do r=1,cdfT_SMAPsm_struc(n)%nr + do c=1,cdfT_SMAPsm_struc(n)%nc + sm_data(t) = sm_field(c,r) + + if(vwc_field(c,r).gt. 5 ) then !MN Aply QC : if VWC > 5 kg/m2 + sm_data(t) = LIS_rc%udef + else + + if(sm_data(t).ne.-9999.0) then + if(ibits(sm_qa(c,r),0,1).eq.0) then + sm_data_b(t) = .true. + else + sm_data(t) = -9999.0 + endif + endif + endif + + t = t+1 + enddo + enddo + + deallocate(sm_qa) + +!-------------------------------------------------------------------------- +! Interpolate to the LIS running domain +!-------------------------------------------------------------------------- + call neighbor_interp(LIS_rc%obs_gridDesc(k,:),& + sm_data_b, sm_data, smobs_b_ip, smobs_ip, & + cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr, & + LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k), & + cdfT_SMAPsm_struc(n)%rlat, cdfT_SMAPsm_struc(n)%rlon,& + cdfT_SMAPsm_struc(n)%n11, LIS_rc%udef, ios) + + deallocate(sm_field) + deallocate(dims) + +#endif + +end subroutine read_NASASMAP_E_data_cdfTransfer + +! MN: the data structure in both 36 km and 9 km products is the same therefore +! read_NASASMAP_E_data_cdfTransfer is similar to read_NASASMAP_data_cdfTransfer + +!BOP +! +! !ROUTINE: read_NASASMAP_data_cdfTransfer +! \label{read_NASASMAP_data_cdfTransfer} +! +! !INTERFACE: +subroutine read_NASASMAP_data_cdfTransfer(n, k, pass, fname, smobs_ip) +! +! !USES: + + use LIS_coreMod, only : LIS_rc, LIS_domain + use LIS_logMod + use LIS_timeMgrMod + use cdfTransfer_NASASMAPsm_Mod, only : cdfT_SMAPsm_struc +#if (defined USE_HDF5) + use hdf5 +#endif + + implicit none +! +! !INPUT PARAMETERS: +! + integer :: n + integer :: k + character (len=*) :: pass + character (len=*) :: fname + real :: smobs_ip(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + + +! !OUTPUT PARAMETERS: +! +! +! !DESCRIPTION: +! This subroutine reads the SMOS NESDIS binary file and applies the data +! quality flags to filter the data. !\normalsize +! +! tb_time_seconds +! Arithmetic average of the same parameters found in the +! fore- and aft-looking groups in the input SPL1CTB granule. +! The resulting parameter thus describes the average of UTC +! acquisition times of SPL1BTB observations whose boresights +! fall within a 36 km EASE-Grid 2.0 cell. The result is then +! expressed in J2000 seconds (the number of seconds since +! 11:58:55.816 on January 1, 2000 UT). +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[fname] name of the RTNASASMAP AMSR-E file +! \item[smobs\_ip] soil moisture data processed to the LIS domain +! \end{description} +! +! +!EOP + +#if (defined USE_HDF5) + + character*100, parameter :: sm_gr_name_D = "Soil_Moisture_Retrieval_Data_AM" + character*100, parameter :: sm_field_name_D = "soil_moisture" + character*100, parameter :: sm_qa_name_D = "retrieval_qual_flag" + character*100, parameter :: sm_gr_name_A = "Soil_Moisture_Retrieval_Data_PM" + character*100, parameter :: sm_field_name_A = "soil_moisture_pm" + character*100, parameter :: sm_qa_name_A = "retrieval_qual_flag_pm" +! MN + character*100, parameter :: vwc_field_name_D = "vegetation_water_content" + character*100, parameter :: vwc_field_name_A = "vegetation_water_content_pm" + + integer(hsize_t), allocatable :: dims(:) + integer(hsize_t), dimension(2) :: dimsm + integer(hsize_t), dimension(2) :: count_file + integer(hsize_t), dimension(2) :: count_mem + integer(hid_t) :: memspace + integer(hid_t) :: dataspace + integer :: memrank = 2 ! scaler--> rank = 0 ; 2D array--> rank = 2 + integer(hsize_t), dimension(2) :: offset_mem = (/0,0/) + integer(hsize_t), dimension(2) :: offset_file = (/0,0/) + integer(hid_t) :: file_id + integer(hid_t) :: sm_gr_id_D,sm_field_id_D,sm_qa_id_D + integer(hid_t) :: sm_gr_id_A,sm_field_id_A,sm_qa_id_A + integer(hid_t) :: vwc_field_id_D ! MN + integer(hid_t) :: vwc_field_id_A ! MN + real, allocatable :: sm_field(:,:) + real, allocatable :: vwc_field(:,:)! MN + integer, allocatable :: sm_qa(:,:) + integer :: c,r,t + logical*1 :: sm_data_b(cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr) + logical*1 :: smobs_b_ip(LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k)) + real :: sm_data(cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr) + integer :: status,ios + + dimsm = (/cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr/) + count_file = (/cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr/) + count_mem = (/cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr/) + + allocate(sm_field(cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr)) + allocate(sm_qa(cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr)) + allocate(vwc_field(cdfT_SMAPsm_struc(n)%nc, cdfT_SMAPsm_struc(n)%nr)) + allocate(dims(2)) + + dims(1) = cdfT_SMAPsm_struc(n)%nc + dims(2) = cdfT_SMAPsm_struc(n)%nr + + call h5open_f(status) + call LIS_verify(status, 'Error opening HDF fortran interface') + + call h5fopen_f(trim(fname),H5F_ACC_RDONLY_F, file_id, status) + call LIS_verify(status, 'Error opening NASASMAP file ') + + if(pass.eq.'D') then + call h5gopen_f(file_id,sm_gr_name_D,sm_gr_id_D, status) + call LIS_verify(status, 'Error opening SM group in NASASMAP file') + + call h5dopen_f(sm_gr_id_D,sm_field_name_D,sm_field_id_D, status) + call LIS_verify(status, 'Error opening SM field in NASASMAP file') + + call h5dget_space_f(sm_field_id_D, dataspace, status) + call LIS_verify(status, 'Error in h5dget_space_f: readNASASMAPObs') + + call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, & + start=offset_file, count=count_file, hdferr=status) + call LIS_verify(status, 'Error setting hyperslab dataspace in readNASASMAPObs') + + call h5screate_simple_f(memrank,dimsm, memspace, status) + call LIS_verify(status, 'Error in h5create_simple_f; read_cdfTransfer_NASASMAPsm') + + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, & + start=offset_mem, count=count_mem, hdferr=status) + call LIS_verify(status, 'Error in h5sselect_hyperslab_f: read_cdfTransfer_NASASMAPsm') + + call h5dread_f(sm_field_id_D, H5T_NATIVE_REAL,sm_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM field from NASASMAPfile') + + call h5dclose_f(sm_field_id_D,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5dopen_f(sm_gr_id_D,sm_qa_name_D,sm_qa_id_D, status) + call LIS_verify(status, 'Error opening SM QA field in NASASMAP file') + + call h5dread_f(sm_qa_id_D, H5T_NATIVE_INTEGER,sm_qa,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM QA field from NASASMAPfile') + + call h5dclose_f(sm_qa_id_D,status) + call LIS_verify(status,'Error in H5DCLOSE call') + +! MN get the vegetation water contnent + call h5dopen_f(sm_gr_id_D,vwc_field_name_D,vwc_field_id_D, status) + call LIS_verify(status, 'Error opening Veg water content field in NASASMAP file') + + call h5dread_f(vwc_field_id_D, H5T_NATIVE_REAL,vwc_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting Veg water content (AM) field from NASASMAPfile') + + call h5dclose_f(vwc_field_id_D,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5gclose_f(sm_gr_id_D,status) + call LIS_verify(status,'Error in H5GCLOSE call') + + else + call h5gopen_f(file_id,sm_gr_name_A,sm_gr_id_A, status) + call LIS_verify(status, 'Error opening SM group in NASASMAP file') + + call h5dopen_f(sm_gr_id_A,sm_field_name_A,sm_field_id_A, status) + call LIS_verify(status, 'Error opening SM field in NASASMAP file') + + call h5dget_space_f(sm_field_id_A, dataspace, status) + call LIS_verify(status, 'Error in h5dget_space_f: readNASASMAPObs') + + call h5sselect_hyperslab_f(dataspace, H5S_SELECT_SET_F, & + start=offset_file, count=count_file, hdferr=status) + call LIS_verify(status, 'Error setting hyperslab dataspace in readNASASMAPObs') + + call h5screate_simple_f(memrank,dimsm, memspace, status) + call LIS_verify(status, 'Error in h5create_simple_f; read_cdfTransfer_NASASMAPsm') + + call h5sselect_hyperslab_f(memspace, H5S_SELECT_SET_F, & + start=offset_mem, count=count_mem, hdferr=status) + call LIS_verify(status, 'Error in h5sselect_hyperslab_f: read_cdfTransfer_NASASMAPsm') + + call h5dread_f(sm_field_id_A, H5T_NATIVE_REAL,sm_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM field from NASASMAPfile') + + call h5dclose_f(sm_field_id_A,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5dopen_f(sm_gr_id_A,sm_qa_name_A,sm_qa_id_A, status) + call LIS_verify(status, 'Error opening SM QA field in NASASMAP file') + + call h5dread_f(sm_qa_id_A, H5T_NATIVE_INTEGER,sm_qa,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting SM QA field from NASASMAPfile') + + call h5dclose_f(sm_qa_id_A,status) + call LIS_verify(status,'Error in H5DCLOSE call') + +! MN get the vegetation water contnent + call h5dopen_f(sm_gr_id_A,vwc_field_name_A,vwc_field_id_A, status) + call LIS_verify(status, 'Error opening Veg water content field in NASASMAP file') + + call h5dread_f(vwc_field_id_A, H5T_NATIVE_REAL,vwc_field,dims,status, & + memspace, dataspace) + call LIS_verify(status, 'Error extracting Veg water content (AM) field from NASASMAPfile') + + call h5dclose_f(vwc_field_id_A,status) + call LIS_verify(status,'Error in H5DCLOSE call') + + call h5gclose_f(sm_gr_id_A,status) + call LIS_verify(status,'Error in H5GCLOSE call') + + endif + + call h5fclose_f(file_id,status) + call LIS_verify(status,'Error in H5FCLOSE call') + + call h5close_f(status) + call LIS_verify(status,'Error in H5CLOSE call') + + sm_data_b = .false. + t = 1 + +! The retrieval_quality_field variable's binary representation consists of bits +! that indicate whether retrieval is performed or not at a given grid cell. +! When retrieval is performed, it contains additional bits to further +! indicate the exit status and quality of the retrieval. The first bit +! indicates the recommended quality (0-means retrieval has recommended quality). + + do r=1,cdfT_SMAPsm_struc(n)%nr + do c=1,cdfT_SMAPsm_struc(n)%nc + sm_data(t) = sm_field(c,r) + + if(vwc_field(c,r).gt. 5 ) then !MN Aply QC : if VWC > 5 kg/m2 + sm_data(t) = LIS_rc%udef + else + + if(sm_data(t).ne.-9999.0) then + if(ibits(sm_qa(c,r),0,1).eq.0) then + sm_data_b(t) = .true. + else + sm_data(t) = -9999.0 + endif + endif + endif + + t = t+1 + enddo + enddo + + deallocate(sm_qa) + +!-------------------------------------------------------------------------- +! Interpolate to the LIS running domain +!-------------------------------------------------------------------------- + call neighbor_interp(LIS_rc%obs_gridDesc(k,:),& + sm_data_b, sm_data, smobs_b_ip, smobs_ip, & + cdfT_SMAPsm_struc(n)%nc*cdfT_SMAPsm_struc(n)%nr, & + LIS_rc%obs_lnc(k)*LIS_rc%obs_lnr(k), & + cdfT_SMAPsm_struc(n)%rlat, cdfT_SMAPsm_struc(n)%rlon,& + cdfT_SMAPsm_struc(n)%n11, LIS_rc%udef, ios) + + deallocate(sm_field) + deallocate(dims) + +#endif + + +end subroutine read_NASASMAP_data_cdfTransfer + + diff --git a/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/write_cdfTransfer_NASASMAPsmobs.F90 b/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/write_cdfTransfer_NASASMAPsmobs.F90 new file mode 100644 index 000000000..48aa3cbd7 --- /dev/null +++ b/lis/dataassim/obs/CDF_Transfer_NASA_SMAPsm/write_cdfTransfer_NASASMAPsmobs.F90 @@ -0,0 +1,122 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.3 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +!BOP +! +! !ROUTINE: write_cdfTransfer_NASASMAPsmobs +! \label{write_cdfTransfer_NASASMAPsmobs} +! +! !REVISION HISTORY: +! 25Jan2008: Sujay Kumar; Initial Specification +! 2 Mar 2022: Mahdi navari; modified for cdf transfer DA +! +! !INTERFACE: +subroutine write_cdfTransfer_NASASMAPsmobs(n, k, OBS_State) +! !USES: + use ESMF + use LIS_coreMod + use LIS_logMod + use LIS_fileIOMod + use LIS_historyMod + use LIS_DAobservationsMod + + implicit none + +! !ARGUMENTS: + + integer, intent(in) :: n + integer, intent(in) :: k + type(ESMF_State) :: OBS_State +! +! !DESCRIPTION: +! +! writes the transformed (interpolated/upscaled/reprojected) +! LPRM AMSRE observations to a file +! +!EOP + type(ESMF_Field) :: smField + logical :: data_update + real, pointer :: smobs(:) + real :: smobs_unsc(LIS_rc%obs_ngrid(k)) + character*100 :: obsname + integer :: ftn + integer :: status + + call ESMF_AttributeGet(OBS_State, "Data Update Status", & + data_update, rc=status) + call LIS_verify(status) + + if(data_update) then + + call ESMF_StateGet(OBS_State, "Observation01",smField, & + rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(smField, localDE=0, farrayPtr=smobs, rc=status) + call LIS_verify(status) + + if(LIS_rc%obs_ngrid(k).gt.0) then + call ESMF_AttributeGet(smfield,"Unscaled Obs",smobs_unsc,& + itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status, 'Error in AttributeGet-Unscaled Obs') + endif + + if(LIS_masterproc) then + ftn = LIS_getNextUnitNumber() + call cdfTransfer_smobsname(n,k,obsname) + + call LIS_create_output_directory('DAOBS') + open(ftn,file=trim(obsname), form='unformatted') + endif + + call LIS_writevar_gridded_obs(ftn,n,k,smobs_unsc) + call LIS_writevar_gridded_obs(ftn,n,k,smobs) + + if(LIS_masterproc) then + call LIS_releaseUnitNumber(ftn) + endif + + endif + +end subroutine write_cdfTransfer_NASASMAPsmobs + +!BOP +! !ROUTINE: cdfTransfer_smobsname +! \label{cdfTransfer_smobsname} +! +! !INTERFACE: +subroutine cdfTransfer_smobsname(n,k,obsname) +! !USES: + use LIS_coreMod, only : LIS_rc + +! !ARGUMENTS: + integer :: n + integer :: k + character(len=*) :: obsname +! +! !DESCRIPTION: +! +!EOP + + character(len=12) :: cdate1 + character(len=12) :: cdate + character(len=10) :: cda + + write(unit=cdate1, fmt='(i4.4, i2.2, i2.2, i2.2, i2.2)') & + LIS_rc%yr, LIS_rc%mo, & + LIS_rc%da, LIS_rc%hr,LIS_rc%mn + + write(unit=cda, fmt='(a2,i2.2)') '.a',k + write(unit=cdate, fmt='(a2,i2.2)') '.d',n + + obsname = trim(LIS_rc%odir)//'/DAOBS/'//cdate1(1:6)//& + '/LISDAOBS_'//cdate1// & + trim(cda)//trim(cdate)//'.1gs4r' + +end subroutine cdfTransfer_smobsname diff --git a/lis/make/default.cfg b/lis/make/default.cfg index 363a09167..87a51f553 100644 --- a/lis/make/default.cfg +++ b/lis/make/default.cfg @@ -866,6 +866,10 @@ enabled: True macro: DA_OBS_GRACE path: dataassim/obs/GRACE +[DA CDF TRANSFER NASA SMAPSM] +enabled: True +macro: DA_CDF_TRANSFER_NASA_SMAPSM +path: dataassim/obs/CDF_Transfer_NASA_SMAPsm #}}} # diff --git a/lis/plugins/LIS_DAobs_pluginMod.F90 b/lis/plugins/LIS_DAobs_pluginMod.F90 index 14491b275..680b12dc3 100644 --- a/lis/plugins/LIS_DAobs_pluginMod.F90 +++ b/lis/plugins/LIS_DAobs_pluginMod.F90 @@ -251,6 +251,10 @@ subroutine LIS_DAobs_plugin use NASASMAPsm_Mod, only : NASASMAPsm_setup #endif +#if ( defined DA_CDF_TRANSFER_NASA_SMAPSM ) + use cdfTransfer_NASASMAPsm_Mod, only : cdfTransfer_NASASMAPsm_setup +#endif + !YK #if ( defined DA_OBS_SMOS_NRT_NN ) use SMOSNRTNNL2sm_Mod, only : SMOSNRTNNL2sm_setup @@ -451,6 +455,10 @@ subroutine LIS_DAobs_plugin external read_NASASMAPsm, write_NASASMAPsmobs #endif +#if ( defined DA_CDF_TRANSFER_NASA_SMAPSM ) + external read_cdfTransfer_NASASMAPsm, write_cdfTransfer_NASASMAPsmobs +#endif + !YK #if ( defined DA_OBS_SMOS_NRT_NN ) external read_SMOSNRTNNL2sm, write_SMOSNRTNNL2smobs @@ -838,6 +846,16 @@ subroutine LIS_DAobs_plugin write_NASASMAPsmobs) #endif +#if ( defined DA_CDF_TRANSFER_NASA_SMAPSM ) + call registerdaobsclass(trim(LIS_CDFTRANSFERNASASMAPsmobsId),"LSM") + call registerdaobssetup(trim(LIS_CDFTRANSFERNASASMAPsmobsId)//char(0),& + cdfTransfer_NASASMAPsm_setup) + call registerreaddaobs(trim(LIS_CDFTRANSFERNASASMAPsmobsId)//char(0),& + read_cdfTransfer_NASASMAPsm) + call registerwritedaobs(trim(LIS_CDFTRANSFERNASASMAPsmobsId)//char(0),& + write_cdfTransfer_NASASMAPsmobs) +#endif + !YK #if ( defined DA_OBS_SMOS_NRT_NN ) call registerdaobsclass(trim(LIS_SMOSNRTNNL2smobsId),"LSM") diff --git a/lis/plugins/LIS_lsmda_pluginMod.F90 b/lis/plugins/LIS_lsmda_pluginMod.F90 index d7915469e..98e892393 100644 --- a/lis/plugins/LIS_lsmda_pluginMod.F90 +++ b/lis/plugins/LIS_lsmda_pluginMod.F90 @@ -2689,6 +2689,26 @@ subroutine LIS_lsmda_plugin trim(LIS_NASASMAPsmobsId )//char(0),NoahMP401_descale_soilm) call registerlsmdaupdatestate(trim(LIS_noahmp401Id)//"+"//& trim(LIS_NASASMAPsmobsId )//char(0),NoahMP401_updatesoilm) +!MN +! Noah-MP.4.0.1 SMAP(NASA) soil moisture with CDF Transfer + call registerlsmdainit(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_dasoilm_init) + call registerlsmdagetstatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_getsoilm) + call registerlsmdasetstatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_setsoilm) + call registerlsmdagetobspred(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_getsmpred) + call registerlsmdaqcstate(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_qcsoilm) + call registerlsmdaqcobsstate(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_qc_soilmobs) + call registerlsmdascalestatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_scale_soilm) + call registerlsmdadescalestatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_descale_soilm) + call registerlsmdaupdatestate(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_CDFTRANSFERNASASMAPsmobsId )//char(0),NoahMP401_updatesoilm) !YK ! Noah-MP.4.0.1 SMOS NRT NN soil moisture diff --git a/lis/plugins/LIS_pluginIndices.F90 b/lis/plugins/LIS_pluginIndices.F90 index 352625ac3..76536bbc3 100644 --- a/lis/plugins/LIS_pluginIndices.F90 +++ b/lis/plugins/LIS_pluginIndices.F90 @@ -272,6 +272,8 @@ module LIS_pluginIndices "SMOS(NESDIS) soil moisture" character*50, public, parameter :: LIS_NASASMAPsmobsId = & "SMAP(NASA) soil moisture" + character*50, public, parameter :: LIS_CDFTRANSFERNASASMAPsmobsId = & + "SMAP(NASA) soil moisture with CDF Transfer" !MN character*50, public, parameter :: LIS_SMOSNRTNNL2smobsId = & "SMOS NRT NN soil moisture" !YK character*50, public, parameter :: LIS_SMAPEOPLsmobsId = & diff --git a/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/NOAHMP36_OUTPUT_LIST.TBL.dyn.full b/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/NOAHMP36_OUTPUT_LIST.TBL.dyn.full new file mode 100755 index 000000000..924c0f29d --- /dev/null +++ b/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/NOAHMP36_OUTPUT_LIST.TBL.dyn.full @@ -0,0 +1,218 @@ +#short_name select? units signconv timeavg? min/max? std? vert.levels grib_id grib_scalefactor longname +# 0 == Instantaneous output +# 1 == Time averaged output +# 2 == Instantaneous and Time averaged output +# 3 == Accumulated output + +#Energy balance components +Swnet: 0 W/m2 DN 1 0 0 1 111 10 # Net shortwave radiation (W/m2) +Lwnet: 0 W/m2 DN 1 0 0 1 112 10 # Net longwave radiation (W/m2) +Qle: 0 W/m2 UP 1 0 0 1 121 10 # Latent heat flux (W/m2) +Qh: 0 W/m2 UP 1 0 0 1 122 10 # Sensible heat flux (W/m2) +Qg: 0 W/m2 DN 1 0 0 1 155 10 # Ground heat flux (W/m2) +Qf: 0 W/m2 S2L 1 0 0 1 229 10 # Energy of fusion (W/m2) +Qv: 0 W/m2 S2V 1 0 0 1 198 10 # Energy of sublimation (W/m2) +Qa: 0 W/m2 DN 1 0 0 1 136 10 # Advective energy (W/m2) +Qtau: 0 N/m2 DN 1 0 0 1 172 10 # Momentum flux (N/m2) +DelSurfHeat: 0 J/m2 INC 1 0 0 1 137 10 # Change in surface heat storage (J/m2) +DelColdCont: 0 J/m2 INC 1 0 0 1 138 10 # Change in snow cold content (J/m2) +BR: 0 - - 1 0 1 1 256 10 # Bowen ratio (-) +EF: 0 - - 1 0 1 1 256 10 # Evaporative fraction (-) +Rnet: 0 W/m2 DN 1 0 1 1 256 10 # Total net radiation (W/m2) + +#Water balance components +Snowf: 0 kg/m2s DN 1 0 0 1 161 10000 # Snowfall rate (kg/m2s) +Rainf: 0 kg/m2s DN 1 0 0 1 162 10000 # Rainfall rate (kg/m2s) +RainfConv: 0 kg/m2s DN 1 0 0 1 63 10000 # Convective rainfall rate (kg/m2s) +TotalPrecip: 1 kg/m2s DN 1 0 0 1 61 10000 # Total precipitation rate (kg/m2s) +Evap: 1 kg/m2s UP 1 0 0 1 57 10000 # Total evapotranspiration (kg/m2s) +Qs: 1 kg/m2s OUT 1 0 0 1 235 10000 # Surface runoff (kg/m2s) +Qrec: 0 kg/m2s IN 1 0 0 1 163 10000 # Recharge (kg/m2s) +Qsb: 1 kg/m2s OUT 1 0 0 1 234 10000 # Subsurface runoff (kg/m2s) +Qsm: 0 kg/m2s S2L 1 0 0 1 99 10000 # Snowmelt (kg/m2s) +Qfz: 0 kg/m2s L2S 1 0 0 1 130 10000 # Refreezing of water in the snowpack (kg/m2s) +Qst: 0 kg/m2s - 1 0 0 1 131 10000 # Snow throughfall (kg/m2s) +DelSoilMoist: 0 kg/m2 INC 1 0 0 1 132 10000 # Change in soil moisture (kg/m2) +DelSWE: 0 kg/m2 INC 1 0 0 1 133 1000 # Change in snow water equivalent (kg/m2) +DelSurfStor: 0 kg/m2 INC 1 0 0 1 134 1000 # Change in surface water storage (kg/m2) +DelIntercept: 0 kg/m2 INC 1 0 0 1 135 1000 # Change in interception storage (kg/m2) +RHMin: 0 - - 1 0 0 1 52 10 # Minimum 2-meter relative humidity (-) +Ch: 0 - - 1 0 0 1 208 10 # Surface exchange coefficient for heat +Cm: 0 - - 1 0 0 1 252 10 # Surface exchange coefficient for momentum + +#Surface state variables +SnowT: 0 K - 1 0 0 1 165 10 # Snow surface temperature (K) +VegT: 0 K - 1 0 0 1 146 10 # Vegetation canopy temperature (K) +BareSoilT: 0 K - 1 0 0 1 147 10 # Temperature of bare soil (K) +AvgSurfT: 0 K - 1 0 0 1 148 10 # Average surface temperature (K) +RadT: 0 K - 1 0 0 1 149 10 # Surface radiative temperature (K) +Albedo: 0 - - 1 0 0 1 84 100 # Surface albedo (-) +SWE: 1 kg/m2 - 1 0 0 1 65 1000 # Snow Water Equivalent (kg/m2) +SWEVeg: 0 kg/m2 - 1 0 0 1 139 1000 # SWE intercepted by vegetation (kg/m2) +SurfStor: 0 kg/m2 - 1 0 0 1 150 1000 # Surface water storage (kg/m2) + +#Subsurface state variables +SoilMoist: 1 m3/m3 - 2 0 0 4 86 1000 # Average layer soil moisture (kg/m2) +SoilTemp: 0 K - 1 0 0 4 85 1000 # Average layer soil temperature (K) +SmLiqFrac: 0 - - 1 0 0 4 160 100 # Average layer fraction of liquid moisture (-) +SmFrozFrac: 0 - - 1 0 0 4 140 100 # Average layer fraction of frozen moisture (-) +SoilWet: 0 - - 1 0 0 1 144 100 # Total soil wetness (-) +RelSMC: 0 m3/m3 - 1 0 0 4 141 1000 # Relative soil moisture (-) +RootTemp: 0 K - 1 0 0 1 142 1000 # Root zone temperature (K) + +#Evaporation components +PotEvap: 0 kg/m2s UP 1 0 0 1 145 1 # Potential evapotranspiration (kg/m2s) +ECanop: 0 kg/m2s UP 1 0 0 1 200 1 # Interception evaporation (kg/m2s) +TVeg: 0 kg/m2s UP 1 0 0 1 210 1 # Vegetation transpiration (kg/m2s) +ESoil: 0 kg/m2s UP 1 0 0 1 199 1 # Bare soil evaporation (kg/m2s) +EWater: 0 kg/m2s UP 1 0 0 1 197 1 # Open water evaporation (kg/m2s) +RootMoist: 0 kg/m2 - 2 0 0 1 171 1 # Root zone soil moisture (kg/m2) +CanopInt: 0 kg/m2 - 1 0 0 1 223 1000 # Total canopy water storage (kg/m2) +EvapSnow: 0 kg/m2s - 1 0 0 1 173 1000 # Snow evaporation (kg/m2s) +SubSnow: 0 kg/m2s - 1 0 0 1 198 1000 # Snow sublimation (kg/m2s) +SubSurf: 0 kg/m2s - 1 0 0 1 143 1000 # Sublimation of the snow free area (kg/m2s) +ACond: 0 m/s - 1 0 0 1 179 100000 # Aerodynamic conductance (m/s) +CCond: 0 m/s - 1 0 0 1 181 100000 # Canopy conductance (m/s) +SoilET: 0 kg/m2 - 1 0 0 1 256 1 # Soil evaporation (kg/m2s) +AResist: 0 s/m - 1 0 0 1 256 1 # Aerodynamic resistance (s/m) + +#Other hydrologic variables +WaterTableD: 0 m - 1 0 0 1 174 1 # Water table depth (m) +TWS: 0 mm - 1 0 0 1 175 1 # Terrestrial Water Storage (mm) +GWS: 0 mm - 1 0 0 1 176 1 # Ground Water Storage (mm) +WT: 0 mm - 1 0 0 1 177 1 # Noah-MP WT variable (mm) + +#Cold season processes +Snowcover: 0 - - 1 0 0 1 238 100 # Snow cover (-) +SAlbedo: 0 - - 0 0 0 1 184 1000 # Albedo of the snow-covered area (-) +SnowTProf: 0 K - 0 0 0 1 239 1000 # Temperature of the snow pack (K) +SnowDepth: 0 m - 1 0 0 1 66 1000 # Snow depth (m) +SLiqFrac: 0 - - 0 0 0 1 185 1000 # Fraction of SWE in the liquid phase (-) + +#Variables to compared against remote sensed data +LWup: 0 W/m2 UP 1 0 0 1 212 1 # Longwave radiation up from the surface (W/m2) + +#Carbon variables +GPP: 0 g/m2s IN 1 0 0 1 256 1 # Gross Primary Production +NPP: 0 g/m2s IN 1 0 0 1 256 1 # Net Primary Production +NEE: 0 g/m2s IN 1 0 0 1 256 1 # Net Ecosystem Exchange +AutoResp: 0 kg/m2s2 UP 1 0 0 1 256 1 # Autotrophic respiration +HeteroResp: 0 kg/m2s2 UP 1 0 0 1 256 1 # Heterotrophic respiration +LeafResp: 0 kg/m2s2 UP 1 0 0 1 256 1 # Leaf respiration +TotSoilCarb: 0 kg/m2 - 1 0 0 1 256 1 # Total soil carbon +TotLivBiom: 0 kg/m2 - 1 0 0 1 256 1 # Total living biomass + +#Forcings +Wind_f: 1 m/s - 1 0 0 1 32 10 # Near surface wind (m/s) +Rainf_f: 1 kg/m2s DN 1 0 0 1 162 1000 # Average rainfall rate +Snowf_f: 0 kg/m2s DN 1 0 0 1 161 1000 # Average snowfall rate +CRainf_f: 0 kg/m2s DN 1 0 0 1 63 1000 # Average convective rainfall rate +Tair_f: 1 K - 1 0 0 1 11 10 # Near surface air temperature +Qair_f: 1 kg/kg - 1 0 0 1 51 1000 # Near surface specific humidity +Psurf_f: 1 Pa - 1 0 0 1 1 10 # Surface pressure +SWdown_f: 1 W/m2 DN 1 0 0 1 204 10 # Surface incident shortwave radiation +LWdown_f: 1 W/m2 DN 1 0 0 1 205 10 # Surface incident longwave radiation +PARDR_f: 0 W/m2 DN 1 0 0 1 256 10 # Surface incident PAR direct +PARDF_f: 0 W/m2 DN 1 0 0 1 256 10 # Surface incident PAR diffuse + +#Additional forcings +DirectSW_f: 0 W/m2 - 1 0 0 1 166 10 # Surface direct incident shortwave radiation +DiffuseSW_f: 0 W/m2 - 1 0 0 1 167 10 # Surface diffuse incident shortwave radiation +NWind_f: 0 m/s N 1 0 0 1 34 10 # Northward wind +EWind_f: 0 m/s E 1 0 0 1 33 10 # Eastward wind +FHeight_f: 0 m - 1 0 0 1 256 10 # Height of forcing variables +Ch_f: 0 - - 1 0 0 1 208 10 # Surface exchange coefficient for heat +Cm_f: 0 - - 1 0 0 1 252 10 # Surface exchange coefficient for momentum +Emiss_f: 0 - - 1 0 0 1 256 10 # Surface emissivity +MixRatio_f: 0 kg/kg - 1 0 0 1 53 10 # Surface mixing ratio +CosZenith_f: 0 - - 1 0 0 1 256 10 # Cosine of zenith angle +Albedo_f: 0 - - 1 0 0 1 84 10 # Surface albedo +CAPE_f: 0 J/kg - 1 0 0 1 157 10 # Convective Available Potential Energy +Z0brd: 0 m - 1 0 0 1 256 1 # Z0brd +T2diag: 0 K - 1 0 0 1 256 1 # Diagnostic t2 +Q2diag: 0 kg/kg - 1 0 0 1 256 1 # Diagnostic q2 +Snowflag_f: 0 - - 1 0 0 1 256 1 # Snowflag +Density_f: 0 kg/m3 - 1 0 0 1 256 1 # Atmospheric density +VaporPress_f: 0 - - 1 0 0 1 256 1 # Vapor pressure +VaporPressDeficit_f: 0 - - 1 0 0 1 256 1 # Vapor pressure deficit + +#Parameters +Landmask: 0 - - 0 0 0 1 81 1 # Land mask (0 - Water, 1 - Land) +Landcover: 0 - - 0 0 0 1 225 1 # Land cover +Soiltype: 0 - - 0 0 0 1 224 1 # Soil type +SandFrac: 0 - - 0 0 0 1 256 1 # Sand fraction +ClayFrac: 0 - - 0 0 0 1 256 1 # Clay fraction +SiltFrac: 0 - - 0 0 0 1 256 1 # Silt fraction +Porosity: 0 - - 0 0 0 1 240 1 # Porosity +Soilcolor: 0 - - 0 0 0 1 256 1 # Soil color +Elevation: 0 m - 0 0 0 1 196 10 # Elevation +Slope: 0 - - 0 0 0 1 222 10 # Slope +LAI: 0 - - 1 0 0 1 182 100 # LAI +SAI: 0 - - 1 0 0 1 256 100 # SAI +Snfralbedo: 0 - - 0 0 0 1 184 100 # Snow fraction albedo +Mxsnalbedo: 0 - - 0 0 0 1 159 100 # Maximum snow albedo +Greenness: 0 - - 1 0 0 1 87 100 # Greenness +Roughness: 0 m - 0 0 0 1 83 10 # Roughness +Tempbot: 0 K - 0 0 0 1 256 10 # Bottom soil temperature + +#Routing +Streamflow: 1 m3/s - 1 0 0 1 256 10 # Streamflow +RiverStor: 0 m3 - 1 0 0 1 256 10 # River storage +RiverDepth: 0 m - 1 0 0 1 256 10 # River depth +RiverVelocity: 0 m/s - 1 0 0 1 256 10 # River velocity +FloodQ: 0 m3/s - 1 0 0 1 256 10 # Flood discharge +FloodEvap: 0 m3 - 1 0 0 1 256 10 # Flood evaporation +FloodStor: 0 m3 - 1 0 0 1 256 10 # Flood storage +FloodDepth: 0 m - 1 0 0 1 256 10 # Flood depth +FloodVelocity: 0 m/s - 1 0 0 1 256 10 # Flood velocity +FloodedFrac: 0 - - 1 0 0 1 256 10 # Flooded fraction +FloodedArea: 0 m2 - 1 0 0 1 256 10 # Flooded area +SurfElev: 0 m - 1 0 0 1 256 10 # Surface elevation +RunoffStor: 0 m3 - 1 0 0 1 256 10 # Runoff storage +BaseflowStor: 0 m3 - 1 0 0 1 256 10 # Baseflow storage + +#Irrigation +Irrigated water: 1 kg/m2s - 1 0 0 1 256 10 # Irrigation amount + +#NoahMP +SnowIce: 0 mm - 1 0 0 3 85 1000 # snow ice contents in layers +SnowLiq: 0 mm - 1 0 0 3 85 1000 # snow water contents in layers +VegCanopT: 0 K - 1 0 0 1 146 10 # Vegetation canopy temperature (K) +CanopVP: 0 Pa - 1 0 0 1 146 10 # Canopy air vapor pressure +CanopWetFrac: 0 - - 1 0 0 1 146 10 # Canopy wet fraction +CanopIntLiq: 0 mm - 1 0 0 1 146 10 # Canopy intercepted liquid water +SnowAge: 0 - - 1 0 0 1 146 10 # snow age factor +AvgGrndT: 0 K - 1 0 0 1 148 10 # Average surface temperature (K) +ActSnowNL: 0 - - 1 0 0 1 148 10 # actual number of snow layers +LeafMass: 0 g/m2 - 1 0 0 1 148 10 # leaf mass +RootMass: 0 g/m2 - 1 0 0 1 148 10 # stem mass +StemMass: 0 g/m2 - 1 0 0 1 148 10 # wood mass +WoodMass: 0 g/m2 - 1 0 0 1 148 10 # mass of wood including woody roots [g/m2] +DeepSoilCarbon: 0 g/m2 - 1 0 0 1 148 10 # stable carbon in deep soil [g/m2] +ShallowSoilCarbon: 0 g/m2 - 1 0 0 1 148 10 # short-lived carbon in shallow soil [g/m2] +RechToGW: 0 m - 1 0 0 1 148 10 # recharge to water table when groundwater is deep +RechFromGW: 0 m - 1 0 0 1 148 10 # recharge from water table when groundwater shallow +SwReflect: 0 W/m2 UP 1 0 0 1 111 10 # Net shortwave radiation (W/m2) +VegT2m: 0 K - 1 0 0 1 146 10 # Vegetation temperature at 2m height (K) +QairT2m: 0 kg/kg - 1 0 0 1 51 1000 # 2-m over vegetation specific humidity +APAR: 0 W/m2 IN 1 0 0 1 111 10 # absorbed PAR energy by canopy (W/m2) +SAV: 0 W/m2 IN 1 0 0 1 111 10 # solar radiation absorbed by vegetation (W/m2) +SAG: 0 W/m2 IN 1 0 0 1 111 10 # solar radiation absorbed by ground (W/m2) +PSCO2: 0 umol/m2s IN 1 0 0 1 111 10 # total photosynthesis of CO2 (umol/m2s) +RsSunlit: 0 s/m - 0 0 0 1 238 100 # sunlit stomatal resistance +RsShaded: 0 s/m - 0 0 0 1 238 100 # shaded stomatal resistance +BCanoGapFrac: 0 - - 0 0 0 1 238 100 # between-canopy gap fraction for beam +WCanoGapFrac: 0 - - 0 0 0 1 238 100 # within-canopy gap fraction for beam +ChVeg: 0 s/m - 1 0 0 1 208 10 # Surface heat exchange coeff. over vegetated fraction +ChBare: 0 s/m - 1 0 0 1 208 10 # Surface heat exchange coeff. over bare soil fraction +QhCano: 0 W/m2 UP 1 0 0 1 122 10 # canopy sensible heat flux (W/m2) +QhBare: 0 W/m2 UP 1 0 0 1 122 10 # bare soil sensible heat flux (W/m2) +EvapHBare: 0 W/m2 UP 1 0 0 1 122 10 # bare soil evaporation heat (W/m2) +EvapHGrnd: 0 W/m2 UP 1 0 0 1 122 10 # ground evaporation heat (W/m2) +GrndHBare: 0 W/m2 UP 1 0 0 1 122 10 # bare ground heat flux (W/m2) +GrndHVeg: 0 W/m2 UP 1 0 0 1 122 10 # vegetated ground heat flux (W/m2) +ChLeaf: 0 s/m - 1 0 0 1 208 10 # Surface exchange coefficient for heat over leaf +ChV2: 0 s/m - 1 0 0 1 208 10 # sensible heat exchange coef. over vegetated fraction +ChB2: 0 s/m - 1 0 0 1 208 10 # sensible heat exchange coeff. over bare ground +fpice: 0 - - 1 0 0 1 208 10 # snow fraction in precipitation + diff --git a/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/README b/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/README new file mode 100644 index 000000000..304b65a33 --- /dev/null +++ b/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/README @@ -0,0 +1,25 @@ +Noah-MP.4.0.1 SMAP Assimilation with Relocated CDF + +This case demonstrates EnKF-based assimilation of NASA SMAP data into +Noah-MP.4.0.1 LSM using relocated CDFs + +It uses: +* MERRA2 +* Retrospective run mode +* Noah-MP.4.0.1 LSM +* Regional domain over Europe +* Time period from 1 June 2020 to 1 Jan 2022. + +This directory contains: +* This README file +* The lis.config file used for this test case +* The NOAHMP36_OUTPUT_LIST.TBL.dyn.full file used by lis.config to select + the output variables. + +Note that lis.config should be edited to make sure the locations of the +parameter and forcing files are specified correctly. + +To run this testcase: +* Compile LIS +* Run the LIS executable using the lis.config file and the testcase input + data. diff --git a/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/lis.config b/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/lis.config new file mode 100755 index 000000000..ea1a74cb1 --- /dev/null +++ b/lis/testcases/dataassim/SMAP_sm_pclimo_transfer/lis.config @@ -0,0 +1,260 @@ +#Overall driver options +Running mode: retrospective +Map projection of the LIS domain: "lambert" # latlon # MN: why latlon should be lambert +Number of nests: 1 +Number of surface model types: 1 +Surface model types: LSM +Surface model output interval: 6hr +Land surface model: "Noah-MP.4.0.1" +Number of met forcing sources: 1 +Blending method for forcings: overlay +Met forcing sources: "MERRA2" # "COAMPSout" +Met forcing chosen ensemble member: 1 +Topographic correction method (met forcing): "none" +Enable spatial downscaling of precipitation: 0 +Spatial upscaling method (met forcing): average +Spatial interpolation method (met forcing): bilinear +Temporal interpolation method (met forcing): linear + +#Runtime options +Forcing variables list file: ./forcing_variables.txt +Output methodology: "2d gridspace" +Output model restart files: 1 +Output data format: netcdf +Output naming style: "3 level hierarchy" +Start mode: restart #coldstart +Starting year: 2020 #2015 ##2000 +Starting month: 05 #03 ##01 +Starting day: 31 #01 ##01 +Starting hour: 23 #0 +Starting minute: 45 #0 +Starting second: 0 +Ending year: 2022 #2015 +Ending month: 01 #03 +Ending day: 01 #01 +Ending hour: 0 +Ending minute: 0 +Ending second: 0 +Undefined value: -9999 +Output directory: "output.da.cdf.transfer.diff.forcing" +Diagnostic output file: "output.da.cdf.transfer.diff.forcing/logs/lislog" +Number of ensembles per tile: 12 + +#The following options are used for subgrid tiling based on vegetation +Maximum number of surface type tiles per grid: 1 +Minimum cutoff percentage (surface type tiles): 0.05 +Maximum number of soil texture tiles per grid: 1 +Minimum cutoff percentage (soil texture tiles): 0.05 +Maximum number of soil fraction tiles per grid: 1 +Minimum cutoff percentage (soil fraction tiles): 0.05 +Maximum number of elevation bands per grid: 1 +Minimum cutoff percentage (elevation bands): 0.05 +Maximum number of slope bands per grid: 1 +Minimum cutoff percentage (slope bands): 0.05 +Maximum number of aspect bands per grid: 1 +Minimum cutoff percentage (aspect bands): 0.05 + +#Processor layout +#Should match the total number of processors used +Number of processors along x: 256 #28 # 16 +Number of processors along y: 1 #16 +Halo size along x: 0 +Halo size along y: 0 + +#Sub-models +Routing model: "none" +Radiative transfer model: none +Number of application models: 0 + +#HYMAP router +HYMAP routing model time step: 15mn +HYMAP routing model output interval: 1da +HYMAP routing model restart interval: 1mo +HYMAP run in ensemble mode: 0 +# method: enter 1 - kinematic; 2 - diffusive +# linear reservoir flag: enter 1 - use; or 2 - do not use linear reservoirs +# evaporation option: enter 1 - compute; or 2 - do not compute evapotation in floodplains +HYMAP routing method: kinematic +HYMAP routing model linear reservoir flag: 1 +HYMAP routing model evaporation option: 2 +HYMAP routing model restart file: LIS_RST_HYMAP_router_201701312345.d01.bin +HYMAP routing model start mode: restart +HYMAP routing LIS output directory: HYMAPTEST + +#---------------------DATA ASSIMILATION ---------------------------------- +#Data Assimilation Options +Number of data assimilation instances: 1 + +Data assimilation algorithm: "EnKF" +Data assimilation set: "SMAP(NASA) soil moisture with CDF Transfer" # "SMAP(NASA) soil moisture" +Number of state variables: 4 +Data assimilation use a trained forward model: 0 +Data assimilation trained forward model output file: none +Data assimilation exclude analysis increments: 0 +Data assimilation output interval for diagnostics: "1da" +Data assimilation number of observation types: 1 +Data assimilation output ensemble spread: 1 +Data assimilation output processed observations: 1 +Data assimilation output innovations: 1 + +Data assimilation scaling strategy: "CDF matching" +Data assimilation observation domain file: ../procLSM_test2/lis_input.noahmp401.EU.nc + +Bias estimation algorithm: "none" +Bias estimation attributes file: "none" +Bias estimation restart output frequency: +Bias estimation start mode: +Bias estimation restart file: + +#Perturbation options +Perturbations start mode: "restart" #"coldstart" +Perturbations restart output interval: "1mo" +#Perturbations restart filename: output.da.cdf.transfer.diff.forcing/DAPERT/202005/LIS_DAPERT_202005312345.d01.bin # "none" +Perturbations restart filename: rst/LIS_DAPERT_202005312345.d01.bin # "none" + +Apply perturbation bias correction: 1 + +Forcing perturbation algorithm: "GMAO scheme" +Forcing perturbation frequency: "1hr" +Forcing attributes file: ../../attribs/forcing_attribs.txt +Forcing perturbation attributes file: ../../attribs/forcing_pertattribs.txt + +State perturbation algorithm: "GMAO scheme" +State perturbation frequency: "6hr" +State attributes file: ../../attribs/noah_sm_attribs.txt +State perturbation attributes file: ../../attribs/noah_sm_pertattribs.txt + +Observation perturbation algorithm: "GMAO scheme" +Observation perturbation frequency: "6hr" +Observation attributes file: ../../attribs/smap_attribs.txt +Observation perturbation attributes file: ../../attribs/smap_pertattribs.txt + +# SMAP DATA ENTRIES: +SMAP(NASA) soil moisture data designation: SPL3SMP +SMAP(NASA) soil moisture data directory: ../input/RS_DATA/SMAP/SPL3SMP.007 +SMAP(NASA) soil moisture use scaled standard deviation model: 0 +SMAP(NASA) soil moisture apply SMAP QC flags: 1 +SMAP(NASA) model CDF file: # ../DA_proc_LSM/OUTPUT.strat.p.climo.15bins/cdf_noahmp401.nc # ../DA_proc_LSM/OUTPUT/cdf_noahmp401.nc +SMAP(NASA) observation CDF file: # ../DA_proc_SMAP/OUTPUT.strat.p.climo.15bins/cdf_smapobs.nc #../DA_proc_SMAP/OUTPUT/cdf_smapobs.nc +SMAP(NASA) soil moisture number of bins in the CDF: 100 +SMAP(NASA) soil moisture use scaled standard deviation model: 0 +SMAP(NASA) CDF read option: 0 +SMAP(NASA) soil moisture Composite Release ID: "R17" + + +#----------------------- Soil moisture CDF transfer ------------- +Use CDF transfer for soil moisture data assimilation: 1 # defualt =0 +Reference domain model CDF file: ref_model_cdf_forcing_diff/stratified_cdf_noahmp401.nc # ref_model_cdf/stratified_cdf_noahmp401.nc +Reference domain obs CDF file: ref_obs_cdf_forcing_diff/stratified_cdf_smapobs.nc # ref_obs_cdf/stratified_cdf_smapobs.nc +#Number of bins in the soil moisture CDF: 100 +Reference domain precipitation climatology data source: Precip.climo.us.nldas2/LVT_MEAN_FINAL.202201010000.d01.nc #Precip.climo.us.merra/LVT_MEAN_FINAL.202201010000.d01.nc +Target domain precipitation climatology data source: Precip.climo.eu.merra/LVT_MEAN_FINAL.202201010000.d01.nc +#---------------------------------------------------------------- + +GLASS LAI data directory: ../GLASSLAIdata + +Irrigation scheme: "none" +Irrigation output interval: "1da" +Irrigation threshold: 0.5 +Irrigation GVF parameter 1: 0.0 +Irrigation GVF parameter 2: 0.4 + +# UMD-landcover and "CROPMAP" crop types (Leff et al. 2004): +Sprinkler irrigation max root depth file: ../input/LS_PARAMETERS/irrigation/conus_modis/maxrootdepth32.txt +#------------------------DOMAIN SPECIFICATION-------------------------- + +#The following options list the choice of parameter maps to be used +LIS domain and parameter data file: ../procLSM_test2/lis_input.noahmp401.EU.nc +Landmask data source: LDT +Landcover data source: LDT +Soil texture data source: LDT +Soil fraction data source: none +Soil color data source: none +Elevation data source: LDT +Slope data source: LDT +Aspect data source: LDT +Curvature data source: none +LAI data source: none +SAI data source: none +Albedo data source: LDT +Max snow albedo data source: LDT +Greenness data source: LDT +Roughness data source: none +Porosity data source: none +Ksat data source: none +B parameter data source: none +Quartz data source: none +Emissivity data source: none + +#--------------------------------FORCINGS---------------------------------- + +MERRA2 forcing directory: ../input/MET_FORCING/MERRA2 +MERRA2 use lowest model level forcing: 1 +MERRA2 use 2m wind fields: 0 +MERRA2 use corrected total precipitation: 1 + +IMERG forcing directory: ../input/MET_FORCING/IMERG +IMERG product: 'final' +IMERG version: V06B + +COAMPS output forcing directory: ../../COAMPS/EU/COAMPS_forcing # ../../COAMPS/COAMPS_forcing # input/MET_FORCING/COAMPS_OUT/bak +COAMPS nest id: 1 # 1 # 1: 45 km ; 2: 15km + +#-----------------------LAND SURFACE MODELS-------------------------- + + +Noah-MP.4.0.1 model timestep: 15mn +Noah-MP.4.0.1 restart output interval: 1mo +#Noah-MP.4.0.1 restart file: output.da.cdf.transfer.diff.forcing/SURFACEMODEL/202005/LIS_RST_NOAHMP401_202005312345.d01.nc +Noah-MP.4.0.1 restart file: rst/LIS_RST_NOAHMP401_202005312345.d01.nc + +#Noah-MP.4.0.1 restart file: ../DA_ensrst/LIS_EnRST_NOAH36_201502282345.d01.nc +Noah-MP.4.0.1 restart file format: netcdf +Noah-MP.4.0.1 soil parameter table: ../input/LS_PARAMETERS/noahmp401_parms/SOILPARM.TBL +Noah-MP.4.0.1 general parameter table: ../input/LS_PARAMETERS/noahmp401_parms/GENPARM.TBL +Noah-MP.4.0.1 MP parameter table: ../input/LS_PARAMETERS/noahmp401_parms/MPTABLE.TBL +Noah-MP.4.0.1 number of soil layers: 4 +Noah-MP.4.0.1 thickness of soil layers: 0.1 0.3 0.6 1.0 +Noah-MP.4.0.1 dynamic vegetation option: 4 # Up to 10 different options +Noah-MP.4.0.1 canopy stomatal resistance option: 1 # 1=Ball-Berry; 2=Jarvis +Noah-MP.4.0.1 soil moisture factor for stomatal resistance: 1 # 1=Noah; 2=CLM; 3=SSiB +Noah-MP.4.0.1 runoff and groundwater option: 1 # 1=SIMGM; 2=SIMTOP; 3=Schaake96; 4=BATS; 5=Miguez-Macho&Fan +Noah-MP.4.0.1 surface layer drag coefficient option: 1 # 1=M-O; 2=Chen97 +Noah-MP.4.0.1 supercooled liquid water option: 1 # 1=NY06; 2=Koren99 +Noah-MP.4.0.1 frozen soil permeability option: 1 # 1=NY06; 2=Koren99 +Noah-MP.4.0.1 radiation transfer option: 3 # 1=gap=F(3D;cosz); 2=gap=0; 3=gap=1-Fveg +Noah-MP.4.0.1 snow surface albedo option: 2 # 1=BATS; 2=CLASS +Noah-MP.4.0.1 rainfall & snowfall option: 1 # 1=Jordan91; 2=BATS; 3=Noah +Noah-MP.4.0.1 lower boundary of soil temperature option: 2 # 1=zero-flux; 2=Noah +Noah-MP.4.0.1 snow&soil temperature time scheme option: 1 # 1=semi-implicit; 2=fully implicit; 3=FSNO for TS +Noah-MP.4.0.1 glacier option: 1 # 1=include phase change; 2=slab ice (Noah) +Noah-MP.4.0.1 surface resistance option: 1 # 1=Sakaguchi and Zeng 2009; 2=Sellers (1992); 3=adjusted Sellers; 4=option1 for non-snow and rsurf_snow for snow +Noah-MP.4.0.1 soil configuration option: 1 # 1=input dominant soil texture; 2=input soil texture varies that varies with depth; 3=soil composition and pedotransfer; 4=input soil properties +Noah-MP.4.0.1 soil pedotransfer function option: 1 # 1=Saxton and Rawls (2006) (used when soil_opt=3) +Noah-MP.4.0.1 crop model option: 0 # 0=No crop model; 1=Liu et al. 2016; 2=Gecros +Noah-MP.4.0.1 urban physics option: 0 # 0=No; 1=Single-layer; 2=Multi-layer BEP scheme; 3=Multi-layer BEM scheme +Noah-MP.4.0.1 reference height of temperature and humidity: 10.0 +Noah-MP.4.0.1 initial surface skin temperature: 288.0 +Noah-MP.4.0.1 initial snow water equivalent: 0.0 +Noah-MP.4.0.1 initial snow depth: 0.0 +Noah-MP.4.0.1 initial total canopy surface water: 0.0 +Noah-MP.4.0.1 initial soil temperatures: 288.0 288.0 288.0 288.0 +Noah-MP.4.0.1 initial total soil moistures: 0.20 0.20 0.20 0.20 +Noah-MP.4.0.1 initial leaf area index: 0.5 +Noah-MP.4.0.1 initial water table depth: 2.5 +Noah-MP.4.0.1 initial water in the aquifer: 4900.0 +Noah-MP.4.0.1 initial water in aquifer and saturated soil: 4900.0 + + +#---------------------------MODEL OUTPUT CONFIGURATION----------------------- +#Specify the list of ALMA variables that need to be featured in the +#LSM model output +Output start year: +Output start month: +Output start day: +Output start hour: +Output start minutes: +Output start seconds: + +Model output attributes file: './NOAHMP36_OUTPUT_LIST.TBL.dyn.full' +