diff --git a/ldt/DAobs/LISlsmTEFFobs/LISlsmTEFF_obsMod.F90 b/ldt/DAobs/LISlsmTEFFobs/LISlsmTEFF_obsMod.F90 index 5d171a908..33e0742c1 100644 --- a/ldt/DAobs/LISlsmTEFFobs/LISlsmTEFF_obsMod.F90 +++ b/ldt/DAobs/LISlsmTEFFobs/LISlsmTEFF_obsMod.F90 @@ -52,6 +52,10 @@ module LISlsmTEFF_obsMod character*20 :: data_category character*20 :: area_of_data character*20 :: write_interval + ! EMK + integer :: num_ens + integer :: num_tiles + integer :: ntiles_pergrid !-------------------------------------------------------- ! interpolation/upscaling weights !-------------------------------------------------------- @@ -127,6 +131,38 @@ subroutine LISlsmTEFF_obsInit() label="LIS soil temperature output directory:",rc=rc) call LDT_verify(rc,'LIS soil temperature output directory: not defined') + ! EMK Add LIS ensemble information + call ESMF_ConfigGetAttribute(LDT_config, lsmteffobs%num_ens, & + label="LIS ensemble size:", rc=rc) + call LDT_verify(rc,'LIS ensemble size: not defined') + if (lsmteffobs%num_ens < 1) then + write(LDT_logunit,*)'[ERR] LIS ensemble size must be at least 1!' + write(LDT_logunit,*)'[ERR] Read in ', lsmteffobs%num_ens + call LDT_endrun() + end if + + call ESMF_ConfigGetAttribute(LDT_config, lsmteffobs%num_tiles, & + label="LIS total number of tiles (including ensembles):", rc=rc) + call LDT_verify(rc,'LIS total number of tiles (including ensembles): ' // & + 'not defined') + if (lsmteffobs%num_tiles < 1) then + write(LDT_logunit,*) & + '[ERR] LIS total number of tiles (including ensembles) must be' & + //'at least 1!' + write(LDT_logunit,*)'[ERR] Read in ', lsmteffobs%num_tiles + call LDT_endrun() + end if + + call ESMF_ConfigGetAttribute(LDT_config, lsmteffobs%ntiles_pergrid, & + label="LIS number of tiles per grid point:", rc=rc) + call LDT_verify(rc,'LIS number of tiles per grid point: not defined') + if (lsmteffobs%num_tiles < 1) then + write(LDT_logunit,*) & + '[ERR] LIS number of tiles per grid point must be at least 1!' + write(LDT_logunit,*)'[ERR] Read in ', lsmteffobs%ntiles_pergrid + call LDT_endrun() + end if + ! WMO-convention specific identifiers if ( lsmteffobs%wstyle == "WMO convention") then call ESMF_ConfigGetAttribute(LDT_config,lsmteffobs%security_class, & diff --git a/ldt/DAobs/LISlsmTEFFobs/readLISlsmTEFFObs.F90 b/ldt/DAobs/LISlsmTEFFobs/readLISlsmTEFFObs.F90 index 4cd00b183..d24976fa8 100644 --- a/ldt/DAobs/LISlsmTEFFobs/readLISlsmTEFFObs.F90 +++ b/ldt/DAobs/LISlsmTEFFobs/readLISlsmTEFFObs.F90 @@ -9,11 +9,11 @@ !-------------------------END NOTICE -- DO NOT EDIT----------------------- #include "LDT_misc.h" !BOP -! +! ! !ROUTINE: readLISlsmTEFFobs ! \label{readLISlsmTEFFobs} ! -! !INTERFACE: +! !INTERFACE: subroutine readLISlsmTEFFobs(n) ! !USES: #if(defined USE_NETCDF3 || defined USE_NETCDF4) @@ -29,19 +29,19 @@ subroutine readLISlsmTEFFobs(n) use LISlsmTEFF_obsMod, only : lsmteffobs use LDT_timeMgrMod ! -! !DESCRIPTION: -! This routine reads the soil temperature fields from a LIS model -! simulation. +! !DESCRIPTION: +! This routine reads the soil temperature fields from a LIS model +! simulation. ! !EOP implicit none integer, intent(in) :: n - character*200 :: fname + character*200 :: fname logical :: file_exists real :: teff_data(LDT_rc%lnc(n),LDT_rc%lnr(n)) - + integer :: t, index integer :: ftn integer :: iret @@ -50,6 +50,8 @@ subroutine readLISlsmTEFFobs(n) real :: teffvalue1d(lsmteffobs%nc*lsmteffobs%nr) real :: Tsoil01value2d(lsmteffobs%nc, lsmteffobs%nr) real :: Tsoil02value2d(lsmteffobs%nc, lsmteffobs%nr) + real :: tsoil(lsmteffobs%nc, lsmteffobs%nr, 4) + integer :: rc1 integer :: c,r character*20 :: vname integer :: varid @@ -59,6 +61,9 @@ subroutine readLISlsmTEFFobs(n) real :: gmt real :: lhour integer :: zone + integer, allocatable :: ntiles_pergrid(:) + integer, allocatable :: str_tind(:) + integer :: gid interface subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & @@ -101,80 +106,67 @@ end subroutine create_lsm_teff_output_filename inquire(file=trim(fname),exist=file_exists) - if(file_exists) then + if(file_exists) then write(LDT_logunit,*) '[INFO] reading LSM output ',trim(fname) - if(lsmteffobs%format.eq."binary") then + if(lsmteffobs%format.eq."binary") then write(LDT_logunit,*) '[ERR] DA preprocessing on the binary format is not ' write(LDT_logunit,*) '[ERR] currently supported. Program stopping....' call LDT_endrun() - -#if 0 - ftn = LDT_getNextUnitNumber() - open(ftn,file=trim(fname), form='unformatted') - if(file_exists) then - do index=1,LDT_MOC_COUNT - call LDT_readLISSingleBinaryVar(n,ftn,LDT_DAobsDataPtr(n,index)%dataEntryPtr) - enddo - else - print*, 'LSM file ',trim(fname),' does not exist' - print*, 'Program stopping.. ' - stop - endif - call LDT_releaseUnitNumber(ftn) -#endif - elseif(lsmteffobs%format.eq."grib1") then + + elseif(lsmteffobs%format.eq."grib1") then write(LDT_logunit,*) '[ERR] DA preprocessing on the grib1 format is not ' write(LDT_logunit,*) '[ERR] currently supported. Program stopping....' call LDT_endrun() -#if 0 - ftn = LDT_getNextUnitNumber() - open(ftn,file=trim(fname), form='unformatted') - if(file_exists) then - do index=1,LDT_MOC_COUNT - call LDT_readLISSingleBinaryVar(n,ftn,LDT_DAobsDataPtr(n,index)%dataEntryPtr) - enddo - else - print*, 'LSM file ',trim(fname),' does not exist' - print*, 'Program stopping.. ' - stop - endif - call LDT_releaseUnitNumber(ftn) -#endif + elseif(lsmteffobs%format.eq."netcdf") then +#if(defined USE_NETCDF3 || defined USE_NETCDF4) - elseif(lsmteffobs%format.eq."netcdf") then -#if(defined USE_NETCDF3 || defined USE_NETCDF4) - iret = nf90_open(path=trim(fname),mode=nf90_nowrite, ncid=ftn) call LDT_verify(iret, 'Error opening file '//trim(fname)) - -! The code looks for instantaneous variables. If it doesn't exist, -! the time averaged data fields will be read in. + +! The code looks for instantaneous variables. If it doesn't exist, +! the time averaged data fields will be read in. iret = nf90_inq_varid(ftn, 'SoilTemp_inst', varid) vname = 'SoilTemp_inst' - if(iret.ne.0) then + if(iret.ne.0) then vname = 'SoilTemp_tavg' endif - - if(lsmteffobs%datares.eq.LDT_rc%gridDesc(n,10)) then - call LDT_readLISSingleNetcdfVar(n,ftn, vname,& - 1,lsmteffobs%nc, lsmteffobs%nr, Tsoil01value2d) - call LDT_readLISSingleNetcdfVar(n,ftn, vname,& - 2,lsmteffobs%nc, lsmteffobs%nr, Tsoil02value2d) - else - call LDT_readLISSingleNetcdfVar(n,ftn, vname,& - 1,lsmteffobs%nc, lsmteffobs%nr, Tsoil01value2d) - call LDT_readLISSingleNetcdfVar(n,ftn, vname,& - 2,lsmteffobs%nc, lsmteffobs%nr, Tsoil02value2d) - endif - iret = nf90_close(ftn) call LDT_verify(iret,'Error in nf90_close') + if ((LDT_rc%lnc(n) .ne. lsmteffobs%nc) .or. & + (LDT_rc%lnr(n) .ne. lsmteffobs%nr)) then + write(LDT_logunit,*)'[ERR] Dimension mismatch for LIS data!' + write(LDT_logunit,*)'[ERR] LDT_rc%lnc, LDT_rc%lnr = ', & + LDT_rc%lnc(n), LDT_rc%lnr(n) + write(LDT_Logunit,*)'[ERR] lsmteffobs%nc, lsmteffobs%nr = ', & + lsmteffobs%nc, lsmteffobs%nr + call LDT_endrun() + end if + tsoil = 0 + allocate(ntiles_pergrid(lsmteffobs%nc * lsmteffobs%nr)) + ntiles_pergrid = lsmteffobs%ntiles_pergrid ! Copy scalar to array + allocate(str_tind(lsmteffobs%nc * lsmteffobs%nr)) + do gid = 1, lsmteffobs%nc * lsmteffobs%nr + str_tind(gid) = ((gid - 1) * lsmteffobs%num_ens) + 1 + end do + call read_LIStsoil_data_usaf(n, lsmteffobs%num_tiles, str_tind, & + ntiles_pergrid, & + lsmteffobs%num_ens, & + fname, tsoil, rc1) + if (rc1 .ne. 0) then + write(LDT_logunit,*) '[ERR] Cannot read from ', trim(fname) + call LDT_endrun() + end if + tsoil01value2d = tsoil(:,:,1) + tsoil02value2d = tsoil(:,:,2) + deallocate(ntiles_pergrid) + deallocate(str_tind) + kk = 1.007 - cc_6am = 0.246; !Descending - cc_6pm = 1.000; !Ascending + cc_6am = 0.246; !Descending + cc_6pm = 1.000; !Ascending do r=1,lsmteffobs%nr do c=1, lsmteffobs%nc @@ -183,14 +175,18 @@ end subroutine create_lsm_teff_output_filename gmt = LDT_rc%hr call LDT_gmt2localtime(gmt, lon, lhour, zone) - if(lhour.eq.6) then + !if(lhour.eq.6) then ! Orig + if(lhour > 4 .and. lhour < 8) then ! EMK for 3-hrly data + if (Tsoil01value2d(c,r).gt.273.15.and.Tsoil02value2d(c,r).gt.273.15) then teffvalue1d(c+(r-1)*lsmteffobs%nc) = & kk * (cc_6am * Tsoil01value2d(c,r) + (1 - cc_6am) * Tsoil02value2d(c,r)) else teffvalue1d(c+(r-1)*lsmteffobs%nc) = LDT_rc%udef endif - elseif(lhour.eq.18) then + !elseif(lhour.eq.18) then ! Orig + elseif (lhour > 16 .and. lhour < 20) then ! EMK for 3-hrly data + if (Tsoil01value2d(c,r).gt.273.15.and.Tsoil02value2d(c,r).gt.273.15) then teffvalue1d(c+(r-1)*lsmteffobs%nc) = & kk * (cc_6pm * Tsoil01value2d(c,r) + (1 - cc_6pm) * Tsoil02value2d(c,r)) @@ -200,9 +196,10 @@ end subroutine create_lsm_teff_output_filename else teffvalue1d(c+(r-1)*lsmteffobs%nc) = LDT_rc%udef endif + enddo enddo - + call transformDataToLDTgrid_teff(n,teffvalue1d,teff_data) #endif @@ -234,7 +231,7 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & use LDT_coreMod, only : LDT_rc use LDT_logMod - implicit none + implicit none ! !ARGUMENTS: integer, intent(IN) :: n character(len=*) :: fname @@ -249,11 +246,11 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & character(len=*), optional :: data_category character(len=*), optional :: area_of_data character(len=*), optional :: write_interval -! -! !DESCRIPTION: +! +! !DESCRIPTION: ! Create the file name for the output data files. It creates both the GSWP ! style of output filenames and the standard LIS style. The convention used -! in LIS creates a filename in a hierarchical style (output directory, +! in LIS creates a filename in a hierarchical style (output directory, ! model name, date, file extention) ! ! 2 level hierarchy @@ -288,15 +285,15 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & ! DD = date \newline ! DT = data time \newline ! DF = data format \newline -! -! The arguments are: +! +! The arguments are: ! \begin{description} ! \item [n] ! index of the domain or nest ! \item [fname] -! the created file name. +! the created file name. ! \item [model\_name] -! string describing the name of the model +! string describing the name of the model ! \item [writeint] ! output writing interval of the model ! \item [style] @@ -307,7 +304,7 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & character(len=10) :: time character(len=5) :: zone integer, dimension(8) :: values - character(len=20) :: mname + character(len=20) :: mname character(len=10) :: cdate character(len=14) :: cdate1 character(len=2) :: fint @@ -322,32 +319,32 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & integer :: i, c mname = 'SURFACEMODEL' - if(wstyle.eq."4 level hierarchy") then + if(wstyle.eq."4 level hierarchy") then write(unit=cdate1, fmt='(i4.4, i2.2, i2.2, i2.2, i2.2)') & LDT_rc%yr, LDT_rc%mo, & LDT_rc%da, LDT_rc%hr,LDT_rc%mn - + dname = trim(odir)//'/' dname = trim(dname)//trim(mname)//'/' - + write(unit=cdate, fmt='(i4.4)') LDT_rc%yr dname = trim(dname)//trim(cdate)//'/' - + write(unit=cdate, fmt='(i4.4, i2.2, i2.2)') LDT_rc%yr, LDT_rc%mo, LDT_rc%da dname = trim(dname)//trim(cdate) - + out_fname = trim(dname)//'/LIS_HIST_'//trim(cdate1) - - write(unit=cdate, fmt='(a2,i2.2)') '.d',n + + write(unit=cdate, fmt='(a2,i2.2)') '.d',n out_fname = trim(out_fname)//trim(cdate) - + select case ( form ) case ( "binary" ) - if(wopt.eq."1d tilespace") then + if(wopt.eq."1d tilespace") then out_fname = trim(out_fname)//'.ts4r' - elseif(wopt.eq."2d gridspace") then + elseif(wopt.eq."2d gridspace") then out_fname = trim(out_fname)//'.gs4r' - elseif(wopt.eq."1d gridspace") then + elseif(wopt.eq."1d gridspace") then out_fname = trim(out_fname)//'.gs4r' endif case ("grib1") @@ -359,31 +356,31 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & case default call ldt_log_msg('ERR: create_lsm_teff_output_filename -- '// & 'Unrecognized output format') - call LDT_endrun + call LDT_endrun endselect elseif(wstyle.eq."3 level hierarchy") then write(unit=cdate1, fmt='(i4.4, i2.2, i2.2, i2.2, i2.2)') & LDT_rc%yr, LDT_rc%mo, & LDT_rc%da, LDT_rc%hr,LDT_rc%mn - + dname = trim(odir)//'/' dname = trim(dname)//trim(mname)//'/' - + write(unit=cdate, fmt='(i4.4, i2.2)') LDT_rc%yr, LDT_rc%mo dname = trim(dname)//trim(cdate)//'/' out_fname = trim(dname)//'LIS_HIST_'//trim(cdate1) - - write(unit=cdate, fmt='(a2,i2.2)') '.d',n + + write(unit=cdate, fmt='(a2,i2.2)') '.d',n out_fname = trim(out_fname)//trim(cdate) - + select case ( form ) case ("binary") - if(wopt.eq."1d tilespace") then + if(wopt.eq."1d tilespace") then out_fname = trim(out_fname)//'.ts4r' - elseif(wopt.eq."2d gridspace") then + elseif(wopt.eq."2d gridspace") then out_fname = trim(out_fname)//'.gs4r' - elseif(wopt.eq."1d gridspace") then + elseif(wopt.eq."1d gridspace") then out_fname = trim(out_fname)//'.gs4r' endif case ("grib1") @@ -395,28 +392,28 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & case default call ldt_log_msg('ERR: create_lsm_teff_output_filename -- '// & 'Unrecognized form value') - call LDT_endrun + call LDT_endrun endselect elseif(wstyle.eq."2 level hierarchy") then write(unit=cdate1, fmt='(i4.4, i2.2, i2.2, i2.2, i2.2)') & LDT_rc%yr, LDT_rc%mo, & LDT_rc%da, LDT_rc%hr,LDT_rc%mn - + dname = trim(odir)//'/' dname = trim(dname)//trim(mname)//'/' - + out_fname = trim(dname)//'LIS_HIST_'//trim(cdate1) - - write(unit=cdate, fmt='(a2,i2.2)') '.d',n + + write(unit=cdate, fmt='(a2,i2.2)') '.d',n out_fname = trim(out_fname)//trim(cdate) - + select case ( form ) case ("binary") - if(wopt.eq."1d tilespace") then + if(wopt.eq."1d tilespace") then out_fname = trim(out_fname)//'.ts4r' - elseif(wopt.eq."2d gridspace") then + elseif(wopt.eq."2d gridspace") then out_fname = trim(out_fname)//'.gs4r' - elseif(wopt.eq."1d gridspace") then + elseif(wopt.eq."1d gridspace") then out_fname = trim(out_fname)//'.gs4r' endif case ("grib1") @@ -428,9 +425,9 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & case default call ldt_log_msg('ERR: create_lsm_teff_output_filename -- '// & 'Unrecognized form value') - call LDT_endrun + call LDT_endrun endselect - elseif(wstyle.eq."WMO convention") then + elseif(wstyle.eq."WMO convention") then if ( .not. present(run_dd) .or. & .not. present(security_class) .or. & .not. present(distribution_class) .or. & @@ -439,15 +436,15 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & .not. present(write_interval) ) then call ldt_log_msg('ERR: create_lsm_teff_output_filename -- '// & 'missing WMO convention identifiers') - call LDT_endrun + call LDT_endrun endif write(unit=cdate1, fmt='(i4.4, i2.2, i2.2)') & LDT_rc%yr, LDT_rc%mo, LDT_rc%da - + write(unit=cdate, fmt='(i2.2, i2.2)') LDT_rc%hr, LDT_rc%mn - - if(map_proj.eq."polar") then + + if(map_proj.eq."polar") then fproj = 'P' print *,"fres ",run_dd(6) if (run_dd(6) .ge. 10.) then @@ -456,10 +453,9 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & write(unit=fres, fmt='(i1)') nint(run_dd(6)) endif fres2 = trim(fres)//'KM' - elseif(map_proj.eq."lambert") then + elseif(map_proj.eq."lambert") then fproj = 'L' print *,"fres ",run_dd(6) -! write(unit=fres, fmt='(f2.0)') run_dd(6) write(unit=fres, fmt='(f3.0)') run_dd(6) if (run_dd(6) .ge. 10.) then write(unit=fres, fmt='(i2)') nint(run_dd(6)) @@ -467,19 +463,19 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & write(unit=fres, fmt='(i1)') nint(run_dd(6)) endif fres2 = trim(fres)//'KM' - elseif(map_proj.eq."mercator") then + elseif(map_proj.eq."mercator") then fproj = 'M' write(unit=fres, fmt='(i2.2)') run_dd(6) fres = trim(fres)//'KM' - elseif(map_proj.eq."gaussian") then + elseif(map_proj.eq."gaussian") then fproj = 'G' - write(unit=fres, fmt='(i2.2)') run_dd(5)*100 + write(unit=fres, fmt='(i2.2)') run_dd(5)*100 fres2 = '0P'//trim(fres)//'DEG' else fproj = 'C' write(unit=fres, fmt='(i10)') nint(run_dd(6)*100) read(unit=fres,fmt='(10a1)') (fres1(i),i=1,10) - c = 0 + c = 0 do i=1,10 if(fres1(i).ne.' '.and.c==0) c = i enddo @@ -503,33 +499,44 @@ subroutine create_lsm_teff_output_filename(n, form, fname, odir, wstyle, wopt, & trim(fproj)//trim(fres2)//& '_AR.'//trim(area_of_data)//& '_PA.'//trim(write_interval)//'-HR-SUM_DD.'//& - trim(cdate1)//'_DT.'//trim(cdate)//'_DF.GR1' + trim(cdate1)//'_DT.'//trim(cdate)//'_DF' + if (form == "netcdf") then + out_fname = trim(out_fname) // ".nc" + else if (form == "grib1") then + out_fname = trim(out_fname) // ".GR1" + else if (form == "grib2") then + out_fname = trim(out_fname) // ".GR2" + else + write(LDT_logunit,*)'[ERR] Invalid LIS file format ', trim(form) + call LDT_endrun() + end if + endif fname = out_fname end subroutine create_lsm_teff_output_filename !BOP -! +! ! !ROUTINE: transformDataToLDTgrid_teff ! \label{transformDataToLDTgrid_teff} ! -! !INTERFACE: +! !INTERFACE: subroutine transformDataToLDTgrid_teff(n, teff_inp, teff_out) -! !USES: +! !USES: use LDT_coreMod use LISlsmTEFF_obsMod implicit none -! !ARGUMENTS: - integer :: n +! !ARGUMENTS: + integer :: n real :: teff_inp(lsmteffobs%nc*lsmteffobs%nr) real :: teff_out(LDT_rc%lnc(n),LDT_rc%lnr(n)) ! -! !DESCRIPTION: -! This routine interpolates or upscales the input data to +! !DESCRIPTION: +! This routine interpolates or upscales the input data to ! the LDT grid. If the input data is finer than the LDT ! grid, the input data is upscaled. If the input data is -! coarser, then it is interpolated to the LDT grid. +! coarser, then it is interpolated to the LDT grid. ! !EOP integer :: ios @@ -540,23 +547,23 @@ subroutine transformDataToLDTgrid_teff(n, teff_inp, teff_out) do r=1,lsmteffobs%nr do c=1, lsmteffobs%nc - if(teff_inp(c+(r-1)*lsmteffobs%nc).ne.LDT_rc%udef) then - teff_data_b(c+(r-1)*lsmteffobs%nc) = .true. + if(teff_inp(c+(r-1)*lsmteffobs%nc).ne.LDT_rc%udef) then + teff_data_b(c+(r-1)*lsmteffobs%nc) = .true. else teff_data_b(c+(r-1)*lsmteffobs%nc) = .false. endif - if(teff_inp(c+(r-1)*lsmteffobs%nc).lt.0) then + if(teff_inp(c+(r-1)*lsmteffobs%nc).lt.0) then teff_inp(c+(r-1)*lsmteffobs%nc) = LDT_rc%udef teff_data_b(c+(r-1)*lsmteffobs%nc) = .false. endif enddo enddo - if(LDT_isLDTatAfinerResolution(n,lsmteffobs%datares)) then + if(LDT_isLDTatAfinerResolution(n,lsmteffobs%datares)) then !-------------------------------------------------------------------------- ! Interpolate to the LDT running domain -!-------------------------------------------------------------------------- +!-------------------------------------------------------------------------- call bilinear_interp(LDT_rc%gridDesc(n,:),& teff_data_b, teff_inp, teffobs_b_ip, teffobs_ip, & lsmteffobs%nc*lsmteffobs%nr, & @@ -579,17 +586,17 @@ subroutine transformDataToLDTgrid_teff(n, teff_inp, teff_out) lsmteffobs%nc*lsmteffobs%nr,& LDT_rc%lnc(n)*LDT_rc%lnr(n),LDT_rc%udef, & lsmteffobs%n11,teff_data_b, teff_inp, teffobs_b_ip,teffobs_ip) - + endif - + do r=1,LDT_rc%lnr(n) do c=1,LDT_rc%lnc(n) - if(teffobs_b_ip(c+(r-1)*LDT_rc%lnc(n))) then + if(teffobs_b_ip(c+(r-1)*LDT_rc%lnc(n))) then teff_out(c,r) = teffobs_ip(c+(r-1)*LDT_rc%lnc(n)) else teff_out(c,r) = LDT_rc%udef endif enddo enddo - + end subroutine transformDataToLDTgrid_teff diff --git a/ldt/SMAP_E_OPL/ARFSSMRETRIEVAL.F90 b/ldt/SMAP_E_OPL/ARFSSMRETRIEVAL.F90 index 137324b4f..3edef91d9 100644 --- a/ldt/SMAP_E_OPL/ARFSSMRETRIEVAL.F90 +++ b/ldt/SMAP_E_OPL/ARFSSMRETRIEVAL.F90 @@ -7,17 +7,21 @@ ! REVISION HISTORY: ! 22 Feb 2022: P.W.LIU; Initial implemetation ! 22 Feb 2022: Yonghwan Kwon; modified for LDT +! 10 Feb 2023: Eric Kemp, modified to output retrievals in netCDF. +! 21 Feb 2023: Eric Kemp, added third LIS time level. ! ! DESCRIPTION: RETRIEVE SMAP SM FOR ARFS -! INPUT : SMAP - L1B Brightness Temperature +! INPUT : SMAP - L1B Brightness Temperature ! OUTPUT: SMAPTB_ARFSGRIDE_ddmmyyy.dat ! NOTES : Inverse Distance Squared with 0.4 deg serching window !------------------------------------------------------------------------- - subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& - ARFS_SNOW,DOY,UTChr) +subroutine ARFSSMRETRIEVAL(SMAPFILE, & + TS_bfresample_01, TS_bfresample_02, TS_bfresample_03, & + ARFS_SNOW, DOY, UTChr, firsttime, secondtime, thirdtime) - !USE HDF5 + !USE HDF5 + use esmf USE VARIABLES USE DATADOMAIN USE FUNCTIONS @@ -25,15 +29,20 @@ subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& USE invdist_temp2smap USE varsio_m USE algo_vpol_m - use LDT_logMod + use LDT_ARFSSM_netcdfMod, only: LDT_ARFSSM_write_netcdf + use LDT_logMod, only: LDT_logunit USE LDT_smap_e_oplMod - + IMPLICIT NONE ! !ARGUMENTS: CHARACTER (len=100) :: SMAPFILE - REAL*4, DIMENSION(2560,1920) :: TS_bfresample_01, TS_bfresample_02 + REAL*4, DIMENSION(2560,1920), intent(in) :: TS_bfresample_01, & + TS_bfresample_02, TS_bfresample_03 REAL*4, DIMENSION(2560,1920) :: ARFS_SNOW, UTChr INTEGER :: DOY + type(ESMF_Time), intent(in) :: firsttime + type(ESMF_Time), intent(in) :: secondtime + type(ESMF_Time), intent(in) :: thirdtime !EOP INTEGER :: i, j, nrow, mcol CHARACTER (len=100) :: fname_TAU @@ -43,7 +52,7 @@ subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& REAL*4, DIMENSION(2560,1920) :: ARFS_TAU, ARFS_CLAY, ARFS_BD, ARFS_OMEGA, ARFS_H INTEGER*1, DIMENSION(2560,1920) :: ARFS_LC, ARFS_SM_FLAG INTEGER*1 :: retrieval_flag - REAL*4, DIMENSION(2560,1920) :: ARFS_TS_01, ARFS_TS_02, ARFS_SM + REAL*4, DIMENSION(2560,1920) :: ARFS_TS_01, ARFS_TS_02, ARFS_TS_03, ARFS_SM REAL*8 ,DIMENSION(:), ALLOCATABLE :: ARFS_FINE_LAT, ARFS_FINE_LON REAL*8 ,DIMENSION(:), ALLOCATABLE :: ARFS_LAT, ARFS_LON REAL :: T1, T2 @@ -55,6 +64,16 @@ subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& integer :: L1B_dir_len,L1B_fname_len real :: utc_check + ! EMK + character(8) :: yyyymmdd + character(6) :: hhmmss + real :: deltasec, wgt + integer :: firstUTCyr, firstUTCmo, firstUTCdy, firstUTChr + integer :: secondUTCyr, secondUTCmo, secondUTCdy, secondUTChr + integer :: thirdUTCyr, thirdUTCmo, thirdUTCdy, thirdUTChr + + real :: TS_A, TS_B + nrow=2560 mcol=1920 @@ -70,6 +89,7 @@ subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& ARFS_FINE_LON = LON(arfs_geo_lon_lf,arfs_geo_lon_rt,arfs_lon_3km_space) CALL RESAMPLETEMP(TS_bfresample_01,ARFS_LAT,ARFS_LON,ARFS_FINE_LAT,ARFS_FINE_LON,ARFS_TS_01) CALL RESAMPLETEMP(TS_bfresample_02,ARFS_LAT,ARFS_LON,ARFS_FINE_LAT,ARFS_FINE_LON,ARFS_TS_02) + CALL RESAMPLETEMP(TS_bfresample_03,ARFS_LAT,ARFS_LON,ARFS_FINE_LAT,ARFS_FINE_LON,ARFS_TS_03) ! IF EVENTUALLY THE RESAMPLING DOES NOT CHANGE TEFF MUCH WE COULD SIMPLELY USE ARFS_TS=TS_bfresample ! UP TO HERE TAKES 38 SECS write (LDT_logunit,*) '[INFO] Finished resampling effective soil temperature' @@ -121,22 +141,56 @@ subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& write (LDT_logunit,*) '[INFO] Finished reading landcover' !generate soil moisture retrievals - write (LDT_logunit,*) '[INFO] Generating soil moisture retrievals' + write (LDT_logunit,*) '[INFO] Generating soil moisture retrievals' ARFS_SM=-9999 ARFS_SM_FLAG=-1 - DO i=1,nrow !ROW LON - DO j=1,mcol !COL LAT - tbv = ARFS_TB(i,j) - if(UTChr(i,j).ge.0) then - utc_check = UTChr(i,j) - floor(UTChr(i,j)) + call ESMF_TimeGet(firsttime, yy=firstUTCyr, mm=firstUTCmo, dd=firstUTCdy, & + h=firstUTChr) + call ESMF_TimeGet(secondtime, yy=secondUTCyr, mm=secondUTCmo, & + dd=secondUTCdy, & + h=secondUTChr) + call ESMF_TimeGet(thirdtime, yy=thirdUTCyr, mm=thirdUTCmo, dd=thirdUTCdy, & + h=thirdUTChr) - if(utc_check.le.0.5) then - Ts = ARFS_TS_01(i,j) - else - Ts = ARFS_TS_02(i,j) - endif - endif + DO j=1,mcol !COL LAT + DO i=1,nrow !ROW LON + + tbv = ARFS_TB(i,j) + + if (UTChr(i,j) < 0) cycle + + if (UTChr(i,j) == firstUTChr) then + TS_A = ARFS_TS_01(i,j) + TS_B = ARFS_TS_02(i,j) + wgt = 1 + else if (UTChr(i,j) > firstUTChr .and. & + UTChr(i,j) < secondUTChr) then + TS_A = ARFS_TS_01(i,j) + TS_B = ARFS_TS_02(i,j) + deltasec = ( UTChr(i,j) - firstUTChr ) * 3600 + wgt = (10800. - deltasec) / 10800. + else if (UTChr(i,j) > firstUTChr .and. & + firstUTChr == 21 .and. secondUTChr == 0) then + TS_A = ARFS_TS_01(i,j) + TS_B = ARFS_TS_02(i,j) + deltasec = ( UTChr(i,j) - firstUTChr ) * 3600 + wgt = (10800. - deltasec) / 10800. + else if (UTChr(i,j) == secondUTChr) then + TS_A = ARFS_TS_02(i,j) + TS_B = ARFS_TS_03(i,j) + wgt = 1 + else + TS_A = ARFS_TS_02(i,j) + TS_B = ARFS_TS_03(i,j) + deltasec = ( UTChr(i,j) - secondUTChr ) * 3600 + wgt = (10800. - deltasec) / 10800. + end if + if (TS_A > 0 .and. TS_B > 0) then + TS = ((wgt)*TS_A) + ((1. - wgt)*TS_B) + else + cycle + end if IF (tbv.GT.0.0.AND.Ts.GT.0.AND.ARFS_SNOW(i,j).LE.SMAPeOPL%SD_thold.AND.ARFS_BD(i,j).NE.-9999.AND.ARFS_LC(i,j).NE.0.AND.& UTChr(i,j).GE.0) THEN @@ -150,11 +204,10 @@ subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& CALL algo_vpol(real(i),real(j),sm_retrieval, tau_return, retrieval_flag) ARFS_SM(i,j)=sm_retrieval ARFS_SM_FLAG(i,j)=retrieval_flag - ELSE - !PRINT*,i, j, "NO RETRIEVAL" + END IF - END DO !jj=1,mcol !COL LAT - END DO !ii=1,nrow !ROW LON + END DO !ii=1,nrow !ROW LON + END DO !jj=1,mcol !COL LAT !write soil moisture retrieval outputs L1B_dir_len = len_trim(SMAPeOPL%L1Bdir) @@ -162,16 +215,22 @@ subroutine ARFSSMRETRIEVAL(SMAPFILE,TS_bfresample_01,TS_bfresample_02,& if(SMAPeOPL%L1Btype.eq.1) then !NRT retrieval_fname = trim(SMAPeOPL%SMoutdir)//"/"//"ARFS_SM_V_"//& - trim(SMAPFILE(L1B_dir_len+18:L1B_fname_len-3))//".dat" + trim(SMAPFILE(L1B_dir_len+18:L1B_fname_len-3))//".nc" + yyyymmdd = trim(SMAPFILE(L1B_fname_len-28:L1B_fname_len-20)) + hhmmss = trim(SMAPFILE(L1B_fname_len-19:L1B_fname_len-13)) elseif(SMAPeOPL%L1Btype.eq.2) then !Historical retrieval_fname = trim(SMAPeOPL%SMoutdir)//"/"//"ARFS_SM_V_"//& - trim(SMAPFILE(L1B_dir_len+14:L1B_fname_len-3))//".dat" + trim(SMAPFILE(L1B_dir_len+14:L1B_fname_len-3))//".nc" + yyyymmdd = trim(SMAPFILE(L1B_fname_len-28:L1B_fname_len-20)) + hhmmss = trim(SMAPFILE(L1B_fname_len-19:L1B_fname_len-13)) endif write (LDT_logunit,*) '[INFO] Writing soil moisture retrieval file ', trim(retrieval_fname) - OPEN(UNIT=151, FILE=retrieval_fname,FORM='UNFORMATTED',ACCESS='DIRECT', RECL=nrow*mcol*4) - WRITE(UNIT=151, REC = 1) ARFS_SM - CLOSE(151) + + ! NOTE: nrow is actually number of columns, mcol is actually number of + ! rows + call LDT_ARFSSM_write_netcdf(nrow, mcol, arfs_sm, retrieval_fname, & + yyyymmdd, hhmmss) write (LDT_logunit,*) '[INFO] Successfully wrote soil moisture retrieval file ', trim(retrieval_fname) write (LDT_logunit,*) '[INFO] Finished generating soil moisture retrievals' diff --git a/ldt/SMAP_E_OPL/LDT_ARFSSM_netcdfMod.F90 b/ldt/SMAP_E_OPL/LDT_ARFSSM_netcdfMod.F90 new file mode 100644 index 000000000..5163841cf --- /dev/null +++ b/ldt/SMAP_E_OPL/LDT_ARFSSM_netcdfMod.F90 @@ -0,0 +1,318 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.4 +! +! Copyright (c) 2022 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT----------------------- +! +! MODULE: LDT_ARFSSM_netcdfMod +! +! REVISION HISTORY: +! 10 Feb 2023: Eric Kemp. Initial implementation. +! +! DESCRIPTION: Outputs SMAP based soil moisture retrieval in netCDF format. +! +!------------------------------------------------------------------------------ + +#include "LDT_misc.h" +#include "LDT_NetCDF_inc.h" + +module LDT_ARFSSM_netcdfMod + + ! Defaults + implicit none + private + + ! Public routines + public :: LDT_ARFSSM_write_netcdf + +contains + +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + ! Subroutine for writing ARFS SM retrieval to netCDF + subroutine LDT_ARFSSM_write_netcdf(nc, nr, arfs_sm, retrieval_fname, & + yyyymmdd, hhmmss) + + ! Imports + use LDT_coreMod, only: LDT_rc, LDT_masterproc + use LDT_logMod, only: LDT_logunit, LDT_endrun, LDT_verify +#if ( defined SPMD) + use mpi +#endif +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + ! Arguments + integer, intent(in) :: nc + integer, intent(in) :: nr + real*4, intent(in) :: arfs_sm(nc,nr) + character(*), intent(in) :: retrieval_fname + character(8), intent(in) :: yyyymmdd + character(6), intent(in) :: hhmmss + + ! Locals + integer :: shuffle, deflate, deflate_level + integer :: iret + integer :: ncid + integer :: dim_ids(3) + real :: dlat, dlon + real :: swlat, swlon + real :: nelat, nelon + integer :: lon_varid, lat_varid, time_varid, sm_varid + character(120) :: time_units + character(8) :: date + character(10) :: time + character(5) :: zone + integer :: values(8) + real, allocatable :: lats(:) + real, allocatable :: lons(:) + integer :: ierr + integer :: i, j + + ! Only the master process handles the file output + if (LDT_masterproc) then + write(LDT_logunit,*)'[INFO] Creating NETCDF file ', & + trim(retrieval_fname) + + ! Copy netCDF compression settings + shuffle = NETCDF_shuffle + deflate = NETCDF_deflate + deflate_level = NETCDF_deflate_level + + ! Create the output file +#if(defined USE_NETCDF3) + iret=nf90_create(path=trim(retrieval_fname),& + cmode=NF90_CLOBBER, ncid=ncid) + call LDT_verify(iret,& + '[ERR] nf90_create failed') +#endif +#if(defined USE_NETCDF4) + iret=nf90_create(path=trim(retrieval_fname),& + cmode=NF90_NETCDF4, ncid=ncid) + call LDT_verify(iret, & + '[ERR] nf90_create failed') +#endif + + ! Write out dimensions headers + call LDT_verify(nf90_def_dim(ncid, 'time', 1, dim_ids(3)), & + '[ERR] nf90_def_dim failed') + call LDT_verify(nf90_def_dim(ncid, 'lat', nr, dim_ids(2)), & + '[ERR] nf90_def_dim failed') + call LDT_verify(nf90_def_dim(ncid, 'lon', nc, dim_ids(1)), & + '[ERR] nf90_def_dim failed') + + ! Map projection + !FIXME: Allow map projections other than lat/lon + select case("latlon") + case ("latlon") + dlon = LDT_rc%gridDesc(1,9) + dlat = LDT_rc%gridDesc(1,10) + swlat = LDT_rc%gridDesc(1,4) + swlon = LDT_rc%gridDesc(1,5) + nelat = LDT_rc%gridDesc(1,7) + nelon = LDT_rc%gridDesc(1,8) + + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL,& + "MAP_PROJECTION", "EQUIDISTANT CYLINDRICAL"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL,& + "SOUTH_WEST_CORNER_LAT", swlat), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL,& + "SOUTH_WEST_CORNER_LON", swlon), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, & + "NORTH_EAST_CORNER_LAT", nelat), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, & + "NORTH_EAST_CORNER_LON", nelon), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, & + "DX", dlon),& + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, & + "DY", dlat), & + '[ERR] nf90_put_att failed') + + case default + write(LDT_logunit,*) & + '[ERR] Only latlon map projection supported for SMAP_E_OPL' + call LDT_endrun() + end select + + ! Include the water points + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, & + "INC_WATER_PTS", "true"), & + '[ERR] nf90_put_att failed') + + ! Construct the longitudes + ! FIXME: Add support for other map projections + call LDT_verify(nf90_def_var(ncid, "lon", NF90_FLOAT, dim_ids(1), & + lon_varid),'[ERR] nf90_def_var failed') +#if(defined USE_NETCDF4) + call LDT_verify(nf90_def_var_deflate(ncid, & + lon_varid, shuffle, deflate, deflate_level),& + '[ERR] nf90_def_var_deflate') +#endif + call LDT_verify(nf90_put_att(ncid, lon_varid, & + "units", "degrees_east"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, lon_varid, & + "long_name", "longitude"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, lon_varid, & + "standard_name", "longitude"), & + '[ERR] nf90_put_att failed') + + ! Construct the latitudes + ! FIXME: Add support for other map projections + call LDT_verify(nf90_def_var(ncid, "lat", NF90_FLOAT, dim_ids(2), & + lat_varid), '[ERR] nf90_def_var failed') +#if(defined USE_NETCDF4) + call LDT_verify(nf90_def_var_deflate(ncid, & + lat_varid, shuffle, deflate, deflate_level),& + '[ERR] nf90_def_var_deflate') + call LDT_verify(nf90_put_att(ncid, lat_varid, & + "units", "degrees_north"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, lat_varid, & + "long_name", "latitude"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, lat_varid, & + "standard_name", "latitude"), & + '[ERR] nf90_put_att failed') +#endif + + ! Define the time array. The valid time will be written as an + ! attribute. + call LDT_verify(nf90_def_var(ncid, 'time', NF90_DOUBLE, & + dimids=dim_ids(3), varid=time_varid), & + '[ERR] nf90_def_var failed') + write(time_units,'(A)') & + "seconds since "//yyyymmdd(1:4)//"-" & + //yyyymmdd(5:6)//"-" & + //yyyymmdd(7:8)//" " & + //hhmmss(1:2)//":"//hhmmss(3:4)//":"//hhmmss(5:6) + call LDT_verify(nf90_put_att(ncid, time_varid, & + "units", trim(time_units)), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, time_varid, & + "long_name", "time"), & + '[ERR] nf90_put_att failed') + + ! Define the soil moisture retrieval + call LDT_verify(nf90_def_var(ncid, "arfs_sm", NF90_FLOAT, & + dimids=dim_ids, & + varid=sm_varid), '[ERR] nf90_def_var failed') +#if (defined USE_NETCDF4) + call LDT_verify(nf90_def_var_deflate(ncid, & + sm_varid, shuffle, deflate, deflate_level), & + '[ERR] nf90_def_var_deflate failed') +#endif + call LDT_verify(nf90_put_att(ncid, sm_varid, & + "units", "m^3/m^3"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, sm_varid, & + "long_name","soil moisture"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, sm_varid, & + '_FillValue', -9999.), & + '[ERR] nf90_put_att failed') + + ! Miscellaneous header information + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, "Conventions", & + "CF-1.10"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, "title", & + "LDT SMAP_E_OPL retrieval"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, "institution", & + "NASA GSFC Hydrological Sciences Laboratory"), & + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, "source", & + "Land Data Toolkit (LDT)"), & + '[ERR] nf90_put_att failed') + call date_and_time(date, time, zone, values) +#ifndef LDT_SKIP_HISTORY + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, "history", & + "created on date: "//date(1:4)//"-"//date(5:6)//"-"//& + date(7:8)//"T"//time(1:2)//":"//time(3:4)//":"//time(5:10)), & + '[ERR] nf90_put_att failed') +#endif + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, "references", & + "Arsenault_etal_GMD_2018, Kumar_etal_EMS_2006"),& + '[ERR] nf90_put_att failed') + call LDT_verify(nf90_put_att(ncid, NF90_GLOBAL, "comment", & + "website: http://lis.gsfc.nasa.gov/"), & + '[ERR] nf90_put_att failed') + + ! We are ready to write the actual data. This requires taking NETCDF + ! out of define mode. + call LDT_verify(nf90_enddef(ncid), & + '[ERR] ncf90_enddef failed') + + ! Write the latitude data + allocate(lats(nr)) + do j = 1, nr + lats(j) = swlat + (j-1)*(dlat) + end do + call LDT_verify(nf90_put_var(ncid, lat_varid, & + lats(:), (/1/), (/nr/)), & + '[ERR] nf90_put_var failed for lats') + deallocate(lats) + + ! Write the longitude data + allocate(lons(nc)) + do i = 1, nc + lons(i) = swlon + (i-1)*(dlon) + end do + call LDT_verify(nf90_put_var(ncid, lon_varid, lons(:),& + (/1/), (/nc/)), & + '[ERR] nf90_put_var failed for lon') + deallocate(lons) + + ! Write the time data + call LDT_verify(nf90_put_var(ncid, time_varid, 0.0), & + '[ERR] nf90_put_var failed for time') + + ! Write the ARFS SM field + call LDT_verify(nf90_put_var(ncid, sm_varid, & + arfs_sm(:,:), & + (/1,1/), (/nc,nr/)), & + '[ERR] nf90_put_var failed for sm') + + ! Close the file and clean up + call LDT_verify(nf90_close(ncid), & + '[ERR] nf90_close failed!') + endif + + ierr = 0 +#if (defined SPMD) + call mpi_barrier(MPI_COMM_WORLD, ierr) +#endif + end subroutine LDT_ARFSSM_write_netcdf +#else + ! Dummy version + subroutine LDT_ARFSSM_write_netcdf(nc, nr, arfs_sm, retrieval_fname, & + yyyymmdd, hhmmss) + use LDT_logMod, only: LDT_logunit, LDT_endrun + implicit none + integer, intent(in) :: nc + integer, intent(in) :: nr + real*4, intent(in) :: arfs_sm(nc,nr) + character(*), intent(in) :: retrieval_fname + character(8), intent(in) :: yyyymmdd + character(6), intent(in) :: hhmmss + write(LDT_logunit,*)'[ERR] LDT not compiled with NETCDF support!' + write(LDT_logunit,*)'Cannot write ARFS SM retrieval in NETCDF format!' + write(LDT_logunit,*)'Recompile with NETCDF support and try again!' + call LDT_endrun() + end subroutine LDT_ARFSSM_write_netcdf +#endif + +end module LDT_ARFSSM_netcdfMod + diff --git a/ldt/SMAP_E_OPL/LDT_smap_e_oplMod.F90 b/ldt/SMAP_E_OPL/LDT_smap_e_oplMod.F90 index d32c497fe..f20b9ed67 100644 --- a/ldt/SMAP_E_OPL/LDT_smap_e_oplMod.F90 +++ b/ldt/SMAP_E_OPL/LDT_smap_e_oplMod.F90 @@ -15,6 +15,12 @@ ! ! !REVISION HISTORY: ! 14 Dec 2021: Yonghwan Kwon, Initial Specification +! 06 Feb 2023: Eric Kemp, now process subset of SMAP fields. +! 14 Feb 2023: Eric Kemp, now uses USAFSI and USAF LIS output. +! 22 Feb 2023: Eric Kemp, ensemble size now in ldt.config file. +! 01 Jul 2023: Mahdi Navari,This edit generates a separate SMAP_filelist +! for each LDT job based on user input. +! Now we can run several LDT jobs in the same directory. ! #include "LDT_misc.h" #include "LDT_NetCDF_inc.h" @@ -38,7 +44,7 @@ module LDT_smap_e_oplMod CLAYfile, Hfile, LCfile character*100 :: dailystats_ref, dailystats_lis character*10 :: date_curr - integer :: L1BresampWriteOpt, L1Btype + integer :: L1BresampWriteOpt, L1Btype, SMAPfilelistSuffixNumber integer :: Teffscale integer :: ntimes,ngrid real, allocatable :: mu_6am_ref(:), mu_6pm_ref(:) !(ngrid) @@ -48,7 +54,9 @@ module LDT_smap_e_oplMod integer, allocatable :: grid_col(:), grid_row(:) !(ngrid) real, allocatable :: ARFS_TBV_COR(:,:) real :: SD_thold - + integer :: num_ens ! Number of ensemble members in LIS USAF file. + integer :: num_tiles ! Total number of tiles in LIS USAF file. + integer :: ntiles_pergrid ! Number of tiles per grid point end type smap_e_opl_dec type(smap_e_opl_dec), public :: SMAPeOPL @@ -61,7 +69,7 @@ subroutine LDT_smap_e_oplInit() ! Imports use ESMF use LDT_coreMod, only: LDT_config - use LDT_logMod, only: LDT_verify + use LDT_logMod, only: LDT_logunit, LDT_endrun, LDT_verify ! Defaults implicit none @@ -103,6 +111,12 @@ subroutine LDT_smap_e_oplInit() call ESMF_ConfigGetAttribute(LDT_config, SMAPeOPL%L1BresampWriteOpt, rc=rc) call LDT_verify(rc, trim(cfg_entry)//" not specified") + cfg_entry = "SMAP_E_OPL filelist suffix number:" + call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + call ESMF_ConfigGetAttribute(LDT_config, SMAPeOPL%SMAPfilelistSuffixNumber, rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + if(SMAPeOPL%L1BresampWriteOpt.eq.1) then cfg_entry = "SMAP_E_OPL L1B resampled output directory:" call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) @@ -143,6 +157,42 @@ subroutine LDT_smap_e_oplInit() call ESMF_ConfigGetAttribute(LDT_config, SMAPeOPL%LISsnowdir, rc=rc) call LDT_verify(rc, trim(cfg_entry)//" not specified") + cfg_entry = "SMAP_E_OPL LIS ensemble size:" + call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + call ESMF_ConfigGetAttribute(LDT_config, SMAPeOPL%num_ens, rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + if (SMAPeOPL%num_ens < 1) then + write(LDT_logunit,*)'[ERR] LIS ensemble size must be at least 1!' + write(LDT_logunit,*)'[ERR] Read in ', SMAPeOPL%num_ens + call LDT_endrun() + end if + + cfg_entry = "SMAP_E_OPL LIS total number of tiles (including ensembles):" + call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + call ESMF_ConfigGetAttribute(LDT_config, SMAPeOPL%num_tiles, rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + if (SMAPeOPL%num_tiles < 1) then + write(LDT_logunit,*) & + '[ERR] LIS total number of tiles (including ensembles) must be' & + //'at least 1!' + write(LDT_logunit,*)'[ERR] Read in ', SMAPeOPL%num_tiles + call LDT_endrun() + end if + + cfg_entry = "SMAP_E_OPL LIS number of tiles per grid point:" + call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + call ESMF_ConfigGetAttribute(LDT_config, SMAPeOPL%ntiles_pergrid, rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + if (SMAPeOPL%num_tiles < 1) then + write(LDT_logunit,*) & + '[ERR] LIS number of tiles per grid point must be at least 1!' + write(LDT_logunit,*)'[ERR] Read in ', SMAPeOPL%ntiles_pergrid + call LDT_endrun() + end if + cfg_entry = "SMAP_E_OPL snow depth threshold:" call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) call LDT_verify(rc, trim(cfg_entry)//" not specified") @@ -192,8 +242,9 @@ subroutine LDT_smap_e_oplRun(n) ! This calls the actual SMAP_E_OPL driver ! !USES: - use LDT_logMod + use esmf use LDT_coreMod + use LDT_logMod implicit none ! !ARGUMENTS: @@ -205,28 +256,42 @@ subroutine LDT_smap_e_oplRun(n) integer :: ftn, ierr character*100 :: fname character*100 :: smap_L1B_filename(10) - character*8 :: yyyymmdd, yyyymmdd_01, yyyymmdd_02 + character*8 :: yyyymmdd, yyyymmdd_01, yyyymmdd_02, yyyymmdd_03 character*6 :: hhmmss(10) - character*4 :: yyyy, yyyy_01, yyyy_02 + character*4 :: yyyy, yyyy_01, yyyy_02, yyyy_03 character*2 :: hh, mm, dd - character*2 :: hh_01, mm_01, dd_01 + character*2 :: hh_01, mm_01, dd_01 character*2 :: hh_02, mm_02, dd_02 + character*2 :: hh_03, mm_03, dd_03 + character*2 :: tmp character*1 :: Orbit integer :: yr, mo, da, hr - integer :: yr_pre, mo_pre, da_pre + integer :: yr_pre, mo_pre, da_pre, hh_pre integer :: yr_02, mo_02, da_02, hr_02 + integer :: yr_03, mo_03, da_03, hr_03 logical :: dir_exists, read_L1Bdata real :: teff_01(LDT_rc%lnc(n),LDT_rc%lnr(n)) real :: teff_02(LDT_rc%lnc(n),LDT_rc%lnr(n)) + real :: teff_03(LDT_rc%lnc(n),LDT_rc%lnr(n)) real :: SnowDepth(LDT_rc%lnc(n),LDT_rc%lnr(n)) real :: TIMEsec(LDT_rc%lnc(n),LDT_rc%lnr(n)) real :: UTChr(LDT_rc%lnc(n),LDT_rc%lnr(n)) integer :: L1B_dir_len integer :: doy_pre, doy_curr + type(ESMF_Calendar) :: calendar + type(ESMF_Time) :: firsttime, secondtime, thirdtime, curtime, prevdaytime + type(ESMF_TimeInterval) :: deltatime + integer :: deltahr + integer :: rc + integer :: col, row + external :: readUSAFSI + external :: readLIS_Teff_usaf + + allocate(LDT_rc%nensem(LDT_rc%nnest)) - ! Resample SMAP L1B to L1C + ! Resample SMAP L1B to L1C call search_SMAPL1B_files(SMAPeOPL%L1Bdir,SMAPeOPL%date_curr,& - SMAPeOPL%L1Btype) + SMAPeOPL%L1Btype, SMAPeOPL%SMAPfilelistSuffixNumber) yyyymmdd = SMAPeOPL%date_curr(1:8) yyyy = SMAPeOPL%date_curr(1:4) @@ -242,10 +307,11 @@ subroutine LDT_smap_e_oplRun(n) trim(SMAPeOPL%L1Bresampledir_02)) endif + write (tmp,'(I2.2)') SMAPeOPL%SMAPfilelistSuffixNumber + ftn = LDT_getNextUnitNumber() - open(unit=ftn, file='./SMAP_L1B_filelist.dat',& + open(unit=ftn, file='./SMAP_L1B_filelist_'//tmp//'.dat',& status='old', iostat=ierr) - fi = 0 do while (ierr .eq. 0) read (ftn, '(a)', iostat=ierr) fname @@ -268,15 +334,29 @@ subroutine LDT_smap_e_oplRun(n) if(i == fi) then write (LDT_logunit,*) '[INFO] Resampling ', trim(smap_L1B_filename(i)) allocate(SMAPeOPL%ARFS_TBV_COR(LDT_rc%lnc(n),LDT_rc%lnr(n))) - call SMAPL1BRESAMPLE(smap_L1B_filename(i),SMAPeOPL%L1Bdir,Orbit,TIMEsec) - write (LDT_logunit,*) '[INFO] Finished resampling ', trim(smap_L1B_filename(i)) - read_L1Bdata = .true. + ! EMK...Process subset of fields. + call SMAPL1BRESAMPLE_subset(smap_L1B_filename(i), & + SMAPeOPL%L1Bdir, Orbit, TIMEsec, rc) + + if (rc == 0) then + write (LDT_logunit,*) '[INFO] Finished resampling ', trim(smap_L1B_filename(i)) + read_L1Bdata = .true. + else + deallocate(SMAPeOPL%ARFS_TBV_COR) + end if elseif(hhmmss(i) /= hhmmss(i+1)) then write (LDT_logunit,*) '[INFO] Resampling ', trim(smap_L1B_filename(i)) allocate(SMAPeOPL%ARFS_TBV_COR(LDT_rc%lnc(n),LDT_rc%lnr(n))) - call SMAPL1BRESAMPLE(smap_L1B_filename(i),SMAPeOPL%L1Bdir,Orbit,TIMEsec) - write (LDT_logunit,*) '[INFO] Finished resampling ', trim(smap_L1B_filename(i)) - read_L1Bdata = .true. + !EMK Process subset of fields. + call SMAPL1BRESAMPLE_subset(smap_L1B_filename(i), & + SMAPeOPL%L1Bdir, Orbit, TIMEsec, rc) + + if (rc == 0) then + write (LDT_logunit,*) '[INFO] Finished resampling ', trim(smap_L1B_filename(i)) + read_L1Bdata = .true. + else + deallocate(SMAPeOPL%ARFS_TBV_COR) + end if endif if(read_L1Bdata) then @@ -285,105 +365,94 @@ subroutine LDT_smap_e_oplRun(n) ! use LIS outputs from previous day read(yyyy,*,iostat=ierr) yr read(mm,*,iostat=ierr) mo - read(dd,*,iostat=ierr) da + read(dd,*,iostat=ierr) da read(hh,*,iostat=ierr) hr - yr_pre = yr - mo_pre = mo - da_pre = da - 1 - if(da_pre.eq.0) then - mo_pre = mo - 1 - - if(mo_pre.eq.0) then - yr_pre = yr - 1 - mo_pre = 12 - da_pre = 31 - else - if(mo_pre.eq.1.or.& - mo_pre.eq.3.or.& - mo_pre.eq.5.or.& - mo_pre.eq.7.or.& - mo_pre.eq.8.or.& - mo_pre.eq.10.or.& - mo_pre.eq.12) then - da_pre = 31 - elseif(mo_pre.eq.2) then - if(mod(yr,4).eq.0) then - da_pre = 29 - else - da_pre = 28 - endif - else - da_pre = 30 - endif - endif - endif + calendar = ESMF_CalendarCreate(ESMF_CALKIND_GREGORIAN, & + name="Gregorian", & + rc=rc) + + ! Set current time + call ESMF_TimeSet(curtime, yy=yr, mm=mo, dd=da, h=hr, m=0, s=0, & + calendar=calendar, rc=rc) + call LDT_verify(rc, '[ERR] in ESMF_TimeSet in LDT_smap_e_oplRun') + ! Go back 24 hours + call ESMF_TimeIntervalSet(deltatime, d=1, rc=rc) + call LDT_verify(rc, & + '[ERR] in ESMF_TimeIntervalSet in LDT_smap_e_oplRun') + prevdaytime = curtime - deltatime + + ! Now, find the nearest 3-hrly time (00Z, 03Z, ..., 21Z) prior + ! to prevdaytime + if (mod(hr, 3) == 0) then + deltahr = 0 + else + deltahr = hr - ((floor(real(hr)/3.))*3) + end if + call ESMF_TimeIntervalSet(deltatime, h=deltahr, rc=rc) + call LDT_verify(rc, & + '[ERR] in ESMF_TimeIntervalSet in LDT_smap_e_oplRun') + firsttime = prevdaytime - deltatime + + ! Now, find the next 3-hrly time (00Z, 03Z, ..., 21Z) after + ! firsttime + call ESMF_TimeIntervalSet(deltatime, h=3, rc=rc) + call LDT_verify(rc, & + '[ERR] in ESMF_TimeIntervalSet in LDT_smap_e_oplRun') + secondtime = firsttime + deltatime + + ! Now, find the next 3-hrly time (00Z, 03Z, ..., 21Z) after + ! secondtime + call ESMF_TimeIntervalSet(deltatime, h=3, rc=rc) + call LDT_verify(rc, & + '[ERR] in ESMF_TimeIntervalSet in LDT_smap_e_oplRun') + thirdtime = secondtime + deltatime + + ! Now, read the first time. + call ESMF_TimeGet(firsttime, yy=yr_pre, mm=mo_pre, dd=da_pre, & + h=hh_pre) write(unit=yyyy_01, fmt='(i4.4)') yr_pre write(unit=mm_01, fmt='(i2.2)') mo_pre write(unit=dd_01, fmt='(i2.2)') da_pre yyyymmdd_01 = trim(yyyy_01)//trim(mm_01)//trim(dd_01) - hh_01 = hh - - call readLIS_Teff(n,yyyymmdd_01,hh_01,Orbit,teff_01) - - yr_02 = yr_pre - mo_02 = mo_pre - da_02 = da_pre - hr_02 = hr + 1 - - if(hr_02.eq.24) then - hr_02 = 0 - da_02 = da_pre + 1 - - if(mo_pre.eq.1.or.& - mo_pre.eq.3.or.& - mo_pre.eq.5.or.& - mo_pre.eq.7.or.& - mo_pre.eq.8.or.& - mo_pre.eq.10.or.& - mo_pre.eq.12) then - if(da_02.gt.31) then - da_02 = 1 - mo_02 = mo_pre + 1 - endif - elseif(mo_pre.eq.2) then - if(mod(yr_02,4).eq.0) then - if(da_02.gt.29) then - da_02 = 1 - mo_02 = mo_pre + 1 - endif - else - if(da_02.gt.28) then - da_02 = 1 - mo_02 = mo_pre + 1 - endif - endif - else - if(da_02.gt.30) then - da_02 = 1 - mo_02 = mo_pre + 1 - endif - endif - - if(mo_02.gt.12) then - mo_02 = 1 - yr_02 = yr_pre + 1 - endif + write(unit=hh_01, fmt='(i2.2)') hh_pre + call readLIS_Teff_usaf(n, yyyymmdd_01, hh_01, Orbit, teff_01, rc) + if (rc .ne. 0) then + write(LDT_logunit,*)'[WARN] No Teff data available...' endif + ! Now, read the second time. + call ESMF_TimeGet(secondtime, yy=yr_02, mm=mo_02, dd=da_02, & + h=hr_02) + write(unit=yyyy_02, fmt='(i4.4)') yr_02 write(unit=mm_02, fmt='(i2.2)') mo_02 write(unit=dd_02, fmt='(i2.2)') da_02 write(unit=hh_02, fmt='(i2.2)') hr_02 yyyymmdd_02 = trim(yyyy_02)//trim(mm_02)//trim(dd_02) + call readLIS_Teff_usaf(n, yyyymmdd_02, hh_02, Orbit, teff_02, rc) + if (rc .ne. 0) then + write(LDT_logunit,*)'[WARN] No Teff data available...' + endif - call readLIS_Teff(n,yyyymmdd_02,hh_02,Orbit,teff_02) + ! Now read the third time. + call ESMF_TimeGet(thirdtime, yy=yr_03, mm=mo_03, dd=da_03, & + h=hr_03) + + write(unit=yyyy_03, fmt='(i4.4)') yr_03 + write(unit=mm_03, fmt='(i2.2)') mo_03 + write(unit=dd_03, fmt='(i2.2)') da_03 + write(unit=hh_03, fmt='(i2.2)') hr_03 + yyyymmdd_03 = trim(yyyy_03)//trim(mm_03)//trim(dd_03) + call readLIS_Teff_usaf(n, yyyymmdd_03, hh_03, Orbit, teff_03, rc) + if (rc .ne. 0) then + write(LDT_logunit,*)'[WARN] No Teff data available...' + endif - ! Scale LIS teff to GEOS teff climatology + ! Scale LIS teff to GEOS teff climatology ! get DOY call get_doy(mo_pre,da_pre,doy_pre) - if(SMAPeOPL%Teffscale.eq.1) then ! get getattributes call getattributes(SMAPeOPL%dailystats_ref,& @@ -402,10 +471,9 @@ subroutine LDT_smap_e_oplRun(n) allocate(SMAPeOPL%grid_row(SMAPeOPL%ngrid)) call read_DailyTeffStats(doy_pre) - ! scale write (LDT_logunit,*) '[INFO] Scaling LIS effective soil temperature' - call scale_teff(n,Orbit,teff_01,teff_02) + call scale_teff(n, Orbit, teff_01, teff_02, teff_03) write (LDT_logunit,*) '[INFO] Finished scaling LIS effective soil temperature' deallocate(SMAPeOPL%mu_6am_ref) @@ -419,25 +487,27 @@ subroutine LDT_smap_e_oplRun(n) deallocate(SMAPeOPL%grid_col) deallocate(SMAPeOPL%grid_row) endif - read_L1Bdata = .false. ! Get snow information from LIS outputs - call readLIS_snow(n,yyyymmdd,hh,SnowDepth) + call readUSAFSI(n, yyyymmdd, hh, SnowDepth, rc) + if (rc .ne. 0) then + write(LDT_logunit,*)'[WARN] No USAFSI data available!' + endif ! Retrieve SMAP soil moisture ! get DOY call get_doy(mo,da,doy_curr) ! get UTC - call get_UTC(n,TIMEsec,UTChr) + call get_UTC(n,TIMEsec,UTChr) ! retrieve ierr = LDT_create_subdirs(len_trim(SMAPeOPL%SMoutdir), & trim(SMAPeOPL%SMoutdir)) - call ARFSSMRETRIEVAL(smap_L1B_filename(i),teff_01,teff_02,& - SnowDepth,doy_curr,UTChr) - + call ARFSSMRETRIEVAL(smap_L1B_filename(i), & + teff_01, teff_02, teff_03, & + SnowDepth, doy_curr, UTChr, firsttime, secondtime, thirdtime) deallocate(SMAPeOPL%ARFS_TBV_COR) endif enddo @@ -446,34 +516,37 @@ subroutine LDT_smap_e_oplRun(n) end subroutine LDT_smap_e_oplRun - subroutine search_SMAPL1B_files(ndir,date_curr,L1Btype) + subroutine search_SMAPL1B_files(ndir,date_curr,L1Btype,suffix) implicit none ! !ARGUMENTS: character (len=*) :: ndir character (len=*) :: date_curr - integer :: L1Btype + integer :: L1Btype,suffix ! !Local variables character*8 :: yyyymmdd character*2 :: hh + character*2 :: tmp character*200 :: list_files yyyymmdd = date_curr(1:8) hh = date_curr(9:10) + write (tmp,'(I2.2)') suffix if(L1Btype.eq.1) then !NRT list_files = 'ls '//trim(ndir)//'/SMAP_L1B_TB_NRT_*'//& trim(yyyymmdd)//'T'//trim(hh) & - //"*.h5 > SMAP_L1B_filelist.dat" + //'*.h5 > SMAP_L1B_filelist_'//trim(tmp)//'.dat' elseif(L1Btype.eq.2) then !Historical list_files = 'ls '//trim(ndir)//'/SMAP_L1B_TB_*'//& trim(yyyymmdd)//'T'//trim(hh) & - //"*.h5 > SMAP_L1B_filelist.dat" + //'*.h5 > SMAP_L1B_filelist_'//trim(tmp)//'.dat' endif call system(trim(list_files)) end subroutine search_SMAPL1B_files + end module LDT_smap_e_oplMod diff --git a/ldt/SMAP_E_OPL/SMAPL1BTOL1C_ARFS_SUBSET.F90 b/ldt/SMAP_E_OPL/SMAPL1BTOL1C_ARFS_SUBSET.F90 new file mode 100644 index 000000000..6cb1d5b15 --- /dev/null +++ b/ldt/SMAP_E_OPL/SMAPL1BTOL1C_ARFS_SUBSET.F90 @@ -0,0 +1,114 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! 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----------------------- +! +! SUBROUTINE: SMAPL1BRESAMPLE +! +! REVISION HISTORY: +! 22 Oct 2021: P.W.LIU; Initial implemetation +! 18 Dec 2021: Yonghwan Kwon; modified for LDT +! 09 Feb 2023: Eric Kemp; now processes subset of fields, no output of +! data to separate binary files. +! +! DESCRIPTION: RESAMPLE SMAPL1B TB TO AIR FORCE GRID +! INPUT : SMAP - L1B Brightness Temperature +! OUTPUT: SMAPTB_ARFSGRIDE_ddmmyyy.dat +! NOTES : Inverse Distance Squared with 0.4 deg serching window +!------------------------------------------------------------------------- + +subroutine SMAPL1BRESAMPLE_SUBSET(SMAPFILE,L1B_dir,Orbit,ARFS_TIME,rc) + + USE VARIABLES + USE DATADOMAIN + USE FUNCTIONS + USE TOOLSUBS + USE invdist_l1b2arfs + USE LDT_logMod + USE LDT_smap_e_oplMod + + IMPLICIT NONE + + INTEGER :: i, j, nrow, mcol + CHARACTER (len=100) :: SMAPFILE + character (len=100) :: L1B_dir + character (len=20) :: variable_name(13) + character (len=100) :: resample_filename(13) + character (len=1) :: Orbit + integer :: var_i + integer :: L1B_dir_len,L1B_fname_len + integer :: ierr + integer :: rc + + REAL*4,DIMENSION(:,:),ALLOCATABLE :: TIME_L1B, TBV_COR_L1B + REAL*4,DIMENSION(:,:),ALLOCATABLE :: LAT_L1B, LON_L1B, SCNANG_L1B + REAL*4,DIMENSION(:),ALLOCATABLE :: ANTSCN_L1B + INTEGER*4,DIMENSION(:,:),ALLOCATABLE :: TBVFLAG_L1B, TBHFLAG_L1B + REAL*8,DIMENSION(:), ALLOCATABLE :: ARFS_LAT, ARFS_LON + REAL*4,DIMENSION(2560,1920) :: ARFS_TIME, ARFS_COR_TBV + + REAL :: T1, T2 + + rc = 0 + + CALL ARFS_GEO + ALLOCATE(ARFS_LAT(arfs_nrow_lat),ARFS_LON(arfs_mcol_lon)) + ARFS_LAT = LAT(arfs_geo_lat_lo,arfs_geo_lat_up,-arfs_lat_space) + ARFS_LON = LON(arfs_geo_lon_lf,arfs_geo_lon_rt,arfs_lon_space) + + !Input (Path/filename, datatypes, lat, lon) Return(DATA,LAT,LON,length of row and col) + ! READ SMAP_L1B DATA FROM HDF5 + + ! EMK...Try fault tolerant NRT version. + CALL GetSMAP_L1B_NRT_SUBSET(SMAPFILE, TIME_L1B, TBV_COR_L1B, & + LAT_L1B, LON_L1B, TBVFLAG_L1B, TBHFLAG_L1B, & + ANTSCN_L1B, SCNANG_L1B, nrow, mcol, ierr) + if (ierr == 1) then + if (nrow == 0 .and. mcol == 0) then + write(LDT_logunit,*)'[ERR] Problem reading ', trim(SMAPFILE) + rc = 1 + return + else + write(LDT_logunit,*)'[ERR] Unknown internal error!' + write(LDT_logunit,*)'[ERR] Aborting...' + call LDT_endrun() + end if + end if + + !Input (DATA,LAT,LON,length of row and col); Return(TB in ARFS GRID) + !CALL L1BTB2ARFS_INVDIS(TIME_L1B, TBV_COR_L1B, TBH_COR_L1B, TBV_L1B, TBH_L1B, SURWAT_V_L1B, SURWAT_H_L1B, & + ! NETD_V_L1B, NETD_H_L1B, LAT_L1B, LON_L1B, TBVFLAG_L1B, TBHFLAG_L1B, ANTSCN_L1B, SCNANG_L1B, nrow, mcol, & + ! ARFS_LAT, ARFS_LON, ARFS_TIME, ARFS_COR_TBV, ARFS_COR_TBH, ARFS_TBV, ARFS_TBH, ARFS_NEDTV, ARFS_NEDTH, & + ! ARFS_SURWAT_V, ARFS_SURWAT_H, ARFS_WTV, ARFS_WTH, ARFS_SAMPLE_V, ARFS_SAMPLE_H) + CALL L1BTB2ARFS_INVDIS_SUBSET(TIME_L1B, TBV_COR_L1B, & + LAT_L1B, LON_L1B, TBVFLAG_L1B, TBHFLAG_L1B, ANTSCN_L1B, SCNANG_L1B, & + nrow, mcol, & + ARFS_LAT, ARFS_LON, ARFS_TIME, ARFS_COR_TBV) + SMAPeOPL%ARFS_TBV_COR = ARFS_COR_TBV + + L1B_dir_len = len_trim(L1B_dir) + L1B_fname_len = len_trim(SMAPFILE) + + if(SMAPeOPL%L1Btype.eq.1) then !NRT + Orbit = trim(SMAPFILE(L1B_dir_len+24:L1B_dir_len+24)) + elseif(SMAPeOPL%L1Btype.eq.2) then !Historical + Orbit = trim(SMAPFILE(L1B_dir_len+20:L1B_dir_len+20)) + endif + + ! Cleanup + deallocate(TBV_COR_L1B) + deallocate(LAT_L1B) + deallocate(LON_L1B) + deallocate(SCNANG_L1B) + deallocate(ANTSCN_L1B) + deallocate(TBVFLAG_L1B) + deallocate(TBHFLAG_L1B) + deallocate(ARFS_LAT) + deallocate(ARFS_LON) + +end subroutine SMAPL1BRESAMPLE_SUBSET diff --git a/ldt/SMAP_E_OPL/TOOLSUBS.F90 b/ldt/SMAP_E_OPL/TOOLSUBS.F90 index 33bf25ce4..2a31460f8 100644 --- a/ldt/SMAP_E_OPL/TOOLSUBS.F90 +++ b/ldt/SMAP_E_OPL/TOOLSUBS.F90 @@ -181,6 +181,8 @@ END SUBROUTINE GetSMAP_L2 SUBROUTINE GetSMAP_L1B(filename, data1_out, data2_out, data3_out, data4_out, data5_out, data6_out, data7_out, data8_out, & data9_out, data10_out, data11_out, data12_out, data13_out, data14_out, data15_out, n, m) + use LDT_logMod, only: LDT_logunit ! EMK + CHARACTER (len=100) :: filename, dataset1, dataset2, dataset3, dataset4, dataset5, dataset6, dataset7 CHARACTER (len=100) :: dataset8, dataset9, dataset10, dataset11, dataset12, dataset13, dataset14, dataset15 #if (defined USE_HDF5) @@ -216,8 +218,9 @@ SUBROUTINE GetSMAP_L1B(filename, data1_out, data2_out, data3_out, data4_out, dat #if (defined USE_HDF5) CALL h5open_f(hdferr) !Initialize hdf5 CALL h5fopen_f (trim(filename),H5F_ACC_RDONLY_F,file_id,hdferr) !Open file -! PRINT*, 'sd_id, hdferr', file_id, hdferr + ! PRINT*, 'sd_id, hdferr', file_id, hdferr CALL h5dopen_f (file_id, trim(dataset1),dataset_id1, hdferr) !Open dataset + call h5dget_space_f(dataset_id1,dspace_id,hdferr) call h5sget_simple_extent_dims_f(dspace_id, dims, maxdims, hdferr) CALL h5dopen_f (file_id, trim(dataset2),dataset_id2, hdferr) !Open dataset @@ -235,7 +238,6 @@ SUBROUTINE GetSMAP_L1B(filename, data1_out, data2_out, data3_out, data4_out, dat CALL h5dopen_f (file_id, trim(dataset14),dataset_id14, hdferr) !Open dataset CALL h5dopen_f (file_id, trim(dataset15),dataset_id15, hdferr) !Open dataset - !PRINT*,'ds_id hdferr', dataset_id1, hdferr !PRINT*,'dims maxdims', dims, maxdims n=dims(1) @@ -285,6 +287,369 @@ SUBROUTINE GetSMAP_L1B(filename, data1_out, data2_out, data3_out, data4_out, dat #endif END SUBROUTINE GetSMAP_L1B + ! Forked version of GetSMAP_L1B to SMAP files in operations, processing + ! a subset of the fields. Also with fault tolerance and some special + ! logic for older NRT files missing the tb_v_surface_corrected field. + ! Eric Kemp/SSAI. + SUBROUTINE GetSMAP_L1B_NRT_subset(filename, tb_time_seconds, & + tb_v_surface_corrected, & + tb_lat, tb_lon, tb_qual_flag_v, tb_qual_flag_h, sc_nadir_angle, & + antenna_scan_angle, n, m, ierr) + + ! Imports + use LDT_logMod, only: LDT_logunit, LDT_endrun ! EMK + + ! Arguments + character(*), intent(in) :: filename + real*4, allocatable, intent(out) :: tb_time_seconds(:,:) + real*4, allocatable, intent(out) :: tb_v_surface_corrected(:,:) + real*4, allocatable, intent(out) :: tb_lat(:,:), tb_lon(:,:) + integer*4, allocatable, intent(out) :: tb_qual_flag_v(:,:) + integer*4, allocatable, intent(out) :: tb_qual_flag_h(:,:) + real*4, allocatable, intent(out) :: sc_nadir_angle(:) + real*4, allocatable, intent(out) :: antenna_scan_angle(:,:) + integer, intent(out) :: m, n + integer, intent(out) :: ierr + +#if (defined USE_HDF5) + + ! Locals + character(100) :: dataset + integer(HID_T) :: file_id, dataset_id, dspace_id + integer(HSIZE_T) :: dims(2), maxdims(2) + integer :: hdferr + logical :: exists, ishdf5 + + ierr = 0 + m = 0 + n = 0 + + ! Make sure file exists + inquire(file=trim(filename), exist=exists) + if (.not. exists) then + write(LDT_logunit,*)'[ERR] Cannot find file ', trim(filename) + ierr = 1 + return + end if + + ! Initialize HDF5 + call h5open_f(hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot initialize HDF5 Fortran interface!' + call h5close_f(hdferr) + ierr = 1 + return + end if + + ! Make sure the file is HDF5 + call h5fis_hdf5_f(trim(filename), ishdf5, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Problem checking if ', trim(filename), & + ' is HDF5' + call h5close_f(hdferr) + ierr = 1 + return + end if + if (.not. ishdf5) then + write(LDT_logunit,*)'[ERR] File ', trim(filename), ' is not HDF5!' + call h5close_f(hdferr) + ierr = 1 + return + end if + + ! Open the file + call h5fopen_f(trim(filename), H5F_ACC_RDONLY_F, file_id, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot open ', trim(filename) + call h5close_f(hdferr) + ierr = 1 + return + end if + + ! Find Tb_lat, plus dimensions + dataset = "/Brightness_Temperature/tb_lat/" + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + call h5dget_space_f(dataset_id, dspace_id, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot find dimensions for ', & + trim(dataset) + call h5dclose_f(dataset_id, hdferr) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + ierr = 1 + return + end if + call h5sget_simple_extent_dims_f(dspace_id, dims, maxdims, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot get dimensions for ', & + trim(dataset) + call h5dclose_f(dataset_id, hdferr) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + ierr = 1 + return + end if + + ! We have the dimensions for the arrays, so let's allocate here. + n = dims(1) + m = dims(2) + allocate(tb_v_surface_corrected(n,m)); tb_v_surface_corrected = 0 + allocate(tb_time_seconds(n,m)); tb_time_seconds = 0 + allocate(tb_lat(n,m)) ; tb_lat = 0 + allocate(tb_lon(n,m)) ; tb_lon = 0 + allocate(tb_qual_flag_v(n,m)) ; tb_qual_flag_v = 0 + allocate(tb_qual_flag_h(n,m)) ; tb_qual_flag_h = 0 + allocate(sc_nadir_angle(n)) ; sc_nadir_angle = 0 + allocate(antenna_scan_angle(n,m)) ; antenna_scan_angle = 0 + + ! Get Tb_lat + call read_dataset_real(file_id, dataset, dataset_id, tb_lat, & + dims, ierr) + if (ierr == 1) return + + ! Get tb_lon + dataset = "/Brightness_Temperature/tb_lon/" + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + call read_dataset_real(file_id, dataset, dataset_id, tb_lon, & + dims, ierr) + if (ierr == 1) return + + ! Get tb_time_seconds + dataset = "/Brightness_Temperature/tb_time_seconds/" + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + call read_dataset_real(file_id, dataset, dataset_id, tb_time_seconds, & + dims, ierr) + if (ierr == 1) return + + ! Get tb_qual_flag_v + dataset = "/Brightness_Temperature/tb_qual_flag_v/" + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + call read_dataset_integer(file_id, dataset, dataset_id, tb_qual_flag_v, & + dims, ierr) + if (ierr == 1) return + + ! Get tb_qual_flag_h + dataset = "/Brightness_Temperature/tb_qual_flag_h/" + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + call read_dataset_integer(file_id, dataset, dataset_id, tb_qual_flag_h, & + dims, ierr) + if (ierr == 1) return + + ! Get sc_nadir_angle + dataset = "/Spacecraft_Data/sc_nadir_angle/" + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + call read_dataset_real_1d(file_id, dataset, dataset_id, sc_nadir_angle, & + dims, ierr) + if (ierr == 1) return + + ! Get antenna_scan_angle + dataset = "/Brightness_Temperature/antenna_scan_angle/" + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + call read_dataset_real(file_id, dataset, dataset_id, & + antenna_scan_angle, dims, ierr) + if (ierr == 1) return + + ! We will try to get tb_v_surface_corrected, but this is missing + ! in older NRT files. If absent, we will substitute tb_v. + ! NOTE: If a *read* error occurs, the file will assumed to be + ! corrupt all arrays will be nuked. + ! Get tb_v_surface_corrected + dataset = "/Brightness_Temperature/tb_v_surface_corrected/" + call get_dataset_id(file_id, dataset, dataset_id, ierr, & + handle_missing_nrt_field=.true.) + if (ierr == 1) then + ierr = 0 + dataset = "/Brightness_Temperature/tb_v/" + write(LDT_logunit,*)'[WARN] Will try substituting ', trim(dataset) + call get_dataset_id(file_id, dataset, dataset_id, ierr) + if (ierr == 1) return + end if + call read_dataset_real(file_id, dataset, dataset_id, & + tb_v_surface_corrected, & + dims, ierr) + if (ierr == 1) return + + ! Clean up + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + ierr = 0 + + return + + contains + + ! Internal subroutine + subroutine get_dataset_id(file_id, dataset, dataset_id, ierr, & + handle_missing_nrt_field) + implicit none + integer(HID_T), intent(in) :: file_id + character(*), intent(in) :: dataset + integer(HID_T), intent(out) :: dataset_id + integer, intent(out) :: ierr + logical, optional, intent(in) :: handle_missing_nrt_field + logical :: link_exists + integer :: hdferr + call h5lexists_f(file_id, trim(dataset), link_exists, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Problem finding ', trim(dataset) + ierr = 1 + if (handle_missing_nrt_field) return + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + endif + if (.not. link_exists) then + write(LDT_logunit,*)'[ERR] Nonexistent dataset ', trim(dataset) + ierr = 1 + if (handle_missing_nrt_field) return + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + endif + call h5dopen_f(file_id, trim(dataset), dataset_id, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot open dataset ', trim(dataset) + ierr = 1 + if (handle_missing_nrt_field) return + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + end if + end subroutine get_dataset_id + + ! Internal subroutine + subroutine read_dataset_real(file_id, dataset, dataset_id, buf, & + dims, ierr) + implicit none + integer(HID_T), intent(in) :: file_id + character(*), intent(in) :: dataset + integer(HID_T), intent(in) :: dataset_id + real*4, intent(inout) :: buf(:,:) + integer(HSIZE_T), intent(in) :: dims(:) + integer, intent(out) :: ierr + integer :: hdferr + call h5dread_f(dataset_id, H5T_IEEE_F32LE, buf, dims, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot read dataset ', trim(dataset) + call h5dclose_f(dataset_id, hdferr) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + endif + call h5dclose_f(dataset_id, hdferr) + if (hdferr == -1) then + write(LDT_Logunit,*)'[ERR] Problem closing dataset ', & + trim(dataset) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + end if + end subroutine read_dataset_real + + ! Internal subroutine + subroutine read_dataset_real_1d(file_id, dataset, dataset_id, buf, & + dims, ierr) + implicit none + integer(HID_T), intent(in) :: file_id + character(*), intent(in) :: dataset + integer(HID_T), intent(in) :: dataset_id + real*4, intent(inout) :: buf(:) + integer(HSIZE_T), intent(in) :: dims(:) + integer, intent(out) :: ierr + integer :: hdferr + call h5dread_f(dataset_id, H5T_IEEE_F32LE, buf, dims, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot read dataset ', trim(dataset) + call h5dclose_f(dataset_id, hdferr) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + endif + call h5dclose_f(dataset_id, hdferr) + if (hdferr == -1) then + write(LDT_Logunit,*)'[ERR] Problem closing dataset ', & + trim(dataset) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + end if + end subroutine read_dataset_real_1d + + ! Internal subroutine + subroutine read_dataset_integer(file_id, dataset, dataset_id, buf, & + dims, ierr) + implicit none + integer(HID_T), intent(in) :: file_id + character(*), intent(in) :: dataset + integer(HID_T), intent(in) :: dataset_id + integer*4, intent(inout) :: buf(:,:) + integer(HSIZE_T), intent(in) :: dims(:) + integer, intent(out) :: ierr + integer :: hdferr + call h5dread_f(dataset_id, H5T_NATIVE_INTEGER, buf, dims, hdferr) + if (hdferr == -1) then + write(LDT_logunit,*)'[ERR] Cannot read dataset ', trim(dataset) + call h5dclose_f(dataset_id, hdferr) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + endif + call h5dclose_f(dataset_id, hdferr) + if (hdferr == -1) then + write(LDT_Logunit,*)'[ERR] Problem closing dataset ', & + trim(dataset) + call h5fclose_f(file_id, hdferr) + call h5close_f(hdferr) + call freeall(ierr) + return + end if + end subroutine read_dataset_integer + + ! Internal subroutine. Warning -- deallocates memory in + ! parent subroutine and resets two variables. This is intended + ! for gracefully handling errors returned from HDF5. + subroutine freeall(ierr) + implicit none + integer, intent(out) :: ierr + if (allocated(tb_v_surface_corrected)) & + deallocate(tb_v_surface_corrected) + if (allocated(tb_lat)) deallocate(tb_lat) + if (allocated(tb_lon)) deallocate(tb_lon) + if (allocated(tb_time_seconds)) deallocate(tb_time_seconds) + if (allocated(tb_qual_flag_v)) deallocate(tb_qual_flag_v) + if (allocated(tb_qual_flag_h)) deallocate(tb_qual_flag_h) + if (allocated(sc_nadir_angle))deallocate(sc_nadir_angle) + if (allocated(antenna_scan_angle)) deallocate(antenna_scan_angle) + m = 0 + n = 0 + ierr = 1 + return + end subroutine freeall +#else + ! Dummy version if LDT was compiled w/o HDF5 support. + write(LDT_logunit,*) & + '[ERR] GetSMAP_L1B_NRT called without HDF5 support!' + write(LDT_logunit,*) & + '[ERR] Recompile LDT with HDF5 support and try again!' + call LDE_endrun() +#endif + end subroutine GetSMAP_L1B_NRT_Subset + SUBROUTINE NEAREST_1d ( nd, xd, yd, ni, xi, yi ) IMPLICIT NONE INTEGER :: i, j, k, nd, ni @@ -397,5 +762,7 @@ SUBROUTINE LSQ_Linear(n, x, y, a, b, r, rmse) !PRINT*, " y-intercept b = ", b !PRINT*, " Correlation r = ", r END SUBROUTINE LSQ_Linear - + + + END MODULE TOOLSUBS diff --git a/ldt/SMAP_E_OPL/invdist_l1b2arfs.F90 b/ldt/SMAP_E_OPL/invdist_l1b2arfs.F90 index 07f09bac5..717c82a17 100644 --- a/ldt/SMAP_E_OPL/invdist_l1b2arfs.F90 +++ b/ldt/SMAP_E_OPL/invdist_l1b2arfs.F90 @@ -8,197 +8,306 @@ ! 26 Oct 2021: P.W.LIU; Initial implemetation ! 28 Jan 2022: P.W.LIU; CHNAGE OUTPUT COORDINATE TO COMPLY LIS OUTPUT ! 22 Feb 2022: Yonghwan Kwon; modified for LDT +! 10 Feb 2023: Eric Kemp; Added code to resample subset of variables. ! ! DESCRIPTION: SUBROUTINE TO RESAMPLE L1B TB ONTO AIR FORCE GRID !------------------------------------------------------------------------- MODULE invdist_l1b2arfs - IMPLICIT NONE - - CONTAINS - SUBROUTINE L1BTB2ARFS_INVDIS(tim, tbvl1b_cor, tbhl1b_cor, tbvl1b, tbhl1b, surwat_v_l1b, surwat_h_l1b, & - netd_v_l1b, netd_h_l1b, lat_l1b, lon_l1b, tbv_qual_flag, tbh_qual_flag, sc_nadir_angle, antenna_scan_angle, nrows_l1btb, ncols_l1btb, & - ref_lat, ref_lon, arfs_tim, arfs_tbv_cor, arfs_tbh_cor, arfs_tbv, arfs_tbh, arfs_nedtv, arfs_nedth, & - arfs_surwatv, arfs_surwath, arfs_wt_cor_tbv, arfs_wt_cor_tbh, arfs_samplenumv, arfs_samplenumh) - - - INTEGER(4) :: ii, jj, k, r, c, rr, rmin, rmax, cc, cmin, cmax, nrows_l1btb, ncols_l1btb - INTEGER(4), PARAMETER :: qualitybit = 0 - REAL(8), PARAMETER :: RE_KM = 6371.228, search_radius = 20.0, PI = 3.141592653589793238, d2r = PI/180.0 - REAL(8) :: gcdist, lat1, lon1, lat2, lon2 - REAL*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: tim, tbvl1b_cor, tbhl1b_cor, tbvl1b, tbhl1b, surwat_v_l1b, surwat_h_l1b, netd_v_l1b, netd_h_l1b - REAL*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: lat_l1b, lon_l1b, antenna_scan_angle - REAL*4,DIMENSION(ncols_l1btb) :: sc_nadir_angle - INTEGER*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: tbv_qual_flag, tbh_qual_flag - INTEGER(4),DIMENSION(:,:),ALLOCATABLE :: zerodistflag - REAL*8,DIMENSION(:),ALLOCATABLE :: ref_lat, ref_lon - REAL*4,DIMENSION(2560,1920) :: arfs_tim, arfs_tbv_cor, arfs_tbh_cor, arfs_tbv, arfs_tbh, arfs_nedtv, arfs_nedth, arfs_surwatv, arfs_surwath - REAL*4,DIMENSION(2560,1920) :: arfs_wt_tim, arfs_wt_cor_tbv, arfs_wt_cor_tbh, arfs_wt_tbv, arfs_wt_tbh, arfs_wt_nedtv, arfs_wt_nedth, arfs_wt_surwatv, arfs_wt_surwath - INTEGER*4,DIMENSION(2560,1920) :: arfs_samplenumv, arfs_samplenumh - - !ALLOCATE(zerodistflag(size(ref_lat),size(ref_lon))) - ALLOCATE(zerodistflag(size(ref_lon),size(ref_lat))) - !INITIAL THE OUTPUT VARIABLES - arfs_tim=0.0 - arfs_tbv_cor=0.0 - arfs_tbh_cor=0.0 - arfs_tbv=0.0 - arfs_tbh=0.0 - arfs_nedtv=0.0 - arfs_nedth=0.0 - arfs_surwatv=0.0 - arfs_surwath=0.0 - arfs_wt_tim=0.0 - arfs_wt_cor_tbv=0.0 - arfs_wt_cor_tbh=0.0 - arfs_wt_tbv=0.0 - arfs_wt_tbh=0.0 - arfs_wt_nedtv=0.0 - arfs_wt_nedth=0.0 - arfs_wt_surwatv=0.0 - arfs_wt_surwath=0.0 - arfs_samplenumv=0.0 - - DO ii = 1,ncols_l1btb - IF (ABS (sc_nadir_angle(ii)) <= 2.0) THEN - DO jj = 1,nrows_l1btb - IF (ABS (antenna_scan_angle(jj,ii)).LE.360.00) THEN - ! FIND ARFS_GRID (r,c) - c = MINLOC(ABS(lat_l1b(jj,ii)-ref_lat(:)),1) !Lat Direction - r = MINLOC(ABS(lon_l1b(jj,ii)-ref_lon(:)),1) !Lon Direction - rmin=r-5 ; IF (rmin < 1) rmin=1 - rmax=r+5 ; IF (rmax > size(ref_lon)) rmax=size(ref_lon) - cmin=c-5 ; IF (cmin < 1) cmin=1 - cmax=c+5 ; IF (cmax > size(ref_lat)) cmax=size(ref_lat) - IF (IBITS (tbv_qual_flag(jj,ii),qualitybit,1) == 0 .AND. IBITS (tbh_qual_flag(jj,jj),qualitybit,1) == 0) THEN !RESAMPLE ONLY WHEN BOTH V and H MEET QUALITY - k=0 - DO rr = rmin,rmax !Lon direction - DO cc =cmin,cmax !Lat direction - lat1 = DBLE (lat_l1b(jj,ii)*d2r) - lon1 = DBLE (lon_l1b(jj,ii)*d2r) - lat2 = DBLE (ref_lat(cc)*d2r) - lon2 = DBLE (ref_lon(rr)*d2r) - - if(lat1.eq.lat2.and.lon1.eq.lon2) then - gcdist = 0. - else - gcdist = RE_KM * DACOS ( DSIN (lat1) * DSIN (lat2) + DCOS (lat1) * DCOS (lat2) * DCOS (lon1-lon2) ) - endif - - IF (gcdist < search_radius) THEN !RESAMPLE ONLY WITHIN THE SEARCH RANGE - IF (gcdist < 0.0001D0) THEN !The TB is right on the grid center - zerodistflag (rr,cc) = 1 - IF ((ABS (tim(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tim(rr,cc) = tim(jj,ii) ; arfs_wt_tim(rr,cc) = 1.0 - END IF - IF ((ABS (surwat_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_surwatv(rr,cc) = surwat_v_l1b(jj,ii) ; arfs_wt_surwatv(rr,cc) = 1.0 - END IF - IF ((ABS (surwat_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_surwath(rr,cc) = surwat_h_l1b(jj,ii) ; arfs_wt_surwath(rr,cc) = 1.0 - END IF - IF ((ABS (netd_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_nedtv(rr,cc) = netd_v_l1b(jj,ii) ; arfs_wt_nedtv(rr,cc) = 1.0 - END IF - IF ((ABS (netd_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_nedth(rr,cc) = netd_h_l1b(jj,ii) ; arfs_wt_nedth(rr,cc) = 1.0 - END IF - IF ((ABS (tbvl1b_cor(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbv_cor(rr,cc) = tbvl1b_cor(jj,ii) ; arfs_wt_cor_tbv(rr,cc) = 1.0 - arfs_samplenumv(rr,cc)=1 !Sample number only calculate for correct tb - k=k+1; - END IF - IF ((ABS (tbhl1b_cor(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbh_cor(rr,cc) = tbhl1b_cor(jj,ii) ; arfs_wt_cor_tbh(rr,cc) = 1.0 - arfs_samplenumh(rr,cc)=1 !Sample number only calculate for correct tb - END IF - IF ((ABS (tbvl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbv(rr,cc) = tbvl1b(jj,ii) ; arfs_wt_tbv(rr,cc) = 1.0 - END IF - IF ((ABS (tbhl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbh(rr,cc) = tbhl1b(jj,ii) ; arfs_wt_tbh(rr,cc) = 1.0 - END IF - ELSE - IF (zerodistflag (rr,cc).EQ.0) THEN - - IF ((ABS (tim(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tim(rr,cc) = arfs_tim(rr,cc) + tim(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_tim(rr,cc) = arfs_wt_tim(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - END IF - IF ((ABS (surwat_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_surwatv(rr,cc) = arfs_surwatv(rr,cc) + surwat_v_l1b(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_surwatv(rr,cc) = arfs_wt_surwatv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - END IF - IF ((ABS (surwat_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_surwath(rr,cc) = arfs_surwath(rr,cc) + surwat_h_l1b(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_surwath(rr,cc) = arfs_wt_surwath(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - END IF - IF ((ABS (netd_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_nedtv(rr,cc) = arfs_nedtv(rr,cc) + netd_v_l1b(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_nedtv(rr,cc) = arfs_wt_nedtv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - END IF - IF ((ABS (netd_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_nedth(rr,cc) = arfs_nedth(rr,cc) + netd_h_l1b(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_nedth(rr,cc) = arfs_wt_nedth(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - END IF - IF ((ABS (tbvl1b_cor(jj,ii) - (-9999.0))).GT.1.0D-7) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbv_cor(rr,cc) = arfs_tbv_cor(rr,cc) + tbvl1b_cor(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_cor_tbv(rr,cc) = arfs_wt_cor_tbv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - arfs_samplenumv(rr,cc)=arfs_samplenumv(rr,cc)+1.0 !Sample number only calculate for correct tb - k=k+1; - END IF - IF ((ABS (tbhl1b_cor(jj,ii) - (-9999.0))).GT.1.0D-7) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbh_cor(rr,cc) = arfs_tbh_cor(rr,cc) + tbhl1b_cor(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_cor_tbh(rr,cc) = arfs_wt_cor_tbh(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - arfs_samplenumh(rr,cc)=arfs_samplenumh(rr,cc)+1.0 !Sample number only calculate for correct tb - END IF - IF ((ABS (tbvl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbv(rr,cc) = arfs_tbv(rr,cc) + tbvl1b(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_tbv(rr,cc) = arfs_wt_tbv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - END IF - IF ((ABS (tbhl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) - arfs_tbh(rr,cc) = arfs_tbh(rr,cc) + tbhl1b(jj,ii) / SNGL (gcdist*gcdist) - arfs_wt_tbh(rr,cc) = arfs_wt_tbh(rr,cc) + 1.0 / SNGL (gcdist*gcdist) - END IF - END IF !(zerodistflag (rr,cc) = 0) - END IF !(gcdist < 0.0001D0)! - END IF !(gcdist < search_radius) - END DO !cc =cmin,cmax - END DO !rr = rmin,rmax - END IF !(IBITS (tbv_qual_flag(jj,ii),qualitybit,1) == 0 .AND. IBITS (tbh_qual_flag(ii,jj),qualitybit,1) == 0) - END IF !(ABS (antenna_scan_angle(ii,jj)) <= 360.00) - END DO !jj=1,2 - END IF !(ABS (sc_nadir_angle(ii)) <= 2.0) - END DO !ii=1,2 - - !APPLY WEIGHTING FUNCTION FOR RESAMPLING - WHERE(arfs_tim.NE.0.0.AND.arfs_wt_tim.NE.0.0) - arfs_tim = arfs_tim / arfs_wt_tim - END WHERE - WHERE(arfs_surwatv.NE.0.0.AND.arfs_wt_surwatv.NE.0.0) - arfs_surwatv = arfs_surwatv / arfs_wt_surwatv - END WHERE - WHERE(arfs_surwath.NE.0.0.AND.arfs_wt_surwath.NE.0.0) - arfs_surwath = arfs_surwath / arfs_wt_surwath - END WHERE - WHERE(arfs_nedtv.NE.0.0.AND.arfs_wt_nedtv.NE.0.0) - arfs_nedtv = arfs_nedtv / arfs_wt_nedtv - END WHERE - WHERE(arfs_nedth.NE.0.0.AND.arfs_wt_nedth.NE.0.0) - arfs_nedth = arfs_nedth / arfs_wt_nedth - END WHERE - WHERE(arfs_tbv_cor.NE.0.0.AND.arfs_wt_cor_tbv.NE.0.0) - arfs_tbv_cor = arfs_tbv_cor / arfs_wt_cor_tbv - END WHERE - WHERE(arfs_tbh_cor.NE.0.0.AND.arfs_wt_cor_tbh.NE.0.0) - arfs_tbh_cor = arfs_tbh_cor / arfs_wt_cor_tbh - END WHERE - WHERE(arfs_tbv.NE.0.0.AND.arfs_wt_tbv.NE.0.0) - arfs_tbv = arfs_tbv / arfs_wt_tbv - END WHERE - WHERE(arfs_tbh.NE.0.0.AND.arfs_wt_tbh.NE.0.0) - arfs_tbh = arfs_tbh / arfs_wt_tbh - END WHERE - - END SUBROUTINE L1BTB2ARFS_INVDIS + IMPLICIT NONE + + CONTAINS + SUBROUTINE L1BTB2ARFS_INVDIS(tim, tbvl1b_cor, tbhl1b_cor, tbvl1b, tbhl1b, surwat_v_l1b, surwat_h_l1b, & + netd_v_l1b, netd_h_l1b, lat_l1b, lon_l1b, tbv_qual_flag, tbh_qual_flag, sc_nadir_angle, antenna_scan_angle, nrows_l1btb, ncols_l1btb, & + ref_lat, ref_lon, arfs_tim, arfs_tbv_cor, arfs_tbh_cor, arfs_tbv, arfs_tbh, arfs_nedtv, arfs_nedth, & + arfs_surwatv, arfs_surwath, arfs_wt_cor_tbv, arfs_wt_cor_tbh, arfs_samplenumv, arfs_samplenumh) + + + INTEGER(4) :: ii, jj, k, r, c, rr, rmin, rmax, cc, cmin, cmax, nrows_l1btb, ncols_l1btb + INTEGER(4), PARAMETER :: qualitybit = 0 + REAL(8), PARAMETER :: RE_KM = 6371.228, search_radius = 20.0, PI = 3.141592653589793238, d2r = PI/180.0 + REAL(8) :: gcdist, lat1, lon1, lat2, lon2 + REAL*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: tim, tbvl1b_cor, tbhl1b_cor, tbvl1b, tbhl1b, surwat_v_l1b, surwat_h_l1b, netd_v_l1b, netd_h_l1b + REAL*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: lat_l1b, lon_l1b, antenna_scan_angle + REAL*4,DIMENSION(ncols_l1btb) :: sc_nadir_angle + INTEGER*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: tbv_qual_flag, tbh_qual_flag + INTEGER(4),DIMENSION(:,:),ALLOCATABLE :: zerodistflag + REAL*8,DIMENSION(:),ALLOCATABLE :: ref_lat, ref_lon + REAL*4,DIMENSION(2560,1920) :: arfs_tim, arfs_tbv_cor, arfs_tbh_cor, arfs_tbv, arfs_tbh, arfs_nedtv, arfs_nedth, arfs_surwatv, arfs_surwath + REAL*4,DIMENSION(2560,1920) :: arfs_wt_tim, arfs_wt_cor_tbv, arfs_wt_cor_tbh, arfs_wt_tbv, arfs_wt_tbh, arfs_wt_nedtv, arfs_wt_nedth, arfs_wt_surwatv, arfs_wt_surwath + INTEGER*4,DIMENSION(2560,1920) :: arfs_samplenumv, arfs_samplenumh + + !ALLOCATE(zerodistflag(size(ref_lat),size(ref_lon))) + ALLOCATE(zerodistflag(size(ref_lon),size(ref_lat))) + !INITIAL THE OUTPUT VARIABLES + arfs_tim=0.0 + arfs_tbv_cor=0.0 + arfs_tbh_cor=0.0 + arfs_tbv=0.0 + arfs_tbh=0.0 + arfs_nedtv=0.0 + arfs_nedth=0.0 + arfs_surwatv=0.0 + arfs_surwath=0.0 + arfs_wt_tim=0.0 + arfs_wt_cor_tbv=0.0 + arfs_wt_cor_tbh=0.0 + arfs_wt_tbv=0.0 + arfs_wt_tbh=0.0 + arfs_wt_nedtv=0.0 + arfs_wt_nedth=0.0 + arfs_wt_surwatv=0.0 + arfs_wt_surwath=0.0 + arfs_samplenumv=0.0 + + DO ii = 1,ncols_l1btb + IF (ABS (sc_nadir_angle(ii)) <= 2.0) THEN + DO jj = 1,nrows_l1btb + IF (ABS (antenna_scan_angle(jj,ii)).LE.360.00) THEN + ! FIND ARFS_GRID (r,c) + c = MINLOC(ABS(lat_l1b(jj,ii)-ref_lat(:)),1) !Lat Direction + r = MINLOC(ABS(lon_l1b(jj,ii)-ref_lon(:)),1) !Lon Direction + rmin=r-5 ; IF (rmin < 1) rmin=1 + rmax=r+5 ; IF (rmax > size(ref_lon)) rmax=size(ref_lon) + cmin=c-5 ; IF (cmin < 1) cmin=1 + cmax=c+5 ; IF (cmax > size(ref_lat)) cmax=size(ref_lat) + IF (IBITS (tbv_qual_flag(jj,ii),qualitybit,1) == 0 .AND. IBITS (tbh_qual_flag(jj,jj),qualitybit,1) == 0) THEN !RESAMPLE ONLY WHEN BOTH V and H MEET QUALITY + k=0 + DO rr = rmin,rmax !Lon direction + DO cc =cmin,cmax !Lat direction + lat1 = DBLE (lat_l1b(jj,ii)*d2r) + lon1 = DBLE (lon_l1b(jj,ii)*d2r) + lat2 = DBLE (ref_lat(cc)*d2r) + lon2 = DBLE (ref_lon(rr)*d2r) + + if(lat1.eq.lat2.and.lon1.eq.lon2) then + gcdist = 0. + else + gcdist = RE_KM * DACOS ( DSIN (lat1) * DSIN (lat2) + DCOS (lat1) * DCOS (lat2) * DCOS (lon1-lon2) ) + endif + + IF (gcdist < search_radius) THEN !RESAMPLE ONLY WITHIN THE SEARCH RANGE + IF (gcdist < 0.0001D0) THEN !The TB is right on the grid center + zerodistflag (rr,cc) = 1 + IF ((ABS (tim(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tim(rr,cc) = tim(jj,ii) ; arfs_wt_tim(rr,cc) = 1.0 + END IF + IF ((ABS (surwat_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_surwatv(rr,cc) = surwat_v_l1b(jj,ii) ; arfs_wt_surwatv(rr,cc) = 1.0 + END IF + IF ((ABS (surwat_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_surwath(rr,cc) = surwat_h_l1b(jj,ii) ; arfs_wt_surwath(rr,cc) = 1.0 + END IF + IF ((ABS (netd_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_nedtv(rr,cc) = netd_v_l1b(jj,ii) ; arfs_wt_nedtv(rr,cc) = 1.0 + END IF + IF ((ABS (netd_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_nedth(rr,cc) = netd_h_l1b(jj,ii) ; arfs_wt_nedth(rr,cc) = 1.0 + END IF + IF ((ABS (tbvl1b_cor(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbv_cor(rr,cc) = tbvl1b_cor(jj,ii) ; arfs_wt_cor_tbv(rr,cc) = 1.0 + arfs_samplenumv(rr,cc)=1 !Sample number only calculate for correct tb + k=k+1; + END IF + IF ((ABS (tbhl1b_cor(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbh_cor(rr,cc) = tbhl1b_cor(jj,ii) ; arfs_wt_cor_tbh(rr,cc) = 1.0 + arfs_samplenumh(rr,cc)=1 !Sample number only calculate for correct tb + END IF + IF ((ABS (tbvl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbv(rr,cc) = tbvl1b(jj,ii) ; arfs_wt_tbv(rr,cc) = 1.0 + END IF + IF ((ABS (tbhl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbh(rr,cc) = tbhl1b(jj,ii) ; arfs_wt_tbh(rr,cc) = 1.0 + END IF + ELSE + IF (zerodistflag (rr,cc).EQ.0) THEN + + IF ((ABS (tim(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tim(rr,cc) = arfs_tim(rr,cc) + tim(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_tim(rr,cc) = arfs_wt_tim(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + IF ((ABS (surwat_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_surwatv(rr,cc) = arfs_surwatv(rr,cc) + surwat_v_l1b(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_surwatv(rr,cc) = arfs_wt_surwatv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + IF ((ABS (surwat_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_surwath(rr,cc) = arfs_surwath(rr,cc) + surwat_h_l1b(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_surwath(rr,cc) = arfs_wt_surwath(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + IF ((ABS (netd_v_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_nedtv(rr,cc) = arfs_nedtv(rr,cc) + netd_v_l1b(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_nedtv(rr,cc) = arfs_wt_nedtv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + IF ((ABS (netd_h_l1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_nedth(rr,cc) = arfs_nedth(rr,cc) + netd_h_l1b(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_nedth(rr,cc) = arfs_wt_nedth(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + IF ((ABS (tbvl1b_cor(jj,ii) - (-9999.0))).GT.1.0D-7) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbv_cor(rr,cc) = arfs_tbv_cor(rr,cc) + tbvl1b_cor(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_cor_tbv(rr,cc) = arfs_wt_cor_tbv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + arfs_samplenumv(rr,cc)=arfs_samplenumv(rr,cc)+1.0 !Sample number only calculate for correct tb + k=k+1; + END IF + IF ((ABS (tbhl1b_cor(jj,ii) - (-9999.0))).GT.1.0D-7) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbh_cor(rr,cc) = arfs_tbh_cor(rr,cc) + tbhl1b_cor(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_cor_tbh(rr,cc) = arfs_wt_cor_tbh(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + arfs_samplenumh(rr,cc)=arfs_samplenumh(rr,cc)+1.0 !Sample number only calculate for correct tb + END IF + IF ((ABS (tbvl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbv(rr,cc) = arfs_tbv(rr,cc) + tbvl1b(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_tbv(rr,cc) = arfs_wt_tbv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + IF ((ABS (tbhl1b(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbh(rr,cc) = arfs_tbh(rr,cc) + tbhl1b(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_tbh(rr,cc) = arfs_wt_tbh(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + END IF !(zerodistflag (rr,cc) = 0) + END IF !(gcdist < 0.0001D0)! + END IF !(gcdist < search_radius) + END DO !cc =cmin,cmax + END DO !rr = rmin,rmax + END IF !(IBITS (tbv_qual_flag(jj,ii),qualitybit,1) == 0 .AND. IBITS (tbh_qual_flag(ii,jj),qualitybit,1) == 0) + END IF !(ABS (antenna_scan_angle(ii,jj)) <= 360.00) + END DO !jj=1,2 + END IF !(ABS (sc_nadir_angle(ii)) <= 2.0) + END DO !ii=1,2 + + !APPLY WEIGHTING FUNCTION FOR RESAMPLING + WHERE(arfs_tim.NE.0.0.AND.arfs_wt_tim.NE.0.0) + arfs_tim = arfs_tim / arfs_wt_tim + END WHERE + WHERE(arfs_surwatv.NE.0.0.AND.arfs_wt_surwatv.NE.0.0) + arfs_surwatv = arfs_surwatv / arfs_wt_surwatv + END WHERE + WHERE(arfs_surwath.NE.0.0.AND.arfs_wt_surwath.NE.0.0) + arfs_surwath = arfs_surwath / arfs_wt_surwath + END WHERE + WHERE(arfs_nedtv.NE.0.0.AND.arfs_wt_nedtv.NE.0.0) + arfs_nedtv = arfs_nedtv / arfs_wt_nedtv + END WHERE + WHERE(arfs_nedth.NE.0.0.AND.arfs_wt_nedth.NE.0.0) + arfs_nedth = arfs_nedth / arfs_wt_nedth + END WHERE + WHERE(arfs_tbv_cor.NE.0.0.AND.arfs_wt_cor_tbv.NE.0.0) + arfs_tbv_cor = arfs_tbv_cor / arfs_wt_cor_tbv + END WHERE + WHERE(arfs_tbh_cor.NE.0.0.AND.arfs_wt_cor_tbh.NE.0.0) + arfs_tbh_cor = arfs_tbh_cor / arfs_wt_cor_tbh + END WHERE + WHERE(arfs_tbv.NE.0.0.AND.arfs_wt_tbv.NE.0.0) + arfs_tbv = arfs_tbv / arfs_wt_tbv + END WHERE + WHERE(arfs_tbh.NE.0.0.AND.arfs_wt_tbh.NE.0.0) + arfs_tbh = arfs_tbh / arfs_wt_tbh + END WHERE + + END SUBROUTINE L1BTB2ARFS_INVDIS + + ! EMK...Only process subset of SMAP L1B fields for NRT operations. + SUBROUTINE L1BTB2ARFS_INVDIS_SUBSET(tim, tbvl1b_cor, & + lat_l1b, lon_l1b, tbv_qual_flag, tbh_qual_flag, & + sc_nadir_angle, antenna_scan_angle, nrows_l1btb, ncols_l1btb, & + ref_lat, ref_lon, arfs_tim, arfs_tbv_cor) + + INTEGER(4) :: ii, jj, k, r, c, rr, rmin, rmax, cc, cmin, cmax, nrows_l1btb, ncols_l1btb + INTEGER(4), PARAMETER :: qualitybit = 0 + REAL(8), PARAMETER :: RE_KM = 6371.228, search_radius = 20.0, PI = 3.141592653589793238, d2r = PI/180.0 + REAL(8) :: gcdist, lat1, lon1, lat2, lon2 + REAL*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: tim, tbvl1b_cor + REAL*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: lat_l1b, lon_l1b, antenna_scan_angle + REAL*4,DIMENSION(ncols_l1btb) :: sc_nadir_angle + INTEGER*4,DIMENSION(nrows_l1btb,ncols_l1btb) :: tbv_qual_flag, tbh_qual_flag + INTEGER(4),DIMENSION(:,:),ALLOCATABLE :: zerodistflag + REAL*8,DIMENSION(:),ALLOCATABLE :: ref_lat, ref_lon + REAL*4,DIMENSION(2560,1920) :: arfs_tim, arfs_tbv_cor + REAL*4,DIMENSION(2560,1920) :: arfs_wt_tim, arfs_wt_cor_tbv + INTEGER*4,DIMENSION(2560,1920) :: arfs_samplenumv + + ALLOCATE(zerodistflag(size(ref_lon),size(ref_lat))) + !INITIAL THE OUTPUT VARIABLES + arfs_tim=0.0 + arfs_tbv_cor=0.0 + arfs_wt_tim=0.0 + arfs_wt_cor_tbv=0.0 + arfs_samplenumv=0.0 + + DO ii = 1,ncols_l1btb + IF (ABS (sc_nadir_angle(ii)) <= 2.0) THEN + DO jj = 1,nrows_l1btb + IF (ABS (antenna_scan_angle(jj,ii)).LE.360.00) THEN + lat1 = DBLE (lat_l1b(jj,ii)*d2r) + lon1 = DBLE (lon_l1b(jj,ii)*d2r) + ! FIND ARFS_GRID (r,c) + c = MINLOC(ABS(lat_l1b(jj,ii)-ref_lat(:)),1) !Lat Direction + r = MINLOC(ABS(lon_l1b(jj,ii)-ref_lon(:)),1) !Lon Direction + rmin=r-5 ; IF (rmin < 1) rmin=1 + rmax=r+5 ; IF (rmax > size(ref_lon)) rmax=size(ref_lon) + cmin=c-5 ; IF (cmin < 1) cmin=1 + cmax=c+5 ; IF (cmax > size(ref_lat)) cmax=size(ref_lat) +! IF (IBITS (tbv_qual_flag(jj,ii),qualitybit,1) == 0 .AND. IBITS (tbh_qual_flag(jj,jj),qualitybit,1) == 0) THEN !RESAMPLE ONLY WHEN BOTH V and H MEET QUALITY + IF (IBITS (tbv_qual_flag(jj,ii),qualitybit,1) == 0 .AND. IBITS (tbh_qual_flag(jj,ii),qualitybit,1) == 0) THEN !RESAMPLE ONLY WHEN BOTH V and H MEET QUALITY + k=0 + ! DO rr = rmin,rmax !Lon direction + ! DO cc =cmin,cmax !Lat direction + DO cc =cmin,cmax !Lat direction + lat2 = DBLE (ref_lat(cc)*d2r) + DO rr = rmin,rmax !Lon direction + lon2 = DBLE (ref_lon(rr)*d2r) + + if(lat1.eq.lat2.and.lon1.eq.lon2) then + gcdist = 0. + else + gcdist = RE_KM * DACOS ( DSIN (lat1) * DSIN (lat2) + DCOS (lat1) * DCOS (lat2) * DCOS (lon1-lon2) ) + endif + + IF (gcdist < search_radius) THEN !RESAMPLE ONLY WITHIN THE SEARCH RANGE + IF (gcdist < 0.0001D0) THEN !The TB is right on the grid center + zerodistflag (rr,cc) = 1 + IF ((ABS (tim(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tim(rr,cc) = tim(jj,ii) ; arfs_wt_tim(rr,cc) = 1.0 + END IF + IF ((ABS (tbvl1b_cor(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbv_cor(rr,cc) = tbvl1b_cor(jj,ii) ; arfs_wt_cor_tbv(rr,cc) = 1.0 + arfs_samplenumv(rr,cc)=1 !Sample number only calculate for correct tb + k=k+1; + END IF + ELSE + IF (zerodistflag (rr,cc).EQ.0) THEN + + IF ((ABS (tim(jj,ii) - (-9999.0)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tim(rr,cc) = arfs_tim(rr,cc) + tim(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_tim(rr,cc) = arfs_wt_tim(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + END IF + IF ((ABS (tbvl1b_cor(jj,ii) - (-9999.0))).GT.1.0D-7) THEN !DO IF NOT FILLVALUE(-9999) + arfs_tbv_cor(rr,cc) = arfs_tbv_cor(rr,cc) + tbvl1b_cor(jj,ii) / SNGL (gcdist*gcdist) + arfs_wt_cor_tbv(rr,cc) = arfs_wt_cor_tbv(rr,cc) + 1.0 / SNGL (gcdist*gcdist) + arfs_samplenumv(rr,cc)=arfs_samplenumv(rr,cc)+1.0 !Sample number only calculate for correct tb + k=k+1; + END IF + + END IF !(zerodistflag (rr,cc) = 0) + END IF !(gcdist < 0.0001D0)! + END IF !(gcdist < search_radius) + ! END DO !cc =cmin,cmax + ! END DO !rr = rmin,rmax + END DO !rr = rmin,rmax + END DO !cc =cmin,cmax + + END IF !(IBITS (tbv_qual_flag(jj,ii),qualitybit,1) == 0 .AND. IBITS (tbh_qual_flag(ii,jj),qualitybit,1) == 0) + END IF !(ABS (antenna_scan_angle(ii,jj)) <= 360.00) + END DO !jj=1,2 + END IF !(ABS (sc_nadir_angle(ii)) <= 2.0) + END DO !ii=1,2 + + !APPLY WEIGHTING FUNCTION FOR RESAMPLING + WHERE(arfs_tim.NE.0.0.AND.arfs_wt_tim.NE.0.0) + arfs_tim = arfs_tim / arfs_wt_tim + END WHERE + WHERE(arfs_tbv_cor.NE.0.0.AND.arfs_wt_cor_tbv.NE.0.0) + arfs_tbv_cor = arfs_tbv_cor / arfs_wt_cor_tbv + END WHERE + + deallocate(zerodistflag) ! EMK cleanup memory + END SUBROUTINE L1BTB2ARFS_INVDIS_SUBSET + END MODULE invdist_l1b2arfs diff --git a/ldt/SMAP_E_OPL/invdist_temp2smap.F90 b/ldt/SMAP_E_OPL/invdist_temp2smap.F90 index 0f7e9be14..b7283fd46 100644 --- a/ldt/SMAP_E_OPL/invdist_temp2smap.F90 +++ b/ldt/SMAP_E_OPL/invdist_temp2smap.F90 @@ -4,112 +4,114 @@ ! SMAP RESOLUTIN AT ARFS GRID !======================================================================= MODULE invdist_temp2smap - IMPLICIT NONE + IMPLICIT NONE - CONTAINS - !SUBROUTINE RESAMPLETEMP(arfs_temp,arfs_lat,arfs_lon,arfs_fine_lat,arfs_fine_lon,arfs_fine_wt,arfs_fine_temp,arfs_wt,arfs_smap_temp) - SUBROUTINE RESAMPLETEMP(arfs_temp,arfs_lat,arfs_lon,arfs_fine_lat,arfs_fine_lon,arfs_smap_temp) + CONTAINS + !SUBROUTINE RESAMPLETEMP(arfs_temp,arfs_lat,arfs_lon,arfs_fine_lat,arfs_fine_lon,arfs_fine_wt,arfs_fine_temp,arfs_wt,arfs_smap_temp) + SUBROUTINE RESAMPLETEMP(arfs_temp,arfs_lat,arfs_lon,arfs_fine_lat,arfs_fine_lon,arfs_smap_temp) - INTEGER(4) :: ii, jj, k, r, c, rr, rmin, rmax, cc, cmin, cmax, nrows_l1btb, ncols_l1btb - REAL(8), PARAMETER :: RE_KM = 6371.228, search_radius = 20.0, PI = 3.141592653589793238, d2r = PI/180.0 - REAL(8) :: gcdist, lat1, lon1, lat2, lon2 - REAL*8 ,DIMENSION(:), ALLOCATABLE :: arfs_lat, arfs_lon - REAL*8 ,DIMENSION(:), ALLOCATABLE :: arfs_fine_lat, arfs_fine_lon - REAL*4,DIMENSION(2560,1920) :: arfs_temp, arfs_smap_temp, arfs_wt - REAL*4,DIMENSION(7680,5760) :: arfs_fine_temp, arfs_fine_wt - INTEGER(4),DIMENSION(2560,1920):: zerodistflag - INTEGER(4),DIMENSION(7680,5760):: zerodistflag_fine + INTEGER(4) :: ii, jj, k, r, c, rr, rmin, rmax, cc, cmin, cmax, nrows_l1btb, ncols_l1btb + REAL(8), PARAMETER :: RE_KM = 6371.228, search_radius = 20.0, PI = 3.141592653589793238, d2r = PI/180.0 + REAL(8) :: gcdist, lat1, lon1, lat2, lon2 + REAL*8 ,DIMENSION(:), ALLOCATABLE :: arfs_lat, arfs_lon + REAL*8 ,DIMENSION(:), ALLOCATABLE :: arfs_fine_lat, arfs_fine_lon + REAL*4,DIMENSION(2560,1920) :: arfs_temp, arfs_smap_temp, arfs_wt + REAL*4,DIMENSION(7680,5760) :: arfs_fine_temp, arfs_fine_wt + INTEGER(4),DIMENSION(2560,1920):: zerodistflag + INTEGER(4),DIMENSION(7680,5760):: zerodistflag_fine - !INITIAL THE OUTPUT VARIABLES - zerodistflag=0 - zerodistflag_fine=0 - arfs_smap_temp=0 - arfs_wt=0 - arfs_fine_temp=0 - arfs_fine_wt=0 - !DISAGGREGATE ARFS TEMP TO FINEER GRID (~3km) - DO ii=1,1920 + !INITIAL THE OUTPUT VARIABLES + zerodistflag=0 + zerodistflag_fine=0 + arfs_smap_temp=0 + arfs_wt=0 + arfs_fine_temp=0 + arfs_fine_wt=0 + !DISAGGREGATE ARFS TEMP TO FINEER GRID (~3km) + DO ii=1,1920 !DO ii=522,524 - c = MINLOC(ABS(arfs_lat(ii)-arfs_fine_lat(:)),1) !Lat Direction r: row in fine scale - !PRINT*, arfs_fine_lat(r), arfs_lat(ii) - DO jj=1,2560 + lat1 = DBLE (arfs_lat(ii)*d2r) + c = MINLOC(ABS(arfs_lat(ii)-arfs_fine_lat(:)),1) !Lat Direction r: row in fine scale + !PRINT*, arfs_fine_lat(r), arfs_lat(ii) + cmin=c-3 ; IF (cmin < 1) cmin=1 + cmax=c+3 ; IF (cmax > size(arfs_fine_lat)) cmax=size(arfs_fine_lat) + + DO jj=1,2560 !DO jj=82,84 - r = MINLOC(ABS(arfs_lon(jj)-arfs_fine_lon(:)),1) !Lat Direction c: colume in fine scale - !PRINT*, arfs_fine_lon(c), arfs_lon(jj) - cmin=c-3 ; IF (cmin < 1) cmin=1 - cmax=c+3 ; IF (cmax > size(arfs_fine_lat)) cmax=size(arfs_fine_lat) - rmin=r-3 ; IF (rmin < 1) rmin=1 - rmax=r+3 ; IF (rmax > size(arfs_fine_lon)) rmax=size(arfs_fine_lon) - !PRINT*,'cmin cmax', cmin, cmax, size(arfs_fine_lon) - DO rr=rmin,rmax - DO cc=cmin,cmax - lat1 = DBLE (arfs_lat(ii)*d2r) - lon1 = DBLE (arfs_lon(jj)*d2r) - lat2 = DBLE (arfs_fine_lat(cc)*d2r) - lon2 = DBLE (arfs_fine_lon(rr)*d2r) - !PRINT*, lat1, lon1, lat2, lon2 - !gcdist = RE_KM * DACOS ( DSIN (lat1) * DSIN (lat2) + DCOS (lat1) * DCOS (lat2) * DCOS (lon1-lon2) ) !original - !---------------------------------------kyh - if(lat1.eq.lat2.and.lon1.eq.lon2) then - gcdist = 0. - else - gcdist = RE_KM * DACOS ( DSIN (lat1) * DSIN (lat2) + DCOS (lat1) * DCOS (lat2) * DCOS (lon1-lon2) ) - endif - !---------------------------------------kyh - !PRINT*, lat1, lon1, lat2, lon2, gcdist - !PRINT*, arfs_temp(ii,jj) - IF (gcdist < search_radius) THEN !RESAMPLE ONLY WITHIN THE SEARCH RANGE - !PRINT*,'HERE 0' - IF (gcdist < 0.0001D0) THEN !The TB is right on the grid center - zerodistflag_fine(rr,cc) = 1 - !PRINT*,'HERE 01' - IF ((ABS (arfs_temp(jj,ii)-(-9999.000)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(0) - arfs_fine_temp(rr,cc) = arfs_temp(jj,ii) ; arfs_fine_wt(rr,cc) = 1.0 - !PRINT*,'Here 1' - ENDIF - ELSE - IF (zerodistflag_fine(rr,cc).EQ.0) THEN !To maintain the corresponding pixel has the same value - IF ((ABS (arfs_temp(jj,ii)-(-9999.000)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(0) - arfs_fine_temp(rr,cc) = arfs_fine_temp(rr,cc)+arfs_temp(jj,ii) / SNGL(gcdist) - arfs_fine_wt(rr,cc) = arfs_fine_wt(rr,cc) + 1.0 / SNGL(gcdist) - !PRINT*,'Here 2' - ENDIF - ENDIF !(zerodistflag_fine(rr,cc).NE.0) - ENDIF !IF (gcdist < 0.0001D0) - ENDIF !(gcdist < search_radius) - ENDDO !cc=cmin,cmax - ENDDO !rr=rmin,rmax - - ENDDO !jj=1:2560 - ENDDO !ii=1,1920 - WHERE(arfs_fine_temp.NE.0.0.AND.arfs_fine_wt.NE.0.0) - arfs_fine_temp=arfs_fine_temp/arfs_fine_wt - ENDWHERE - !UPSCALL ARFS FINE TEMP TO ARFS GRID (at ~33km) - DO ii=1,1920 - c=3*ii-1 !3x3 FINE GRID EQ TO A ARFS GRID - DO jj=1,2560 - r=3*jj-1 !3x3 FINE GRID EQ TO A ARFS GRID - cmin=c-5 ; IF (cmin < 1) cmin=1 - cmax=c+5 ; IF (cmax > size(arfs_fine_lat)) cmax=size(arfs_fine_lat) - rmin=r-5 ; IF (rmin < 1) rmin=1 - rmax=r+5 ; IF (rmax > size(arfs_fine_lon)) rmax=size(arfs_fine_lon) - DO rr=rmin,rmax - DO cc=cmin,cmax - IF (abs(arfs_fine_temp(rr,cc)).GT.1.0D-7) THEN !DO WHEN T ~0 (change to fillvalue) - arfs_smap_temp(jj,ii)=arfs_smap_temp(jj,ii)+arfs_fine_temp(rr,cc); arfs_wt(jj,ii)=arfs_wt(jj,ii)+1.0 !Weight the point 1 - ENDIF - ENDDO !cc=cmin,cmax - ENDDO !rr=rmin,rmax - ENDDO !jj=1:2560 - ENDDO !ii=1:1920 - WHERE(arfs_smap_temp.NE.0.0.AND.arfs_wt.NE.0.0) - arfs_smap_temp=arfs_smap_temp/arfs_wt - ENDWHERE - WHERE(arfs_smap_temp.EQ.0.0) - arfs_smap_temp=-9999 - ENDWHERE + lon1 = DBLE (arfs_lon(jj)*d2r) + r = MINLOC(ABS(arfs_lon(jj)-arfs_fine_lon(:)),1) !Lat Direction c: colume in fine scale + !PRINT*, arfs_fine_lon(c), arfs_lon(jj) + rmin=r-3 ; IF (rmin < 1) rmin=1 + rmax=r+3 ; IF (rmax > size(arfs_fine_lon)) rmax=size(arfs_fine_lon) + !PRINT*,'cmin cmax', cmin, cmax, size(arfs_fine_lon) + !DO rr=rmin,rmax + ! DO cc=cmin,cmax + DO cc=cmin,cmax + lat2 = DBLE (arfs_fine_lat(cc)*d2r) + DO rr=rmin,rmax + lon2 = DBLE (arfs_fine_lon(rr)*d2r) + !PRINT*, lat1, lon1, lat2, lon2 + !gcdist = RE_KM * DACOS ( DSIN (lat1) * DSIN (lat2) + DCOS (lat1) * DCOS (lat2) * DCOS (lon1-lon2) ) !original + !---------------------------------------kyh + if(lat1.eq.lat2.and.lon1.eq.lon2) then + gcdist = 0. + else + gcdist = RE_KM * DACOS ( DSIN (lat1) * DSIN (lat2) + DCOS (lat1) * DCOS (lat2) * DCOS (lon1-lon2) ) + endif + !---------------------------------------kyh + !PRINT*, lat1, lon1, lat2, lon2, gcdist + !PRINT*, arfs_temp(ii,jj) + IF (gcdist < search_radius) THEN !RESAMPLE ONLY WITHIN THE SEARCH RANGE + !PRINT*,'HERE 0' + IF (gcdist < 0.0001D0) THEN !The TB is right on the grid center + zerodistflag_fine(rr,cc) = 1 + !PRINT*,'HERE 01' + IF ((ABS (arfs_temp(jj,ii)-(-9999.000)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(0) + arfs_fine_temp(rr,cc) = arfs_temp(jj,ii) ; arfs_fine_wt(rr,cc) = 1.0 + !PRINT*,'Here 1' + ENDIF + ELSE + IF (zerodistflag_fine(rr,cc).EQ.0) THEN !To maintain the corresponding pixel has the same value + IF ((ABS (arfs_temp(jj,ii)-(-9999.000)).GT.1.0D-7)) THEN !DO IF NOT FILLVALUE(0) + arfs_fine_temp(rr,cc) = arfs_fine_temp(rr,cc)+arfs_temp(jj,ii) / SNGL(gcdist) + arfs_fine_wt(rr,cc) = arfs_fine_wt(rr,cc) + 1.0 / SNGL(gcdist) + !PRINT*,'Here 2' + ENDIF + ENDIF !(zerodistflag_fine(rr,cc).NE.0) + ENDIF !IF (gcdist < 0.0001D0) + ENDIF !(gcdist < search_radius) + ENDDO !rr=rmin,rmax + ENDDO !cc=cmin,cmax + ENDDO !jj=1:2560 + ENDDO !ii=1,1920 + WHERE(arfs_fine_temp.NE.0.0.AND.arfs_fine_wt.NE.0.0) + arfs_fine_temp=arfs_fine_temp/arfs_fine_wt + ENDWHERE + !UPSCALL ARFS FINE TEMP TO ARFS GRID (at ~33km) + DO ii=1,1920 + c=3*ii-1 !3x3 FINE GRID EQ TO A ARFS GRID + cmin=c-5 ; IF (cmin < 1) cmin=1 + cmax=c+5 ; IF (cmax > size(arfs_fine_lat)) cmax=size(arfs_fine_lat) + DO jj=1,2560 + r=3*jj-1 !3x3 FINE GRID EQ TO A ARFS GRID + rmin=r-5 ; IF (rmin < 1) rmin=1 + rmax=r+5 ; IF (rmax > size(arfs_fine_lon)) rmax=size(arfs_fine_lon) + DO cc=cmin,cmax + DO rr=rmin,rmax + IF (abs(arfs_fine_temp(rr,cc)).GT.1.0D-7) THEN !DO WHEN T ~0 (change to fillvalue) + arfs_smap_temp(jj,ii)=arfs_smap_temp(jj,ii)+arfs_fine_temp(rr,cc); arfs_wt(jj,ii)=arfs_wt(jj,ii)+1.0 !Weight the point 1 + ENDIF + ENDDO !rr=rmin,rmax + ENDDO !cc=cmin,cmax + ENDDO !jj=1:2560 + ENDDO !ii=1:1920 + WHERE(arfs_smap_temp.NE.0.0.AND.arfs_wt.NE.0.0) + arfs_smap_temp=arfs_smap_temp/arfs_wt + ENDWHERE + WHERE(arfs_smap_temp.EQ.0.0) + arfs_smap_temp=-9999 + ENDWHERE - END SUBROUTINE RESAMPLETEMP + END SUBROUTINE RESAMPLETEMP END MODULE invdist_temp2smap diff --git a/ldt/SMAP_E_OPL/readLISTeff_usaf.F90 b/ldt/SMAP_E_OPL/readLISTeff_usaf.F90 new file mode 100644 index 000000000..b26bb0c39 --- /dev/null +++ b/ldt/SMAP_E_OPL/readLISTeff_usaf.F90 @@ -0,0 +1,405 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT---------------------- + +#include "LDT_misc.h" +#include "LDT_NetCDF_inc.h" + +!BOP +! +! !ROUTINE: readLIS_Teff_usaf +! \label{readLISTeff} +! +! !REVISION HISTORY: +! 14 FEB 2023: Eric Kemp, Initial Specification +! 1 July 2023: Mahdi Navari, Modified the code to read the unperturbed +! ensemble. This is a temporary fix to overcome the bias in +! soil temperature from a bug in SMAP soil moisture DA +! PR#1385 +! +! !INTERFACE: +subroutine readLIS_Teff_usaf(n, yyyymmdd, hh, Orbit, teff, rc) +! !USES: + use LDT_coreMod + use LDT_domainMod + use LDT_logMod + use LDT_smap_e_oplMod + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + character(8), intent(in) :: yyyymmdd + character(2), intent(in) :: hh + character(1), intent(in) :: Orbit + real, intent(out) :: teff(LDT_rc%lnc(n),LDT_rc%lnr(n)) + integer, intent(out) :: rc + +!EOP + integer :: c,r + real :: tsoil(LDT_rc%lnc(n),LDT_rc%lnr(n),4) + character(255) :: fname + logical :: file_exists + integer :: rc1 + integer :: gid + integer :: nens + integer, allocatable :: str_tind(:) + integer, allocatable :: ntiles_pergrid(:) + real, parameter :: kk = 1.007 + real, parameter :: cc_6am = 0.246 + real, parameter :: cc_6pm = 1.000 + + external :: create_LISsoilT_filename_usaf + external :: read_LIStsoil_data_usaf + + ! Initializations + teff = LDT_rc%udef + tsoil = LDT_rc%udef + rc = 1 ! Assume error by default, update below + + ! Set up basic info on Air Force product + nens = SMAPeOPL%num_ens + allocate(str_tind(LDT_rc%gnc(n) * LDT_rc%gnr(n))) + if (SMAPeOPL%ntiles_pergrid .ne. 1) then + write(LDT_logunit,*) & + '[ERR] Current SMAP_E_OPL code assumes ntiles_pergrid = 1' + write(LDT_logunit,*) & + '[ERR] Actual value is ', SMAPeOPL%ntiles_pergrid + call LDT_endrun() + end if + do gid = 1, (LDT_rc%gnc(n) * LDT_rc%gnr(n)) + str_tind(gid) = ((gid - 1) * nens) + 1 + end do + allocate(ntiles_pergrid(LDT_rc%gnc(n) * LDT_rc%gnr(n))) + ntiles_pergrid = SMAPeOPL%ntiles_pergrid + + call create_LISsoilT_filename_usaf(SMAPeOPL%LISdir, & + yyyymmdd, hh, fname) + + inquire(file=trim(fname), exist=file_exists) + if (file_exists) then + write(LDT_logunit,*) '[INFO] Reading ', trim(fname) + call read_LIStsoil_data_usaf(n, SMAPeOPL%num_tiles, str_tind, & + ntiles_pergrid, nens, & + fname, tsoil, rc1) + if (rc1 .ne. 0) then + write(LDT_logunit,*) '[ERR] Cannot read from ', trim(fname) + return + endif + write(LDT_logunit,*) '[INFO] Finished reading ', trim(fname) + + ! Calculate effective temperature + do r = 1, LDT_rc%lnr(n) + do c = 1, LDT_rc%lnc(n) + if (Orbit == "D") then + if (tsoil(c,r,1) > 273.15 .and. tsoil(c,r,2) > 273.15) then + teff(c,r) = kk * & + (cc_6am * tsoil(c,r,1) + (1 - cc_6am) * tsoil(c,r,2)) + endif + elseif (Orbit == "A") then + if (tsoil(c,r,1) > 273.15 .and. tsoil(c,r,2) > 273.15) then + teff(c,r) = kk * & + (cc_6pm * tsoil(c,r,1) + (1 - cc_6pm) * tsoil(c,r,2)) + endif + endif + end do + end do + rc = 0 + end if + + if (allocated(str_tind)) deallocate(str_tind) + if (allocated(ntiles_pergrid)) deallocate(ntiles_pergrid) + +end subroutine readLIS_Teff_usaf + +!BOP +! +! !ROUTINE: read_LIStsoil_data_usaf +! \label{read_LIStsoil_data_usaf} +! +! !INTERFACE: +subroutine read_LIStsoil_data_usaf(n, ntiles, str_tind, ntiles_pergrid, nens, & + fname, tsoil, rc) +! +! !USES: + use LDT_logMod + use LDT_coreMod + use LDT_smap_e_oplMod +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none + +! +! !INPUT PARAMETERS: +! + integer, intent(in) :: n + integer, intent(in) :: ntiles + integer, intent(in) :: str_tind(LDT_rc%gnc(n) * LDT_rc%gnc(n)) + integer, intent(in) :: ntiles_pergrid(LDT_rc%gnc(n) * LDT_rc%gnc(n)) + integer, intent(in) :: nens + character(*), intent(in) :: fname + real, intent(inout) :: tsoil(LDT_rc%lnc(n),LDT_rc%lnr(n),4) + integer, intent(out) :: rc +!EOP + + integer :: ios, ncid + integer :: ntiles_dimid, ntiles_local + integer :: SoilTemp_profiles_dimid, SoilTemp_profiles + integer :: SoilTemp_inst_id + real, allocatable :: SoilTemp_inst_tiles(:,:) + real :: SoilTemp_inst_ensmean_1layer(LDT_rc%gnc(n), LDT_rc%gnr(n)) + integer :: k, c, r + + external :: calc_gridded_ensmean_1layer + + rc = 1 ! Initialize as error, reset near bottom + + ! Sanity checks + if (LDT_rc%gnc(n) .ne. LDT_rc%lnc(n)) then + write(LDT_logunit,*)'[ERR] Mismatched dimensions!' + write(LDT_logunit,*)'[ERR] LDT_rc%gnc(n) = ', LDT_rc%gnc(n) + write(LDT_logunit,*)'[ERR] LDT_rc%lnc(n) = ', LDT_rc%lnc(n) + call LDT_endrun() + end if + if (LDT_rc%gnr(n) .ne. LDT_rc%lnr(n)) then + write(LDT_logunit,*)'[ERR] Mismatched dimensions!' + write(LDT_logunit,*)'[ERR] LDT_rc%gnr(n) = ', LDT_rc%gnr(n) + write(LDT_logunit,*)'[ERR] LDT_rc%lnr(n) = ', LDT_rc%lnr(n) + call LDT_endrun() + end if + +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + ios = nf90_open(path=trim(fname), mode=NF90_NOWRITE ,ncid=ncid) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Error opening file' // trim(fname) + return + end if + + ios = nf90_inq_dimid(ncid, "ntiles", ntiles_dimid) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot find ntiles in ' // trim(fname) + ios = nf90_close(ncid=ncid) + return + end if + + ios = nf90_inquire_dimension(ncid, ntiles_dimid, len=ntiles_local) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot get dimension ntiles in ' & + // trim(fname) + ios = nf90_close(ncid=ncid) + return + end if + + ! if (ntiles .ne. LDT_rc%glbntiles_red(n)) then + ! write(LDT_logunit,*)'[ERR] Dimension mismatch!' + ! write(LDT_logunit,*)'[ERR] ntiles = ', ntiles + ! write(LDT_logunit,*)'[ERR] LDT_rc%glbntiles_red(n) = ', & + ! LDT_rc%glbntiles_red(n) + ! call LDT_endrun() + ! end if + if (ntiles_local .ne. ntiles) then + write(LDT_logunit,*)'[ERR] Dimension mismatch!' + write(LDT_logunit,*)'[ERR] ntiles = ', ntiles_local + write(LDT_logunit,*)'[ERR] Expected ', ntiles + call LDT_endrun() + end if + + ios = nf90_inq_dimid(ncid, "SoilTemp_profiles", SoilTemp_profiles_dimid) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot find SoilTemp_profiles in ' & + // trim(fname) + ios = nf90_close(ncid=ncid) + return + end if + + ios = nf90_inquire_dimension(ncid, SoilTemp_profiles_dimid, & + len=SoilTemp_profiles) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot get dimension SoilTemp_profiles in ' & + // trim(fname) + ios = nf90_close(ncid=ncid) + return + end if + if (SoilTemp_profiles .ne. 4) then + write(LDT_logunit,*)'[ERR] Dimension mismatch!' + write(LDT_logunit,*)'[ERR] SoilTemp_profiles should be 4, but is ', & + SoilTemp_profiles + call LDT_endrun() + end if + + ios = nf90_inq_varid(ncid, 'SoilTemp_inst', SoilTemp_inst_id) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot find SoilTemp_inst in ' // trim(fname) + ios = nf90_close(ncid=ncid) + return + end if + + allocate(SoilTemp_inst_tiles(ntiles, SoilTemp_profiles)) + SoilTemp_inst_tiles = 0 + + ios = nf90_get_var(ncid, SoilTemp_inst_id, SoilTemp_inst_tiles, & + start=(/1, 1/), & + count=(/ntiles, SoilTemp_profiles/)) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot read SoilTemp_inst in ' // trim(fname) + deallocate(SoilTemp_inst_tiles) + ios = nf90_close(ncid=ncid) + return + end if + + ! Calculate ensemble mean in 2d grid space, for each soil layer + do k = 1, SoilTemp_profiles + !call calc_gridded_ensmean_1layer(n, ntiles, str_tind, ntiles_pergrid, & + ! nens, & + ! SoilTemp_inst_tiles(:,k), & + ! SoilTemp_inst_ensmean_1layer) + + call calc_gridded_lastens_1layer(n, ntiles, str_tind, ntiles_pergrid, & + nens, & + SoilTemp_inst_tiles(:,k), & + SoilTemp_inst_ensmean_1layer) ! Note to minimize the code chnage: SoilTemp_inst_ensmean_1layer is actually + ! SoilTemp_inst_lastens_1layer (ens #12) + + do r = 1, LDT_rc%lnr(n) + do c = 1, LDT_rc%lnc(n) + if (SoilTemp_inst_ensmean_1layer(c,r) > 0) then + tsoil(c,r,k) = SoilTemp_inst_ensmean_1layer(c,r) + end if + end do + end do + end do + + ! Clean up + deallocate(SoilTemp_inst_tiles) + + ios = nf90_close(ncid=ncid) + if (ios .ne. 0) then + write(LDT_logunit,*) '[WARN] Error closing ' // trim(fname) + return + end if + + rc = 0 ! No error detected + +#endif + +end subroutine read_LIStsoil_data_usaf + +! Subroutine for calculating 2d gridded ensemble mean for a single soil layer, +! from tiled data. +subroutine calc_gridded_ensmean_1layer(n, ntiles, str_tind, ntiles_pergrid, & + nens, gvar_tile, gvar) + + ! Imports + use LDT_coreMod, only: LDT_rc, LDT_domain, LDT_masterproc + use LDT_logMod, only: LDT_logunit + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: n + integer, intent(in) :: ntiles + integer, intent(in) :: str_tind(LDT_rc%gnc(n) * LDT_rc%gnr(n)) + integer, intent(in) :: ntiles_pergrid(LDT_rc%gnc(n) * LDT_rc%gnr(n)) + integer, intent(in) :: nens + real, intent(in) :: gvar_tile(ntiles) + real, intent(out) :: gvar(LDT_rc%gnc(n), LDT_rc%gnr(n)) + + ! Locals + integer :: m, r, c, gid, stid, tid + + gvar = 0 + if (LDT_masterproc) then + do r = 1, LDT_rc%gnr(n) + do c = 1, LDT_rc%gnc(n) + gid = c + ((r-1) * LDT_rc%gnc(n)) + stid = str_tind(gid) + if (ntiles_pergrid(gid) > 0) then + do m = 1, nens + tid = stid + m - 1 + gvar(c,r) = gvar(c,r) + gvar_tile(tid) + enddo + gvar(c,r) = gvar(c,r) / nens + end if + end do + end do + end if +end subroutine calc_gridded_ensmean_1layer + +! Subroutine for extracting last ensemble member for a single soil layer, +! from tiled data. + subroutine calc_gridded_lastens_1layer(n, ntiles, str_tind, ntiles_pergrid, & + nens, gvar_tile, gvar) + + ! Imports + use LDT_coreMod, only: LDT_rc, LDT_domain, LDT_masterproc + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: n + integer, intent(in) :: ntiles + !real, intent(in) :: gvar_tile(LDT_rc%glbntiles_red(n)) + integer, intent(in) :: str_tind(LDT_rc%gnc(n) * LDT_rc%gnr(n)) + integer, intent(in) :: ntiles_pergrid(LDT_rc%gnc(n) * LDT_rc%gnr(n)) + real, intent(out) :: gvar(LDT_rc%gnc(n), LDT_rc%gnr(n)) + real, intent(in) :: gvar_tile(ntiles) + integer, intent(in) :: nens + + ! Locals + integer :: m, r, c, gid, stid, tid + + if (LDT_masterproc) then + gvar = 0 + do r = 1, LDT_rc%gnr(n) + do c = 1, LDT_rc%gnc(n) + gid = c + ((r-1) * LDT_rc%gnc(n)) + stid = str_tind(gid) + if (ntiles_pergrid(gid) > 0) then + m = nens + tid = stid + m - 1 + gvar(c,r) = gvar_tile(tid) + end if + end do + end do + end if + end subroutine calc_gridded_lastens_1layer + +!BOP +! !ROUTINE: create_LISsoilT_filename_usaf +! \label{create_LISsoilT_filename_usaf} +! +! !INTERFACE: +subroutine create_LISsoilT_filename_usaf(LISdir, yyyymmdd, hh, filename) + !USES: + + implicit none + !ARGUMENTS: + character(*), intent(in) :: LISdir + character(8), intent(in) :: yyyymmdd + character(2), intent(in) :: hh + character(*), intent(out) :: filename +!EOP + + filename = trim(LISdir) & + // '/PS.AFWA' & + // '_SC.U' & + // '_DI.C' & + // '_DC.ANLYS' & + // '_GP.LIS' & + // '_GR.C0P09DEG' & + // '_AR.GLOBAL' & + // '_PA.03-HR-SUM' & + // '_DD.' // yyyymmdd & + // '_DT.' // hh // '00' & + // '_DF.nc' + +end subroutine create_LISsoilT_filename_usaf diff --git a/ldt/SMAP_E_OPL/readUSAFSI.F90 b/ldt/SMAP_E_OPL/readUSAFSI.F90 new file mode 100644 index 000000000..a905b6f2c --- /dev/null +++ b/ldt/SMAP_E_OPL/readUSAFSI.F90 @@ -0,0 +1,215 @@ +!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- +! NASA Goddard Space Flight Center +! Land Information System Framework (LISF) +! Version 7.5 +! +! Copyright (c) 2020 United States Government as represented by the +! Administrator of the National Aeronautics and Space Administration. +! All Rights Reserved. +!-------------------------END NOTICE -- DO NOT EDIT---------------------- + +#include "LDT_misc.h" +#include "LDT_NetCDF_inc.h" +!BOP +! +! !ROUTINE: readUSAFSI +! \label{readUSAFSI} +! +! !REVISION HISTORY: +! 13 Feb 2023: Eric Kemp, Initial Specification +! +! !INTERFACE: +subroutine readUSAFSI(n, yyyymmdd, hh, SnowDepth, rc) + +! !USES: + use LDT_coreMod + use LDT_logMod + use LDT_smap_e_oplMod + + implicit none + +! !ARGUMENTS: + integer, intent(in) :: n + character*8, intent(in) :: yyyymmdd + character*2, intent(in) :: hh + real, intent(out):: SnowDepth(LDT_rc%lnc(n),LDT_rc%lnr(n)) + integer, intent(out) :: rc + +!EOP + character*255 :: fname + character*8 :: yyyymmdd_snow + character*4 :: yyyy_snow + character*2 :: mm_snow, dd_snow, hh_snow + logical :: file_exists + integer :: yr, mo, da, hr + integer :: nn + integer :: ierr + integer :: rc1 + + external :: create_USAFSI_filename + external :: read_USAFSI_data + + rc = 1 ! Default to error, will update below if USAFSI file read in. + SnowDepth = LDT_rc%udef + + yyyymmdd_snow = yyyymmdd + yyyy_snow = yyyymmdd(1:4) + mm_snow = yyyymmdd(5:6) + dd_snow = yyyymmdd(7:8) + hh_snow = hh + + read(yyyy_snow, *, iostat=ierr) yr + read(mm_snow, *, iostat=ierr) mo + read(dd_snow, *, iostat=ierr) da + read(hh_snow, *, iostat=ierr) hr + + nn = 0 + do while (nn <= 24) + + call create_USAFSI_filename(SMAPeOPL%LISsnowdir, & + yyyymmdd_snow, hh_snow, fname) + inquire(file=trim(fname), exist=file_exists) + + if (file_exists) then + write(LDT_logunit,*) '[INFO] Reading snow depth from ', trim(fname) + call read_USAFSI_data(n, fname, SnowDepth, rc1) + if (rc1 == 0) then + write(LDT_logunit,*) '[INFO] Finished reading snow outputs from ', & + trim(fname) + rc = 0 + exit + end if + end if + + ! Go back one hour + hr = hr - 1 + if (hr < 0) then + da = da - 1 + hr = 23 + if (da == 0) then + mo = mo - 1 + if (mo == 0) then + yr = yr - 1 + mo = 12 + da = 31 + else + if (mo == 1 .or. & + mo == 3 .or. & + mo == 5 .or. & + mo == 7 .or. & + mo == 8 .or. & + mo == 10 .or. & + mo == 12) then + da = 31 + elseif (mo == 2) then + if (mod(yr,4) == 0) then + da = 29 + else + da = 28 + endif + else + da = 30 + endif + endif + endif + endif + + write(unit=yyyy_snow, fmt='(i4.4)') yr + write(unit=mm_snow, fmt='(i2.2)') mo + write(unit=dd_snow, fmt='(i2.2)') da + write(unit=hh_snow, fmt='(i2.2)') hr + yyyymmdd_snow = trim(yyyy_snow) // trim(mm_snow) // trim(dd_snow) + + nn = nn + 1 + + end do +end subroutine readUSAFSI + +!BOP +! +! !ROUTINE: read_USAFSI_data +! \label{read_USAFSI_data} +! +! !INTERFACE: +subroutine read_USAFSI_data(n, fname, SnowDepth, rc) +! +! !USES: + use LDT_logMod + use LDT_coreMod + use LDT_smap_e_oplMod +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif + + implicit none +! +! !INPUT PARAMETERS: +! + integer, intent(in) :: n + character(*), intent(in) :: fname + real, intent(inout) :: SnowDepth(LDT_rc%lnc(n),LDT_rc%lnr(n)) + integer, intent(out) :: rc +!EOP + + integer :: ios, nid + integer :: snowdepthid + + rc = 1 ! Initialize as error, reset near bottom + +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + ios = nf90_open(path=trim(fname), mode=NF90_NOWRITE, ncid=nid) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Error opening file' // trim(fname) + return + end if + + ios = nf90_inq_varid(nid, 'snoanl', snowdepthid) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot find snoanl in ' // trim(fname) + ios = nf90_close(ncid=nid) + return + end if + + ios = nf90_get_var(nid, snowdepthid, SnowDepth, & + start=(/1, 1/), & + count=(/LDT_rc%lnc(n), LDT_rc%lnr(n)/)) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Cannot read snoanl in ' // trim(fname) + ios = nf90_close(ncid=nid) + return + end if + + ios = nf90_close(ncid=nid) + if (ios .ne. 0) then + write(LDT_logunit,*)'[WARN] Error closing ' // trim(fname) + return + end if + + rc = 0 ! No error detected! + +#endif + +end subroutine read_USAFSI_data + +!BOP +! !ROUTINE: create_USAFSI_filename +! \label{create_USAFSI_filename} +! +! !INTERFACE: +subroutine create_USAFSI_filename(LISdir, yyyymmdd, hh, filename) +! !USES: + + implicit none + +! !ARGUMENTS: + character(*), intent(in) :: LISdir + character(8), intent(in) :: yyyymmdd + character(2), intent(in) :: hh + character(*), intent(out) :: filename +!EOP + + filename = trim(LISdir) // '/USAFSI_' // trim(yyyymmdd) & + // trim(hh) // '.nc' + +end subroutine create_USAFSI_filename + diff --git a/ldt/SMAP_E_OPL/scale_teff.F90 b/ldt/SMAP_E_OPL/scale_teff.F90 index 1691b59ee..c7c97fd07 100644 --- a/ldt/SMAP_E_OPL/scale_teff.F90 +++ b/ldt/SMAP_E_OPL/scale_teff.F90 @@ -12,12 +12,13 @@ ! ! !ROUTINE: scale_teff ! \label{scale_teff} -! -! !REVISION HISTORY: +! +! !REVISION HISTORY: ! 12 JAN 2022: Yonghwan Kwon, Initial Specification -! -! !INTERFACE: -subroutine scale_teff(n,Orbit,teff_01,teff_02) +! 21 Feb 2023: Eric Kemp, added third time level +! +! !INTERFACE: +subroutine scale_teff(n, Orbit, teff_01, teff_02, teff_03) ! !USES: use LDT_coreMod use LDT_smap_e_oplMod @@ -25,51 +26,52 @@ subroutine scale_teff(n,Orbit,teff_01,teff_02) implicit none ! !ARGUMENTS: integer, intent(in) :: n - character*1 :: Orbit + character*1, intent(in) :: Orbit real :: teff_01(LDT_rc%lnc(n),LDT_rc%lnr(n)) real :: teff_02(LDT_rc%lnc(n),LDT_rc%lnr(n)) + real :: teff_03(LDT_rc%lnc(n),LDT_rc%lnr(n)) !EOP integer :: t integer :: col, row real :: mu_ref, mu_lis real :: sigma_ref, sigma_lis - real :: teff_01_tmp, teff_02_tmp + real :: teff_01_tmp, teff_02_tmp, teff_03_tmp - do t=1,SMAPeOPL%ngrid + do t = 1, SMAPeOPL%ngrid col = SMAPeOPL%grid_col(t) row = SMAPeOPL%grid_row(t) - if(Orbit.eq."D") then + if (Orbit.eq."D") then mu_ref = SMAPeOPL%mu_6am_ref(t) mu_lis = SMAPeOPL%mu_6am_lis(t) sigma_ref = SMAPeOPL%sigma_6am_ref(t) - sigma_lis = SMAPeOPL%sigma_6am_lis(t) - elseif(Orbit.eq."A") then + sigma_lis = SMAPeOPL%sigma_6am_lis(t) + elseif (Orbit.eq."A") then mu_ref = SMAPeOPL%mu_6pm_ref(t) mu_lis = SMAPeOPL%mu_6pm_lis(t) sigma_ref = SMAPeOPL%sigma_6pm_ref(t) sigma_lis = SMAPeOPL%sigma_6pm_lis(t) endif - if(teff_01(col,row).ne.-9999.0) then - if(mu_ref.gt.0.and.& - mu_lis.gt.0.and.& - sigma_ref.gt.0.and.& - sigma_lis.gt.0) then + if (teff_01(col,row).ne.-9999.0) then + if (mu_ref.gt.0.and.& + mu_lis.gt.0.and.& + sigma_ref.gt.0.and.& + sigma_lis.gt.0) then teff_01_tmp = (teff_01(col,row) - mu_lis) * & sigma_ref/sigma_lis + mu_ref - teff_01(col,row) = teff_01_tmp + teff_01(col,row) = teff_01_tmp endif endif - if(teff_02(col,row).ne.-9999.0) then - if(mu_ref.gt.0.and.& - mu_lis.gt.0.and.& - sigma_ref.gt.0.and.& - sigma_lis.gt.0) then + if (teff_02(col,row).ne.-9999.0) then + if (mu_ref.gt.0.and.& + mu_lis.gt.0.and.& + sigma_ref.gt.0.and.& + sigma_lis.gt.0) then teff_02_tmp = (teff_02(col,row) - mu_lis) * & sigma_ref/sigma_lis + mu_ref @@ -78,6 +80,19 @@ subroutine scale_teff(n,Orbit,teff_01,teff_02) endif endif + if (teff_03(col,row).ne.-9999.0) then + if (mu_ref.gt.0.and.& + mu_lis.gt.0.and.& + sigma_ref.gt.0.and.& + sigma_lis.gt.0) then + + teff_03_tmp = (teff_03(col,row) - mu_lis) * & + sigma_ref/sigma_lis + mu_ref + + teff_03(col,row) = teff_03_tmp + endif + endif + enddo end subroutine scale_teff diff --git a/ldt/core/LDT_SigmaMod.F90 b/ldt/core/LDT_SigmaMod.F90 index 4b9ac9df3..f806b507b 100644 --- a/ldt/core/LDT_SigmaMod.F90 +++ b/ldt/core/LDT_SigmaMod.F90 @@ -178,14 +178,16 @@ subroutine diagnoseSingleSigma(n,obs, metrics) gmt = LDT_rc%hr call LDT_gmt2localtime(gmt, lon, lhour, zone) - if(lhour.eq.6) then + !if(lhour.eq.6) then ! Orig + if(lhour > 4 .and. lhour < 8) then ! EMK for 3-hrly data metrics%sxx_sigma_6am(t,j,k) = metrics%sxx_sigma_6am(t,j,k) + & obs%value(t1,k)*obs%value(t1,k) metrics%sx_sigma_6am(t,j,k) = metrics%sx_sigma_6am(t,j,k) + & obs%value(t1,k) metrics%count_sigma_6am(t,j,k) = & metrics%count_sigma_6am(t,j,k) + 1 - elseif(lhour.eq.18) then + !elseif(lhour.eq.18) then ! Orig + elseif (lhour > 16 .and. lhour < 20) then ! EMK for 3-hrly data metrics%sxx_sigma_6pm(t,j,k) = metrics%sxx_sigma_6pm(t,j,k) + & obs%value(t1,k)*obs%value(t1,k) metrics%sx_sigma_6pm(t,j,k) = metrics%sx_sigma_6pm(t,j,k) + & diff --git a/lis/dataassim/obs/SMAP_E_OPLsm/read_SMAPEOPLsm.F90 b/lis/dataassim/obs/SMAP_E_OPLsm/read_SMAPEOPLsm.F90 index 049de1ce1..ee45b14d1 100644 --- a/lis/dataassim/obs/SMAP_E_OPLsm/read_SMAPEOPLsm.F90 +++ b/lis/dataassim/obs/SMAP_E_OPLsm/read_SMAPEOPLsm.F90 @@ -14,6 +14,7 @@ ! ! !REVISION HISTORY: ! 6 Jun 2022: Yonghwan Kwon; Updated for use with SMAP_E_OPL soil moisture +! 23 Feb 2023: Eric Kemp; Updated to read netCDF file. ! ! !INTERFACE: subroutine read_SMAPEOPLsm(n, k, OBS_State, OBS_Pert_State) @@ -143,8 +144,11 @@ subroutine read_SMAPEOPLsm(n, k, OBS_State, OBS_Pert_State) write(hh,'(i2.2)') LIS_rc%hr if(LIS_masterproc) then + !list_files = trim(smobsdir)//'/ARFS_SM_*' & + ! //trim(yyyymmdd)//'T'//trim(hh)//'*.dat' + !EMK...Use netCDF list_files = trim(smobsdir)//'/ARFS_SM_*' & - //trim(yyyymmdd)//'T'//trim(hh)//'*.dat' + //trim(yyyymmdd)//'T'//trim(hh)//'*.nc' write(LIS_logunit,*) & '[INFO] Searching for ',trim(list_files) rc = create_filelist(trim(list_files)//char(0), & @@ -426,7 +430,9 @@ end subroutine read_SMAPEOPLsm subroutine read_SMAPEOPLsm_data(n, k,fname, smobs_inp, time) ! ! !USES: - +#if (defined USE_NETCDF3 || defined USE_NETCDF4) + use netcdf +#endif use LIS_coreMod use LIS_logMod use LIS_timeMgrMod @@ -452,28 +458,205 @@ subroutine read_SMAPEOPLsm_data(n, k,fname, smobs_inp, time) integer :: ios, nid integer :: c,r integer :: ftn1 - - ftn1 = LIS_getNextUnitNumber() - open(unit=ftn1,file=fname,form='unformatted',access='direct',recl=4*nc*nr,status='old') - read(ftn1, rec=1) sm_raw - close(1) - call LIS_releaseUnitNumber(ftn1) - + ! EMK + logical :: file_exists + character(255) :: map_projection + integer :: ncid, dim_ids(3), var_id + integer :: ntime, nlat, nlon + real, allocatable :: tmp(:,:,:) + integer :: rc + + ! Old code to use binary data + !ftn1 = LIS_getNextUnitNumber() + !open(unit=ftn1,file=fname,form='unformatted',access='direct',recl=4*nc*nr,status='old') + + ! read(ftn1, rec=1) sm_raw + ! close(1) + ! call LIS_releaseUnitNumber(ftn1) + + ! sm_data_b = .false. + + ! do r=1,SMAPEOPLsm_struc(n)%nr + ! do c=1,SMAPEOPLsm_struc(n)%nc + ! if (sm_raw(c,r)>=0.and.& + ! sm_raw(c,r)<=100) then + + ! sm_in(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = sm_raw(c,r) + ! sm_data_b(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = .true. + ! else + ! sm_in(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = LIS_rc%udef + ! sm_data_b(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = .false. + ! endif + ! enddo + ! enddo + + ! EMK...Read netCDF data + sm_in = LIS_rc%udef sm_data_b = .false. - do r=1,SMAPEOPLsm_struc(n)%nr - do c=1,SMAPEOPLsm_struc(n)%nc - if (sm_raw(c,r)>=0.and.& - sm_raw(c,r)<=100) then +#if(defined USE_NETCDF3 || defined USE_NETCDF4) + + ! See if the file exists. + inquire(file=trim(fname), exist=file_exists) + if (.not. file_exists) then + write(LIS_logunit,*)'[ERR] Cannot find ', trim(fname) + return + end if + + ! Open the file + rc = nf90_open(path=trim(fname), & + mode=NF90_NOWRITE, & + ncid=ncid) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot open ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + return + end if + + ! Read the map projection + rc = nf90_get_att(ncid=ncid, & + varid=NF90_GLOBAL, & + name='MAP_PROJECTION', & + values=map_projection) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read MAP_PROJECTION from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + + ! Sanity check map projection + ! TODO: Support other map projections + if (trim(map_projection) .ne. 'EQUIDISTANT CYLINDRICAL') then + write(LIS_logunit,*) & + '[ERR] Unrecognized map projection found in SMAP file!' + write(LIS_logunit,*) '[ERR] Expected EQUIDISTANT CYLINDRICAL' + write(LIS_logunit,*) '[ERR] Found ',trim(map_projection) + write(LIS_logunit,*) '[ERR] LIS will skip file ', trim(fname) + rc = nf90_close(ncid) + return + end if + + ! Get dimension IDs + rc = nf90_inq_dimid(ncid=ncid, & + name='time', & + dimid=dim_ids(3)) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read time dimension from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + rc = nf90_inq_dimid(ncid=ncid, & + name='lat', & + dimid=dim_ids(2)) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read lat dimension from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + rc = nf90_inq_dimid(ncid=ncid, & + name='lon', & + dimid=dim_ids(1)) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read lon dimension from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + + ! Get actual dimension sizes + rc = nf90_inquire_dimension(ncid=ncid, & + dimid=dim_ids(3), & + len=ntime) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read time dimension from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + rc = nf90_inquire_dimension(ncid=ncid, & + dimid=dim_ids(2), & + len=nlat) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read lat dimension from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + rc = nf90_inquire_dimension(ncid=ncid, & + dimid=dim_ids(1), & + len=nlon) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read lon dimension from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + + ! Sanity check the dimensions + if (ntime .ne. 1) then + write(LIS_logunit,*)'[ERR] Expected time dimension to be 1' + write(LIS_logunit,*)'[ERR] Found ', ntime, ' from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + if (nlat .ne. SMAPEOPLsm_struc(n)%nr) then + write(LIS_logunit,*)'[ERR] Expected lat dimension to be ', & + SMAPEOPLsm_struc(n)%nr + write(LIS_logunit,*)'[ERR] Found ', nlat, ' from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + if (nlon .ne. SMAPEOPLsm_struc(n)%nc) then + write(LIS_logunit,*)'[ERR] Expected lon dimension to be ', & + SMAPEOPLsm_struc(n)%nc + write(LIS_logunit,*)'[ERR] Found ', nlon, ' from ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + + ! Fetch the variable id + rc = nf90_inq_varid(ncid=ncid, & + name='arfs_sm', & + varid=var_id) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read arfs_sm ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + return + end if + + ! Read the retrievals + allocate(tmp(nlon, nlat, ntime)) + rc = nf90_get_var(ncid=ncid, & + varid=var_id, & + values=tmp) + if (rc .ne. 0) then + write(LIS_logunit,*)'[ERR] Cannot read arfs_sm ', trim(fname) + write(LIS_logunit,*)'[ERR] LIS will continue...' + rc = nf90_close(ncid) + deallocate(tmp) + return + end if + rc = nf90_close(ncid) + + do r = 1, nlat + do c = 1, nlon + if (tmp(c,r,1) >= 0 .and. & + tmp(c,r,1) <= 1) then + sm_in(c + (r-1)*nc) = tmp(c,r,1)*100 + sm_data_b(c + (r-1)*nc) = .true. + end if + end do + end do + deallocate(tmp) - sm_in(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = sm_raw(c,r) - sm_data_b(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = .true. - else - sm_in(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = LIS_rc%udef - sm_data_b(c+(r-1)*SMAPEOPLsm_struc(n)%nc) = .false. - endif - enddo - enddo +#endif !-------------------------------------------------------------------------- ! Interpolate to the LIS running domain @@ -495,7 +678,7 @@ subroutine read_SMAPEOPLsm_data(n, k,fname, smobs_inp, time) sm_data_b,sm_in, smobs_b_ip, smobs_ip) endif -!overwrite the input data +!overwrite the input data do r=1,LIS_rc%obs_lnr(k) do c=1,LIS_rc%obs_lnc(k) if(smobs_ip(c+(r-1)*LIS_rc%obs_lnc(k)).ne.-9999.0) then