diff --git a/lis/dataassim/obs/ML_SNWD/ML_SNWD_Mod.F90 b/lis/dataassim/obs/ML_SNWD/ML_SNWD_Mod.F90 new file mode 100755 index 000000000..6d65718f9 --- /dev/null +++ b/lis/dataassim/obs/ML_SNWD/ML_SNWD_Mod.F90 @@ -0,0 +1,292 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center Land Information System (LIS) v7.2 +! +! Copyright (c) 2015 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: ML_SNWD_Mod +! +! !DESCRIPTION: +! This module contains interfaces and subroutines to +! handle machine leaning ML snow depth retrievals. +! +! !REVISION HISTORY: +! 29 Jul 2024 Devon Dunmire; Initial Specification +! +! This reader can be used for netcdf data in the following format: +! filename: 'SD_yyyymmdd_.nc' with variables 'SD' [m],'lat' and 'lon' + + +module ML_SNWD_Mod +! !USES: + use ESMF +!EOP + implicit none + + PRIVATE + +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: ML_SNWD_setup +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + PUBLIC :: ML_SNWD_struc + + type, public :: ML_SNWD_dec + logical :: startMode + real :: ssdev + integer :: snowfield + integer :: nc,nr + real, allocatable :: snwd(:,:) + real, allocatable :: snwdtime(:,:) + end type ML_SNWD_dec + + type(ML_SNWD_dec), allocatable :: ML_SNWD_struc(:) + +contains +!BOP +! +! !ROUTINE: ML_SNWD_setup +! \label{ML_SNWD_setup} +! +! !INTERFACE: + subroutine ML_SNWD_setup(k, OBS_State, OBS_Pert_State) +! !USES: + + use LIS_coreMod + use LIS_timeMgrMod + use LIS_historyMod + use LIS_perturbMod + use LIS_DAobservationsMod + use LIS_logMod + use netcdf + + 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 ML snow depth data. +! +! The arguments are: +! \begin{description} +! \item[OBS\_State] observation state +! \item[OBS\_Pert\_State] observation perturbations state +! \end{description} +!EOP + integer :: n + integer :: ftn + integer :: i + 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 :: MLsnowobsdir + 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, allocatable :: ssdev(:) + + allocate(ML_SNWD_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,"ML snow depth data directory:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetattribute(LIS_config,MLsnowobsdir,& + rc=status) + call LIS_verify(status,'ML snow depth data directory: not defined') + call ESMF_AttributeSet(OBS_State(n),"Data Directory",& + MLsnowobsdir, rc=status) + call LIS_verify(status) + 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 ML snow depth data specifications' + +!---------------------------------------------------------------------------- +! Create the array containers that will contain the observations and +! the perturbations. Since there is only one layer +! being assimilated, the array size is LIS_rc%ngrid(n). +!---------------------------------------------------------------------------- + + do n=1,LIS_rc%nnest + + write(unit=temp,fmt='(i2.2)') 1 + read(unit=temp,fmt='(2a1)') vid + + obsField(n) = ESMF_FieldCreate(grid=LIS_obsvecGrid(n,k),& + arrayspec=realarrspec,& + 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) + + call ESMF_StateAdd(OBS_State(n),(/obsField(n)/),rc=status) + call LIS_verify(status) + + 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) + + ssdev = obs_pert%ssdev(1) + ML_SNWD_struc(n)%ssdev =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 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) + + call ESMF_StateAdd(OBS_Pert_State(n),(/pertField(n)/),rc=status) + call LIS_verify(status) + endif + + deallocate(vname) + deallocate(varmax) + deallocate(varmin) + deallocate(ssdev) + enddo +!------------------------------------------------------------- +! set up the ML domain %and interpolation weights. +!------------------------------------------------------------- + + do n=1,LIS_rc%nnest + call ESMF_ConfigFindLabel(LIS_config,& + "ML snow depth domain x-dimension size:",rc=status) + call ESMF_ConfigGetAttribute(LIS_config,ML_SNWD_struc(n)%nc,rc=status) + call LIS_verify(status,'[ERR] ML snow depth domain x-dimension size: not defined') + + call ESMF_ConfigFindLabel(LIS_config,& + "ML snow depth domain y-dimension size:",rc=status) + call ESMF_ConfigGetAttribute(LIS_config,ML_SNWD_struc(n)%nr,rc=status) + call LIS_verify(status,'[ERR] ML snow depth domain y-dimension size: not defined') + + allocate(ML_SNWD_struc(n)%snwd(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k))) + allocate(ML_SNWD_struc(n)%snwdtime(& + LIS_rc%obs_lnc(k), LIS_rc%obs_lnr(k))) + ML_SNWD_struc(n)%snwd = LIS_rc%udef + ML_SNWD_struc(n)%snwdtime = -1 + + call LIS_registerAlarm("ML snow depth read alarm",& + 86400.0, 86400.0) + + ML_SNWD_struc(n)%startMode = .true. + enddo + + write(LIS_logunit,*) '[INFO] Created ESMF States to hold ML observations data' + + end subroutine ML_SNWD_setup + +end module ML_SNWD_Mod + diff --git a/lis/dataassim/obs/ML_SNWD/read_ML_SNWD.F90 b/lis/dataassim/obs/ML_SNWD/read_ML_SNWD.F90 new file mode 100755 index 000000000..f0db26f2c --- /dev/null +++ b/lis/dataassim/obs/ML_SNWD/read_ML_SNWD.F90 @@ -0,0 +1,425 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center Land Information System (LIS) v7.2 +! +! Copyright (c) 2015 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +! +! +! This reader can be used for netcdf files in the following format: +! name: 'SD_yyyymmdd_.nc' with variables 'SD' (m), 'lat' and 'lon' +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +#include "LIS_misc.h" +!BOP +! !ROUTINE: read_ML_SNWD +! \label{read_ML_SNWD} +! +! !REVISION HISTORY: +! 29 Aug 2019: Hans Lievens; Initial Specification +! +! !INTERFACE: +subroutine read_ML_SNWD(n,k,OBS_State,OBS_Pert_State) +! !USES: + use ESMF + use LIS_mpiMod + use LIS_coreMod + use LIS_timeMgrMod + use LIS_logMod + use map_utils + use LIS_DAobservationsMod + use LIS_pluginIndices, only : LIS_ML_SNWD_obsId + use ML_SNWD_Mod, only : ML_SNWD_struc + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + type(ESMF_State) :: OBS_State + type(ESMF_State) :: OBS_Pert_State +! +! !DESCRIPTION: +! +! This routine reads Sentinel-1 snow depth observations in netcdf4 format +! The data is read at 0z every day and is kept in memory. At +! each timestep, a subset of data is chosen for use in DA if +! the local time of the grid point is 6AM. +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[k] index of the assimilation instance +! \item[OBS\_State] observations state +! \item[OBS\_Pert\_State] observations perturbation state +! \end{description} +! +!EOP + type(ESMF_Field) :: snowField,pertfield + logical :: alarmCheck + logical :: data_upd, file_exists + logical :: dataflag(LIS_npes) + logical :: dataflag_local + integer :: c,r, p, t + character*100 :: obsdir + character*80 :: ML_filename + integer :: ftn + real :: lon, lhour + real :: gmt + real :: dt + integer :: zone + integer :: grid_index + real :: ssdev(LIS_rc%obs_ngrid(k)) + real, pointer :: obsl(:) + integer :: gid(LIS_rc%obs_ngrid(k)) + integer :: assimflag(LIS_rc%obs_ngrid(k)) + integer :: status, iret, ierr + integer :: fnd + real :: snwd_current(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + integer :: days(12) + data days /31,28,31,30,31,30,31,31,30,31,30,31/ !BZ + + + + call ESMF_AttributeGet(OBS_State,"Data Directory",& + obsdir, rc=status) + call LIS_verify(status,'Error in AttributeGet: Data Directory') + +!------------------------------------------------------------------------- +! Read the data at 0z daily. +!------------------------------------------------------------------------- + alarmCheck = LIS_isAlarmRinging(LIS_rc, "ML snow depth read alarm") + + if(alarmCheck.or.ML_SNWD_struc(n)%startMode) then + ML_SNWD_struc(n)%startMode = .false. + + ML_SNWD_struc(n)%snwd = LIS_rc%udef + ML_SNWD_struc(n)%snwdtime = -1 + + call ML_SNWD_filename(ML_filename,obsdir,& + LIS_rc%yr,LIS_rc%mo,LIS_rc%da) + + inquire(file=ML_filename,exist=file_exists) + if(file_exists) then + + write(LIS_logunit,*) '[INFO] Reading ML SNWD data ',ML_filename + call read_ML_SNWD_data(n,k, ML_filename, ML_SNWD_struc(n)%snwd) + + +!------------------------------------------------------------------------- +! Store the GMT corresponding to 6AM localtime at each grid point +!------------------------------------------------------------------------- + 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) + lon = LIS_obs_domain(n,k)%lon(grid_index) + + lhour = 6.0 + call LIS_localtime2gmt(gmt,lon,lhour,zone) + ML_SNWD_struc(n)%snwdtime(c,r) = gmt + + endif + enddo + enddo + + endif + endif + +!------------------------------------------------------------------------- +! Update the OBS_State +!------------------------------------------------------------------------- + + call ESMF_StateGet(OBS_State,"Observation01",snowfield,& + rc=status) + call LIS_verify(status, 'Error: StateGet Observation01') + + call ESMF_FieldGet(snowfield,localDE=0,farrayPtr=obsl,rc=status) + call LIS_verify(status, 'Error: FieldGet') + + obsl = LIS_rc%udef + +!------------------------------------------------------------------------- +! Update the OBS_State by subsetting to the local grid time +!------------------------------------------------------------------------- + + 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 - ML_SNWD_struc(n)%snwdtime(c,r))*3600.0 + lon = LIS_obs_domain(n,k)%lon(grid_index) + + if(dt.ge.0.and.dt.lt.LIS_rc%ts) then + obsl(LIS_obs_domain(n,k)%gindex(c,r)) = & + ML_SNWD_struc(n)%snwd(c,r) + endif + + endif + enddo + enddo + + dataflag_local = .false. + +!------------------------------------------------------------------------- +! Apply LSM based quality control and screening of observations +!------------------------------------------------------------------------- + + call lsmdaqcobsstate(trim(LIS_rc%lsm)//"+"& + //trim(LIS_ML_SNWD_obsId)//char(0), & + n, k,OBS_state) + + snwd_current = LIS_rc%udef + call LIS_checkForValidObs(n,k,obsl,fnd,snwd_current) + + if(fnd.eq.0) then + dataflag_local = .false. + else + dataflag_local = .true. + endif + +#if (defined SPMD) + call MPI_ALLGATHER(dataflag_local,1, MPI_LOGICAL, dataflag(:),& + 1, MPI_LOGICAL, LIS_mpi_comm, ierr) +#endif + data_upd = .false. + + do p=1,LIS_npes + data_upd = data_upd.or.dataflag(p) + enddo + + 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, 'Error: AttributeSet in Data Update Status') + + + call ESMF_StateGet(OBS_Pert_State,"Observation01",pertfield,& + rc=status) + call LIS_verify(status, 'ESMF_StateGet for Observation01 for OBS_Pert_State failed in read_ML_SNWD') + + if(LIS_rc%obs_ngrid(k).gt.0) then + +!linearly scale the observation err + ssdev = ML_SNWD_struc(n)%ssdev + do t=1,LIS_rc%obs_ngrid(k) + if(obsl(t).ne.-9999.0) then + ssdev(t) = ML_SNWD_struc(n)%ssdev !+ 0.05*obsl(t) + endif + enddo + + call ESMF_AttributeSet(pertField,"Standard Deviation",& + ssdev,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status) + +! To be added (for EnKS test): +! call ESMF_AttributeSet(OBS_State,& +! name="Data averaging factor",& +! value=float(days(LIS_rc%mo)),rc=status) +! call LIS_verify(status) + + call ESMF_AttributeSet(snowfield,"Grid Number",& + gid,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status,'Error: AttributeSet in Grid Number') + + call ESMF_AttributeSet(snowfield,"Assimilation Flag",& + assimflag,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status, 'Error: AttributeSet in Assimilation Flag') + + endif + else + call ESMF_AttributeSet(OBS_State,"Data Update Status",& + .false., rc=status) + call LIS_verify(status, "Error: AttributeSet Data Update Status") + return + end if + + +end subroutine read_ML_SNWD + + + +!BOP +! +! !ROUTINE: read_ML_SNWD_data +! \label{read_ML_SNWD_data} +! +! !INTERFACE: +subroutine read_ML_SNWD_data(n, k, fname, snwd_ip) +! +! !USES: +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + use LIS_coreMod + use LIS_logMod + use map_utils, only : latlon_to_ij + use ML_SNWD_Mod, only : ML_SNWD_struc + + implicit none +! +! !INPUT PARAMETERS: +! + integer :: n + integer :: k + character (len=*) :: fname + +! !OUTPUT PARAMETERS: +! +! !DESCRIPTION: +! This subroutine reads the ML NETCDF files +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[fname] name of the ML SNWD file +! \item[SDobs\_ip] snow depth data processed to the LIS domain +! \end{description} +! +! !FILES USED: +! +! !REVISION HISTORY: +! +!EOP + real :: SD(ML_SNWD_struc(n)%nr,ML_SNWD_struc(n)%nc) + real :: lat_nc(ML_SNWD_struc(n)%nr) + real :: lon_nc(ML_SNWD_struc(n)%nc) + real :: snwd_ip(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + integer :: nsnwd_ip(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + logical :: file_exists + integer :: c,r,i,j + integer :: stn_col,stn_row + real :: col,row + integer :: nid + integer :: SDId,latId,lonId + integer :: ios + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + +! SD = LIS_rc%udef +! lat_nc = LIS_rc%udef +! lon_nc = LIS_rc%udef +! snwd_ip = LIS_rc%udef + + inquire(file=fname, exist=file_exists) + if(file_exists) then + write(LIS_logunit,*) 'Reading ',trim(fname) + ios = nf90_open(path=trim(fname),mode=NF90_NOWRITE,ncid=nid) + call LIS_verify(ios,'Error opening file '//trim(fname)) + + ! variables + ios = nf90_inq_varid(nid, 'SD',SDid) + call LIS_verify(ios, 'Error nf90_inq_varid: snow depth data') + + ios = nf90_inq_varid(nid, 'lat',latid) + call LIS_verify(ios, 'Error nf90_inq_varid: latitude data') + + ios = nf90_inq_varid(nid, 'lon',lonid) + call LIS_verify(ios, 'Error nf90_inq_varid: longitude data') + + !values + ios = nf90_get_var(nid, SDid, SD) + call LIS_verify(ios, 'Error nf90_get_var: SD') + + ios = nf90_get_var(nid, latid, lat_nc) + call LIS_verify(ios, 'Error nf90_get_var: lat') + + ios = nf90_get_var(nid, lonid, lon_nc) + call LIS_verify(ios, 'Error nf90_get_var: lon') + + ! close file + ios = nf90_close(ncid=nid) + call LIS_verify(ios,'Error closing file '//trim(fname)) + + + + snwd_ip = 0 + nsnwd_ip = 0 + + ! Map the data by averaging + do i=1,ML_SNWD_struc(n)%nr + do j=1,ML_SNWD_struc(n)%nc + + call latlon_to_ij(LIS_domain(n)%lisproj,& + lat_nc(i),lon_nc(j),col,row) + stn_col = nint(col) + stn_row = nint(row) + + if(SD(i,j).ge.-999.and.& + stn_col.gt.0.and.stn_col.le.LIS_rc%obs_lnc(k).and.& + stn_row.gt.0.and.stn_row.le.LIS_rc%obs_lnr(k)) then + snwd_ip(stn_col,stn_row) = snwd_ip(stn_col,stn_row) + SD(i,j) + nsnwd_ip(stn_col,stn_row) = nsnwd_ip(stn_col,stn_row) + 1 + endif + enddo + enddo + + do r=1,LIS_rc%obs_lnr(k) + do c=1,LIS_rc%obs_lnc(k) + if(nsnwd_ip(c,r).ne.0) then + !Sentinel-1 SD bias correction factor 1.0/0.85 if wanted + ! is implemented here, now disabled + snwd_ip(c,r) = 1.0*snwd_ip(c,r)/nsnwd_ip(c,r) ! NoahMP4 snow DA wants it in meter + else + snwd_ip(c,r) = LIS_rc%udef + endif + enddo + enddo + + + endif + +#endif + +end subroutine read_ML_SNWD_data + + + +!BOP +! +! !ROUTINE: ML_SNWD_filename +! \label{ML_SWND_filename} +! +! !INTERFACE: +subroutine ML_SNWD_filename(filename, ndir, yr, mo, da) + + implicit none +! !ARGUMENTS: + character*80 :: filename + integer :: yr, mo, da + character (len=*) :: ndir +! +! !DESCRIPTION: +! This subroutine creates a timestamped ML filename +! +! The arguments are: +! \begin{description} +! \item[name] name of the ML filename +! \item[ndir] name of the ML root directory +! \item[yr] current year +! \item[mo] current month +! \item[da] current day +! \end{description} +!EOP + + character (len=4) :: fyr + character (len=2) :: fmo,fda + + write(unit=fyr, fmt='(i4.4)') yr + write(unit=fmo, fmt='(i2.2)') mo + write(unit=fda, fmt='(i2.2)') da + + filename = trim(ndir)//'/SD_'//trim(fyr)//trim(fmo)//trim(fda)//'_.nc' + +end subroutine ML_SNWD_filename + + + diff --git a/lis/dataassim/obs/ML_SNWD/write_ML_SNWDobs.F90 b/lis/dataassim/obs/ML_SNWD/write_ML_SNWDobs.F90 new file mode 100755 index 000000000..0419ac514 --- /dev/null +++ b/lis/dataassim/obs/ML_SNWD/write_ML_SNWDobs.F90 @@ -0,0 +1,130 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center Land Information System (LIS) v7.2 +! +! Copyright (c) 2015 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_ML_SNWDobs +! \label{write_ML_SNWDobs} +! +! !REVISION HISTORY: +! 30 Aug 2019: Hans Lievens; Initial Specification +! +! !INTERFACE: +subroutine write_ML_SNWDobs(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) +! ML observations to a file +! +!EOP + type(ESMF_Field) :: snowField + logical :: data_update + real, pointer :: snowobs(:) + character*100 :: obsname + integer :: ftn + integer :: status + +! integer :: c,r +! real :: testdata(LIS_rc%obs_lnc(k), LIS_rc%obs_lnr(k)) + + 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",snowField, & + rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(snowField, localDE=0, farrayPtr=snowobs, rc=status) + call LIS_verify(status) + + if(LIS_masterproc) then + ftn = LIS_getNextUnitNumber() + call ML_SNWD_obsname(n,k,obsname) + + call LIS_create_output_directory('DAOBS') + open(ftn,file=trim(obsname), form='unformatted') + endif + +#if 0 + testdata = 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 + testdata(c,r) = snowobs(LIS_obs_domain(n,k)%gindex(c,r)) + endif + enddo + enddo + if(LIS_localPet.eq.241) then + open(100,file='test_out.bin',form='unformatted') + write(100) testdata + close(100) + stop + endif +#endif + call LIS_writevar_gridded_obs(ftn,n,k,snowobs) + + if(LIS_masterproc) then + call LIS_releaseUnitNumber(ftn) + endif + + endif + +end subroutine write_ML_SNWDobs + +!BOP +! !ROUTINE: ML_SNWD_obsname +! \label{ML_SNWD_obsname} +! +! !INTERFACE: +subroutine ML_SNWD_obsname(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=cda, fmt='(a2,i2.2)') '.a',k + write(unit=cdate, fmt='(a2,i2.2)') '.d',n + + 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 + + obsname = trim(LIS_rc%odir)//'/DAOBS/'//cdate1(1:6)//& + '/LISDAOBS_'//cdate1// & + trim(cda)//trim(cdate)//'.1gs4r' + +end subroutine ML_SNWD_obsname diff --git a/lis/plugins/LIS_DAobs_pluginMod.F90 b/lis/plugins/LIS_DAobs_pluginMod.F90 old mode 100755 new mode 100644 index 2c32f31ab..bf025c615 --- a/lis/plugins/LIS_DAobs_pluginMod.F90 +++ b/lis/plugins/LIS_DAobs_pluginMod.F90 @@ -222,6 +222,11 @@ subroutine LIS_DAobs_plugin use S1_SNWD_Mod, only : S1_SNWD_setup #endif +! Devon Dunmire added ML snow depth DA +#if ( defined DA_OBS_ML_SNWD ) + use ML_SNWD_Mod, only : ML_SNWD_setup +#endif + #if ( defined DA_OBS_S1_sigmaVVSM ) use S1_sigmaVVSM_Mod, only : S1_sigmaVVSM_setup #endif @@ -448,6 +453,11 @@ subroutine LIS_DAobs_plugin external read_S1_SNWD, write_S1_SNWDobs #endif +! Devon Dunmire added ML snow depth obs +#if ( defined DA_OBS_ML_SNWD ) + external read_ML_SNWD, write_ML_SNWDobs +#endif + #if ( defined DA_OBS_S1_sigmaVVSM) external read_S1_sigmaVVSM, write_S1_sigmaVVSMobs #endif @@ -755,6 +765,19 @@ subroutine LIS_DAobs_plugin write_S1_SNWDobs) #endif +! Devon Dunmire added ML snow depth case +#if ( defined DA_OBS_ML_SNWD ) +!ML SNWD snow obs + call registerdaobsclass(trim(LIS_ML_SNWD_obsId),"LSM") + call registerdaobssetup(trim(LIS_ML_SNWD_obsId)//char(0), & + ML_SNWD_setup) + call registerreaddaobs(trim(LIS_ML_SNWD_obsId)//char(0), & + read_ML_SNWD) + call registerwritedaobs(trim(LIS_ML_SNWD_obsId)//char(0), & + write_ML_SNWDobs) +#endif + + #if ( defined DA_OBS_S1_sigmaVVSM ) !S1 backscatter obs VVSM call registerdaobsclass(trim(LIS_S1_sigmaVVSM_obsId),"LSM") diff --git a/lis/plugins/LIS_lsmda_pluginMod.F90 b/lis/plugins/LIS_lsmda_pluginMod.F90 old mode 100755 new mode 100644 index 41826edd9..f0ec39a7a --- a/lis/plugins/LIS_lsmda_pluginMod.F90 +++ b/lis/plugins/LIS_lsmda_pluginMod.F90 @@ -3340,6 +3340,28 @@ subroutine LIS_lsmda_plugin trim(LIS_S1_SNWD_obsId)//char(0),noahmp401_updatesnowvars) call registerlsmdaqcobsstate(trim(LIS_noahmp401Id)//"+"//& trim(LIS_S1_SNWD_obsId)//char(0),noahmp401_qc_snowobs) + +! Devon Dunmire added ML snow depth + call registerlsmdainit(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_dasnow_init) + call registerlsmdagetstatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_getsnowvars) + call registerlsmdasetstatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_setsnowsimplevars) + call registerlsmdagetobspred(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_getsnowpred) + call registerlsmdaqcstate(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_qcsnow) + call registerlsmdaqcobsstate(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_qc_snowobs) + call registerlsmdascalestatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_scale_snow) + call registerlsmdadescalestatevar(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_descale_snow) + call registerlsmdaupdatestate(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_updatesnowvars) + call registerlsmdaqcobsstate(trim(LIS_noahmp401Id)//"+"//& + trim(LIS_ML_SNWD_obsId)//char(0),noahmp401_qc_snowobs) #if ( defined DA_OBS_GCOMW_AMSR2L3SND ) call registerlsmdainit(trim(LIS_noahmp401Id)//"+"//& diff --git a/lis/plugins/LIS_pluginIndices.F90 b/lis/plugins/LIS_pluginIndices.F90 old mode 100755 new mode 100644 index 13e223d8f..0c0932801 --- a/lis/plugins/LIS_pluginIndices.F90 +++ b/lis/plugins/LIS_pluginIndices.F90 @@ -246,6 +246,7 @@ module LIS_pluginIndices character*50, public, parameter :: LIS_SMMRSNWDsnowobsId = "SMMR snow depth" ! Hans Lievens added S1 character*50, public, parameter :: LIS_S1_SNWD_obsId = "S1 snow depth" + character*50, public, parameter :: LIS_ML_SNWD_obsId = "ML snow depth" character*50, public, parameter :: LIS_S1_sigmaVVSM_obsId = "S1 backscatter VVSM" character*50, public, parameter :: LIS_S1_sigmaVVSMLAI_obsId = "S1 backscatter VVSMLAI" character*50, public, parameter :: LIS_S1_sigmaVHSM_obsId = "S1 backscatter VHSM"