diff --git a/docs/LIS_users_guide/user_cfg_table.adoc b/docs/LIS_users_guide/user_cfg_table.adoc index e6a0acb33..fe4e73a02 100644 --- a/docs/LIS_users_guide/user_cfg_table.adoc +++ b/docs/LIS_users_guide/user_cfg_table.adoc @@ -97,6 +97,7 @@ | CRTM2EM | On | CMEM | On | Tau Omega | On +| WCM | On |==== .Applications @@ -149,6 +150,8 @@ | DA OBS THYSM | On | DA OBS LPRM_AMSREsm | On | DA OBS SMMR_SNWD | On +| DA OBS S1_sigmaVVSM | On +| DA_OBS_S1_sigmaVVVHSMLAI| On | DA OBS SSMI_SNWD | On | DA OBS ANSA_SNWD | On | DA OBS GCOMW_AMSR2L3SND | On diff --git a/lis/configs/lis.config.adoc b/lis/configs/lis.config.adoc index 82382f346..1b90412c7 100644 --- a/lis/configs/lis.config.adoc +++ b/lis/configs/lis.config.adoc @@ -1229,6 +1229,8 @@ Acceptable values are: |"`ANSA snow depth`" | ANSA snow depth |"`SMMR snow depth`" | SMMR snow depth |"`SMMI snow depth`" | SMMI snow depth +|"`S1 backscatter VVSM`" | S1 backscatter VV for SM +|"`S1 backscatter VVVHSMLAI`" | S1 backscatter VV and VH joint assim for SM and LAI |"`AMSR-E SWE`" | AMSR-E SWE |"`PMW snow`" | PMW-based SWE or snow depth |"`MODIS SCF`" | MODIS SCF @@ -2315,6 +2317,18 @@ SMMR snow depth data directory: ./SMMR .... +[[sssec_s1backscatterda,S1 backscatter Assimilation]] +==== S1 backscatter Assimilation + +`S1 backscatter data directory:` specifies the directory for the +S1 backscatter data. + +.Example _lis.config_ entry +.... +S1 backscatter data directory: ./S1 +.... + + [[sssec_ssmisnowdepthda,SSMI snow depth Assimilation]] ==== SSMI snow depth Assimilation @@ -3311,6 +3325,7 @@ endif::devonly[] |CRTM2EM | CRTM2EM |CMEM | CMEM |"`Tau Omega`" | "`Tau Omega`" +|"WCM" | World Cloud Model. Note: "LIS_CRTM" and "LIS-CMEM" must be enabled. |==== `RTM invocation frequency:` specifies the invocation frequency diff --git a/lis/core/LIS_DAobservationsMod.F90 b/lis/core/LIS_DAobservationsMod.F90 index 49b61cde5..b0552f730 100644 --- a/lis/core/LIS_DAobservationsMod.F90 +++ b/lis/core/LIS_DAobservationsMod.F90 @@ -36,6 +36,7 @@ module LIS_DAobservationsMod ! !REVISION HISTORY: ! ! 21 Jun 2006: Sujay Kumar; Initial implementation +! 10 Mar 2022: Michel Bechtold (MB); Bug Fix for multiple obstypes ! use ESMF use LIS_coreMod @@ -324,7 +325,7 @@ subroutine create_obsdomain_objects() integer, allocatable :: deblklist(:,:,:) type(ESMF_DistGrid) :: obsgridDG(LIS_rc%ndas), gridEnsDG(LIS_rc%ndas) integer :: stid, enid - integer :: n, i, k + integer :: n, i, k, m integer :: status, ierr call readObsDomainInput() @@ -337,8 +338,9 @@ subroutine create_obsdomain_objects() do n=1,LIS_rc%nnest do k=1,LIS_rc%ndas -! odeltas(k) = LIS_rc%obs_ngrid(k)*LIS_rc%nobtypes(k) - odeltas(k) = LIS_rc%obs_ngrid(k) + !MB: bug fix multiple obstypes + odeltas(k) = LIS_rc%obs_ngrid(k)*LIS_rc%nobtypes(k) +! odeltas(k) = LIS_rc%obs_ngrid(k) enddo !----------------------------------------------------------------------------- ! The grid, tile space sizes of the decomposed domains are gathered @@ -375,13 +377,15 @@ subroutine create_obsdomain_objects() allocate(deblklist(1,2,LIS_npes)) do i=0,LIS_npes-1 stid = LIS_ooffsets(k,i)+1 - enid = stid + LIS_obs_ngrids(k,i)-1 + !MB: bug fix multiple obstypes + enid = stid + LIS_obs_ngrids(k,i)*LIS_rc%nobtypes(k)-1 deblklist(:,1,i+1) = (/stid/) deblklist(:,2,i+1) = (/enid/) enddo + !MB: bug fix multiple obstypes obsgridDG(k) = ESMF_DistGridCreate(minIndex=(/1/), & - maxIndex=(/LIS_rc%obs_glbngrid(k)/),& + maxIndex=(/LIS_rc%obs_glbngrid(k)*LIS_rc%nobtypes(k)/),& deBlockList=deblklist,rc=status) call LIS_verify(status, 'ESMF_DistGridCreate failed: obsgridDG') deallocate(deblklist) @@ -389,13 +393,15 @@ subroutine create_obsdomain_objects() allocate(deblklist(2,2,LIS_npes)) do i=0,LIS_npes-1 stid = LIS_ooffsets(k,i)+1 - enid = stid+LIS_obs_ngrids(k,i)-1 + !MB: bug fix multiple obstypes + enid = stid+LIS_obs_ngrids(k,i)*LIS_rc%nobtypes(k)-1 deblklist(:,1,i+1) = (/stid,1/) deblklist(:,2,i+1) = (/enid,LIS_rc%nensem(n)/) enddo + !MB: bug fix multiple obstypes gridEnsDG(k) = ESMF_DistGridCreate(minIndex=(/1,1/),& - maxIndex=(/LIS_rc%obs_glbngrid(k),LIS_rc%nensem(n)/),& + maxIndex=(/LIS_rc%obs_glbngrid(k)*LIS_rc%nobtypes(k),LIS_rc%nensem(n)/),& deBlockList=deblklist,rc=status) call LIS_verify(status, 'ESMF_DistGridCreate failed: gridEnsDG') @@ -2035,12 +2041,16 @@ subroutine LIS_writevar_innov(ftn, n, k, varid, var) #else gtmp = var #endif - if(LIS_masterproc) then + if(LIS_masterproc) then + !MB: bug fix multiple obstypes do i=1,LIS_rc%nobtypes(k) + count1=1 + (i-1)*SUM(LIS_obs_ngrids(k,0:i-2)) gtmp1 = LIS_rc%udef - count1=1 - do l=1,LIS_npes + count1=count1+MIN0(1,l-1)*& + (LIS_rc%nobtypes(k)-i)*LIS_obs_ngrids(k,l-2) +& + MIN0(1,l-1)*& + (i-1)*LIS_obs_ngrids(k,l-1) do t =1, LIS_obs_ngrids(k,l-1) c = LIS_obs_domain(n,k)%glb_col(l,t) r = LIS_obs_domain(n,k)%glb_row(l,t) @@ -2053,8 +2063,8 @@ subroutine LIS_writevar_innov(ftn, n, k, varid, var) enddo enddo #if ( defined USE_NETCDF3 || defined USE_NETCDF4 ) - ierr = nf90_put_var(ftn,varid,gtmp1,(/1,1,1/),& - (/LIS_rc%obs_gnc(k),LIS_rc%obs_gnr(k),i/)) + ierr = nf90_put_var(ftn,varid,gtmp1,(/1,1,i/),& + (/LIS_rc%obs_gnc(k),LIS_rc%obs_gnr(k),1/)) call LIS_verify(ierr,'nf90_put_var failed in LIS_historyMod') #endif diff --git a/lis/core/LIS_histDataMod.F90 b/lis/core/LIS_histDataMod.F90 index f5e931293..283102083 100644 --- a/lis/core/LIS_histDataMod.F90 +++ b/lis/core/LIS_histDataMod.F90 @@ -317,6 +317,8 @@ module LIS_histDataMod public :: LIS_MOC_RTM_EMISSIVITY public :: LIS_MOC_RTM_TB public :: LIS_MOC_RTM_SM + public :: LIS_MOC_RTM_Sig0VV + public :: LIS_MOC_RTM_Sig0VH ! Irrigation public :: LIS_MOC_IRRIGATEDWATER @@ -822,6 +824,8 @@ module LIS_histDataMod integer :: LIS_MOC_RTM_EMISSIVITY = -9999 integer :: LIS_MOC_RTM_TB = -9999 integer :: LIS_MOC_RTM_SM = -9999 + integer :: LIS_MOC_RTM_Sig0VV = -9999 + integer :: LIS_MOC_RTM_Sig0VH = -9999 integer :: LIS_MOC_IRRIGATEDWATER @@ -4324,6 +4328,28 @@ subroutine LIS_histDataInit(n, ntiles) n,1,ntiles,(/"m3/m3"/),1,(/"-"/),1,1,1) endif + call ESMF_ConfigFindLabel(modelSpecConfig,"RTM Sig0VV:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_rtm_list,& + "RTM_Sig0VV",& + "rtm_Sigma0_VV",& + "rtm Sigma0 VV",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_RTM_COUNT,LIS_MOC_RTM_Sig0VV,& + LIS_histData(n)%head_rtm_list,& + n,1,ntiles,(/"dB"/),1,(/"-"/),1,1,1) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"RTM Sig0VH:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_rtm_list,& + "RTM_Sig0VH",& + "rtm_Sigma0_VH",& + "rtm Sigma0 VH",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_RTM_COUNT,LIS_MOC_RTM_Sig0VH,& + LIS_histData(n)%head_rtm_list,& + n,1,ntiles,(/"dB"/),1,(/"-"/),1,1,1) + endif + call ESMF_ConfigFindLabel(modelSpecConfig,"Irrigated water:",rc=rc) call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_irrig_list, & "IrrigatedWater",& diff --git a/lis/dataassim/algorithm/enkf/enkf_Mod.F90 b/lis/dataassim/algorithm/enkf/enkf_Mod.F90 index 5e025fd96..f0bb57e92 100644 --- a/lis/dataassim/algorithm/enkf/enkf_Mod.F90 +++ b/lis/dataassim/algorithm/enkf/enkf_Mod.F90 @@ -26,7 +26,9 @@ module enkf_Mod ! ! !REVISION HISTORY: ! 27 Feb 2005: Sujay Kumar; Initial Specification -! +! 30 Jun 2022: Michel Bechtold (MB); Enabling multiple observation types +! 30 Jun 2022: Zdenko Heyaert, Michel Bechtold; bug fix for multiple DA instances + ! !USES: use ESMF use enkf_types @@ -736,9 +738,10 @@ subroutine writeInnovationOutput(n,k) call LIS_verify(nf90_def_dim(ftn,& vardimname,LIS_rc%nobtypes(k),dimId(3)),& 'nf90_def_dim failed for ninnov_'//trim(finst)) + !MB add third dimension for multiple obstypes call LIS_verify(nf90_def_var(ftn,varname,& nf90_float,& - dimids = dimID(1:2), varID=ninnov_Id),& + dimids = dimID(1:3), varID=ninnov_Id),& 'nf90_def_var failed for '//trim(varname)//' in enkf_mod') #if(defined USE_NETCDF4) @@ -763,9 +766,10 @@ subroutine writeInnovationOutput(n,k) vardimname,LIS_rc%nobtypes(k),dimId(3)),& 'nf90_def_dim failed for innov_'//trim(finst)) + !MB add third dimension for multiple obstypes call LIS_verify(nf90_def_var(ftn,varname,& nf90_float,& - dimids = dimID(1:2), varID=innov_Id),& + dimids = dimID(1:3), varID=innov_Id),& 'nf90_def_var failed for innov') #if(defined USE_NETCDF4) @@ -791,9 +795,10 @@ subroutine writeInnovationOutput(n,k) vardimname,LIS_rc%nobtypes(k),dimId(3)),& 'nf90_def_dim failed for analysis_residual_'//trim(finst)) + !MB add third dimension for multiple obstypes call LIS_verify(nf90_def_var(ftn,varname,& nf90_float,& - dimids = dimID(1:2), varID=ares_Id),& + dimids = dimID(1:3), varID=ares_Id),& 'nf90_def_var failed for '//trim(varname)//' in enkf_mod') #if(defined USE_NETCDF4) @@ -818,9 +823,10 @@ subroutine writeInnovationOutput(n,k) vardimname,LIS_rc%nobtypes(k),dimId(3)),& 'nf90_def_dim failed for forecast_sigma_'//trim(finst)) + !MB add third dimension for multiple obstypes call LIS_verify(nf90_def_var(ftn,varname,& nf90_float,& - dimids = dimID(1:2), varID=forecast_sigma_Id),& + dimids = dimID(1:3), varID=forecast_sigma_Id),& 'nf90_def_var for forecast_sigma failed in enkf_mod') #if(defined USE_NETCDF4) @@ -985,7 +991,8 @@ subroutine writeEnsembleSpread(n,k) LIS_rc%nstvars(k),state_size, stvar) do v=1,LIS_rc%nstvars(k) - call LIS_writevar_spread(ftn,n,k,ensspread_id(v), & + ! MB: bug fix for multiple DA instances + call LIS_writevar_spread(ftn,n,LIS_rc%lsm_index,ensspread_id(v), & stvar(v,:),v) enddo @@ -1117,7 +1124,8 @@ subroutine writeAnalysisIncr(n,k) endif do v=1,LIS_rc%nstvars(k) - call LIS_writevar_incr(ftn,n,k,incr_id(v), & + ! MB: bug fix for multiple DA instances + call LIS_writevar_incr(ftn,n,LIS_rc%lsm_index,incr_id(v), & enkf_struc(n,k)%anlys_incr(v,:),v) enddo diff --git a/lis/dataassim/obs/S1_sigmaVVSM/S1_sigmaVVSM_Mod.F90 b/lis/dataassim/obs/S1_sigmaVVSM/S1_sigmaVVSM_Mod.F90 new file mode 100644 index 000000000..3af54ae5f --- /dev/null +++ b/lis/dataassim/obs/S1_sigmaVVSM/S1_sigmaVVSM_Mod.F90 @@ -0,0 +1,309 @@ +!-----------------------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: S1_sigmaVVSM_Mod +! +! !DESCRIPTION: +! This module contains interfaces and subroutines to +! handle S1 backscatter retrievals. +! +! !REVISION HISTORY: +! 29 Aug 2019 Hans Lievens; Initial Specification for sigma depth +! 9 Mar 2021 Isis Brangers, Michel Bechtold; Adaptation for backscatter +! 29 Jun 2022: Louise Busschaert; getting domain dims from files +! +module S1_sigmaVVSM_Mod +! !USES: + use ESMF +!EOP + implicit none + + PRIVATE + +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: S1_sigmaVVSM_setup +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + PUBLIC :: S1_sigma_struc + + type, public :: S1_sigma_dec + logical :: startMode + real :: ssdev + integer :: sigmafield + integer :: nc,nr + real, allocatable :: sigma(:,:) + real, allocatable :: sigmatime(:,:) + end type S1_sigma_dec + + type(S1_sigma_dec), allocatable :: S1_sigma_struc(:) + +contains +!BOP +! +! !ROUTINE: S1_sigmaVVSM_setup +! \label{S1_sigmaVVSM_setup} +! +! !INTERFACE: + subroutine S1_sigmaVVSM_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 structures required for S1 backscatter 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 :: S1sigmaobsdir + character*80 :: S1_firstfile + character*80 :: S1_filename + 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(:) + real :: dx, dy + integer :: NX, NY + integer :: ncid + integer :: flist + integer :: ios + character*100 :: infile + character*100 :: xname, yname + + + allocate(S1_sigma_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,"S1 backscatter data directory:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetattribute(LIS_config,S1sigmaobsdir,& + rc=status) + call LIS_verify(status,'S1 backscatter data directory: not defined') + call ESMF_AttributeSet(OBS_State(n),"Data Directory",& + S1sigmaobsdir, 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 S1 backscatter 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) + S1_sigma_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 S1 domain %and interpolation weights. +!------------------------------------------------------------- + + ! Open first file to get nx ny dimensions + call system('ls ./' // trim(S1sigmaobsdir) // ' > ./S1_listfiles.txt') + + open(flist, file=trim('./S1_listfiles.txt'), & + status='old', iostat=status) + + read(flist, '(a)', iostat=status) S1_firstfile + S1_filename = trim(S1sigmaobsdir) // '/' // trim(S1_firstfile) + + ios = nf90_open(path=S1_filename,& + mode=NF90_NOWRITE,ncid=ncid) + call LIS_verify(ios,& + 'Error reading in S1 data dimensions: Error opening file '// S1_filename) + ios = nf90_inquire_dimension(ncid,1,yname,NY) + ios = nf90_inquire_dimension(ncid,2,xname,NX) + ios = nf90_close(ncid) + + do n=1,LIS_rc%nnest + S1_sigma_struc(n)%nc = NX + S1_sigma_struc(n)%nr = NY + + allocate(S1_sigma_struc(n)%sigma(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k))) + allocate(S1_sigma_struc(n)%sigmatime(& + LIS_rc%obs_lnc(k), LIS_rc%obs_lnr(k))) + S1_sigma_struc(n)%sigma = LIS_rc%udef + S1_sigma_struc(n)%sigmatime = -1 + + call LIS_registerAlarm("S1 backscatter read alarm",& + 86400.0, 86400.0) + + S1_sigma_struc(n)%startMode = .true. + enddo + + write(LIS_logunit,*) '[INFO] Created ESMF States to hold S1 observations data' + + end subroutine S1_sigmaVVSM_setup + +end module S1_sigmaVVSM_Mod diff --git a/lis/dataassim/obs/S1_sigmaVVSM/read_S1_sigmaVVSM.F90 b/lis/dataassim/obs/S1_sigmaVVSM/read_S1_sigmaVVSM.F90 new file mode 100644 index 000000000..7a397543f --- /dev/null +++ b/lis/dataassim/obs/S1_sigmaVVSM/read_S1_sigmaVVSM.F90 @@ -0,0 +1,466 @@ +!-----------------------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----------------------- +#include "LIS_misc.h" +!BOP +! !ROUTINE: read_S1_sigmaVVSM +! \label{read_S1_sigmaVVSM} +! +! !REVISION HISTORY: +! 29 Aug 2019: Hans Lievens; Initial Specification for snow depth +! 9 Mar 2021: Isis Brangers, Michel Bechtold; adaptation to backscatter +! 29 Jun 2022: Louise Busschaert; getting domain dims from files +! +! !INTERFACE: +subroutine read_S1_sigmaVVSM(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_S1_sigmaVVSM_obsId + use S1_sigmaVVSM_Mod, only : S1_sigma_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 backscatter 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) :: sigmaField,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 :: S1_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 :: sigma_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, "S1 backscatter read alarm") + + if(alarmCheck.or.S1_sigma_struc(n)%startMode) then + S1_sigma_struc(n)%startMode = .false. + + S1_sigma_struc(n)%sigma = LIS_rc%udef + S1_sigma_struc(n)%sigmatime = -1 + + call S1_sigmaVVSM_filename(S1_filename,obsdir,& + LIS_rc%yr,LIS_rc%mo,LIS_rc%da) + + inquire(file=S1_filename,exist=file_exists) + if(file_exists) then + + write(LIS_logunit,*) '[INFO] Reading S1 sigma data ',S1_filename + call read_S1_sigmaVVSM_data(n,k, S1_filename, S1_sigma_struc(n)%sigma) + + +!------------------------------------------------------------------------- +! 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) + S1_sigma_struc(n)%sigmatime(c,r) = gmt + + endif + enddo + enddo + + endif + endif + +!------------------------------------------------------------------------- +! Update the OBS_State +!------------------------------------------------------------------------- + + call ESMF_StateGet(OBS_State,"Observation01",sigmafield,& + rc=status) + call LIS_verify(status, 'Error: StateGet Observation01') + + call ESMF_FieldGet(sigmafield,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 - S1_sigma_struc(n)%sigmatime(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)) = & + S1_sigma_struc(n)%sigma(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_S1_sigmaVVSM_obsId)//char(0), & + n, k,OBS_state) + + sigma_current = LIS_rc%udef + call LIS_checkForValidObs(n,k,obsl,fnd,sigma_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_S1_sigma') + + if(LIS_rc%obs_ngrid(k).gt.0) then + +!linearly scale the observation err + ssdev = S1_sigma_struc(n)%ssdev + do t=1,LIS_rc%obs_ngrid(k) + if(obsl(t).ne.-9999.0) then + ssdev(t) = S1_sigma_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) + + call ESMF_AttributeSet(sigmafield,"Grid Number",& + gid,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status,'Error: AttributeSet in Grid Number') + + call ESMF_AttributeSet(sigmafield,"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_S1_sigmaVVSM + + + +!BOP +! +! !ROUTINE: read_S1_sigmaVVSM_data +! \label{read_S1_sigmaVVSM_data} +! +! !INTERFACE: +subroutine read_S1_sigmaVVSM_data(n, k, fname, sigma_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 S1_sigmaVVSM_Mod, only : S1_sigma_struc + + implicit none +! +! !INPUT PARAMETERS: +! + integer :: n + integer :: k + character (len=*) :: fname + +! !OUTPUT PARAMETERS: +! +! !DESCRIPTION: +! This subroutine reads the S1 NETCDF files +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[fname] name of the S1 sigma file +! \item[s0vvobs\_ip] sigma depth data processed to the LIS domain +! \end{description} +! +! !FILES USED: +! +! !REVISION HISTORY: +! +!EOP + real :: sigma(S1_sigma_struc(n)%nr,S1_sigma_struc(n)%nc) + real :: lat_nc(S1_sigma_struc(n)%nr) + real :: lat_nc_fold(S1_sigma_struc(n)%nr) + real :: lon_nc(S1_sigma_struc(n)%nc) + real :: sigma_ip(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + real :: sigma_fill(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + integer :: nsigma_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 :: sigmaId,latId,lonId + integer :: ios + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + +! sigma = LIS_rc%udef +! lat_nc = LIS_rc%udef +! lon_nc = LIS_rc%udef +! sigma_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, 'g0vv',sigmaid) + if(ios.lt.0) then + ios = nf90_inq_varid(nid, 's0vv',sigmaid) + endif + call LIS_verify(ios, 'Error nf90_inq_varid: backscatter 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, sigmaid, sigma) + call LIS_verify(ios, 'Error nf90_get_var: s0vv') + + ios = nf90_get_var(nid, latid, lat_nc_fold) + call LIS_verify(ios, 'Error nf90_get_var: lat') + !new g0 dataset of Hans with flipped lat variable + do i=1,size(lat_nc_fold) + lat_nc(i)=lat_nc_fold(size(lat_nc_fold)+1-i) + end do + 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)) + + + +! ! Initialize + sigma_ip = 0 + nsigma_ip = 0 + + ! Interpolate the data by averaging + do i=1,S1_sigma_struc(n)%nr + do j=1,S1_sigma_struc(n)%nc + + call latlon_to_ij(LIS_domain(n)%lisproj,& + !lat_nc(i),lon_nc(j),col,row) + lat_nc(S1_sigma_struc(n)%nr-(i-1)),lon_nc(j),col,row) + stn_col = nint(col) + stn_row = nint(row) + + if(sigma(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 + ! add in linear scale for later averaging + sigma_ip(stn_col,stn_row) = sigma_ip(stn_col,stn_row) + 10.**(sigma(i,j)/10.) + nsigma_ip(stn_col,stn_row) = nsigma_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(nsigma_ip(c,r).ne.0) then + ! average in linear scale + sigma_ip(c,r) = sigma_ip(c,r)/nsigma_ip(c,r) + else + sigma_ip(c,r) = LIS_rc%udef + endif + enddo + enddo + + + ! Fix stripes in obs caused by regridding (missing obs for LIS grid cell) + sigma_fill = 0 + nsigma_ip = 0 + do r=2,LIS_rc%obs_lnr(k)-1 + do c=2,LIS_rc%obs_lnc(k)-1 + if (sigma_ip(c,r).ne.LIS_rc%udef) then + sigma_fill(c-1,r-1)= sigma_fill(c-1,r-1) + sigma_ip(c,r) + sigma_fill(c-1,r)= sigma_fill(c-1,r) + sigma_ip(c,r) + sigma_fill(c-1,r+1)= sigma_fill(c-1,r+1) + sigma_ip(c,r) + sigma_fill(c,r-1)= sigma_fill(c,r-1) + sigma_ip(c,r) + sigma_fill(c,r+1)= sigma_fill(c,r+1) + sigma_ip(c,r) + sigma_fill(c+1,r-1)= sigma_fill(c+1,r-1) + sigma_ip(c,r) + sigma_fill(c+1,r)= sigma_fill(c+1,r) + sigma_ip(c,r) + sigma_fill(c+1,r+1)= sigma_fill(c+1,r+1) + sigma_ip(c,r) + + nsigma_ip(c-1,r-1)= nsigma_ip(c-1,r-1) + 1 + nsigma_ip(c-1,r)= nsigma_ip(c-1,r) + 1 + nsigma_ip(c-1,r+1)= nsigma_ip(c-1,r+1) + 1 + nsigma_ip(c,r-1)= nsigma_ip(c,r-1) + 1 + nsigma_ip(c,r+1)= nsigma_ip(c,r+1) + 1 + nsigma_ip(c+1,r-1)= nsigma_ip(c+1,r-1) + 1 + nsigma_ip(c+1,r)= nsigma_ip(c+1,r) + 1 + nsigma_ip(c+1,r+1)= nsigma_ip(c+1,r+1) + 1 + endif + enddo + enddo + + do r=1,LIS_rc%obs_lnr(k) + do c=1,LIS_rc%obs_lnc(k) + if (sigma_ip(c,r).eq.LIS_rc%udef.and.nsigma_ip(c,r).gt.0) then + ! average in linear scale + sigma_ip(c,r) = sigma_fill(c,r)/nsigma_ip(c,r) + endif + if (sigma_ip(c,r).ne.LIS_rc%udef) then + ! back to log scale, dB + sigma_ip(c,r) = 10.*LOG10(sigma_ip(c,r)) + endif + enddo + enddo + + endif + +#endif + +end subroutine read_S1_sigmaVVSM_data + + + +!BOP +! +! !ROUTINE: S1_sigmaVVSM_filename +! \label{S1_SWND_filename} +! +! !INTERFACE: +subroutine S1_sigmaVVSM_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 S1 filename +! +! The arguments are: +! \begin{description} +! \item[name] name of the S1 filename +! \item[ndir] name of the S1 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)//'/S1_g0_'//trim(fyr)//trim(fmo)//trim(fda)//'.nc' + +end subroutine S1_sigmaVVSM_filename + + + diff --git a/lis/dataassim/obs/S1_sigmaVVSM/write_S1_sigmaVVSMobs.F90 b/lis/dataassim/obs/S1_sigmaVVSM/write_S1_sigmaVVSMobs.F90 new file mode 100644 index 000000000..db1d08fa3 --- /dev/null +++ b/lis/dataassim/obs/S1_sigmaVVSM/write_S1_sigmaVVSMobs.F90 @@ -0,0 +1,127 @@ +!-----------------------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_S1_sigmaVVSMobs +! \label{write_S1_sigmaVVSMobs} +! +! !REVISION HISTORY: +! 30 Aug 2019: Hans Lievens; Initial Specification +! +! !INTERFACE: +subroutine write_S1_sigmaVVSMobs(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) +! S1 observations to a file +! +!EOP + type(ESMF_Field) :: sigmaField + logical :: data_update + real, pointer :: sigmaobs(:) + 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",sigmaField, & + rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(sigmaField, localDE=0, farrayPtr=sigmaobs, rc=status) + call LIS_verify(status) + + if(LIS_masterproc) then + ftn = LIS_getNextUnitNumber() + call S1_sigmaVVSM_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) = sigmaobs(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,sigmaobs) + + if(LIS_masterproc) then + call LIS_releaseUnitNumber(ftn) + endif + + endif + +end subroutine write_S1_sigmaVVSMobs + +!BOP +! !ROUTINE: S1_sigmaVVSM_obsname +! \label{S1_sigmaVVSM_obsname} +! +! !INTERFACE: +subroutine S1_sigmaVVSM_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 S1_sigmaVVSM_obsname diff --git a/lis/dataassim/obs/S1_sigmaVVVHSMLAI/S1_sigmaVVVHSMLAI_Mod.F90 b/lis/dataassim/obs/S1_sigmaVVVHSMLAI/S1_sigmaVVVHSMLAI_Mod.F90 new file mode 100644 index 000000000..d3b9eb7b9 --- /dev/null +++ b/lis/dataassim/obs/S1_sigmaVVVHSMLAI/S1_sigmaVVVHSMLAI_Mod.F90 @@ -0,0 +1,328 @@ +!-----------------------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: S1_sigmaVVVHSMLAI_Mod +! +! !DESCRIPTION: +! This module contains interfaces and subroutines to +! handle S1 backscatter retrievals. +! +! !REVISION HISTORY: +! 29 Aug 2019 Hans Lievens; Initial Specification for sigma depth +! 9 Mar 2021 Isis Brangers, Michel Bechtold; Adaptation for backscatter +! 29 Jun 2022: Louise Busschaert; getting domain dims from files +! +module S1_sigmaVVVHSMLAI_Mod +! !USES: + use ESMF +!EOP + implicit none + + PRIVATE + +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: S1_sigmaVVVHSMLAI_setup +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + PUBLIC :: S1_sigma_struc + + type, public :: S1_sigma_dec + logical :: startMode + real :: ssdev_inp + integer :: s_vvField + integer :: s_vhField + integer :: nc,nr + real, allocatable :: s_vv(:,:) + real, allocatable :: s_vh(:,:) + real, allocatable :: sigmatime(:,:) + end type S1_sigma_dec + + type(S1_sigma_dec), allocatable :: S1_sigma_struc(:) + +contains +!BOP +! +! !ROUTINE: S1_sigmaVVVHSMLAI_setup +! \label{S1_sigmaVVVHSMLAI_setup} +! +! !INTERFACE: + subroutine S1_sigmaVVVHSMLAI_setup(k, OBS_State, OBS_Pert_State) +! !USES: + use ESMF + 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 structures required for S1 backscatter 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,m + 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 :: S1sigmaobsdir + character*80 :: S1_firstfile + character*80 :: S1_filename + 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(:) + real :: dx, dy + integer :: NX, NY + integer :: ncid + integer :: flist + integer :: ios + character*100 :: infile + character*100 :: xname, yname + + allocate(S1_sigma_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,"S1 backscatter data directory:",& + rc=status) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetattribute(LIS_config,S1sigmaobsdir,& + rc=status) + call LIS_verify(status,'S1 backscatter data directory: not defined') + + call ESMF_AttributeSet(OBS_State(n),"Data Directory",& + S1sigmaobsdir, 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 S1 backscatter 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) + + do m=1,LIS_rc%nobtypes(k) + write(unit=temp,fmt='(i2.2)') m + read(unit=temp,fmt='(2a1)') vid + + obsField = ESMF_FieldCreate(arrayspec=realarrspec,& + grid=LIS_vecGrid(n),& + name="Observation"//vid(1)//vid(2),rc=status) + call LIS_verify(status) + + call ESMF_StateAdd(OBS_State(n),(/obsField/),rc=status) + call LIS_verify(status) + enddo + + allocate(ssdev(LIS_rc%obs_ngrid(n))) + + if(trim(LIS_rc%perturb_obs(k)).ne."none") then + allocate(obs_pert%vname(LIS_rc%nobtypes(k))) + allocate(obs_pert%perttype(LIS_rc%nobtypes(k))) + allocate(obs_pert%ssdev(LIS_rc%nobtypes(k))) + allocate(obs_pert%stdmax(LIS_rc%nobtypes(k))) + allocate(obs_pert%zeromean(LIS_rc%nobtypes(k))) + allocate(obs_pert%tcorr(LIS_rc%nobtypes(k))) + allocate(obs_pert%xcorr(LIS_rc%nobtypes(k))) + allocate(obs_pert%ycorr(LIS_rc%nobtypes(k))) + allocate(obs_pert%ccorr(LIS_rc%nobtypes(k),LIS_rc%nobtypes(k))) + + call LIS_readPertAttributes(LIS_rc%nobtypes(k),& + LIS_rc%obspertAttribfile(k),obs_pert) + + do m=1,LIS_rc%nobtypes(k) + ssdev = obs_pert%ssdev(m) + S1_sigma_struc(n)%ssdev_inp =obs_pert%ssdev(m) + + write(unit=temp,fmt='(i2.2)') m + read(unit=temp,fmt='(2a1)') vid + + 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(m), 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(m), rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(pertField(n),"Ensure Zero Mean",& + obs_pert%zeromean(m),rc=status) + call LIS_verify(status) + + call ESMF_AttributeSet(pertField(n),"Temporal Correlation Scale",& + obs_pert%tcorr(m),rc=status) + + call ESMF_AttributeSet(pertField(n),"X Correlation Scale",& + obs_pert%xcorr(m),rc=status) + + call ESMF_AttributeSet(pertField(n),"Y Correlation Scale",& + obs_pert%ycorr(m),rc=status) + + call ESMF_AttributeSet(pertField(n),"Cross Correlation Strength",& + obs_pert%ccorr(m,:),itemCount=LIS_rc%nobtypes(k),rc=status) + + call ESMF_StateAdd(OBS_Pert_State(n),(/pertField(n)/),rc=status) + call LIS_verify(status) + enddo + deallocate(ssdev) + endif + + deallocate(vname) + deallocate(varmax) + deallocate(varmin) + enddo +!------------------------------------------------------------- +! set up the S1 domain %and interpolation weights. +!------------------------------------------------------------- + + ! Open first file to get nx ny dimensions + call system('ls ./' // trim(S1sigmaobsdir) // ' > ./S1_listfiles.txt') + + open(flist, file=trim('./S1_listfiles.txt'), & + status='old', iostat=status) + + read(flist, '(a)', iostat=status) S1_firstfile + S1_filename = trim(S1sigmaobsdir) // '/' // trim(S1_firstfile) + + ios = nf90_open(path=S1_filename,& + mode=NF90_NOWRITE,ncid=ncid) + call LIS_verify(ios,& + 'Error reading in S1 data dimensions: Error opening file '// S1_filename) + ios = nf90_inquire_dimension(ncid,1,yname,NY) + ios = nf90_inquire_dimension(ncid,2,xname,NX) + ios = nf90_close(ncid) + + do n=1,LIS_rc%nnest + S1_sigma_struc(n)%nc = NX + S1_sigma_struc(n)%nr = NY + + allocate(S1_sigma_struc(n)%s_vv(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k))) + allocate(S1_sigma_struc(n)%s_vh(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k))) + allocate(S1_sigma_struc(n)%sigmatime(& + LIS_rc%obs_lnc(k), LIS_rc%obs_lnr(k))) + S1_sigma_struc(n)%s_vv = LIS_rc%udef + S1_sigma_struc(n)%s_vh = LIS_rc%udef + S1_sigma_struc(n)%sigmatime = -1 + + call LIS_registerAlarm("S1 backscatter read alarm",& + 86400.0, 86400.0) + + S1_sigma_struc(n)%startMode = .true. + enddo + + write(LIS_logunit,*) '[INFO] Created ESMF States to hold S1 observations data' + + end subroutine S1_sigmaVVVHSMLAI_setup + +end module S1_sigmaVVVHSMLAI_Mod + diff --git a/lis/dataassim/obs/S1_sigmaVVVHSMLAI/read_S1_sigmaVVVHSMLAI.F90 b/lis/dataassim/obs/S1_sigmaVVVHSMLAI/read_S1_sigmaVVVHSMLAI.F90 new file mode 100644 index 000000000..4a275696e --- /dev/null +++ b/lis/dataassim/obs/S1_sigmaVVVHSMLAI/read_S1_sigmaVVVHSMLAI.F90 @@ -0,0 +1,488 @@ +!-----------------------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----------------------- +#include "LIS_misc.h" +!BOP +! !ROUTINE: read_S1_sigmaVVVHSMLAI +! \label{read_S1_sigmaVVVHSMLAI} +! +! !REVISION HISTORY: +! 29 Aug 2019: Hans Lievens; Initial Specification for snow depth +! 9 Mar 2021: Isis Brangers, Michel Bechtold; adaptation to backscatter +! 29 Jun 2022: Louise Busschaert; getting domain dims from files +! !INTERFACE: +subroutine read_S1_sigmaVVVHSMLAI(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_S1_sigmaVVVHSMLAI_obsId + use S1_sigmaVVVHSMLAI_Mod, only : S1_sigma_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 backscatter 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) :: s_vvField,s_vhField,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 :: S1_filename + integer :: ftn + real :: lon, lhour + real :: gmt + real :: dt + integer :: zone + integer :: grid_index + real :: ssdev(LIS_rc%obs_ngrid(k)) + real, pointer :: s_vv(:), s_vh(:) + integer :: gid(LIS_rc%obs_ngrid(k)) + integer :: assimflag_vv(LIS_rc%obs_ngrid(k)) + integer :: assimflag_vh(LIS_rc%obs_ngrid(k)) + integer :: status, iret, ierr + integer :: fnd + real :: svv_current(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + real :: svh_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, "S1 backscatter read alarm") + + if(alarmCheck.or.S1_sigma_struc(n)%startMode) then + S1_sigma_struc(n)%startMode = .false. + + S1_sigma_struc(n)%s_vv = LIS_rc%udef + S1_sigma_struc(n)%s_vh = LIS_rc%udef + S1_sigma_struc(n)%sigmatime = -1 + + call S1_sigmaVVVHSMLAI_filename(S1_filename,obsdir,& + LIS_rc%yr,LIS_rc%mo,LIS_rc%da) + + inquire(file=S1_filename,exist=file_exists) + if(file_exists) then + + write(LIS_logunit,*) '[INFO] Reading S1 sigma data ',S1_filename + call read_S1_sigmaVVVHSMLAI_data(n,k, S1_filename, S1_sigma_struc(n)%s_vv, & + S1_sigma_struc(n)%s_vh) + + +!------------------------------------------------------------------------- +! 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) + S1_sigma_struc(n)%sigmatime(c,r) = gmt + + endif + enddo + enddo + + endif + endif + +!------------------------------------------------------------------------- +! Update the OBS_State +!------------------------------------------------------------------------- + + call ESMF_StateGet(OBS_State,"Observation01",s_vvField,& + rc=status) + call LIS_verify(status, 'Error: StateGet Observation01') + + call ESMF_StateGet(OBS_State,"Observation02",s_vhField,& + rc=status) + call LIS_verify(status, 'Error: StateGet Observation02') + + call ESMF_FieldGet(s_vvField,localDE=0,farrayPtr=s_vv,rc=status) + call LIS_verify(status, 'Error: FieldGet') + + call ESMF_FieldGet(s_vhField,localDE=0,farrayPtr=s_vh,rc=status) + call LIS_verify(status, 'Error: FieldGet') + + s_vv = LIS_rc%udef + s_vh = 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 - S1_sigma_struc(n)%sigmatime(c,r))*3600.0 + lon = LIS_obs_domain(n,k)%lon(grid_index) + + if(dt.ge.0.and.dt.lt.LIS_rc%ts) then + s_vv(LIS_obs_domain(n,k)%gindex(c,r)) = & + S1_sigma_struc(n)%s_vv(c,r) + s_vh(LIS_obs_domain(n,k)%gindex(c,r)) = & + S1_sigma_struc(n)%s_vh(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_S1_sigmaVVVHSMLAI_obsId)//char(0), & + n, k,OBS_state) + + svv_current = LIS_rc%udef + svh_current = LIS_rc%udef + call LIS_checkForValidObs(n,k,s_vv,fnd,svv_current) + call LIS_checkForValidObs(n,k,s_vh,fnd,svh_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(s_vv(t).ne.-9999.0) then + assimflag_vv(t) = 1 + else + assimflag_vv(t) = 0 + endif + enddo + do t=1,LIS_rc%obs_ngrid(k) + gid(t) = t + if(s_vh(t).ne.-9999.0) then + assimflag_vh(t) = 1 + else + assimflag_vh(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_S1_sigma') + + call ESMF_StateGet(OBS_Pert_State,"Observation02",pertfield,& + rc=status) + call LIS_verify(status, 'ESMF_StateGet for Observation02 for OBS_Pert_State failed in read_S1_sigma') + + if(LIS_rc%obs_ngrid(k).gt.0) then + + !linearly scale the observation err + ssdev = S1_sigma_struc(n)%ssdev_inp + do t=1,LIS_rc%obs_ngrid(k) + if(s_vv(t).ne.-9999.0) then + ssdev(t) = S1_sigma_struc(n)%ssdev_inp !+ 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) + + call ESMF_AttributeSet(s_vvField,"Grid Number",& + gid,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status,'Error: AttributeSet in Grid Number') + + call ESMF_AttributeSet(s_vvField,"Assimilation Flag",& + assimflag_vv,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status, 'Error: AttributeSet in Assimilation Flag') + + call ESMF_AttributeSet(s_vhField,"Grid Number",& + gid,itemCount=LIS_rc%obs_ngrid(k),rc=status) + call LIS_verify(status,'Error: AttributeSet in Grid Number') + + call ESMF_AttributeSet(s_vhField,"Assimilation Flag",& + assimflag_vh,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_S1_sigmaVVVHSMLAI + + + +!BOP +! +! !ROUTINE: read_S1_sigmaVVVHSMLAI_data +! \label{read_S1_sigmaVVVHSMLAI_data} +! +! !INTERFACE: +subroutine read_S1_sigmaVVVHSMLAI_data(n, k, fname, svv_ip, svh_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 S1_sigmaVVVHSMLAI_Mod, only : S1_sigma_struc + + implicit none +! +! !INPUT PARAMETERS: +! + integer :: n + integer :: k + character (len=*) :: fname + +! !OUTPUT PARAMETERS: +! +! !DESCRIPTION: +! This subroutine reads the S1 NETCDF files +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest +! \item[fname] name of the S1 sigma file +! \item[svvobs\_ip] sigma depth data processed to the LIS domain +! \end{description} +! +! !FILES USED: +! +! !REVISION HISTORY: +! +!EOP + real :: s_vv(S1_sigma_struc(n)%nr,S1_sigma_struc(n)%nc) + real :: s_vh(S1_sigma_struc(n)%nr,S1_sigma_struc(n)%nc) + real :: lat_nc(S1_sigma_struc(n)%nr) + real :: lat_nc_fold(S1_sigma_struc(n)%nr) + real :: lon_nc(S1_sigma_struc(n)%nc) + real :: svv_ip(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + real :: svh_ip(LIS_rc%obs_lnc(k),LIS_rc%obs_lnr(k)) + integer :: ns_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 :: svvId,svhId,latId,lonId + integer :: ios + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + +! s_vv = LIS_rc%udef +! lat_nc = LIS_rc%udef +! lon_nc = LIS_rc%udef +! svv_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, 'g0vv',svvid) + if(ios.lt.0) then + ios = nf90_inq_varid(nid, 's0vv',svvid) + endif + call LIS_verify(ios, 'Error nf90_inq_varid: backscatter data') + + ios = nf90_inq_varid(nid, 'g0vh',svhid) + if(ios.lt.0) then + ios = nf90_inq_varid(nid, 's0vh',svhid) + endif + call LIS_verify(ios, 'Error nf90_inq_varid: backscatter 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, svvid, s_vv) + call LIS_verify(ios, 'Error nf90_get_var: s0vv') + + ios = nf90_get_var(nid, svhid, s_vh) + call LIS_verify(ios, 'Error nf90_get_var: s0vh') + + ios = nf90_get_var(nid, latid, lat_nc_fold) + call LIS_verify(ios, 'Error nf90_get_var: lat') + !new g0 dataset of Hans with flipped lat variable + do i=1,size(lat_nc_fold) + lat_nc(i)=lat_nc_fold(size(lat_nc_fold)+1-i) + end do + 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)) + + + +! ! Initialize + svv_ip = 0 + svh_ip = 0 + ns_ip = 0 + + ! Interpolate the data by averaging + do i=1,S1_sigma_struc(n)%nr + do j=1,S1_sigma_struc(n)%nc + + call latlon_to_ij(LIS_domain(n)%lisproj,& + !lat_nc(i),lon_nc(j),col,row) + lat_nc(S1_sigma_struc(n)%nr-(i-1)),lon_nc(j),col,row) + stn_col = nint(col) + stn_row = nint(row) + + if(s_vv(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 + ! add in linear scale for later averaging + svv_ip(stn_col,stn_row) = svv_ip(stn_col,stn_row) + 10.**(s_vv(i,j)/10.) + svh_ip(stn_col,stn_row) = svh_ip(stn_col,stn_row) + 10.**(s_vh(i,j)/10.) + ns_ip(stn_col,stn_row) = ns_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(ns_ip(c,r).ne.0.and.svv_ip(c,r).gt.0.0.and.svh_ip(c,r).gt.0.0) then + ! average in linear scale + svv_ip(c,r) = svv_ip(c,r)/ns_ip(c,r) + svh_ip(c,r) = svh_ip(c,r)/ns_ip(c,r) + else + svv_ip(c,r) = LIS_rc%udef + svh_ip(c,r) = LIS_rc%udef + endif + enddo + enddo + + + do r=1,LIS_rc%obs_lnr(k) + do c=1,LIS_rc%obs_lnc(k) + if(svv_ip(c,r).ne.LIS_rc%udef) then + ! back to log scale, dB + svv_ip(c,r) = 10.*LOG10(svv_ip(c,r)); + endif + if(svh_ip(c,r).ne.LIS_rc%udef) then + ! back to log scale, dB + svh_ip(c,r) = 10.*LOG10(svh_ip(c,r)); + endif + enddo + enddo + + endif + +#endif + +end subroutine read_S1_sigmaVVVHSMLAI_data + + + +!BOP +! +! !ROUTINE: S1_sigmaVVVHSMLAI_filename +! \label{S1_SWND_filename} +! +! !INTERFACE: +subroutine S1_sigmaVVVHSMLAI_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 S1 filename +! +! The arguments are: +! \begin{description} +! \item[name] name of the S1 filename +! \item[ndir] name of the S1 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)//'/S1_g0_'//trim(fyr)//trim(fmo)//trim(fda)//'.nc' + +end subroutine S1_sigmaVVVHSMLAI_filename + + + diff --git a/lis/dataassim/obs/S1_sigmaVVVHSMLAI/write_S1_sigmaVVVHSMLAIobs.F90 b/lis/dataassim/obs/S1_sigmaVVVHSMLAI/write_S1_sigmaVVVHSMLAIobs.F90 new file mode 100644 index 000000000..826a1dcfe --- /dev/null +++ b/lis/dataassim/obs/S1_sigmaVVVHSMLAI/write_S1_sigmaVVVHSMLAIobs.F90 @@ -0,0 +1,133 @@ +!-----------------------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_S1_sigmaVVVHSMLAIobs +! \label{write_S1_sigmaVVVHSMLAIobs} +! +! !REVISION HISTORY: +! 30 Aug 2019: Hans Lievens; Initial Specification +! +! !INTERFACE: +subroutine write_S1_sigmaVVVHSMLAIobs(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) +! S1 observations to a file +! +!EOP + type(ESMF_Field) :: s_vvField, s_vhField + logical :: data_update + real, pointer :: obs_vv(:), obs_vh(:) + character*100 :: obsname_vv, obsname_vh + 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",s_vvField, & + rc=status) + call LIS_verify(status) + + call ESMF_StateGet(OBS_State, "Observation02",s_vhField, & + rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(s_vvField, localDE=0, farrayPtr=obs_vv, rc=status) + call LIS_verify(status) + + call ESMF_FieldGet(s_vhField, localDE=0, farrayPtr=obs_vh, rc=status) + call LIS_verify(status) + +! ! save VV + if(LIS_masterproc) then + ftn = LIS_getNextUnitNumber() + call S1_sigmaVVVHSMLAI_obsname(n,k,obsname_vv,obsname_vh) + + call LIS_create_output_directory('DAOBS') + + open(ftn,file=trim(obsname_vv), form='unformatted') + endif + call LIS_writevar_gridded_obs(ftn,n,k,obs_vv) + if(LIS_masterproc) then + call LIS_releaseUnitNumber(ftn) + endif +! ! save VH + if(LIS_masterproc) then + ftn = LIS_getNextUnitNumber() + call S1_sigmaVVVHSMLAI_obsname(n,k,obsname_vv,obsname_vh) + + call LIS_create_output_directory('DAOBS') + + open(ftn,file=trim(obsname_vh), form='unformatted') + endif + call LIS_writevar_gridded_obs(ftn,n,k,obs_vh) + if(LIS_masterproc) then + call LIS_releaseUnitNumber(ftn) + endif + + endif + +end subroutine write_S1_sigmaVVVHSMLAIobs + +!BOP +! !ROUTINE: S1_sigmaVVVHSMLAI_obsname +! \label{S1_sigmaVVVHSMLAI_obsname} +! +! !INTERFACE: +subroutine S1_sigmaVVVHSMLAI_obsname(n,k,obsname_vv,obsname_vh) +! !USES: + use LIS_coreMod, only : LIS_rc + +! !ARGUMENTS: + integer :: n + integer :: k + character(len=*) :: obsname_vv, obsname_vh +! +! !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_vv = trim(LIS_rc%odir)//'/DAOBS/'//cdate1(1:6)//& + '/LISDAOBS_VV_'//cdate1// & + trim(cda)//trim(cdate)//'.1gs4r' + obsname_vh = trim(LIS_rc%odir)//'/DAOBS/'//cdate1(1:6)//& + '/LISDAOBS_VH_'//cdate1// & + trim(cda)//trim(cdate)//'.1gs4r' +end subroutine S1_sigmaVVVHSMLAI_obsname diff --git a/lis/interp/map_utils.F90 b/lis/interp/map_utils.F90 index 43df2ac20..02b263ab0 100644 --- a/lis/interp/map_utils.F90 +++ b/lis/interp/map_utils.F90 @@ -145,6 +145,8 @@ module map_utils ! Brent L. Shaw, NOAA/FSL (CSU/CIRA) ! 09 Apr 2001 - Added compare\_projections routine to compare two ! sets of projection parameters. +! 30 Jun 2021 - Michel Bechtold, Samuel Scherrer, bug fix for latlon projection +! occurred at edges of subdomains in parallel mode implicit none @@ -681,6 +683,7 @@ SUBROUTINE llij_gauss(lat, lon, proj, i, j) slon360 = proj%lon1 ENDIF deltalon = lon360 - slon360 + IF (deltalon .LT. 0) deltalon = deltalon + 360. ! Compute i/j i = deltalon/proj%dlon + 1. @@ -1131,7 +1134,7 @@ SUBROUTINE llij_latlon(lat, lon, proj, i, j) REAL :: deltalat REAL :: deltalon - REAL :: lon360 + REAL :: lon360,slon360 ! Compute deltalat and deltalon as the difference between the input @@ -1146,9 +1149,19 @@ SUBROUTINE llij_latlon(lat, lon, proj, i, j) ELSE lon360 = lon ENDIF - deltalon = lon360 - proj%lon1 - IF (deltalon .LT. 0) deltalon = deltalon + 360. - + ! deltalon = lon360 - proj%lon1 + ! IF (deltalon .LT. 0) deltalon = deltalon + 360. + ! MB: June,2021: deltalon+360 causes + ! bug at edges of subdomains in parallel mode. + ! Correct handling is like done in llij_gauss + IF (proj%lon1 .LT. 0) THEN + slon360 = proj%lon1 + 360. + ELSE + slon360 = proj%lon1 + ENDIF + deltalon = lon360 - slon360 + IF (deltalon .LT. 0) deltalon = deltalon + 360. + ! Compute i/j i = deltalon/proj%dlon + 1. j = deltalat/proj%dlat + 1. diff --git a/lis/make/default.cfg b/lis/make/default.cfg index 363a09167..53a1839ef 100644 --- a/lis/make/default.cfg +++ b/lis/make/default.cfg @@ -569,6 +569,11 @@ enabled: True macro: RTMS_TAU_OMEGA path: rtms/TauOmegaRTM +[WCM] +enabled: True +macro: RTMS_WCM +path: rtms/WCMRTM + #}}} # @@ -746,6 +751,16 @@ enabled: True macro: DA_OBS_SMMR_SNWD path: dataassim/obs/SMMR_SNWD +[DA OBS S1_sigmaVVSM] +enabled: True +macro: DA_OBS_S1_sigmaVVSM +path: dataassim/obs/S1_sigmaVVSM + +[DA OBS S1_sigmaVVVHSMLAI] +enabled: True +macro: DA_OBS_S1_sigmaVVVHSMLAI +path: dataassim/obs/S1_sigmaVVVHSMLAI + [DA OBS SSMI_SNWD] enabled: True macro: DA_OBS_SSMI_SNWD @@ -1141,13 +1156,18 @@ macro: SM_NOAHMP_3_6 path: surfacemodels/land/noahmp.3.6 dependent_comps: WRF coupling, CMEM, + WCM, virtual_da, virtual_optue, virtual_routing, virtual_irrigation WRF coupling path: surfacemodels/land/noahmp.3.6/cpl_wrf_noesmf CMEM path: surfacemodels/land/noahmp.3.6/sfc_cmem3 +WCM path: surfacemodels/land/noahmp.3.6/sfc_wcm virtual_da path: surfacemodels/land/noahmp.3.6/da_soilm, + surfacemodels/land/noahmp.3.6/da_soilmLAI, + surfacemodels/land/noahmp.3.6/da_sig0VV, + surfacemodels/land/noahmp.3.6/da_sig0VVVH, surfacemodels/land/noahmp.3.6/da_snow, surfacemodels/land/noahmp.3.6/da_snowSVM, surfacemodels/land/noahmp.3.6/da_LAI, diff --git a/lis/plugins/LIS_DAobs_pluginMod.F90 b/lis/plugins/LIS_DAobs_pluginMod.F90 index 14491b275..1cabf5202 100644 --- a/lis/plugins/LIS_DAobs_pluginMod.F90 +++ b/lis/plugins/LIS_DAobs_pluginMod.F90 @@ -22,6 +22,7 @@ module LIS_DAobs_pluginMod ! !REVISION HISTORY: ! 27 Feb 2005; Sujay Kumar Initial Specification ! 11 Aug 2016: Mahdi Navari, PILDAS added +! 30 Jun 2021: Sara Modanesi, Michel Bechtold; Added Sentinel 1 observations ! !EOP implicit none @@ -217,6 +218,14 @@ subroutine LIS_DAobs_plugin use SMMRSNWDsnow_Mod, only : SMMRSNWDsnow_setup #endif +#if ( defined DA_OBS_S1_sigmaVVSM ) + use S1_sigmaVVSM_Mod, only : S1_sigmaVVSM_setup +#endif + +#if ( defined DA_OBS_S1_sigmaVVVHSMLAI ) + use S1_sigmaVVVHSMLAI_Mod, only : S1_sigmaVVVHSMLAI_setup +#endif + #if ( defined DA_OBS_SSMI_SNWD ) use SSMISNWDsnow_Mod, only : SSMISNWDsnow_setup #endif @@ -403,6 +412,14 @@ subroutine LIS_DAobs_plugin external read_SMMRSNWDsnow, write_SMMRSNWDsnowobs #endif +#if ( defined DA_OBS_S1_sigmaVVSM) + external read_S1_sigmaVVSM, write_S1_sigmaVVSMobs +#endif + +#if ( defined DA_OBS_S1_sigmaVVVHSMLAI) + external read_S1_sigmaVVVHSMLAI, write_S1_sigmaVVVHSMLAIobs +#endif + #if ( defined DA_OBS_SSMI_SNWD ) external read_SSMISNWDsnow, write_SSMISNWDsnowobs #endif @@ -670,6 +687,28 @@ subroutine LIS_DAobs_plugin write_SMMRSNWDsnowobs) #endif +#if ( defined DA_OBS_S1_sigmaVVSM ) +!S1 backscatter obs VVSM + call registerdaobsclass(trim(LIS_S1_sigmaVVSM_obsId),"LSM") + call registerdaobssetup(trim(LIS_S1_sigmaVVSM_obsId)//char(0), & + S1_sigmaVVSM_setup) + call registerreaddaobs(trim(LIS_S1_sigmaVVSM_obsId)//char(0), & + read_S1_sigmaVVSM) + call registerwritedaobs(trim(LIS_S1_sigmaVVSM_obsId)//char(0), & + write_S1_sigmaVVSMobs) +#endif + +#if ( defined DA_OBS_S1_sigmaVVVHSMLAI ) +!S1 backscatter obs VVVHSMLAI + call registerdaobsclass(trim(LIS_S1_sigmaVVVHSMLAI_obsId),"LSM") + call registerdaobssetup(trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0), & + S1_sigmaVVVHSMLAI_setup) + call registerreaddaobs(trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0), & + read_S1_sigmaVVVHSMLAI) + call registerwritedaobs(trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0), & + write_S1_sigmaVVVHSMLAIobs) +#endif + #if ( defined DA_OBS_SSMI_SNWD ) !SSMI SNWD snow obs call registerdaobsclass(trim(LIS_SSMISNWDsnowobsId),"LSM") diff --git a/lis/plugins/LIS_RTM_pluginMod.F90 b/lis/plugins/LIS_RTM_pluginMod.F90 index d7fcff1d7..e32fc6732 100644 --- a/lis/plugins/LIS_RTM_pluginMod.F90 +++ b/lis/plugins/LIS_RTM_pluginMod.F90 @@ -22,6 +22,7 @@ module LIS_RTM_pluginMod ! 21 Mar 09 Sujay Kumar Initial Specification ! 20 Oct 10 Yudong Tian Added CRTM2 and CRTM2 EMonly ! 08 Feb 11 Yudong Tian Added CMEM3 +! 04 Sep 20 Sara Modanesi Added WCM (Water Cloud Model) ! !EOP implicit none @@ -75,6 +76,12 @@ subroutine LIS_RTM_plugin TauOmegaRTM_run #endif +#if ( defined RTMS_WCM ) + use WCMRTM_Mod, only : WCMRTM_initialize, WCMRTM_f2t,& + WCMRTM_geometry, & + WCMRTM_run +#endif + #if ( defined RTMS_CRTM ) ! CRTM 1.x ! call registerrtminit(trim(LIS_crtmId)//char(0), CRTM_Kmatrix_initialize) @@ -117,6 +124,15 @@ subroutine LIS_RTM_plugin TauOmegaRTM_geometry) call registerrtmrun(trim(LIS_tauomegartmId)//char(0),TauOmegaRTM_run) #endif + +#if ( defined RTMS_WCM ) +!WCMRTM + call registerrtminit(trim(LIS_wcmrtmId)//char(0),WCMRTM_initialize) + call registerrtmf2t(trim(LIS_wcmrtmId)//char(0),WCMRTM_f2t) + call registergeometry2rtm(trim(LIS_wcmrtmId)//char(0), & + WCMRTM_geometry) + call registerrtmrun(trim(LIS_wcmrtmId)//char(0),WCMRTM_run) +#endif #endif end subroutine LIS_RTM_plugin end module LIS_RTM_pluginMod diff --git a/lis/plugins/LIS_lsmda_pluginMod.F90 b/lis/plugins/LIS_lsmda_pluginMod.F90 index d7915469e..3d7a7353e 100644 --- a/lis/plugins/LIS_lsmda_pluginMod.F90 +++ b/lis/plugins/LIS_lsmda_pluginMod.F90 @@ -32,6 +32,7 @@ module LIS_lsmda_pluginMod ! enabled the compilation of JULES 5.3 DA ! 13 Dec 2019: Eric Kemp, replaced LDTSI with USAFSI ! 17 Feb 2020: Yeosang Yoon, added SNODEP & USAFSI Assimilation for Jules 5.x +! 30 Jun 2021: Sara Modanesi, Michel Bechtold, added Sentinel-1 data assimilation ! 06 Jun 2022: Yonghwan Kwon, added SMAP_E_OPL soil moisture Assimilation ! added GVF data assimilation ! @@ -189,6 +190,7 @@ subroutine LIS_lsmda_plugin #if ( defined SM_NOAHMP_3_6 ) use noahmp36_dasoilm_Mod + use noahmp36_dasoilmLAI_Mod use noahmp36_dasnow_Mod use noahmp36_dasnodep_Mod use noahmp36_tws_DAlogMod, only : noahmp36_tws_DAlog @@ -409,6 +411,19 @@ subroutine LIS_lsmda_plugin external noahmp36_scale_soilm external noahmp36_descale_soilm external noahmp36_updatesoilm + external noahmp36_getSig0VVpred + external noahmp36_qc_Sig0VVobs + external noahmp36_getSig0VHpred + external noahmp36_qc_Sig0VHobs + external noahmp36_getSig0VVVHpred + external noahmp36_qc_Sig0VVVHobs + + external noahmp36_getsoilmLAI + external noahmp36_setsoilmLAI + external noahmp36_qcsoilmLAI + external noahmp36_scale_soilmLAI + external noahmp36_descale_soilmLAI + external noahmp36_updatesoilmLAI external noahmp36_getsnowvars external noahmp36_setsnowvars @@ -1968,6 +1983,53 @@ subroutine LIS_lsmda_plugin #if ( defined SM_NOAHMP_3_6 ) + +! S1 backscatter DA VVSM + call registerlsmdainit(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_dasoilm_init) + call registerlsmdagetstatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_getsoilm) + call registerlsmdasetstatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_setsoilm) + call registerlsmdagetobspred(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_getSig0VVpred) + call registerlsmdaqcstate(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_qcsoilm) + call registerlsmdaqcobsstate(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_qc_Sig0VVobs) + call registerlsmdascalestatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_scale_soilm) + call registerlsmdadescalestatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_descale_soilm) + call registerlsmdaupdatestate(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVSM_obsId)//char(0),noahmp36_updatesoilm) + + +! S1 backscatter DA VVVHSMLAI + call registerlsmdainit(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_dasoilmLAI_init) + call registerlsmdagetstatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_getsoilmLAI) + call registerlsmdasetstatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_setsoilmLAI) + call registerlsmdagetobspred(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_getSig0VVVHpred) + call registerlsmdaqcstate(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_qcsoilmLAI) + call registerlsmdaqcobsstate(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_qc_Sig0VVVHobs) + call registerlsmdascalestatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_scale_soilmLAI) + call registerlsmdadescalestatevar(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_descale_soilmLAI) + call registerlsmdaupdatestate(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_updatesoilmLAI) + call registerlsmdaobstransform(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_transform_veg) + call registerlsmdamapobstolsm(trim(LIS_noahmp36Id)//"+"//& + trim(LIS_S1_sigmaVVVHSMLAI_obsId)//char(0),noahmp36_map_veg) + + ! Noahmp-3.6 synthetic soil moisture call registerlsmdainit(trim(LIS_noahmp36Id)//"+"//& trim(LIS_synsmId)//char(0),noahmp36_dasoilm_init) diff --git a/lis/plugins/LIS_lsmrtm_pluginMod.F90 b/lis/plugins/LIS_lsmrtm_pluginMod.F90 index 421c62db6..cf2c11ae4 100644 --- a/lis/plugins/LIS_lsmrtm_pluginMod.F90 +++ b/lis/plugins/LIS_lsmrtm_pluginMod.F90 @@ -23,6 +23,7 @@ module LIS_lsmrtm_pluginMod ! 08 Feb 11 Yudong Tian; Reistered for CMEM3 ! 18 Dec 18 Peter Shellito; Registered for NoahMp36 CMEM3 ! 21 Feb 19 Peter Shellito; Registered for Noah36 CMEM3 +! 04 Sep 20 Sara Modanesi; Registered for NoahMp36 WCM ! !EOP implicit none @@ -89,6 +90,9 @@ subroutine LIS_lsmrtm_plugin #if ( defined RTMS_CMEM ) external noahmp36_sfc2cmem3 #endif +#if ( defined RTMS_WCM ) + external noahmp36_sfc2wcm +#endif #endif #if ( defined SM_MOSAIC ) @@ -174,6 +178,19 @@ subroutine LIS_lsmrtm_plugin #endif #endif +#if ( defined SM_NOAHMP_3_6 ) +#if ( defined RM_RTM_FORWARD ) + call registerlsmf2t(trim(LIS_noahmp36Id)//"+"& + //trim(LIS_RTMforwardId)//char(0),NoahMP36_f2t) +#endif + +#if ( defined RTMS_WCM) + call registerlsm2rtm(trim(LIS_wcmRTMId)//"+"//& + trim(LIS_noahmp36Id)//char(0), noahmp36_sfc2wcm) +#endif +#endif + + #if ( defined SM_MOSAIC ) #if ( defined RM_RTM_FORWARD ) call registerlsmf2t(trim(LIS_mosaicId)//"+"//& diff --git a/lis/plugins/LIS_pluginIndices.F90 b/lis/plugins/LIS_pluginIndices.F90 index 352625ac3..a5e4f6f17 100644 --- a/lis/plugins/LIS_pluginIndices.F90 +++ b/lis/plugins/LIS_pluginIndices.F90 @@ -30,7 +30,8 @@ module LIS_pluginIndices ! 27 Jan 2014: Shugong Wang, added HRAP projection ! 4 Nov 2014: Jonathan Case, added support for daily NESDIS/VIIRS GVF for Noah ! 16 Aug 2016: Mahdi Navari, added PILDAS -! +! 28 Aug 2020: Sara Modanesi, added WCM in Radiative Transfer Models +! 30 Jun 2021: Sara Modanesi, Michel Bechtold, added Sentinel-1 DA cases for SM and LAI updating !EOP PRIVATE @@ -243,6 +244,8 @@ module LIS_pluginIndices character*50, public, parameter :: LIS_ANSASCFsnowobsId = "ANSA SCF" character*50, public, parameter :: LIS_ANSASNWDsnowobsId = "ANSA snow depth" character*50, public, parameter :: LIS_SMMRSNWDsnowobsId = "SMMR snow depth" + character*50, public, parameter :: LIS_S1_sigmaVVSM_obsId = "S1 backscatter VVSM" + character*50, public, parameter :: LIS_S1_sigmaVVVHSMLAI_obsId = "S1 backscatter VVVHSMLAI" character*50, public, parameter :: LIS_SSMISNWDsnowobsId = "SSMI snow depth" character*50, public, parameter :: LIS_AMSREsweobsId = "AMSR-E SWE" ! character*50, public, parameter :: LIS_AMSREsnowobsId = "AMSR-E snow" !yliu @@ -355,6 +358,7 @@ module LIS_pluginIndices character*50, public, parameter :: LIS_crtm2EMId = "CRTM2EM" character*50, public, parameter :: LIS_cmem3Id = "CMEM" character*50, public, parameter :: LIS_tauomegaRTMId = "Tau Omega" + character*50, public, parameter :: LIS_wcmRTMId = "WCM" !------------------------------------------------------------------------- ! Land Slide Models !------------------------------------------------------------------------- diff --git a/lis/rtms/WCMRTM/WCMRTM_Mod.F90 b/lis/rtms/WCMRTM/WCMRTM_Mod.F90 new file mode 100644 index 000000000..bcdc69aba --- /dev/null +++ b/lis/rtms/WCMRTM/WCMRTM_Mod.F90 @@ -0,0 +1,531 @@ +!-----------------------BEGIN--------------------------------------------- +! 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----------------------- +#include "LIS_misc.h" +module WCMRTM_Mod +!BOP +! +! !MODULE: WCMRTM_Mod +! +! !DESCRIPTION: +! This module provides the routines to control the execution of +! the Water Cloud model (Ulaby et al., 1990). The routine allows to compute +! both the backscatter in VV and the backscatter in VH. It is written +! considering a static incidence angle of 37°. This can be improved in cas +! a local incidence angle for each observation is available. Table of +! calibrated parameters (different for the two polarizations are added also +! in the configuration file). Each pixel of the study area needs to have +! assigned A,B,C,D parameters for each polarization, associated with +! longitude and latitude of the same pixel (pixel where the soil moisture is +! equa to nan need to be removed) +! +! !HISTORY: +! 28 Aug 2020: Sara Modanesi +! 26 Mar 2021 Sara Modanesi: Added specifications for forward states +! !USES: + + +#if (defined RTMS) + + use ESMF + use LIS_coreMod + use LIS_RTMMod + use LIS_logMod + + implicit none + + PRIVATE +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: WCMRTM_initialize + public :: WCMRTM_f2t + public :: WCMRTM_run + public :: WCMRTM_output + public :: WCMRTM_geometry +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + public :: wcm_struc +!EOP + type, public :: wcm_type_dec + + character(len=256) :: AA_VV_tbl_name !! from NoahMP36_lsmMod.F90 + character(len=256) :: BB_VV_tbl_name + character(len=256) :: CC_VV_tbl_name + character(len=256) :: DD_VV_tbl_name + + character(len=256) :: AA_VH_tbl_name + character(len=256) :: BB_VH_tbl_name + character(len=256) :: CC_VH_tbl_name + character(len=256) :: DD_VH_tbl_name + + real, allocatable :: AA_VV(:) + real, allocatable :: BB_VV(:) + real, allocatable :: CC_VV(:) + real, allocatable :: DD_VV(:) + + real, allocatable :: AA_VH(:) + real, allocatable :: BB_VH(:) + real, allocatable :: CC_VH(:) + real, allocatable :: DD_VH(:) + + real, allocatable :: lone(:) + real, allocatable :: late(:) + !-------output------------! + real, allocatable :: Sig0VV(:) + real, allocatable :: Sig0VH(:) + end type wcm_type_dec + + type(wcm_type_dec), allocatable :: wcm_struc(:) + + SAVE + +contains +!BOP +! +! !ROUTINE: WCMRTM_initialize +! \label{WCMRTM_initialize} +! +! !INTERFACE: + subroutine WCMRTM_initialize() +! !USES: + +! !DESCRIPTION: +! +! This routine creates the datatypes and allocates memory for noahMP3.6-specific +! variables. It also invokes the routine to read the runtime specific options +! for noahMP3.6 from the configuration file. +! +! The routines invoked are: +! \begin{description} +! \item[readWCMRTMcrd](\ref{readWCMRTMcrd}) \newline +! +!EOP + implicit none + + integer :: rc + integer :: n,t + integer :: ierr + integer , parameter :: OPEN_OK = 0 + character*128 :: message + +!allocate memory for nest + allocate(wcm_struc(LIS_rc%nnest)) + + do n=1,LIS_rc%nnest +!allocate memory for all tile in current nest + + + allocate(wcm_struc(n)%AA_VV(LIS_rc%glbngrid(n))) + allocate(wcm_struc(n)%BB_VV(LIS_rc%glbngrid(n))) + allocate(wcm_struc(n)%CC_VV(LIS_rc%glbngrid(n))) + allocate(wcm_struc(n)%DD_VV(LIS_rc%glbngrid(n))) + + allocate(wcm_struc(n)%AA_VH(LIS_rc%glbngrid(n))) + allocate(wcm_struc(n)%BB_VH(LIS_rc%glbngrid(n))) + allocate(wcm_struc(n)%CC_VH(LIS_rc%glbngrid(n))) + allocate(wcm_struc(n)%DD_VH(LIS_rc%glbngrid(n))) + + allocate(wcm_struc(n)%lone(LIS_rc%glbngrid(n))) + allocate(wcm_struc(n)%late(LIS_rc%glbngrid(n))) + + allocate(wcm_struc(n)%Sig0VV(LIS_rc%npatch(n,LIS_rc%lsm_index))) + allocate(wcm_struc(n)%Sig0VH(LIS_rc%npatch(n,LIS_rc%lsm_index))) + + call add_sfc_fields(n,LIS_sfcState(n), "Soil Moisture Content") + call add_sfc_fields(n,LIS_sfcState(n), "Leaf Area Index") + + enddo + +!----------------A,B,C and D parameter tables for VV pol---------------------! + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM AA_VV parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%AA_VV_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM AA_VV parameter table: not defined") + enddo + + + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM BB_VV parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%BB_VV_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM BB_VV parameter table: not defined") + enddo + + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM CC_VV parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%CC_VV_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM CC_VV parameter table: not defined") + enddo + + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM DD_VV parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%DD_VV_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM DD_VV parameter table: not defined") + enddo + + + +!----------------A,B,C and D parameter tables for VH pol---------------------! + + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM AA_VH parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%AA_VH_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM AA_VH parameter table: not defined") + enddo + + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM BB_VH parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%BB_VH_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM BB_VH parameter table: not defined") + enddo + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM CC_VH parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%CC_VH_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM CC_VH parameter table: not defined") + enddo + call ESMF_ConfigFindLabel(LIS_config, "WCMRTM DD_VH parameter table:",rc = rc) + do n=1, LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config,wcm_struc(n)%DD_VH_tbl_name, rc=rc) + call LIS_verify(rc, "WCMRTM DD_VH parameter table: not defined") + enddo +!------------------------------------------------------------------------- +!!!!READ IN A PARAMETER TABLE for VV pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%AA_VV_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening AA_VV_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*)wcm_struc(n)%AA_VV(t),wcm_struc(n)%lone(t),& + wcm_struc(n)%late(t) + enddo + CLOSE (19) + enddo +!---------------------------------------------------------------------------- +!-----READ IN B PARAMETER TABLE for VV pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%BB_VV_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening BB_VV_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*)wcm_struc(n)%BB_VV(t) + enddo + + ! READ (19,*) B_val + CLOSE (19) + enddo +!----------------------------------------------------------------------------- +!-----READ IN C PARAMETER TABLE for VV pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%CC_VV_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening CC_VV_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*) wcm_struc(n)%CC_VV(t) + enddo + +! READ (19,*) C_val + CLOSE (19) + enddo +!------------------------------------------------------------------------------- +!-----READ IN D PARAMETER TABLE for VV pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%DD_VV_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening DD_VV_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*) wcm_struc(n)%DD_VV(t) + enddo + +! READ (19,*) D_val + CLOSE (19) + enddo + + +!------------------------------------------------------------------------- +!!!!READ IN A PARAMETER TABLE for VH pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%AA_VH_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening AA_VH_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*)wcm_struc(n)%AA_VH(t) + enddo + CLOSE (19) + enddo +!---------------------------------------------------------------------------- +!-----READ IN B PARAMETER TABLE for VH pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%BB_VH_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening BB_VH_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*)wcm_struc(n)%BB_VH(t) + enddo + + ! READ (19,*) B_val + CLOSE (19) + enddo +!----------------------------------------------------------------------------- +!-----READ IN C PARAMETER TABLE for VH pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%CC_VH_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening CC_VH_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*) wcm_struc(n)%CC_VH(t) + enddo + +! READ (19,*) C_val + CLOSE (19) + enddo +!------------------------------------------------------------------------------- +!-----READ IN D PARAMETER TABLE for VH pol + do n=1, LIS_rc%nnest + OPEN(19, FILE=trim(wcm_struc(n)%DD_VH_tbl_name),FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) + IF(ierr .NE. OPEN_OK ) THEN + WRITE(message,FMT='(A)') & + 'failure opening DD_VH_PARM.TBL' + CALL wrf_error_fatal ( message ) + END IF + do t=1,LIS_rc%glbngrid(n) + READ (19,*) wcm_struc(n)%DD_VH(t) + enddo + +! READ (19,*) D_val + CLOSE (19) + enddo + + do n=1,LIS_rc%nnest !added fields to State 26032021 + call add_fields_toState(n,LIS_forwardState(n),"WCM_Sig0VV") + call add_fields_toState(n,LIS_forwardState(n),"WCM_Sig0VH") + enddo + + end subroutine WCMRTM_initialize + + subroutine add_fields_toState(n, inState,varname) !added subroutine add-fields_toState 26032021 + + use LIS_logMod, only : LIS_verify + use LIS_coreMod, only : LIS_vecTile + + implicit none + + integer :: n + type(ESMF_State) :: inState + character(len=*) :: varname + + type(ESMF_Field) :: varField + type(ESMF_ArraySpec) :: arrspec + integer :: status + real :: sum + call ESMF_ArraySpecSet(arrspec,rank=1,typekind=ESMF_TYPEKIND_R4,& + rc=status) + call LIS_verify(status) + + varField = ESMF_FieldCreate(arrayspec=arrSpec, & + grid=LIS_vecTile(n), name=trim(varname), & + rc=status) + call LIS_verify(status, 'Error in field_create of '//trim(varname)) + + call ESMF_StateAdd(inState, (/varField/), rc=status) + call LIS_verify(status, 'Error in StateAdd of '//trim(varname)) + + end subroutine add_fields_toState +!!-------------------------------------------------------------------------------- + + subroutine add_sfc_fields(n, sfcState,varname) + + implicit none + + integer :: n + type(ESMF_State) :: sfcState + character(len=*) :: varname + + type(ESMF_Field) :: varField + type(ESMF_ArraySpec) :: arrspec + integer :: status + real :: sum + call ESMF_ArraySpecSet(arrspec,rank=1,typekind=ESMF_TYPEKIND_R4,& + rc=status) + call LIS_verify(status) + + varField = ESMF_FieldCreate(arrayspec=arrSpec, & + grid=LIS_vecTile(n), name=trim(varname), & + rc=status) + call LIS_verify(status, 'Error in field_create of '//trim(varname)) + + call ESMF_StateAdd(sfcState, (/varField/), rc=status) + call LIS_verify(status, 'Error in StateAdd of '//trim(varname)) + + end subroutine add_sfc_fields + + + subroutine WCMRTM_f2t(n) + + implicit none + + integer, intent(in) :: n + + end subroutine WCMRTM_f2t + + + subroutine WCMRTM_geometry(n) + implicit none + integer, intent(in) :: n + + end subroutine WCMRTM_geometry + !Do nothing for now + subroutine WCMRTM_run(n) + use LIS_histDataMod +! !USES: + implicit none + + integer, intent(in) :: n + + integer :: t,p + integer :: status + integer :: col,row + real :: A_VV_cal,B_VV_cal,C_VV_cal,D_VV_cal,A_VH_cal,& + B_VH_cal, C_VH_cal, D_VH_cal,lon,lat,lon1,lat1 + real, pointer :: sm(:), lai(:) + real :: sigmabare_VV,sigmabare_VH,s0VV_s_db, s0VH_s_dB, & + sigmacan_VV, sigmacan_VH,sigmasoil_VV,sigmasoil_VH,& + tt_VV, tt_VH + + real :: theta, ctheta + real, pointer :: sig0val(:) !added for forward states 26032021 + + + theta = 0. !incidence angle in radians (i.e., rad(37° for backscatter or rad(0°) for gamma0) + ctheta = cos(theta) + + +! map surface properties to SFC + call getsfcvar(LIS_sfcState(n), "Soil Moisture Content",& + sm) + call getsfcvar(LIS_sfcState(n), "Leaf Area Index", & + lai) + +!--------------------------------------------- +! Tile loop +!-------------------------------------------- + do t=1, LIS_rc%npatch(n,LIS_rc%lsm_index) + row = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%row + col = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%col + lat = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lat + lon = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lon + do p=1,LIS_rc%glbngrid(n) + lon1= wcm_struc(n)%lone(p) + lat1= wcm_struc(n)%late(p) + if (lon1 .eq. lon .and. lat1 .eq. lat) then + A_VV_cal=wcm_struc(n)%AA_VV(p) + B_VV_cal=wcm_struc(n)%BB_VV(p) + C_VV_cal=wcm_struc(n)%CC_VV(p) + D_VV_cal=wcm_struc(n)%DD_VV(p) + A_VH_cal=wcm_struc(n)%AA_VH(p) + B_VH_cal=wcm_struc(n)%BB_VH(p) + C_VH_cal=wcm_struc(n)%CC_VH(p) + D_VH_cal=wcm_struc(n)%DD_VH(p) + endif + enddo + + if(.not.isNaN(sm(t)).and. sm(t).ne.LIS_rc%udef) then + !bare soil backscatter in db + s0VV_s_db=C_VV_cal+D_VV_cal*sm(t) + s0VH_s_db=C_VH_cal+D_VH_cal*sm(t) + !bare soil backscatter in linear units + sigmabare_VV=10.**(s0VV_s_db/10.) + sigmabare_VH=10.**(s0VH_s_db/10.) + !attenuation + tt_VV=exp(-2.*B_VV_cal*lai(t)/ctheta) + tt_VH=exp(-2.*B_VH_cal*lai(t)/ctheta) + !attenuated soil backscatter + sigmasoil_VV=tt_VV*sigmabare_VV + sigmasoil_VH=tt_VH*sigmabare_VH + !vegetation backscatter in linear units + sigmacan_VV=(1.-tt_VV)*ctheta*(A_VV_cal*lai(t)) + sigmacan_VH=(1.-tt_VH)*ctheta*(A_VH_cal*lai(t)) + !total backscatter + wcm_struc(n)%Sig0VV(t)=10.*log10(sigmacan_VV+sigmasoil_VV) + wcm_struc(n)%Sig0VH(t)=10.*log10(sigmacan_VH+sigmasoil_VH) + else + wcm_struc(n)%Sig0VV(t)=LIS_rc%udef + wcm_struc(n)%Sig0VH(t)=LIS_rc%udef + + endif + + call LIS_diagnoseRTMOutputVar(n, t,LIS_MOC_RTM_Sig0VV,value= & + wcm_struc(n)%Sig0VV(t), & + vlevel=1, unit="dB",direction="-") + + call LIS_diagnoseRTMOutputVar(n, t,LIS_MOC_RTM_Sig0VH,value= & + wcm_struc(n)%Sig0VH(t), & + vlevel=1, unit="dB",direction="-") + enddo + + call getsfcvar(LIS_forwardState(n), "WCM_Sig0VV", sig0val) !added for forward states 26032021 + sig0val = wcm_struc(n)%Sig0VV + + call getsfcvar(LIS_forwardState(n),"WCM_Sig0VH", sig0val) + sig0val = wcm_struc(n)%Sig0VH + + end subroutine WCMRTM_run + + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + subroutine getsfcvar(sfcState, varname, var) +! !USES: + + implicit none + + type(ESMF_State) :: sfcState + type(ESMF_Field) :: varField + character(len=*) :: varname + real, pointer :: var(:) + integer :: status + + call ESMF_StateGet(sfcState, trim(varname), varField, rc=status) + call LIS_verify(status, 'Error in StateGet: CMEM3_handlerMod '//trim(varname)) + call ESMF_FieldGet(varField, localDE=0,farrayPtr=var, rc=status) + call LIS_verify(status, 'Error in FieldGet: CMEM3_handlerMod '//trim(varname)) + + end subroutine getsfcvar + +!!!!BOP +!!!! !ROUTINE: WCMRTM_output +!!!! \label{WCMRTM_output} +!!!! +!!!! !INTERFACE: + subroutine WCMRTM_output(n) + integer, intent(in) :: n + end subroutine WCMRTM_output +#endif +end module WCMRTM_Mod + + + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_sig0VV/noahmp36_getSig0VVpred.F90 b/lis/surfacemodels/land/noahmp.3.6/da_sig0VV/noahmp36_getSig0VVpred.F90 new file mode 100644 index 000000000..893505683 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_sig0VV/noahmp36_getSig0VVpred.F90 @@ -0,0 +1,88 @@ +!-----------------------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: noahmp36_getSig0VVpred +! \label{noahmp36_getSig0VVpred} +! +! !REVISION HISTORY: +! 13 May 2021: Sara Modanesi, Michel Bechtold; Initiated specifications for VV polarization +! 18 June 2021: Michel Bechtold; bug fix for obs_ngrid that was problematic for for multiple obs +! +! !INTERFACE: + +subroutine noahmp36_getSig0VVpred(n,k,obs_pred) + +! !USES: + + use ESMF + use LIS_coreMod, only : LIS_rc,LIS_surface + use LIS_RTMMod, only : LIS_forwardState, LIS_RTM_run + use LIS_logMod, only : LIS_verify + + use noahmp36_lsmMod + +!EOP + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + real :: obs_pred(LIS_rc%obs_ngrid(k),LIS_rc%nensem(n)) + integer :: count1(LIS_rc%obs_ngrid(k),LIS_rc%nensem(n)) +! +! !DESCRIPTION: +! +! Returns the Backscatter VV obs pred (model's estimate of +! observations) for data assimilation +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[obs\_pred] model's estimate of observations \newline +! \end{description} +!EOP + + integer :: i,t,m,gid + real, pointer :: s0VV(:) + type(ESMF_Field) :: varField + integer :: status + +!!call the forward model +call LIS_RTM_run(n) + + call ESMF_StateGet(LIS_forwardState(n), "WCM_Sig0VV", varField, rc=status) +!LIS_histDAtaMod + call LIS_verify(status, & + "Error in StateGet in noahmp36_getSig0VVPred for WCM_Sig0VV") + + call ESMF_FieldGet(varField, localDE=0,farrayPtr=s0VV, rc=status)!variable +!s0VV previously defined + call LIS_verify(status, & + 'Error in FieldGet in noahmp36_getSig0VVPred for WCM_Sig0VV') + + obs_pred = 0.0 + count1 = 0.0 + do i=1,LIS_rc%npatch(n,LIS_rc%lsm_index),LIS_rc%nensem(n) + do m=1,LIS_rc%nensem(n) + t = i+m-1 + gid = LIS_surface(n,LIS_rc%lsm_index)%tile(t)%index + obs_pred(gid,m)= obs_pred(gid,m) + s0VV(t) + count1(gid,m) = count1(gid,m) +1 + enddo + enddo + + do i=1,LIS_rc%obs_ngrid(k) + do m=1,LIS_rc%nensem(n) + obs_pred(i,m) = obs_pred(i,m)/(count1(i,m)) + enddo + enddo + +end subroutine noahmp36_getSig0VVpred + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_sig0VV/noahmp36_qc_Sig0VVobs.F90 b/lis/surfacemodels/land/noahmp.3.6/da_sig0VV/noahmp36_qc_Sig0VVobs.F90 new file mode 100644 index 000000000..0f191d4e1 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_sig0VV/noahmp36_qc_Sig0VVobs.F90 @@ -0,0 +1,260 @@ +!-----------------------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: noahmp36_qc_Sig0VVobs +! \label{noahmp36_qc_Sig0VVobs} +! +! !REVISION HISTORY: +! 26/03/2021 Sara Modanesi: Initial specifications +! 12/05/2021 Sara Modanesi: added specifications for Sig0VV S1 obs only and removed flag for veg. cover +! 18/06/2021 Michel Bechtold: assimilation flag for urban and water tiles +! +! !INTERFACE: +subroutine noahmp36_qc_Sig0VVobs(n,k,OBS_State) +! !USES: + use ESMF + use LIS_coreMod + use LIS_logMod, only : LIS_verify + use LIS_constantsMod, only : LIS_CONST_TKFRZ + use LIS_DAobservationsMod + use noahmp36_lsmMod + use module_sf_noahlsm_36 !, only: MAXSMC !MN + + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + type(ESMF_State) :: OBS_State +! +! !DESCRIPTION: +! +! This subroutine performs any model-based QC of the observation +! prior to data assimilation. Here the backscatter observations +! are flagged when LSM indicates that (1) rain is falling (2) +! soil is frozen or (3) ground is fully or partially covered +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[OBS\_State] ESMF state container for observations \newline +! \end{description} +! +!EOP + type(ESMF_Field) :: sigmaField + + real, pointer :: obsl(:) + integer :: t + integer :: gid + integer :: status + real :: lat,lon + +! mn + integer :: SOILTYP ! soil type index [-] + real :: smc1(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: smc2(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: smc3(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: smc4(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o1(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o2(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o3(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o4(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc1(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc2(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc3(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc4(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: vegt(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: SMCMAX(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: SMCWLT(LIS_rc%npatch(n,LIS_rc%lsm_index)) + + real :: rainf_obs(LIS_rc%obs_ngrid(k)) + real :: sneqv_obs(LIS_rc%obs_ngrid(k)) + real :: sca_obs(LIS_rc%obs_ngrid(k)) + real :: t1_obs(LIS_rc%obs_ngrid(k)) + real :: smcwlt_obs(LIS_rc%obs_ngrid(k)) + real :: smcmax_obs(LIS_rc%obs_ngrid(k)) + real :: smc1_obs(LIS_rc%obs_ngrid(k)) + real :: smc2_obs(LIS_rc%obs_ngrid(k)) + real :: smc3_obs(LIS_rc%obs_ngrid(k)) + real :: smc4_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o1_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o2_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o3_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o4_obs(LIS_rc%obs_ngrid(k)) + real :: stc1_obs(LIS_rc%obs_ngrid(k)) + real :: stc2_obs(LIS_rc%obs_ngrid(k)) + real :: stc3_obs(LIS_rc%obs_ngrid(k)) + real :: stc4_obs(LIS_rc%obs_ngrid(k)) + real :: vegt_obs(LIS_rc%obs_ngrid(k)) + +!-----this part is derived from ./lis/dataassim/obs/s1_sigma/read_S1_sigma.F90 + call ESMF_StateGet(OBS_State,"Observation01",sigmaField,& + rc=status) ! + call LIS_verify(status,& + "ESMF_StateGet failed in noahmp36_qc_Sig0VVobs sigmaField") + + call ESMF_FieldGet(sigmaField,localDE=0,farrayPtr=obsl,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet failed in noahmp36_qc_Sigobs obsl") + +!--------------------------------------------------------------------------- + + do t=1, LIS_rc%npatch(n,LIS_rc%lsm_index) + smc1(t) = noahmp36_struc(n)%noahmp36(t)%smc(1) + smc2(t) = noahmp36_struc(n)%noahmp36(t)%smc(2) + smc3(t) = noahmp36_struc(n)%noahmp36(t)%smc(3) + smc4(t) = noahmp36_struc(n)%noahmp36(t)%smc(4) + + sh2o1(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(1) + sh2o2(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(2) + sh2o3(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(3) + sh2o4(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(4) +!--------------------------------------------------------------------------------------------------------- + ! MN NOTE:sstc contain soil and snow temperature first snow + ! temperature and then soil temeprature. + ! But the number of snow layers changes from 0 to 3 +!--------------------------------------------------------------------------------------------------------- + stc1(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+1) + stc2(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+2) + stc3(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+3) + stc4(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+4) + + vegt(t) = noahmp36_struc(n)%noahmp36(t)%vegetype + + SOILTYP = NOAHMP36_struc(n)%noahmp36(t)%soiltype + SMCMAX(t) = MAXSMC (SOILTYP) + SMCWLT(t) = WLTSMC (SOILTYP) + enddo + + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%prcp,& + rainf_obs) ! MN prcp is total precip + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%sneqv,& + sneqv_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%fsno,& + sca_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%tg,& + t1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smcmax, & + smcmax_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smcwlt,& + smcwlt_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc1,& + smc1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc2,& + smc2_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc3,& + smc3_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc4,& + smc4_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o1,& + sh2o1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o2,& + sh2o2_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o3,& + sh2o3_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o4,& + sh2o4_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc1,& + stc1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc2,& + stc2_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc3,& + stc3_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc4,& + stc4_obs) + + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + vegt,& + vegt_obs) + + do t = 1,LIS_rc%obs_ngrid(k) +!------------------start loop considering one obs-------------------------- + if(obsl(t).ne.LIS_rc%udef) then +! MN: check for rain + if(rainf_obs(t).gt.3E-6) then ! Var name Noah36 --> rainf + obsl(t) = LIS_rc%udef +! print*, 'rainf ',gid,t,NOAHMP36_struc(n)%noahmp36(t)%prcp +! MN: check for frozen soil + elseif(abs(smc1_obs(t)- & + sh2o1_obs(t)).gt.0.0001) then + obsl(t) = LIS_rc%udef + elseif(abs(smc2_obs(t)- & + sh2o2_obs(t)).gt.0.0001) then + obsl(t) = LIS_rc%udef + elseif(abs(smc3_obs(t)- & + sh2o3_obs(t)).gt.0.0001) then + obsl(t) = LIS_rc%udef + elseif(abs(smc4_obs(t)- & + sh2o4_obs(t)).gt.0.0001) then + obsl(t) = LIS_rc%udef + elseif(stc1_obs(t).le.LIS_CONST_TKFRZ) then + obsl(t) = LIS_rc%udef + elseif(stc2_obs(t).le.LIS_CONST_TKFRZ) then + obsl(t) = LIS_rc%udef + elseif(stc3_obs(t).le.LIS_CONST_TKFRZ) then + obsl(t) = LIS_rc%udef + elseif(stc4_obs(t).le.LIS_CONST_TKFRZ) then + obsl(t) = LIS_rc%udef + elseif(t1_obs(t).le.LIS_CONST_TKFRZ) then ! Var name Noah36 --> t1 + obsl(t) = LIS_rc%udef +! elseif(vegt_obs(t).le.4) then !forest types ! Var name Noah36 --> vegt +! obsl(t) = LIS_rc%udef + elseif(vegt_obs(t).eq.13) then !urban ! Var name Noah36 --> vegt + obsl(t) = LIS_rc%udef + elseif(vegt_obs(t).eq.17) then !urban ! Var name Noah36 --> vegt + obsl(t) = LIS_rc%udef + ! MN: check for snow + elseif(sneqv_obs(t).gt.0.001) then + obsl(t) = LIS_rc%udef + elseif(sca_obs(t).gt.0.0001) then ! Var name sca + obsl(t) = LIS_rc%udef + endif + endif + enddo + +end subroutine noahmp36_qc_Sig0VVobs + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_sig0VVVH/noahmp36_getSig0VVVHpred.F90 b/lis/surfacemodels/land/noahmp.3.6/da_sig0VVVH/noahmp36_getSig0VVVHpred.F90 new file mode 100644 index 000000000..1ee040bfb --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_sig0VVVH/noahmp36_getSig0VVVHpred.F90 @@ -0,0 +1,127 @@ +!-----------------------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: noahmp36_getSig0VVVHpred +! \label{noahmp36_getSig0VVVHpred} +! +! !REVISION HISTORY: +! 10 Mar 2022: Michel Bechtold; module for joint assim of VV and VH +! +! !INTERFACE: + +subroutine noahmp36_getSig0VVVHpred(n,k,obs_pred) + +! !USES: + + use ESMF + !use LIS_coreMod, only : LIS_rc,LIS_surface + use LIS_RTMMod, only : LIS_forwardState, LIS_RTM_run + !use LIS_logMod, only : LIS_verify + + use noahmp36_lsmMod + + ! added + use LIS_constantsMod + use LIS_coreMod + use LIS_logMod + use LIS_dataAssimMod + use LIS_DAobservationsMod + +!EOP + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + real :: obs_pred(LIS_rc%obs_ngrid(k)*2,LIS_rc%nensem(n)) + real :: obs_pred_VV(LIS_rc%obs_ngrid(k),LIS_rc%nensem(n)) + real :: obs_pred_VH(LIS_rc%obs_ngrid(k),LIS_rc%nensem(n)) + integer :: count1(LIS_rc%obs_ngrid(k)*2,LIS_rc%nensem(n)) +! +! !DESCRIPTION: +! +! Returns the Backscatter VV and VH obs pred (model's estimate of +! observations) for data assimilation +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[obs\_pred] model's estimate of observations \newline +! \end{description} +!EOP + + integer :: i,t,m,gid + real, pointer :: s0VV(:), s0VH(:) + type(ESMF_Field) :: varField + integer :: status + +!!call the forward model +call LIS_RTM_run(n) + + call ESMF_StateGet(LIS_forwardState(n), "WCM_Sig0VV", varField, rc=status) +!LIS_histDAtaMod + call LIS_verify(status, & + "Error in StateGet in noahmp36_getSig0Pred for WCM_Sig0VV") + + call ESMF_FieldGet(varField, localDE=0,farrayPtr=s0VV, rc=status)!variable +!s0VV previously defined + call LIS_verify(status, & + 'Error in FieldGet in noahmp36_getSig0Pred for WCM_Sig0VV') + + call ESMF_StateGet(LIS_forwardState(n), "WCM_Sig0VH", varField, rc=status) + call LIS_verify(status, & + "Error in StateGet in noahmp36_getSig0Pred for WCM_Sig0VH") + + call ESMF_FieldGet(varField, localDE=0,farrayPtr=s0VH, rc=status) + call LIS_verify(status, & + 'Error in FieldGet in noahmp36_getSig0Pred for WCM_Sig0VH') + obs_pred = 0.0 + count1 = 0.0 + do i=1,LIS_rc%npatch(n,LIS_rc%lsm_index),LIS_rc%nensem(n) + do m=1,LIS_rc%nensem(n) + t = i+m-1 + gid = LIS_surface(n,LIS_rc%lsm_index)%tile(t)%index + obs_pred(gid,m)= obs_pred(gid,m) + s0VV(t) + count1(gid,m) = count1(gid,m) +1 + enddo + enddo + + !!! new + + + call LIS_convertPatchSpaceToObsEnsSpace(n,k,& + LIS_rc%lsm_index, & + s0VV,& + obs_pred_VV) + + call LIS_convertPatchSpaceToObsEnsSpace(n,k,& + LIS_rc%lsm_index, & + s0VH,& + obs_pred_VH) + + + do i=1,LIS_rc%npatch(n,LIS_rc%lsm_index),LIS_rc%nensem(n) + do m=1,LIS_rc%nensem(n) + t = i+m-1 + gid = LIS_surface(n,LIS_rc%lsm_index)%tile(t)%index + & + LIS_rc%ngrid(n) + obs_pred(gid,m)= obs_pred(gid,m) + s0VH(t) + count1(gid,m) = count1(gid,m) +1 + enddo + enddo + + do i=1,LIS_rc%obs_ngrid(k)*2 + do m=1,LIS_rc%nensem(n) + obs_pred(i,m) = obs_pred(i,m)/(count1(i,m)) + enddo + enddo + +end subroutine noahmp36_getSig0VVVHpred + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_sig0VVVH/noahmp36_qc_Sig0VVVHobs.F90 b/lis/surfacemodels/land/noahmp.3.6/da_sig0VVVH/noahmp36_qc_Sig0VVVHobs.F90 new file mode 100644 index 000000000..a4657b16f --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_sig0VVVH/noahmp36_qc_Sig0VVVHobs.F90 @@ -0,0 +1,304 @@ +!-----------------------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: noahmp36_qc_Sig0obs +! \label{noahmp36_qc_Sig0obs} +! +! !REVISION HISTORY: +! 26/03/2021 Sara Modanesi: Initial specifications +! 02/04/2021 Sara Modanesi: added specifications for both Sig0VV and Sig0VH S1 obs and removed flag for veg. cover +! 10/03/2022 Michel Bechtold: module for joint assim of VV and VH +! +! !INTERFACE: +subroutine noahmp36_qc_Sig0VVVHobs(n,k,OBS_State) +! !USES: + use ESMF + use LIS_coreMod + use LIS_logMod, only : LIS_verify + use LIS_constantsMod, only : LIS_CONST_TKFRZ + use LIS_DAobservationsMod + use noahmp36_lsmMod + use module_sf_noahlsm_36 !, only: MAXSMC !MN + + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: k + type(ESMF_State) :: OBS_State +! +! !DESCRIPTION: +! +! This subroutine performs any model-based QC of the observation +! prior to data assimilation. Here the backscatter observations +! are flagged when LSM indicates that (1) rain is falling (2) +! soil is frozen or (3) ground is fully or partially covered with snow +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[OBS\_State] ESMF state container for observations \newline +! \end{description} +! +!EOP + type(ESMF_Field) :: s_vvField,s_vhField + + real, pointer :: s_vv(:),s_vh(:) + integer :: t + integer :: gid + integer :: status + real :: lat,lon + +! mn + integer :: SOILTYP ! soil type index [-] + real :: smc1(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: smc2(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: smc3(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: smc4(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o1(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o2(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o3(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: sh2o4(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc1(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc2(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc3(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: stc4(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: vegt(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: SMCMAX(LIS_rc%npatch(n,LIS_rc%lsm_index)) + real :: SMCWLT(LIS_rc%npatch(n,LIS_rc%lsm_index)) + + real :: rainf_obs(LIS_rc%obs_ngrid(k)) + real :: sneqv_obs(LIS_rc%obs_ngrid(k)) + real :: sca_obs(LIS_rc%obs_ngrid(k)) + real :: t1_obs(LIS_rc%obs_ngrid(k)) + real :: smcwlt_obs(LIS_rc%obs_ngrid(k)) + real :: smcmax_obs(LIS_rc%obs_ngrid(k)) + real :: smc1_obs(LIS_rc%obs_ngrid(k)) + real :: smc2_obs(LIS_rc%obs_ngrid(k)) + real :: smc3_obs(LIS_rc%obs_ngrid(k)) + real :: smc4_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o1_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o2_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o3_obs(LIS_rc%obs_ngrid(k)) + real :: sh2o4_obs(LIS_rc%obs_ngrid(k)) + real :: stc1_obs(LIS_rc%obs_ngrid(k)) + real :: stc2_obs(LIS_rc%obs_ngrid(k)) + real :: stc3_obs(LIS_rc%obs_ngrid(k)) + real :: stc4_obs(LIS_rc%obs_ngrid(k)) + real :: vegt_obs(LIS_rc%obs_ngrid(k)) + + call ESMF_StateGet(OBS_State,"Observation01",s_vvField,& + rc=status) ! + call LIS_verify(status,& + "ESMF_StateGet failed in noahmp36_qc_Sig0obs s_vv") + + call ESMF_StateGet(OBS_State,"Observation02",s_vhField,& + rc=status) ! + call LIS_verify(status,& + "ESMF_StateGet failed in noahmp36_qc_Sig0obs s_vh") + + call ESMF_FieldGet(s_vvField,localDE=0,farrayPtr=s_vv,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet failed in noahmp36_qc_Sigobs s_vv") + + call ESMF_FieldGet(s_vhField,localDE=0,farrayPtr=s_vh,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet failed in noahmp36_qc_Sigobs s_vh") + +!--------------------------------------------------------------------------- + + do t=1, LIS_rc%npatch(n,LIS_rc%lsm_index) + smc1(t) = noahmp36_struc(n)%noahmp36(t)%smc(1) + smc2(t) = noahmp36_struc(n)%noahmp36(t)%smc(2) + smc3(t) = noahmp36_struc(n)%noahmp36(t)%smc(3) + smc4(t) = noahmp36_struc(n)%noahmp36(t)%smc(4) + + sh2o1(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(1) + sh2o2(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(2) + sh2o3(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(3) + sh2o4(t) = noahmp36_struc(n)%noahmp36(t)%sh2o(4) +!--------------------------------------------------------------------------------------------------------- + ! MN NOTE:sstc contain soil and snow temperature first snow + ! temperature and then soil temeprature. + ! But the number of snow layers changes from 0 to 3 +!--------------------------------------------------------------------------------------------------------- + stc1(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+1) + stc2(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+2) + stc3(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+3) + stc4(t) = noahmp36_struc(n)%noahmp36(t)%sstc(NOAHMP36_struc(n)%nsnow+4) + + vegt(t) = noahmp36_struc(n)%noahmp36(t)%vegetype + + SOILTYP = NOAHMP36_struc(n)%noahmp36(t)%soiltype + SMCMAX(t) = MAXSMC (SOILTYP) + SMCWLT(t) = WLTSMC (SOILTYP) + enddo + + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%prcp,& + rainf_obs) ! MN prcp is total precip + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%sneqv,& + sneqv_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%fsno,& + sca_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + noahmp36_struc(n)%noahmp36(:)%tg,& + t1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smcmax, & + smcmax_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smcwlt,& + smcwlt_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc1,& + smc1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc2,& + smc2_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc3,& + smc3_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + smc4,& + smc4_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o1,& + sh2o1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o2,& + sh2o2_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o3,& + sh2o3_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + sh2o4,& + sh2o4_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc1,& + stc1_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc2,& + stc2_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc3,& + stc3_obs) + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + stc4,& + stc4_obs) + + call LIS_convertPatchSpaceToObsSpace(n,k,& + LIS_rc%lsm_index, & + vegt,& + vegt_obs) + + do t = 1,LIS_rc%obs_ngrid(k) +!------------------start loop considering s_vv-------------------------- + if(s_vv(t).ne.LIS_rc%udef) then +! MN: check for rain + if(rainf_obs(t).gt.3E-6) then ! Var name Noah36 --> rainf + s_vv(t) = LIS_rc%udef +! MN: check for frozen soil + elseif(abs(smc1_obs(t)- & + sh2o1_obs(t)).gt.0.0001) then + s_vv(t) = LIS_rc%udef + elseif(abs(smc2_obs(t)- & + sh2o2_obs(t)).gt.0.0001) then + s_vv(t) = LIS_rc%udef + elseif(abs(smc3_obs(t)- & + sh2o3_obs(t)).gt.0.0001) then + s_vv(t) = LIS_rc%udef + elseif(abs(smc4_obs(t)- & + sh2o4_obs(t)).gt.0.0001) then + s_vv(t) = LIS_rc%udef + elseif(stc1_obs(t).le.LIS_CONST_TKFRZ) then + s_vv(t) = LIS_rc%udef + elseif(stc2_obs(t).le.LIS_CONST_TKFRZ) then + s_vv(t) = LIS_rc%udef + elseif(stc3_obs(t).le.LIS_CONST_TKFRZ) then + s_vv(t) = LIS_rc%udef + elseif(stc4_obs(t).le.LIS_CONST_TKFRZ) then + s_vv(t) = LIS_rc%udef + elseif(t1_obs(t).le.LIS_CONST_TKFRZ) then ! Var name Noah36 --> t1 + s_vv(t) = LIS_rc%udef + elseif(vegt_obs(t).eq.13) then !urban ! Var name Noah36 --> vegt + s_vv(t) = LIS_rc%udef + elseif(vegt_obs(t).eq.17) then !water ! Var name Noah36 --> vegt + s_vv(t) = LIS_rc%udef + ! MN: check for snow + elseif(sneqv_obs(t).gt.0.001) then + s_vv(t) = LIS_rc%udef + elseif(sca_obs(t).gt.0.0001) then ! Var name sca + s_vv(t) = LIS_rc%udef + endif + endif +!------------------start loop considering s_vh-------------------------- + if(s_vh(t).ne.LIS_rc%udef) then +! MN: check for rain + if(rainf_obs(t).gt.3E-6) then ! Var name Noah36 --> rainf + s_vh(t) = LIS_rc%udef +! MN: check for frozen soil + elseif(abs(smc1_obs(t)- & + sh2o1_obs(t)).gt.0.0001) then + s_vh(t) = LIS_rc%udef + elseif(abs(smc2_obs(t)- & + sh2o2_obs(t)).gt.0.0001) then + s_vh(t) = LIS_rc%udef + elseif(abs(smc3_obs(t)- & + sh2o3_obs(t)).gt.0.0001) then + s_vh(t) = LIS_rc%udef + elseif(abs(smc4_obs(t)- & + sh2o4_obs(t)).gt.0.0001) then + s_vh(t) = LIS_rc%udef + elseif(stc1_obs(t).le.LIS_CONST_TKFRZ) then + s_vh(t) = LIS_rc%udef + elseif(stc2_obs(t).le.LIS_CONST_TKFRZ) then + s_vh(t) = LIS_rc%udef + elseif(stc3_obs(t).le.LIS_CONST_TKFRZ) then + s_vh(t) = LIS_rc%udef + elseif(stc4_obs(t).le.LIS_CONST_TKFRZ) then + s_vh(t) = LIS_rc%udef + elseif(t1_obs(t).le.LIS_CONST_TKFRZ) then ! Var name Noah36 --> t1 + s_vh(t) = LIS_rc%udef + elseif(vegt_obs(t).eq.13) then !urban ! Var name Noah36 --> vegt + s_vh(t) = LIS_rc%udef + elseif(vegt_obs(t).eq.17) then !water ! Var name Noah36 --> vegt + s_vh(t) = LIS_rc%udef + ! MN: check for snow + elseif(sneqv_obs(t).gt.0.001) then + s_vh(t) = LIS_rc%udef + elseif(sca_obs(t).gt.0.0001) then ! Var name sca + s_vh(t) = LIS_rc%udef + endif + endif + enddo + +end subroutine noahmp36_qc_Sig0VVVHobs + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_dasoilmLAI_Mod.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_dasoilmLAI_Mod.F90 new file mode 100644 index 000000000..c11d591ee --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_dasoilmLAI_Mod.F90 @@ -0,0 +1,83 @@ +!-----------------------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" +module noahmp36_dasoilmLAI_Mod +!BOP +! +! !MODULE: noahmp36_dasoilmLAI_Mod +! +! !DESCRIPTION: +! +! !REVISION HISTORY: + +! Sujay Kumar; Initial Code +! 9 Sep 2016: Mahdi Navari; Modified for NoahMP36 ! +! 30 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! !USES: + use ESMF + use LIS_coreMod + use LIS_dataAssimMod + use LIS_logMod + + implicit none + + PRIVATE +!----------------------------------------------------------------------------- +! !PUBLIC MEMBER FUNCTIONS: +!----------------------------------------------------------------------------- + public :: noahmp36_dasoilmLAI_init +!----------------------------------------------------------------------------- +! !PUBLIC TYPES: +!----------------------------------------------------------------------------- + public :: noahmp36_dasm_struc +!EOP + + type, public :: dasm_dec + real, allocatable :: model_xrange(:,:,:) + real, allocatable :: model_cdf(:,:,:) + real, allocatable :: model_mu(:) + + integer :: nbins + integer :: ntimes + integer :: scal + + end type dasm_dec + + type(dasm_dec), allocatable :: noahmp36_dasm_struc(:) + +contains +!BOP +! +! !ROUTINE: noahmp36_dasoilmLAI_init +! \label{noahmp36_dasoilmLAI_init} +! +! !INTERFACE: + subroutine noahmp36_dasoilmLAI_init(k) +! !USES: +! !DESCRIPTION: +! +!EOP + + + implicit none + integer :: k + integer :: n + character*100 :: modelcdffile(LIS_rc%nnest) + integer :: status + integer :: ngrid + + if(.not.allocated(noahmp36_dasm_struc)) then + allocate(noahmp36_dasm_struc(LIS_rc%nnest)) + endif + + ! scaling not active for SM and LAI updating with S1 DA + + end subroutine noahmp36_dasoilmLAI_init +end module noahmp36_dasoilmLAI_Mod diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_descale_soilmLAI.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_descale_soilmLAI.F90 new file mode 100644 index 000000000..7179a945a --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_descale_soilmLAI.F90 @@ -0,0 +1,50 @@ +!-----------------------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: noahmp36_descale_soilmLAI +! \label{noahmp36_descale_soilmLAI} +! +! !REVISION HISTORY: +! 27Feb2005: Sujay Kumar; Initial Specification +! 25Jun2006: Sujay Kumar: Updated for the ESMF design +! 1 Aug 2016: Mahdi Navari; Modified for NoahMP36 +! 18 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! +! !INTERFACE: +subroutine noahmp36_descale_soilmLAI(n, LSM_State, LSM_Incr_State) + +! !USES: + use ESMF + use LIS_coreMod, only : LIS_rc + use LIS_logMod, only : LIS_verify + use noahmp36_lsmMod + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + type(ESMF_State) :: LSM_State + type(ESMF_State) :: LSM_Incr_State +! +! !DESCRIPTION: +! +! Descales soilmoisture and LAI related state prognostic variables for +! data assimilation +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[LSM\_State] ESMF State container for LSM state variables \newline +! \end{description} +!EOP + +! Scaling not active for SM and LAI updating through S1 data assimilation + +end subroutine noahmp36_descale_soilmLAI + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_getsoilmLAI.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_getsoilmLAI.F90 new file mode 100644 index 000000000..519aecf8b --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_getsoilmLAI.F90 @@ -0,0 +1,92 @@ +!-----------------------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: noahmp36_getsoilmLAI +! \label{noahmp36_getsoilmLAI} +! +! !REVISION HISTORY: +! 27Feb2005: Sujay Kumar; Initial Specification +! 25Jun2006: Sujay Kumar: Updated for the ESMF design +! 1 Aug 2016: Mahdi Navari; Modified for NoahMP36 +! To do: makes it general for x layers (currently hard coded for 4 layers) +! 18 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! !INTERFACE: +subroutine noahmp36_getsoilmLAI(n, LSM_State) + +! !USES: + use ESMF + use LIS_coreMod, only : LIS_rc + use LIS_logMod, only : LIS_verify + use noahmp36_lsmMod + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + type(ESMF_State) :: LSM_State +! +! !DESCRIPTION: +! +! Returns the soilmoisture and LAI related state prognostic variables for +! data assimilation +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[LSM\_State] ESMF State container for LSM state variables \newline +! \end{description} +!EOP + type(ESMF_Field) :: sm1Field + type(ESMF_Field) :: sm2Field + type(ESMF_Field) :: sm3Field + type(ESMF_Field) :: sm4Field + type(ESMF_Field) :: laiField,lfmassField + integer :: t + integer :: status + real, pointer :: soilm1(:) + real, pointer :: soilm2(:) + real, pointer :: soilm3(:) + real, pointer :: soilm4(:) + real, pointer :: lai(:) + character*100 :: lsm_state_objs(5) + + + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 1",sm1Field,rc=status) + call LIS_verify(status,'ESMF_StateGet failed for sm1 in noahmp36_getsoilmLAI') + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 2",sm2Field,rc=status) + call LIS_verify(status,'ESMF_StateGet failed for sm2 in noahmp36_getsoilmLAI') + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 3",sm3Field,rc=status) + call LIS_verify(status,'ESMF_StateGet failed for sm3 in noahmp36_getsoilmLAI') + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 4",sm4Field,rc=status) + call LIS_verify(status,'ESMF_StateGet failed for sm4 in noahmp36_getsoilmLAI') + call ESMF_StateGet(LSM_State,"LAI",laiField,rc=status) + call LIS_verify(status,'ESMF_StateGet failed for LAI in noahmp36_getsoilmLAI') + + call ESMF_FieldGet(sm1Field,localDE=0,farrayPtr=soilm1,rc=status) + call LIS_verify(status,'ESMF_FieldGet failed for sm1 in noahmp36_getsoilmLAI') + call ESMF_FieldGet(sm2Field,localDE=0,farrayPtr=soilm2,rc=status) + call LIS_verify(status,'ESMF_FieldGet failed for sm2 in noahmp36_getsoilmLAI') + call ESMF_FieldGet(sm3Field,localDE=0,farrayPtr=soilm3,rc=status) + call LIS_verify(status,'ESMF_FieldGet failed for sm3 in noahmp36_getsoilmLAI') + call ESMF_FieldGet(sm4Field,localDE=0,farrayPtr=soilm4,rc=status) + call LIS_verify(status,'ESMF_FieldGet failed for sm4 in noahmp36_getsoilmLAI') + call ESMF_FieldGet(laiField,localDE=0,farrayPtr=lai,rc=status) + call LIS_verify(status,'ESMF_FieldGet failed for LAI in noahmp36_getsoilmLAI') + + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + soilm1(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(1) + soilm2(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(2) + soilm3(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(3) + soilm4(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(4) + lai(t) = NOAHMP36_struc(n)%noahmp36(t)%lai + enddo + +end subroutine noahmp36_getsoilmLAI + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_qcsoilmLAI.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_qcsoilmLAI.F90 new file mode 100644 index 000000000..90619a352 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_qcsoilmLAI.F90 @@ -0,0 +1,192 @@ +!-----------------------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: noahmp36_qcsoilmLAIveg +! \label{noahmp36_qcsoilmLAIveg} +! +! !REVISION HISTORY: +! 27Feb2005: Sujay Kumar; Initial Specification +! 25Jun2006: Sujay Kumar: Updated for the ESMF design +! 1 Aug 2016: Mahdi Navari; Modified for NoahMP36 +! 18 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! +! !INTERFACE: +subroutine noahmp36_qcsoilmLAI(n, LSM_State) + +! !USES: + use ESMF + use LIS_coreMod + use LIS_logMod + use noahmp36_lsmMod + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + type(ESMF_State) :: LSM_State +! +! !DESCRIPTION: +! +! Returns the soilmoisture related state prognostic variables for +! data assimilation +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[LSM\_State] ESMF State container for LSM state variables \newline +! \end{description} +!EOP + type(ESMF_Field) :: sm1Field + type(ESMF_Field) :: laiField + integer :: t + integer :: status + real, pointer :: soilm1(:) + real, pointer :: lai(:) + real :: smmax1!,smmax2,smmax3,smmax4 + real :: smmin1!,smmin2,smmin3,smmin4 + real :: laimax + real :: laimin + integer :: gid + real :: laitmp + + logical :: update_flag(LIS_rc%ngrid(n)) + real :: perc_violation(LIS_rc%ngrid(n)) + real :: laimean(LIS_rc%ngrid(n)) + integer :: nlaimean(LIS_rc%ngrid(n)) + integer :: N_ens + real :: state_tmp(LIS_rc%nensem(n)),state_mean + + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 1",sm1Field,rc=status) + call LIS_verify(status,& + "ESMF_StateGet for Soil Moisture Layer 1 failed in noahmp36_qcsoilmLAI") + + call ESMF_FieldGet(sm1Field,localDE=0,farrayPtr=soilm1,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet for Soil Moisture Layer 1 failed in noahmp36_qcsoilmLAI") + + call ESMF_AttributeGet(sm1Field,"Max Value",smmax1,rc=status) + call LIS_verify(status,& + "ESMF_AttributeGet: Max Value failed in noahmp36_qcsoilmLAI") + + call ESMF_AttributeGet(sm1Field,"Min Value",smmin1,rc=status) + call LIS_verify(status,& + "ESMF_AttributeGet: Min Value failed in noahmp36_qcsoilmLAI") + + call ESMF_StateGet(LSM_State,"LAI",laiField,rc=status) + call LIS_verify(status,& + "ESMF_StateGet for LAI failed in noahmp36_qcsoilmLAI") + call ESMF_FieldGet(laiField,localDE=0,farrayPtr=lai,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet for LAI failed in noahmp36_qcsoilmLAI") + + call ESMF_AttributeGet(laiField,"Max Value",laimax,rc=status) + call LIS_verify(status,& + "ESMF_AttributeGet for LAI Max Value failed in noahmp36_qcsoilmLAI") + call ESMF_AttributeGet(laiField,"Min Value",laimin,rc=status) + call LIS_verify(status,& + "ESMF_AttributeGet for LAI Min Value failed in noahmp36_qcsoilmLAI") + + + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + + if(soilm1(t).gt.smmax1) soilm1(t) = smmax1 + if(soilm1(t).lt.smmin1) soilm1(t) = smmin1 + enddo + + update_flag = .true. + perc_violation = 0.0 + laimean = 0.0 + nlaimean = 0 + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%row) + + laitmp = lai(t) + + if(laitmp.lt.laimin.or.laitmp.gt.laimax) then + update_flag(gid) = .false. + perc_violation(gid) = perc_violation(gid) +1 + endif + + enddo + + do gid=1,LIS_rc%ngrid(n) + perc_violation(gid) = perc_violation(gid)/LIS_rc%nensem(n) + enddo + +! For ensembles that are unphysical, compute the +! ensemble average after excluding them. This +! is done only if the majority of the ensemble +! members are good (>60%) + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%row) + if(.not.update_flag(gid)) then + if(perc_violation(gid).lt.0.8) then + if((lai(t).gt.laimin).and.& + (lai(t).lt.laimax)) then + laimean(gid) = laimean(gid) + & + lai(t) + nlaimean(gid) = nlaimean(gid) + 1 + endif + endif + endif + enddo + + do gid=1,LIS_rc%ngrid(n) + if(nlaimean(gid).gt.0) then + laimean(gid) = laimean(gid)/nlaimean(gid) + endif + enddo + + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%row) + + laitmp = lai(t) + +! If the update is unphysical, simply set to the average of +! the good ensemble members. If all else fails, do not +! update. + + if(update_flag(gid)) then + lai(t) = laitmp + elseif(perc_violation(gid).lt.0.8) then + if(laitmp.lt.laimin.or.laitmp.gt.laimax) then + lai(t) = laimean(gid) + else + lai(t) = lai(t) + endif + endif + enddo + +#if 0 + N_ens = LIS_rc%nensem(n) + do t=1,N_ens + state_tmp(t) = lai(t) + enddo + state_mean =sum(state_tmp)/N_ens + + write(113,fmt='(i4.4,i2.2,i2.2,i2.2,i2.2,i2.2,21F8.3)') & + LIS_rc%yr, LIS_rc%mo, LIS_rc%da, LIS_rc%hr, & + LIS_rc%mn, LIS_rc%ss, & + state_mean, & + state_tmp +#endif +end subroutine noahmp36_qcsoilmLAI + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_scale_soilmLAI.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_scale_soilmLAI.F90 new file mode 100644 index 000000000..d1ff47911 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_scale_soilmLAI.F90 @@ -0,0 +1,49 @@ +!-----------------------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: noahmp36_scale_soilmLAI +! \label{noahmp36_scale_soilmLAI} +! +! !REVISION HISTORY: +! 27Feb2005: Sujay Kumar; Initial Specification +! 25Jun2006: Sujay Kumar: Updated for the ESMF design +! 1 Aug 2016: Mahdi Navari; Modified for NoahMP36 +! 18 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! +! !INTERFACE: +subroutine noahmp36_scale_soilmLAI(n, LSM_State) + +! !USES: + use ESMF + use LIS_coreMod, only : LIS_rc + use LIS_logMod, only : LIS_verify + use noahmp36_lsmMod + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + type(ESMF_State) :: LSM_State +! +! !DESCRIPTION: +! +! Scales soilmoisture and LAI related state prognostic variables for +! data assimilation +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[LSM\_State] ESMF State container for LSM state variables \newline +! \end{description} +!EOP + +! Scaling not active for SM and LAI updating through S1 obs DA + +end subroutine noahmp36_scale_soilmLAI + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_setsoilmLAI.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_setsoilmLAI.F90 new file mode 100644 index 000000000..118cbd795 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_setsoilmLAI.F90 @@ -0,0 +1,531 @@ +!-----------------------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: noahmp36_setsoilmLAI +! \label{noahmp36_setsoilmLAI} +! +! !REVISION HISTORY: +! 27Feb2005: Sujay Kumar; Initial Specification +! 25Jun2006: Sujay Kumar: Updated for the ESMF design +! 9 Sep 2016: Mahdi Navari: Modified for NoahMP36 +! 18Apr 2018: Mahdi Navari: Bug fixed +! 18 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! +! Apply the update if it met the update conditions +! Update conditions: +! 1- Prior SM(sh2o) + increment > MIN_THRESHOLD +! 2- Prior SM(sh2o) + increment < sm_threshold +! There are 3 cases +! 1- If all the ensemble members met the update conditions --> apply the update +! 2- If more than 50% of the ensemble members met the update condition --> +! apply the update for that members and set the other member to the mean +! value of the ensemble (i.e. mean of the members that met the conditions) +! 3- If less then 50% of the ensemble members met the update conditions --> +! adjust the states + + +! !INTERFACE: +subroutine noahmp36_setsoilmLAI(n, LSM_State) +! !USES: + use ESMF + use LIS_coreMod + use LIS_logMod + use noahmp36_lsmMod + use module_sf_noahlsm_36 !MN + use NOAHMP_VEG_PARAMETERS_36 + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + type(ESMF_State) :: LSM_State +! +! !DESCRIPTION: +! +! This routine assigns the soil moisture and LAI prognostic variables to noah's +! model space. +! +!EOP + real, parameter :: MIN_THRESHOLD = 0.02 + real :: MAX_threshold + real :: sm_threshold + type(ESMF_Field) :: sm1Field + type(ESMF_Field) :: sm2Field + type(ESMF_Field) :: sm3Field + type(ESMF_Field) :: sm4Field + type(ESMF_Field) :: laiField,lfmassField + real, pointer :: soilm1(:) + real, pointer :: soilm2(:) + real, pointer :: soilm3(:) + real, pointer :: soilm4(:) + real, pointer :: lai(:) + real :: lfmass + integer :: t, j,i, gid, m, t_unpert + integer :: status + real :: delta(4) + real :: delta1,delta2,delta3,delta4 + real :: tmpval + logical :: bounds_violation + integer :: nIter + logical :: update_flag(LIS_rc%ngrid(n)) + logical :: ens_flag(LIS_rc%nensem(n)) +! mn + integer :: SOILTYP ! soil type index [-] + real :: SMCMAX , SMCWLT + real :: tmp(LIS_rc%nensem(n)), tmp0(LIS_rc%nensem(n)) + real :: tmp1(LIS_rc%nensem(n)),tmp2(LIS_rc%nensem(n)),tmp3(LIS_rc%nensem(n)),tmp4(LIS_rc%nensem(n)) + logical :: update_flag_tile(LIS_rc%npatch(n,LIS_rc%lsm_index)) + logical :: flag_ens(LIS_rc%ngrid(n)) + logical :: flag_tmp(LIS_rc%nensem(n)) + logical :: update_flag_ens(LIS_rc%ngrid(n)) + logical :: update_flag_new(LIS_rc%ngrid(n)) + integer :: RESULT, pcount, icount + real :: MaxEnsSM1 ,MaxEnsSM2 ,MaxEnsSM3 ,MaxEnsSM4 + real :: MinEnsSM1 ,MinEnsSM2 ,MinEnsSM3 ,MinEnsSM4 + real :: MaxEns_sh2o1, MaxEns_sh2o2, MaxEns_sh2o3, MaxEns_sh2o4 + real :: smc_rnd, smc_tmp + real :: sh2o_tmp, sh2o_rnd + INTEGER, DIMENSION (1) :: seed + + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 1",sm1Field,rc=status) + call LIS_verify(status,& + "ESMF_StateSet: Soil Moisture Layer 1 failed in noahmp36_setsoilmLAI") + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 2",sm2Field,rc=status) + call LIS_verify(status,& + "ESMF_StateSet: Soil Moisture Layer 2 failed in noahmp36_setsoilmLAI") + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 3",sm3Field,rc=status) + call LIS_verify(status,& + "ESMF_StateSet: Soil Moisture Layer 3 failed in noahmp36_setsoilmLAI") + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 4",sm4Field,rc=status) + call LIS_verify(status,& + "ESMF_StateSet: Soil Moisture Layer 4 failed in noahmp36_setsoilmLAI") + call ESMF_StateGet(LSM_State,"LAI",laiField,rc=status) + call LIS_verify(status,& + "ESMF_StateSet: LAI failed in noahmp36_setsoilmLAI") + + call ESMF_FieldGet(sm1Field,localDE=0,farrayPtr=soilm1,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 1 failed in noahmp36_setsoilmLAI") + call ESMF_FieldGet(sm2Field,localDE=0,farrayPtr=soilm2,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 2 failed in noahmp36_setsoilmLAI") + call ESMF_FieldGet(sm3Field,localDE=0,farrayPtr=soilm3,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 3 failed in noahmp36_setsoilmLAI") + call ESMF_FieldGet(sm4Field,localDE=0,farrayPtr=soilm4,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 4 failed in noahmp36_setsoilmLAI") + call ESMF_FieldGet(laiField,localDE=0,farrayPtr=lai,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: LAI failed in noahmp36_setsoilmLAI") + + update_flag = .true. + update_flag_tile= .true. + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + + ! MN: NOTE: SMCMAX and SMCWLT are not stored in the data structure but we + ! can get module variables MAXSMC and WLTSMC from the module_sf_noahlsm_36 + SOILTYP = NOAHMP36_struc(n)%noahmp36(t)%soiltype + MAX_THRESHOLD = MAXSMC (SOILTYP) + sm_threshold = MAXSMC (SOILTYP) - 0.02 + + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%row) + + !MN: delta = X(+) - X(-) + !NOTE: "noahmp36_updatesoilm.F90" updates the soilm_(t) + delta1 = soilm1(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(1) + delta2 = soilm2(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(2) + delta3 = soilm3(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(3) + delta4 = soilm4(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(4) + + ! MN: check MIN_THRESHOLD < volumetric liquid soil moisture < threshold + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(1)+delta1.gt.MIN_THRESHOLD .and.& + NOAHMP36_struc(n)%noahmp36(t)%sh2o(1)+delta1.lt.& + sm_threshold) then + update_flag(gid) = update_flag(gid).and.(.true.) + ! MN save the flag for each tile (col*row*ens) (64*44)*20 + update_flag_tile(t) = update_flag_tile(t).and.(.true.) + else + update_flag(gid) = update_flag(gid).and.(.false.) + update_flag_tile(t) = update_flag_tile(t).and.(.false.) + endif + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(2)+delta2.gt.MIN_THRESHOLD .and.& + NOAHMP36_struc(n)%noahmp36(t)%sh2o(2)+delta2.lt.sm_threshold) then + update_flag(gid) = update_flag(gid).and.(.true.) + update_flag_tile(t) = update_flag_tile(t).and.(.true.) + else + update_flag(gid) = update_flag(gid).and.(.false.) + update_flag_tile(t) = update_flag_tile(t).and.(.false.) + endif + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(3)+delta3.gt.MIN_THRESHOLD .and.& + NOAHMP36_struc(n)%noahmp36(t)%sh2o(3)+delta3.lt.sm_threshold) then + update_flag(gid) = update_flag(gid).and.(.true.) + update_flag_tile(t) = update_flag_tile(t).and.(.true.) + else + update_flag(gid) = update_flag(gid).and.(.false.) + update_flag_tile(t) = update_flag_tile(t).and.(.false.) + endif + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(4)+delta4.gt.MIN_THRESHOLD .and.& + NOAHMP36_struc(n)%noahmp36(t)%sh2o(4)+delta4.lt.sm_threshold) then + update_flag(gid) = update_flag(gid).and.(.true.) + update_flag_tile(t) = update_flag_tile(t).and.(.true.) + else + update_flag(gid) = update_flag(gid).and.(.false.) + update_flag_tile(t) = update_flag_tile(t).and.(.false.) + endif + enddo + +!----------------------------------------------------------------------------------------- +! MN create new falg: if update falg for 50% of the ensemble members is true +! then update the state variabels +!----------------------------------------------------------------------------------------- + update_flag_ens = .True. + do i=1,LIS_rc%npatch(n,LIS_rc%lsm_index),LIS_rc%nensem(n) + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(i)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(i)%row) + flag_tmp=update_flag_tile(i:i+LIS_rc%nensem(n)-1) + !flag_tmp=update_flag_tile((i-1)*LIS_rc%nensem(n)+1:(i)*LIS_rc%nensem(n)) + pcount = COUNT(flag_tmp) ! Counts the number of .TRUE. elements + if (pcount.lt.LIS_rc%nensem(n)*0.5) then ! 50% + update_flag_ens(gid)= .False. + endif + update_flag_new(gid)= update_flag(gid).or.update_flag_ens(gid) ! new flag + enddo + + ! MN print + if(i.eq.66) then !i=66 ! --> domain's center 1376 + if(LIS_rc%hr.eq.12) then + write(2001,'(I4, 2x, 3(I2,x), 2x, 23(L1,2x))'),& + i, LIS_rc%mo, LIS_rc%da, LIS_rc%hr,update_flag_tile& + ((i-1)*LIS_rc%nensem(n)+1:(i)*LIS_rc%nensem(n)),& + update_flag_ens(i), update_flag_new(i), update_flag(i) + endif !mn + endif + + ! update step + ! loop over grid points + do i=1,LIS_rc%npatch(n,LIS_rc%lsm_index),LIS_rc%nensem(n) + + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(i)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(i)%row) + + !if(update_flag(gid)) then + if(update_flag_new(gid)) then +!----------------------------------------------------------------------------------------- + ! 1- update the states + ! 1-1- if update flag for tile is TRUE --> apply the DA update + ! 1-2- if update flag for tile is FALSE --> resample the states + ! 2- adjust the states +!----------------------------------------------------------------------------------------- + ! store update value for cases that update_flag_tile & update_flag_new are TRUE + ! update_flag_tile = TRUE --> means met the both min and max threshold + + tmp1 = LIS_rc%udef + tmp2 = LIS_rc%udef + tmp3 = LIS_rc%udef + tmp4 = LIS_rc%udef + !icount = 1 + do m=1,LIS_rc%nensem(n) + t = i+m-1 + !t = (i-1)*LIS_rc%nensem(n)+m + + if(update_flag_tile(t)) then + + tmp1(m) = soilm1(t) + tmp2(m) = soilm2(t) + tmp3(m) = soilm3(t) + tmp4(m) = soilm4(t) + !icount = icount + 1 + endif + enddo + + MaxEnsSM1 = -10000 + MaxEnsSM2 = -10000 + MaxEnsSM3 = -10000 + MaxEnsSM4 = -10000 + + MinEnsSM1 = 10000 + MinEnsSM2 = 10000 + MinEnsSM3 = 10000 + MinEnsSM4 = 10000 + + do m=1,LIS_rc%nensem(n) + if(tmp1(m).ne.LIS_rc%udef) then + MaxEnsSM1 = max(MaxEnsSM1, tmp1(m)) + MaxEnsSM2 = max(MaxEnsSM2, tmp2(m)) + MaxEnsSM3 = max(MaxEnsSM3, tmp3(m)) + MaxEnsSM4 = max(MaxEnsSM4, tmp4(m)) + + MinEnsSM1 = min(MinEnsSM1, tmp1(m)) + MinEnsSM2 = min(MinEnsSM2, tmp2(m)) + MinEnsSM3 = min(MinEnsSM3, tmp3(m)) + MinEnsSM4 = min(MinEnsSM4, tmp4(m)) + + endif + enddo + + + ! loop over tile + do m=1,LIS_rc%nensem(n) + t = i+m-1 + !t = (i-1)*LIS_rc%nensem(n)+m + + ! MN check update status for each tile + if(update_flag_tile(t)) then + + delta1 = soilm1(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(1) + delta2 = soilm2(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(2) + delta3 = soilm3(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(3) + delta4 = soilm4(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(4) + + NOAHMP36_struc(n)%noahmp36(t)%sh2o(1) = NOAHMP36_struc(n)%noahmp36(t)%sh2o(1)+& + delta1 + NOAHMP36_struc(n)%noahmp36(t)%smc(1) = soilm1(t) + if(soilm1(t).lt.0) then + print*, 'setsoilm1 ',t,soilm1(t) + stop + endif + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(2)+delta2.gt.MIN_THRESHOLD .and.& + NOAHMP36_struc(n)%noahmp36(t)%sh2o(2)+delta2.lt.sm_threshold) then + NOAHMP36_struc(n)%noahmp36(t)%sh2o(2) = NOAHMP36_struc(n)%noahmp36(t)%sh2o(2)+& + soilm2(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(2) + NOAHMP36_struc(n)%noahmp36(t)%smc(2) = soilm2(t) + if(soilm2(t).lt.0) then + print*, 'setsoilm2 ',t,soilm2(t) + stop + endif + endif + + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(3)+delta3.gt.MIN_THRESHOLD .and.& + NOAHMP36_struc(n)%noahmp36(t)%sh2o(3)+delta3.lt.sm_threshold) then + NOAHMP36_struc(n)%noahmp36(t)%sh2o(3) = NOAHMP36_struc(n)%noahmp36(t)%sh2o(3)+& + soilm3(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(3) + NOAHMP36_struc(n)%noahmp36(t)%smc(3) = soilm3(t) + if(soilm3(t).lt.0) then + print*, 'setsoilm3 ',t,soilm3(t) + stop + endif + endif + + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(4)+delta4.gt.MIN_THRESHOLD .and.& + NOAHMP36_struc(n)%noahmp36(t)%sh2o(4)+delta4.lt.sm_threshold) then + NOAHMP36_struc(n)%noahmp36(t)%sh2o(4) = NOAHMP36_struc(n)%noahmp36(t)%sh2o(4)+& + soilm4(t)-NOAHMP36_struc(n)%noahmp36(t)%smc(4) + NOAHMP36_struc(n)%noahmp36(t)%smc(4) = soilm4(t) + + if(soilm4(t).lt.0) then + print*, 'setsoilm4 ',t,soilm4(t) + stop + endif + endif + + +!----------------------------------------------------------------------------------------- + ! randomly resample the smc from [MIN_THRESHOLD, Max value from DA @ that tiem step] +!----------------------------------------------------------------------------------------- + else + +!----------------------------------------------------------------------------------------- +! set the soil moisture to the ensemble mean +!----------------------------------------------------------------------------------------- + + ! use mean value + ! Assume sh2o = smc (i.e. ice content=0) + smc_tmp = (MaxEnsSM1 - MinEnsSM1)/2 + MinEnsSM1 + NOAHMP36_struc(n)%noahmp36(t)%sh2o(1) = smc_tmp + NOAHMP36_struc(n)%noahmp36(t)%smc(1) = smc_tmp + + smc_tmp = (MaxEnsSM2 - MinEnsSM2)/2 + MinEnsSM2 + NOAHMP36_struc(n)%noahmp36(t)%sh2o(2) = smc_tmp + NOAHMP36_struc(n)%noahmp36(t)%smc(2) = smc_tmp + + smc_tmp = (MaxEnsSM3 - MinEnsSM3)/2 + MinEnsSM3 + NOAHMP36_struc(n)%noahmp36(t)%sh2o(3) = smc_tmp + NOAHMP36_struc(n)%noahmp36(t)%smc(3) = smc_tmp + + smc_tmp = (MaxEnsSM4 - MinEnsSM4)/2 + MinEnsSM4 + NOAHMP36_struc(n)%noahmp36(t)%sh2o(4) = smc_tmp + NOAHMP36_struc(n)%noahmp36(t)%smc(4) = smc_tmp + + + endif ! flag for each tile + + enddo ! loop over tile + + else ! if update_flag_new(gid) is FALSE + if(LIS_rc%pert_bias_corr.eq.1) then + !-------------------------------------------------------------------------- + ! if no update is made, then we need to readjust the ensemble if pert bias + ! correction is turned on because the forcing perturbations may cause + ! biases to persist. + !-------------------------------------------------------------------------- + bounds_violation = .true. + nIter = 0 + ens_flag = .true. + + do while(bounds_violation) + niter = niter + 1 + !t_unpert = i*LIS_rc%nensem(n) + t_unpert = i+LIS_rc%nensem(n)-1 + do j=1,4 + delta(j) = 0.0 + do m=1,LIS_rc%nensem(n)-1 + t = i+m-1 + !t = (i-1)*LIS_rc%nensem(n)+m + + if(m.ne.LIS_rc%nensem(n)) then + delta(j) = delta(j) + & + (NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) - & + NOAHMP36_struc(n)%noahmp36(t_unpert)%sh2o(j)) + endif + + enddo + enddo + + do j=1,4 + delta(j) =delta(j)/(LIS_rc%nensem(n)-1) + do m=1,LIS_rc%nensem(n)-1 + t = i+m-1 + !t = (i-1)*LIS_rc%nensem(n)+m + SOILTYP = NOAHMP36_struc(n)%noahmp36(t)%soiltype + MAX_THRESHOLD = MAXSMC (SOILTYP) + sm_threshold = MAXSMC (SOILTYP) - 0.02 + + tmpval = NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) - & + delta(j) + if(tmpval.le.MIN_THRESHOLD) then + NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) = & + max(NOAHMP36_struc(n)%noahmp36(t_unpert)%sh2o(j),& + MIN_THRESHOLD) + NOAHMP36_struc(n)%noahmp36(t)%smc(j) = & + max(NOAHMP36_struc(n)%noahmp36(t_unpert)%smc(j),& + MIN_THRESHOLD) + ens_flag(m) = .false. + elseif(tmpval.ge.sm_threshold) then + NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) = & + min(NOAHMP36_struc(n)%noahmp36(t_unpert)%sh2o(j),& + sm_threshold) + NOAHMP36_struc(n)%noahmp36(t)%smc(j) = & + min(NOAHMP36_struc(n)%noahmp36(t_unpert)%smc(j),& + sm_threshold) + ens_flag(m) = .false. + endif + enddo + enddo + + !-------------------------------------------------------------------------- + ! Recalculate the deltas and adjust the ensemble + !-------------------------------------------------------------------------- + do j=1,4 + delta(j) = 0.0 + do m=1,LIS_rc%nensem(n)-1 + t = i+m-1 + !t = (i-1)*LIS_rc%nensem(n)+m + if(m.ne.LIS_rc%nensem(n)) then + delta(j) = delta(j) + & + (NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) - & + NOAHMP36_struc(n)%noahmp36(t_unpert)%sh2o(j)) + endif + enddo + enddo + + do j=1,4 + delta(j) =delta(j)/(LIS_rc%nensem(n)-1) + do m=1,LIS_rc%nensem(n)-1 + t = i+m-1 + !t = (i-1)*LIS_rc%nensem(n)+m + + if(ens_flag(m)) then + tmpval = NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) - & + delta(j) + SOILTYP = NOAHMP36_struc(n)%noahmp36(t)%soiltype + MAX_THRESHOLD = MAXSMC (SOILTYP) + + if(.not.(tmpval.le.0.0 .or.& + tmpval.gt.(MAX_THRESHOLD))) then + + NOAHMP36_struc(n)%noahmp36(t)%smc(j) = & + NOAHMP36_struc(n)%noahmp36(t)%smc(j) - delta(j) + NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) = & + NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) - delta(j) + bounds_violation = .false. + endif + endif + + tmpval = NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) + SOILTYP = NOAHMP36_struc(n)%noahmp36(t)%soiltype + MAX_THRESHOLD = MAXSMC (SOILTYP) + + if(tmpval.le.0.0 .or.& + tmpval.gt.(MAX_THRESHOLD)) then + bounds_violation = .true. + else + bounds_violation = .false. + endif + enddo + enddo + + if(nIter.gt.10.and.bounds_violation) then + !-------------------------------------------------------------------------- + ! All else fails, set to the bounds + !-------------------------------------------------------------------------- + + write(LIS_logunit,*) '[ERR] Ensemble structure violates physical bounds ' + write(LIS_logunit,*) '[ERR] Please adjust the perturbation settings ..' + + do j=1,4 + do m=1,LIS_rc%nensem(n) + t = i+m-1 + !t = (i-1)*LIS_rc%nensem(n)+m + + SOILTYP = NOAHMP36_struc(n)%noahmp36(t)%soiltype + MAX_THRESHOLD = MAXSMC (SOILTYP) + + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(j).gt.MAX_THRESHOLD.or.& + NOAHMP36_struc(n)%noahmp36(t)%smc(j).gt.MAX_THRESHOLD) then + NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) = MAX_THRESHOLD + NOAHMP36_struc(n)%noahmp36(t)%smc(j) = MAX_THRESHOLD + endif + + if(NOAHMP36_struc(n)%noahmp36(t)%sh2o(j).lt.MIN_THRESHOLD.or.& + NOAHMP36_struc(n)%noahmp36(t)%smc(j).lt.MIN_THRESHOLD) then + NOAHMP36_struc(n)%noahmp36(t)%sh2o(j) = MIN_THRESHOLD + NOAHMP36_struc(n)%noahmp36(t)%smc(j) = MIN_THRESHOLD + endif +! print*, i, m +! print*, 'smc',t, NOAHMP36_struc(n)%noahmp36(t)%smc(:) +! print*, 'sh2o ',t,NOAHMP36_struc(n)%noahmp36(t)%sh2o(:) +! print*, 'max ',t,MAX_THRESHOLD !NOAHMP36_struc(n)%noahmp36(t)%smcmax + enddo +! call LIS_endrun() + enddo + endif + end do + endif + endif + enddo + + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) +! XLAI = max(LFMASS*LAPM,laimin) + + if(sla(NOAHMP36_struc(n)%noahmp36(t)%vegetype).ne.0) then + NOAHMP36_struc(n)%noahmp36(t)%lai = lai(t) + lfmass = lai(t)/(sla(NOAHMP36_struc(n)%noahmp36(t)%vegetype)/1000.0) + NOAHMP36_struc(n)%noahmp36(t)%lfmass = lfmass + endif + enddo + +end subroutine noahmp36_setsoilmLAI + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_updatesoilmLAI.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_updatesoilmLAI.F90 new file mode 100644 index 000000000..120367e81 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_updatesoilmLAI.F90 @@ -0,0 +1,236 @@ +!-----------------------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: noahmp36_updatesoilmLAI +! \label{noahmp36_updatesoilmLAI} +! +! !REVISION HISTORY: +! 27Feb2005: Sujay Kumar; Initial Specification +! 25Jun2006: Sujay Kumar: Updated for the ESMF design +! 9 Sep 2016: Mahdi Navari; Modified for NoahMP36 +! To do: makes it general for x layers (currently hard coded for 4 layers) +! 18 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! +! !INTERFACE: +subroutine noahmp36_updatesoilmLAI(n, LSM_State, LSM_Incr_State) +! !USES: + use ESMF + use LIS_coreMod + use LIS_logMod + use noahmp36_lsmMod + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + type(ESMF_State) :: LSM_State + type(ESMF_State) :: LSM_Incr_State +! +! !DESCRIPTION: +! +! This routine assigns the soil moisture prognostic variables to noah's +! model space. +! +!EOP + + type(ESMF_Field) :: sm1Field + type(ESMF_Field) :: sm2Field + type(ESMF_Field) :: sm3Field + type(ESMF_Field) :: sm4Field + type(ESMF_Field) :: laiField + type(ESMF_Field) :: sm1IncrField + type(ESMF_Field) :: sm2IncrField + type(ESMF_Field) :: sm3IncrField + type(ESMF_Field) :: sm4IncrField + type(ESMF_Field) :: laiIncrField + + real, pointer :: soilm1(:) + real, pointer :: soilm2(:) + real, pointer :: soilm3(:) + real, pointer :: soilm4(:) + real, pointer :: lai(:) + real, pointer :: soilmIncr1(:) + real, pointer :: soilmIncr2(:) + real, pointer :: soilmIncr3(:) + real, pointer :: soilmIncr4(:) + real, pointer :: laiincr(:) + integer :: t,i,m,gid + integer :: status + + real :: laitmp,laimax,laimin + + logical :: update_flag(LIS_rc%ngrid(n)) + real :: perc_violation(LIS_rc%ngrid(n)) + + real :: laimean(LIS_rc%ngrid(n)) + integer :: nlaimean(LIS_rc%ngrid(n)) + + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 1",sm1Field,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 1 failed in noahmp36_updatesoilmLAI") + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 2",sm2Field,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 2 failed in noahmp36_updatesoilmLAI") + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 3",sm3Field,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 3 failed in noahmp36_updatesoilmLAI") + call ESMF_StateGet(LSM_State,"Soil Moisture Layer 4",sm4Field,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 4 failed in noahmp36_updatesoilmLAI") + + call ESMF_FieldGet(sm1Field,localDE=0,farrayPtr=soilm1,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 1 failed in noahmp36_updatesoilmLAI") + call ESMF_FieldGet(sm2Field,localDE=0,farrayPtr=soilm2,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 2 failed in noahmp36_updatesoilmLAI") + call ESMF_FieldGet(sm3Field,localDE=0,farrayPtr=soilm3,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 3 failed in noahmp36_updatesoilmLAI") + call ESMF_FieldGet(sm4Field,localDE=0,farrayPtr=soilm4,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 4 failed in noahmp36_updatesoilmLAI") + + call ESMF_StateGet(LSM_Incr_State,"Soil Moisture Layer 1",sm1IncrField,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 1 failed in noahmp36_updatesoilmLAI") + call ESMF_StateGet(LSM_Incr_State,"Soil Moisture Layer 2",sm2IncrField,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 2 failed in noahmp36_updatesoilmLAI") + call ESMF_StateGet(LSM_Incr_State,"Soil Moisture Layer 3",sm3IncrField,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 3 failed in noahmp36_updatesoilmLAI") + call ESMF_StateGet(LSM_Incr_State,"Soil Moisture Layer 4",sm4IncrField,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: Soil Moisture Layer 4 failed in noahmp36_updatesoilmLAI") + + call ESMF_FieldGet(sm1IncrField,localDE=0,farrayPtr=soilmIncr1,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 1 failed in noahmp36_updatesoilmLAI") + call ESMF_FieldGet(sm2IncrField,localDE=0,farrayPtr=soilmIncr2,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 2 failed in noahmp36_updatesoilmLAI") + call ESMF_FieldGet(sm3IncrField,localDE=0,farrayPtr=soilmIncr3,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 3 failed in noahmp36_updatesoilmLAI") + call ESMF_FieldGet(sm4IncrField,localDE=0,farrayPtr=soilmIncr4,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: Soil Moisture Layer 4 failed in noahmp36_updatesoilmLAI") + + + call ESMF_StateGet(LSM_State,"LAI",laiField,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: LSM_State, failed in noahmp36_updatesoilmLAI") + + call ESMF_StateGet(LSM_Incr_State,"LAI",laiIncrField,rc=status) + call LIS_verify(status,& + "ESMF_StateGet: LSM_Incr_State LAI failed in noahmp36_updatesoilmLAI") + + call ESMF_FieldGet(laiField,localDE=0,farrayPtr=lai,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: laiField failed in noahmp36_updatesoilmLAI") + + call ESMF_FieldGet(laiIncrField,localDE=0,farrayPtr=laiincr,rc=status) + call LIS_verify(status,& + "ESMF_FieldGet: laiIncrField failed in noahmp36_updatesoilmLAI") + + call ESMF_AttributeGet(laiField,"Max Value",laimax,rc=status) + call LIS_verify(status,& + "ESMF_AttributeGet: laiField Max Value failed in noahmp36_updatesoilmLAI") + call ESMF_AttributeGet(laiField,"Min Value",laimin,rc=status) + call LIS_verify(status,& + "ESMF_AttributeGet: laiField Min Value failed in noahmp36_updatesoilmLAI") + + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + soilm1(t) = soilm1(t) + soilmIncr1(t) + soilm2(t) = soilm2(t) + soilmIncr2(t) + soilm3(t) = soilm3(t) + soilmIncr3(t) + soilm4(t) = soilm4(t) + soilmIncr4(t) + enddo + + + update_flag = .true. + perc_violation = 0.0 + laimean = 0.0 + nlaimean = 0 + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%row) + + laitmp = lai(t) + laiincr(t) + + + if(laitmp.lt.laimin.or.laitmp.gt.laimax) then + update_flag(gid) = .false. + perc_violation(gid) = perc_violation(gid) +1 + endif + + enddo + + do gid=1,LIS_rc%ngrid(n) + perc_violation(gid) = perc_violation(gid)/LIS_rc%nensem(n) + enddo + +! For ensembles that are unphysical, compute the +! ensemble average after excluding them. This +! is done only if the majority of the ensemble +! members are good (>60%) + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%row) + if(.not.update_flag(gid)) then + if(perc_violation(gid).lt.0.8) then + if((lai(t)+laiincr(t).gt.laimin).and.& + (lai(t)+laiincr(t).lt.laimax)) then + laimean(gid) = laimean(gid) + & + lai(t) + laiincr(t) + nlaimean(gid) = nlaimean(gid) + 1 + endif + endif + endif + enddo + + do gid=1,LIS_rc%ngrid(n) + if(nlaimean(gid).gt.0) then + laimean(gid) = laimean(gid)/nlaimean(gid) + endif + enddo + + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + gid = LIS_domain(n)%gindex(& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%col,& + LIS_surface(n,LIS_rc%lsm_index)%tile(t)%row) + + laitmp = lai(t) + laiincr(t) + +! If the update is unphysical, simply set to the average of +! the good ensemble members. If all else fails, do not +! update. + + if(update_flag(gid)) then + lai(t) = laitmp + elseif(perc_violation(gid).lt.0.8) then + if(laitmp.lt.laimin.or.laitmp.gt.laimax) then + lai(t) = laimean(gid) + else + lai(t) = lai(t) + laiincr(t) + endif + endif + enddo + +end subroutine noahmp36_updatesoilmLAI + diff --git a/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_write_soilmLAI.F90 b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_write_soilmLAI.F90 new file mode 100644 index 000000000..1c91beef0 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/da_soilmLAI/noahmp36_write_soilmLAI.F90 @@ -0,0 +1,73 @@ +!-----------------------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: noahmp36_write_soilmLAI +! \label{noahmp36_write_soilmLAI} +! +! !REVISION HISTORY: +! 27Feb2005: Sujay Kumar; Initial Specification +! 25Jun2006: Sujay Kumar: Updated for the ESMF design +! 1 Aug 2016: Mahdi Navari; Modified for NoahMP36 +! To do: make it general for x layers (currently hard coded for 4 layers) +! 18 Jun 2021: Michel Bechtold: SM and LAI updating with S1 backscatter w/ WCM +! +! !INTERFACE: +subroutine noahmp36_write_soilmLAI(ftn,n, LSM_State) + +! !USES: + use ESMF + use LIS_coreMod, only : LIS_rc + use noahmp36_lsmMod + use LIS_historyMod, only : LIS_writevar_restart + implicit none +! !ARGUMENTS: + integer, intent(in) :: ftn + integer, intent(in) :: n + type(ESMF_State) :: LSM_State +! +! !DESCRIPTION: +! +! Returns the soilmoisture and LAI related state prognostic variables for +! data assimilation +! +! The arguments are: +! \begin{description} +! \item[n] index of the nest \newline +! \item[LSM\_State] ESMF State container for LSM state variables \newline +! \end{description} +!EOP + integer :: t + real, allocatable :: tmp(:) + + allocate(tmp(LIS_rc%npatch(n,LIS_rc%lsm_index))) + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + tmp(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(1) + enddo + call LIS_writevar_restart(ftn,n,1,tmp) + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + tmp(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(2) + enddo + call LIS_writevar_restart(ftn,n,1,tmp) + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + tmp(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(3) + enddo + call LIS_writevar_restart(ftn,n,1,tmp) + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + tmp(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(4) + enddo + call LIS_writevar_restart(ftn,n,1,tmp) + deallocate(tmp) + +end subroutine noahmp36_write_soilmLAI + diff --git a/lis/surfacemodels/land/noahmp.3.6/sfc_wcm/noahmp36_sfc2wcm.F90 b/lis/surfacemodels/land/noahmp.3.6/sfc_wcm/noahmp36_sfc2wcm.F90 new file mode 100644 index 000000000..6b5d18816 --- /dev/null +++ b/lis/surfacemodels/land/noahmp.3.6/sfc_wcm/noahmp36_sfc2wcm.F90 @@ -0,0 +1,59 @@ +!-----------------------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----------------------- +#include "LIS_misc.h" +!BOP +! !ROUTINE: Noahmp36_sfc2wcm +! \label{Noahmp36_sfc2wcm} +! +! !REVISION HISTORY: +! 4 Sep 2020: Sara Modanesi; Initial Specification + +! !INTERFACE: +subroutine noahmp36_sfc2wcm(n, sfcState) +! !USES: + use ESMF + use LIS_coreMod + use LIS_logMod, only : LIS_verify + use LIS_constantsMod, only : LIS_CONST_RHOFW + + use Noahmp36_lsmMod + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + type(ESMF_State) :: sfcState + +! FUNCTIONS + +! +! !DESCRIPTION: +! This subroutine assigns the noahmp36 specific surface variables +! to WCM. +! +!EOP + type(ESMF_Field) :: smcField, laiField + real, pointer :: soil_moisture_content(:), leaf_area_index(:) + integer :: t,status + + call ESMF_StateGet(sfcState,"Soil Moisture Content",smcField,rc=status) + call LIS_verify(status) + call ESMF_FieldGet(smcField,localDE=0,farrayPtr=soil_moisture_content,rc=status) + call LIS_verify(status) + + call ESMF_StateGet(sfcState,"Leaf Area Index",laiField,rc=status) + call LIS_verify(status) + call ESMF_FieldGet(laiField,localDE=0,farrayPtr=leaf_area_index, rc=status) + call LIS_verify(status) + + do t=1,LIS_rc%npatch(n,LIS_rc%lsm_index) + soil_moisture_content(t) = NOAHMP36_struc(n)%noahmp36(t)%smc(1) + Leaf_Area_Index(t) = NOAHMP36_struc(n)%noahmp36(t)%lai + enddo + +end subroutine noahmp36_sfc2wcm