From 27467c49a0d02dd42da302951c237ad236928594 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 2 Oct 2024 07:33:10 -0400 Subject: [PATCH 1/2] Fixes #55. ifdef out HDF4 reader --- CHANGELOG.md | 6 +- .../clsm_ensupd_read_obs.F90 | 7159 +++++++++-------- 2 files changed, 3587 insertions(+), 3578 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 41bf48d..9032d5a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - More optimal distribution of tiles on processors for cubed-sphere tile space. +- Ifdef'd out HDF4 reader code in `GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90` due to lack of HDF4 Fortran + support. ### Fixed @@ -43,9 +45,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Moved external `GEOSgcm_GridComp` repository to under `GEOSldas/src/Components` for consistency with directory structure of GEOSgcm and GEOSadas ([PR #27](https://github.com/GEOS-ESM/GEOSldas_GridComp/pull/27), [PR #30](https://github.com/GEOS-ESM/GEOSldas_GridComp/pull/30)). -- Changed lenkf.j.template to python string ([PR #16](https://github.com/GEOS-ESM/GEOSldas_GridComp/pull/16)). +- Changed lenkf.j.template to python string ([PR #16](https://github.com/GEOS-ESM/GEOSldas_GridComp/pull/16)). + - ----------------------------- ## [v1.0.1] - 2024-04-10 diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index e87b3df..f7eade5 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -1,10 +1,10 @@ -! -! this file contains subroutines for reading and processing observations +! +! this file contains subroutines for reading and processing observations ! for the GEOS5 land EnKF update algorithm ! ! reichle, 27 Jan 2005 ! reichle, 10 Jan 2011 - replaced "UVA" with "LPRM" -! reichle, 1 Jul 2015 - clarified definition of obs time stamp +! reichle, 1 Jul 2015 - clarified definition of obs time stamp ! (J2000 seconds w/ 'TT12' epoch) ! lcandre2, 10 Jul 2021 - confirmed and cleaned up MODIS obs @@ -15,7 +15,7 @@ ! ********************************************************************* module clsm_ensupd_read_obs - + use MAPL_BaseMod, ONLY: & MAPL_UNDEF @@ -23,17 +23,17 @@ module clsm_ensupd_read_obs MAPL_TICE use io_hdf5, ONLY: & - hdf5read + hdf5read use EASE_conv, ONLY: & ease_convert, & ease_extent - + use LDAS_ensdrv_globals, ONLY: & logit, & logunit, & nodata_tolfrac_generic - + use clsm_ensupd_glob_param, ONLY: & unitnum_obslog @@ -58,10 +58,10 @@ module clsm_ensupd_read_obs f2l_real, & f2l_real8, & f2l_logical - + use LDAS_TilecoordRoutines, ONLY: & get_tile_num_from_latlon - + use LDAS_ensdrv_mpi, ONLY: & root_proc, & numprocs, & @@ -77,26 +77,31 @@ module clsm_ensupd_read_obs use clsm_ensupd_upd_routines, ONLY: & dist_km2deg - + implicit none - include 'mpif.h' - + include 'mpif.h' + private - - public :: collect_obs + + public :: collect_obs contains - - ! ***************************************************************** + ! ***************************************************************** + +! NOTE: This code uses HDF4. HDF4 currently does not support +! building Fortran interfaces due to bugs. For now we +! comment out the code that uses HDF4 but keep it in +! the code base for future reference. +#if 0 subroutine read_ae_l2_sm_hdf( & N_files, fnames, N_data, lon, lat, ae_l2_sm, ease_col, ease_row ) - + ! read soil moisture data from one or more AMSR-E Land hdf files ! ! return ONLY valid data points (ie. excluding no-data-values) - ! that also pass initial QC (based on "Surface Type" QC flag and on + ! that also pass initial QC (based on "Surface Type" QC flag and on ! Heterogeneity_Index ) ! ! reichle, 20 Sep 2005 - added "Surface Type" QC flag @@ -106,21 +111,21 @@ subroutine read_ae_l2_sm_hdf( & ! reichle, 10 Jan 2011 - revised "Surface Type" QC flag implicit none - + integer, intent(in) :: N_files - + character(*), dimension(N_files), intent(in) :: fnames - + integer, intent(out) :: N_data - + real, dimension(:), pointer :: lon, lat, ae_l2_sm ! output - + integer, dimension(:), pointer, optional :: ease_col, ease_row ! output - + ! local parameters - + integer, parameter :: N_fields = 7 - + character(19), dimension(N_fields), parameter :: field_names = (/ & 'Longitude ', & ! 1 'Latitude ', & ! 2 @@ -129,22 +134,22 @@ subroutine read_ae_l2_sm_hdf( & 'Row_Index ', & ! 5 'Column_Index ', & ! 6 'Heterogeneity_Index' /) ! 7 - + real, parameter :: scale_fac_Soil_Moisture = 1000.0; - + integer, parameter :: nodata = -9999 ! NOTE: integer - + ! Initial QC: ! - ! "Surface_Type": Only "low vegetation" or "moderate vegetation" (and no precip, - ! frozen ground, etc) passes. + ! "Surface_Type": Only "low vegetation" or "moderate vegetation" (and no precip, + ! frozen ground, etc) passes. ! - ! "Heterogeneity_Index": Discard all pixels with - ! Heterogeneity_Index>max_Heterogeneity_Index + ! "Heterogeneity_Index": Discard all pixels with + ! Heterogeneity_Index>max_Heterogeneity_Index ! because they are likely mixed land/water pixels. See matlab code ! "detect_coast_in_AMSRE_Land.m" ! in land01:/home/reichle/NSIPP/AMSR/AMSR_E_Land/matlab/ - ! In this subroutine, Heterogeneity_Index is used on its raw form and + ! In this subroutine, Heterogeneity_Index is used on its raw form and ! never scaled into units of Kelvin. ! ! Further info: @@ -152,148 +157,148 @@ subroutine read_ae_l2_sm_hdf( & ! http://nsidc.org/data/amsre/versions.html ! -reichle, 20 Sep 2005 ! -reichle, 8 Feb 2006 - + integer, parameter :: max_Heterogeneity_Index = 500 ! = 5 Kelvin - - ! declarations of hdf functions - + + ! declarations of hdf functions + integer :: hopen, vfstart, vsfatch, vsqfnelt, vsfseek, vsfsfld, vsfread integer :: vsfdtch, vfend, hclose - - + + ! declarations of hdf-related parameters and variables - - integer, dimension(N_files) :: file_id, vdata_id - + + integer, dimension(N_files) :: file_id, vdata_id + integer :: status, n_read, n_records, record_pos - + integer, parameter :: num_dds_block = 0 ! only important for writing hdf - + integer, parameter :: vdata_ref = 7 ! works for AMSRE_L2_Land - + integer, parameter :: DFACC_READ = 1 ! from hdf.inc integer, parameter :: FULL_INTERLACE = 0 ! from hdf.inc - + ! local variables logical :: must_stop - + integer, dimension(N_files) :: N_data_tmp - + integer :: i, j, k, k_off - + integer, dimension(:), allocatable :: surface_type_qc_flag integer, dimension(:), allocatable :: Heterogeneity_Index - + integer*2, dimension(:), allocatable :: tmpint2vec real, dimension(:), allocatable :: tmprealvec - + character(len=*), parameter :: Iam = 'read_ae_l2_sm_hdf' character(len=400) :: err_msg ! ------------------------------------------------------------- - + ! determine number of data to be read from each file - + do j=1,N_files - - ! open and "start" hdf file - file_id(j) = hopen( fnames(j), DFACC_READ, num_dds_block ) - + ! open and "start" hdf file + + file_id(j) = hopen( fnames(j), DFACC_READ, num_dds_block ) + status = vfstart(file_id(j)) - + ! select vdata block that contains fields of interest - - vdata_id(j) = vsfatch(file_id(j), vdata_ref, 'r') - + + vdata_id(j) = vsfatch(file_id(j), vdata_ref, 'r') + ! determine number of records in vdata - + status = vsqfnelt(vdata_id(j), n_records) - + N_data_tmp(j) = n_records - + end do - + ! allocate pointers (must be deallocated outside this subroutine) - + must_stop = .false. - + if ( associated(lon) .or. associated(lat) .or. associated(ae_l2_sm) ) then must_stop = .true. end if - - if ( present(ease_col) ) then + + if ( present(ease_col) ) then if (associated(ease_col)) must_stop = .true. end if - + if ( present(ease_row) ) then if (associated(ease_row)) must_stop = .true. end if - - if (must_stop) then + + if (must_stop) then err_msg = 'output pointers must not be associated/allocated on input.' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + N_data = sum(N_data_tmp) - + allocate(lon(N_data)) allocate(lat(N_data)) if (present(ease_col)) allocate(ease_col(N_data)) if (present(ease_row)) allocate(ease_row(N_data)) allocate(ae_l2_sm(N_data)) - + allocate(surface_type_qc_flag(N_data)) allocate(Heterogeneity_Index(N_data)) ! read hdf data into arrays, concatenate data from N_files files - + k_off = 0 - + do j=1,N_files - + allocate(tmprealvec(N_data_tmp(j))) allocate(tmpint2vec(N_data_tmp(j))) - + do i=1,N_fields - + ! go to start of record (zero-based count) - + record_pos = vsfseek(vdata_id(j), 0) - + ! pick the field to be read - + status = vsfsfld(vdata_id(j), field_names(i)) - + ! read data - + select case (i) - + case (1) - + n_read = vsfread( vdata_id(j), tmprealvec, & N_data_tmp(j), FULL_INTERLACE) - + lon(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - + case (2) - + n_read = vsfread( vdata_id(j), tmprealvec, & N_data_tmp(j), FULL_INTERLACE) - + lat(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - + case (3) - + n_read = vsfread( vdata_id(j) ,tmpint2vec, & N_data_tmp(j), FULL_INTERLACE) surface_type_qc_flag(k_off+1:k_off+N_data_tmp(j)) = tmpint2vec - + ! overwrite surface_type_qc_flag with common pass/fail ! (negative number -1 means fail) - + do k=1,N_data_tmp(j) ! keep only data with @@ -302,19 +307,19 @@ subroutine read_ae_l2_sm_hdf( & ! surface_type_qc_flag=256 ("low veg", and no other bits set) ! ! http://nsidc.org/data/docs/daac/ae_land_l2b_soil_moisture.gd.html - + if (.not. ( & (surface_type_qc_flag(k+k_off)==128) .or. & (surface_type_qc_flag(k+k_off)==256) ) ) & surface_type_qc_flag(k+k_off) = -1 - + end do - + case (4) - + n_read = vsfread(vdata_id(j), tmpint2vec, & N_data_tmp(j), FULL_INTERLACE) - + do k=1,N_data_tmp(j) if (tmpint2vec(k)/=nodata) then ae_l2_sm(k+k_off) = real(tmpint2vec(k)) / & @@ -323,164 +328,164 @@ subroutine read_ae_l2_sm_hdf( & ae_l2_sm(k+k_off) = real(tmpint2vec(k)) end if end do - + case (5) - + if (present(ease_row)) then - + n_read = vsfread( vdata_id(j) ,tmpint2vec, & N_data_tmp(j), FULL_INTERLACE) - + ease_row(k_off+1:k_off+N_data_tmp(j)) = tmpint2vec - + end if - + case (6) - + if (present(ease_col)) then - + n_read = vsfread( vdata_id(j) ,tmpint2vec, & N_data_tmp(j), FULL_INTERLACE) - + ease_col(k_off+1:k_off+N_data_tmp(j)) = tmpint2vec - + end if - + case (7) - + n_read = vsfread( vdata_id(j) ,tmpint2vec, & N_data_tmp(j), FULL_INTERLACE) - + Heterogeneity_Index(k_off+1:k_off+N_data_tmp(j)) = tmpint2vec - + case default - + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown case') end select - + if (n_read/=N_data_tmp(j)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'ERROR reading hdf') end if - + end do - + ! clean up deallocate(tmprealvec) deallocate(tmpint2vec) - + ! close hdf files - + status = vsfdtch(vdata_id(j)) - + status = vfend(file_id(j)) - + status = hclose(file_id(j)) - + ! prepare next j - + k_off = k_off + N_data_tmp(j) - + end do - + ! ------------------------------------- ! ! eliminate no-data-values and data that fail initial QC - + j = 0 - + do i=1,N_data - + if ( (ae_l2_sm(i)>0.) .and. & ! any neg is nodata (surface_type_qc_flag(i)>=0) .and. & ! any neg is fail (Heterogeneity_Index(i)<=max_Heterogeneity_Index) & - ) then - + ) then + j=j+1 - + ae_l2_sm(j) = ae_l2_sm(i) lon(j) = lon(i) lat(j) = lat(i) if (present(ease_col)) ease_col(j) = ease_col(i) if (present(ease_row)) ease_row(j) = ease_row(i) - + end if - + end do - + N_data = j - + deallocate(surface_type_qc_flag) deallocate(Heterogeneity_Index) end subroutine read_ae_l2_sm_hdf - + ! ***************************************************************** - + subroutine read_obs_ae_l2_sm( & work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, ae_l2_sm, std_ae_l2_sm ) - + ! Read observations of surface soil moisture from AMSR-E L2 files - ! Set flag "found_obs" to true if observations are available + ! Set flag "found_obs" to true if observations are available ! for assimilation. ! ! If there are N > 1 observations in a given tile, ! a "super-observation" is computed by averaging the N observations ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! reichle, 25 Jul 2005 ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + character(*), intent(in) :: work_path character(*), intent(in) :: exp_id type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim, N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: ae_l2_sm real, intent(out), dimension(N_catd) :: std_ae_l2_sm logical, intent(out) :: found_obs - + ! --------------- - + ! locals - + ! AMSR_E_L2_Land files are available for approximately 50min periods ! covering one ascending or descending swath. The filename indicates ! the start time of the swath ! "ae_time_offset" is used to find the mean time of the the interval ! which is approximately the time of the equator overpass. ! This time is assigned to all observations of the swath. - + integer, parameter :: ae_time_offset = 1500 ! 25 minutes in seconds - + character(4) :: DDHH character(6) :: YYYYMM character(8) :: date_string @@ -489,110 +494,110 @@ subroutine read_obs_ae_l2_sm( & character(400) :: cmd type(date_time_type) :: date_time_tmp - + integer :: i, ind, N_tmp, N_files character(300), dimension(:), allocatable :: fnames - + real, dimension(:), pointer :: tmp_obs, tmp_lat, tmp_lon integer, dimension(:), pointer :: tmp_tile_num - - integer, dimension(N_catd) :: N_obs_in_tile + + integer, dimension(N_catd) :: N_obs_in_tile character(len=*), parameter :: Iam = 'read_obs_ae_l2_sm' ! ------------------------------------------------------------------- - - nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num ) - + + nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num ) + ! --------------- - + ! initialize - + found_obs = .false. - ! find files that are within half-open interval + ! find files that are within half-open interval ! [date_time-dtstep_assim/2,date_time+dtstep_assim/2) - - date_time_tmp = date_time - + + date_time_tmp = date_time + call augment_date_time( -(dtstep_assim/2 + ae_time_offset), date_time_tmp ) - + ! get tmp file name and remove file if it exists - + call date_and_time(date_string, time_string) ! f90 intrinsic function - + tmpfname = trim(work_path) // '/' // 'tmp.' // trim(exp_id) & // '.' // date_string // time_string - cmd = '/bin/rm -f ' // tmpfname - + cmd = '/bin/rm -f ' // tmpfname + call Execute_command_line(trim(cmd)) - + ! identify all files within current assimilation interval ! (list all files within hourly intervals) - + do i=1,(dtstep_assim/3600) - + write (YYYYMM,'(i6.6)') date_time_tmp%year*100 + date_time_tmp%month write (DDHH, '(i4.4)') date_time_tmp%day *100 + date_time_tmp%hour - + cmd = 'ls ' // trim(this_obs_param%path) // '/' // YYYYMM(1:4) // & '/M' // YYYYMM(5:6) // '/' // trim(this_obs_param%name) & // '*' // YYYYMM // DDHH // '*' - + if (trim(this_obs_param%descr)=='ae_l2_sm_a') then - + cmd = trim(cmd) // '_A.hdf' - + elseif (trim(this_obs_param%descr)=='ae_l2_sm_d') then - + cmd = trim(cmd) // '_D.hdf' - + else - + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown descr') end if - + cmd = trim(cmd) // ' >> ' // trim(tmpfname) - + call Execute_command_line(trim(cmd)) - + call augment_date_time( 3600, date_time_tmp ) - + end do - + ! find out how many need to be read - + tmpfname2 = trim(tmpfname) // '.wc' - + cmd = 'wc -w ' // trim(tmpfname) // ' > ' // trim(tmpfname2) - + call Execute_command_line(trim(cmd)) - + open(10, file=tmpfname2, form='formatted', action='read') - + read(10,*) N_files - + close(10,status='delete') - + ! load file names into "fnames" - + open(10, file=tmpfname, form='formatted', action='read') - + if (N_files>0) then - + allocate(fnames(N_files)) - + do i=1,N_files read(10,'(a)') fnames(i) end do - + end if - + close(10,status='delete') - + ! read observations: ! ! 1.) read N_tmp observations and their lat/lon info from file @@ -610,30 +615,30 @@ subroutine read_obs_ae_l2_sm( & call read_ae_l2_sm_hdf( & N_files, fnames, & N_tmp, tmp_lon, tmp_lat, tmp_obs ) - + if (logit) then - + write (logunit,*) 'read_obs_ae_l2_sm: read ', N_tmp, & ' at date_time = ', date_time, ' from ' do i=1,N_files write (logunit,*) trim(fnames(i)) end do write (logunit,*) '----------' - + end if deallocate(fnames) else - + N_tmp = 0 - + end if ! ------------------------------------------------------------------ ! note QC and no-data-value block in subroutine read_ae_l2_sm_hdf() - + ! ------------------------------------------------------------------ ! ! 2.) for each observation @@ -641,237 +646,239 @@ subroutine read_obs_ae_l2_sm( & ! b) determine tile within grid cell that contains lat/lon if (N_tmp>0) then - + allocate(tmp_tile_num(N_tmp)) - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_tmp, tmp_lat, tmp_lon, & this_obs_param, & tmp_tile_num ) - - + + ! ---------------------------------------------------------------- ! ! 3.) compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + ae_l2_sm = 0. N_obs_in_tile = 0 - + do i=1,N_tmp - + ind = tmp_tile_num(i) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this step eliminates obs outside domain - + ae_l2_sm(ind) = ae_l2_sm(ind) + tmp_obs(i) - + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 - + end if - + end do - + ! normalize - + do i=1,N_catd - + if (N_obs_in_tile(i)>1) then - + ae_l2_sm(i) = ae_l2_sm(i)/real(N_obs_in_tile(i)) - + elseif (N_obs_in_tile(i)==0) then - + ae_l2_sm(i) = this_obs_param%nodata - + end if - + end do - + ! clean up - + if (associated(tmp_tile_num)) deallocate(tmp_tile_num) - + ! -------------------------------- - + ! set observation error standard deviation - + do i=1,N_catd std_ae_l2_sm(i) = this_obs_param%errstd end do - + ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if - + ! clean up - + if (associated(tmp_obs)) deallocate(tmp_obs) if (associated(tmp_lon)) deallocate(tmp_lon) if (associated(tmp_lat)) deallocate(tmp_lat) end subroutine read_obs_ae_l2_sm - + +#endif + ! *************************************************************************** subroutine read_ae_sm_LPRM_bin( & this_obs_param, N_files, fnames, N_data, lon, lat, ae_sm_LPRM, ease_col, ease_row ) - + ! read soil moisture data from one or more AMSR-E LPRM bin files ! ! return ONLY valid data points (ie. excluding no-data-values) - ! + ! ! no QC in addition to what was done in matlab-preprocessing ! ! reichle, 20 Feb 2009 - + implicit none - + type(obs_param_type), intent(in) :: this_obs_param - + integer, intent(in) :: N_files - + character(*), dimension(N_files), intent(in) :: fnames - + integer, intent(out) :: N_data - + real, dimension(:), pointer :: lon, lat, ae_sm_LPRM ! output - + integer, dimension(:), pointer, optional :: ease_col, ease_row ! output - + ! local variables - + logical :: must_stop - + integer, dimension(N_files) :: N_data_tmp - + integer :: i, j, k_off - + integer, dimension(:), allocatable :: tmpintvec real, dimension(:), allocatable :: tmprealvec character(len=*), parameter :: Iam = 'read_ae_sm_LPRM_bin' character(len=400) :: err_msg - + ! ------------------------------------------------------------- - + ! make sure pointers are not allocated or associated - + must_stop = .false. - + if ( associated(lon) .or. associated(lat) .or. associated(ae_sm_LPRM) ) then must_stop = .true. end if - - if ( present(ease_col) ) then + + if ( present(ease_col) ) then if (associated(ease_col)) must_stop = .true. end if - + if ( present(ease_row) ) then if (associated(ease_row)) must_stop = .true. end if - + if (must_stop) then err_msg = 'output pointers must not be associated/allocated on input' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if ! determine number of data to be read from each file - + N_data = 0 - + do j=1,N_files - - ! open file - + + ! open file + open( 10, file=trim(fnames(j)), form='unformatted',convert='big_endian', action='read' ) - + read( 10) N_data_tmp(j) - + close(10, status='keep') - + end do - + ! allocate pointers (must be deallocated outside this subroutine!) - + N_data = sum(N_data_tmp) - + allocate(lon(N_data)) allocate(lat(N_data)) if (present(ease_col)) allocate(ease_col(N_data)) if (present(ease_row)) allocate(ease_row(N_data)) allocate(ae_sm_LPRM(N_data)) - + ! read data into arrays, concatenate data from N_files files - + ! format of AMSR_sm_LPRM_EASE_bin files: ! ! record 1 -- N_data int*4 - ! record 2 -- lon( 1:N_data) real*4 - ! record 3 -- lat( 1:N_data) real*4 - ! record 4 -- sm_C( 1:N_data) real*4 - ! record 5 -- sm_X( 1:N_data) real*4 - ! record 6 -- od_C( 1:N_data) real*4 - ! record 7 -- od_X( 1:N_data) real*4 - ! record 8 -- res_C(1:N_data) real*4 - ! record 9 -- res_X(1:N_data) real*4 - ! record 10 -- ts_C( 1:N_data) real*4 - ! record 11 -- ts_X( 1:N_data) real*4 - ! record 12 -- rfi_C(1:N_data) int*4 - ! record 13 -- rfi_X(1:N_data) int*4 + ! record 2 -- lon( 1:N_data) real*4 + ! record 3 -- lat( 1:N_data) real*4 + ! record 4 -- sm_C( 1:N_data) real*4 + ! record 5 -- sm_X( 1:N_data) real*4 + ! record 6 -- od_C( 1:N_data) real*4 + ! record 7 -- od_X( 1:N_data) real*4 + ! record 8 -- res_C(1:N_data) real*4 + ! record 9 -- res_X(1:N_data) real*4 + ! record 10 -- ts_C( 1:N_data) real*4 + ! record 11 -- ts_X( 1:N_data) real*4 + ! record 12 -- rfi_C(1:N_data) int*4 + ! record 13 -- rfi_X(1:N_data) int*4 ! record 14 -- ind_i(1:N_data) int*4 zero-based (!) EASE row index ! record 15 -- ind_j(1:N_data) int*4 zero-based (!) EASE col index ! record 16 -- time( 1:N_data) real*4 minutes since beginning of half-orbit - + k_off = 0 - + do j=1,N_files - + allocate(tmprealvec(N_data_tmp(j))) - + if (present(ease_col)) allocate(tmpintvec(N_data_tmp(j))) - + open (10, file=trim(fnames(j)), form='unformatted',convert='big_endian', action='read' ) - + ! re-read N_data - - read (10) N_data_tmp(j) - - ! read data as needed - - read (10) tmprealvec; lon(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - read (10) tmprealvec; lat(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - + + read (10) N_data_tmp(j) + + ! read data as needed + + read (10) tmprealvec; lon(k_off+1:k_off+N_data_tmp(j)) = tmprealvec + read (10) tmprealvec; lat(k_off+1:k_off+N_data_tmp(j)) = tmprealvec + if (this_obs_param%descr(13:14)=='_C') then - + read (10) tmprealvec; ae_sm_LPRM(k_off+1:k_off+N_data_tmp(j)) = tmprealvec read (10) ! skip sm_X - + else if (this_obs_param%descr(13:14)=='_X') then - + read (10) ! skip sm_C read (10) tmprealvec; ae_sm_LPRM(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - - else - + + else + err_msg = 'unknown descr, ' // this_obs_param%descr call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + if (present(ease_col) .and. present(ease_row)) then - + read (10) ! skip od_C read (10) ! skip od_X read (10) ! skip res_C @@ -880,119 +887,119 @@ subroutine read_ae_sm_LPRM_bin( & read (10) ! skip ts_X read (10) ! skip rfi_C read (10) ! skip rfi_X - + read (10) tmpintvec; ease_col(k_off+1:k_off+N_data_tmp(j)) = tmpintvec read (10) tmpintvec; ease_row(k_off+1:k_off+N_data_tmp(j)) = tmpintvec - + ! read (10) ! skip time - + end if - + ! clean up - + close(10, status='keep') - + deallocate(tmprealvec) if (allocated(tmpintvec)) deallocate(tmpintvec) - + ! prepare next j - + k_off = k_off + N_data_tmp(j) - + end do - + ! ------------------------------------- ! - ! eliminate no-data-values - + ! eliminate no-data-values + j = 0 - + do i=1,N_data - + if (ae_sm_LPRM(i)>0.) then ! any neg is nodata - + j=j+1 - + ae_sm_LPRM(j) = ae_sm_LPRM(i) lon(j) = lon(i) lat(j) = lat(i) if (present(ease_col)) ease_col(j) = ease_col(i) if (present(ease_row)) ease_row(j) = ease_row(i) - + end if - + end do - + N_data = j - + end subroutine read_ae_sm_LPRM_bin - + ! ***************************************************************** - + subroutine read_obs_ae_sm_LPRM( & work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, ae_sm_LPRM, std_ae_sm_LPRM ) - + ! Read observations of surface soil moisture from AMSR-E sm LPRM files ! (Richard de Jeu, Vrije Universiteit Amsterdam) - ! Set flag "found_obs" to true if observations are available + ! Set flag "found_obs" to true if observations are available ! for assimilation. ! ! If there are N > 1 observations in a given tile, ! a "super-observation" is computed by averaging the N observations ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! reichle, 20 Feb 2009 ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + character(*), intent(in) :: work_path character(*), intent(in) :: exp_id type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim, N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: ae_sm_LPRM real, intent(out), dimension(N_catd) :: std_ae_sm_LPRM logical, intent(out) :: found_obs - + ! --------------- - + ! locals - + ! AMSR_E_L2_Land files are available for approximately 50min periods ! covering one ascending or descending swath. The filename indicates ! the start time of the swath ! "ae_time_offset" is used to find the mean time of the the interval ! which is approximately the time of the equator overpass. ! This time is assigned to all observations of the swath. - + integer, parameter :: ae_time_offset = 1500 ! 25 minutes in seconds - + character(4) :: DDHH character(6) :: YYYYMM character(8) :: date_string @@ -1001,113 +1008,113 @@ subroutine read_obs_ae_sm_LPRM( & character(400) :: cmd type(date_time_type) :: date_time_tmp - + integer :: i, ind, N_tmp, N_files character(300), dimension(:), allocatable :: fnames - + real, dimension(:), pointer :: tmp_obs, tmp_lat, tmp_lon integer, dimension(:), pointer :: tmp_tile_num - - integer, dimension(N_catd) :: N_obs_in_tile + + integer, dimension(N_catd) :: N_obs_in_tile character(len=*), parameter :: Iam = 'read_obs_ae_sm_LPRM' ! ------------------------------------------------------------------- - - nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num ) - + + nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num ) + ! --------------- - + ! initialize - + found_obs = .false. - - ! find files that are within half-open interval + + ! find files that are within half-open interval ! [date_time-dtstep_assim/2,date_time+dtstep_assim/2) - - date_time_tmp = date_time - + + date_time_tmp = date_time + call augment_date_time( -(dtstep_assim/2 + ae_time_offset), date_time_tmp ) - + ! get tmp file name and remove file if it exists - + call date_and_time(date_string, time_string) ! f90 intrinsic function - + tmpfname = trim(work_path) // '/' // 'tmp.' // trim(exp_id) & // '.' // date_string // time_string - cmd = '/bin/rm -f ' // tmpfname - + cmd = '/bin/rm -f ' // tmpfname + call Execute_command_line(trim(cmd)) - + ! identify all files within current assimilation interval ! (list all files within hourly intervals) - + do i=1,(dtstep_assim/3600) - + write (YYYYMM,'(i6.6)') date_time_tmp%year*100 + date_time_tmp%month write (DDHH, '(i4.4)') date_time_tmp%day *100 + date_time_tmp%hour - + cmd = 'ls ' // trim(this_obs_param%path) // '/' // YYYYMM(1:4) // & '.' // YYYYMM(5:6) // '/' // trim(this_obs_param%name) & - // YYYYMM // DDHH(1:2) // '.' // DDHH(3:4) // '??' - + // YYYYMM // DDHH(1:2) // '.' // DDHH(3:4) // '??' + if (this_obs_param%descr(1:12)=='ae_sm_LPRM_a') then - + cmd = trim(cmd) // '_A.bin' - + elseif (this_obs_param%descr(1:12)=='ae_sm_LPRM_d') then - + cmd = trim(cmd) // '_D.bin' - + else - + !write (logunit,*) 'read_obs_ae_sm_LPRM(): unknown descr, STOPPING.' !stop - + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown descr') - + end if - + cmd = trim(cmd) // ' >> ' // trim(tmpfname) - + call Execute_command_line(trim(cmd)) - + call augment_date_time( 3600, date_time_tmp ) - + end do - + ! find out how many need to be read - + tmpfname2 = trim(tmpfname) // '.wc' - + cmd = 'wc -w ' // trim(tmpfname) // ' > ' // trim(tmpfname2) - + call Execute_command_line(trim(cmd)) - + open(10, file=tmpfname2, form='formatted', action='read') - + read(10,*) N_files - + close(10,status='delete') - + ! load file names into "fnames" - + open(10, file=tmpfname, form='formatted', action='read') - + if (N_files>0) then - + allocate(fnames(N_files)) - + do i=1,N_files read(10,'(a)') fnames(i) end do - + end if - + close(10,status='delete') - + ! read observations: ! ! 1.) read N_tmp observations and their lat/lon info from file @@ -1121,13 +1128,13 @@ subroutine read_obs_ae_sm_LPRM( & ! 1.) read N_tmp observations and their lat/lon info from file if (N_files>0) then - + call read_ae_sm_LPRM_bin( & this_obs_param, N_files, fnames, & N_tmp, tmp_lon, tmp_lat, tmp_obs ) - + if (logit) then - + write (logunit,*) 'read_obs_ae_sm_LPRM: read ', N_tmp, & ' at date_time = ', date_time, ' from ' do i=1,N_files @@ -1138,17 +1145,17 @@ subroutine read_obs_ae_sm_LPRM( & end if deallocate(fnames) - + else - + N_tmp = 0 - + end if - + ! ------------------------------------------------------------------ ! QC is done in matlab pre-processing to bin files - + ! ------------------------------------------------------------------ ! ! 2.) for each observation @@ -1156,82 +1163,82 @@ subroutine read_obs_ae_sm_LPRM( & ! b) determine tile within grid cell that contains lat/lon if (N_tmp>0) then - + allocate(tmp_tile_num(N_tmp)) - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_tmp, tmp_lat, tmp_lon, & this_obs_param, & tmp_tile_num ) - + ! ---------------------------------------------------------------- ! ! 3.) compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + ae_sm_LPRM = 0. N_obs_in_tile = 0 - + do i=1,N_tmp - + ind = tmp_tile_num(i) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this step eliminates obs outside domain - + ae_sm_LPRM(ind) = ae_sm_LPRM(ind) + tmp_obs(i) - + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 - + end if - + end do - + ! normalize - + do i=1,N_catd - + if (N_obs_in_tile(i)>1) then - + ae_sm_LPRM(i) = ae_sm_LPRM(i)/real(N_obs_in_tile(i)) - + elseif (N_obs_in_tile(i)==0) then - + ae_sm_LPRM(i) = this_obs_param%nodata - + end if - + end do - + ! clean up - + if (associated(tmp_tile_num)) deallocate(tmp_tile_num) - + ! -------------------------------- - + ! set observation error standard deviation - + do i=1,N_catd std_ae_sm_LPRM(i) = this_obs_param%errstd end do - + ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if - + ! clean up - + if (associated(tmp_obs)) deallocate(tmp_obs) if (associated(tmp_lon)) deallocate(tmp_lon) if (associated(tmp_lat)) deallocate(tmp_lat) @@ -1246,61 +1253,61 @@ subroutine read_obs_sm_ASCAT( & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, sm_ASCAT, std_sm_ASCAT ) - + !--------------------------------------------------------------------- - ! - ! Routine to read in ASCAT surface degree of saturation obs. - ! Output is found_obs, sm_ASCAT, std_sm_ASCAT + ! + ! Routine to read in ASCAT surface degree of saturation obs. + ! Output is found_obs, sm_ASCAT, std_sm_ASCAT ! ! Reads in obs provided by Wolfgang Wagner (TUW), converted to - ! once hourly binary files (and projected onto EASE grid). + ! once hourly binary files (and projected onto EASE grid). ! Updating to the EUMETSAT BUFR (DGG) files will require a new ! reader, and changes to file-name / time-stamping ! - ! Draper, May 2011. + ! Draper, May 2011. ! Based on read_obs_ae_sm_LPRM ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + character(*), intent(in) :: work_path character(*), intent(in) :: exp_id type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim, N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: sm_ASCAT real, intent(out), dimension(N_catd) :: std_sm_ASCAT logical, intent(out) :: found_obs - + ! --------------- - - ! For the files generated from the TUW timseries, each + + ! For the files generated from the TUW timseries, each ! file is time-stamped with the time of the observations, rounded - ! down to the nearest obs + ! down to the nearest obs ! E(minutes past hour)=30 -> use this as the offset ! Will need to be updated if using EUMETSAT BUFR files integer, parameter :: ae_time_offset = 1800 ! 30 minutes in seconds - + character(4) :: DDHH character(6) :: YYYYMM character(8) :: date_string @@ -1309,114 +1316,114 @@ subroutine read_obs_sm_ASCAT( & character(400) :: cmd type(date_time_type) :: date_time_tmp - + integer :: i, ind, N_tmp, N_files character(300), dimension(:), allocatable :: fnames - + real, dimension(:), pointer :: tmp_obs, tmp_lat, tmp_lon integer, dimension(:), pointer :: tmp_tile_num - - integer, dimension(N_catd) :: N_obs_in_tile + + integer, dimension(N_catd) :: N_obs_in_tile real, parameter :: tol = 1e-2 character(len=*), parameter :: Iam = 'read_obs_sm_ASCAT' ! ------------------------------------------------------------------- - - nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num ) - + + nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num ) + ! --------------- - + ! initialize - + found_obs = .false. - ! find files that are within half-open interval + ! find files that are within half-open interval ! [date_time-dtstep_assim/2,date_time+dtstep_assim/2) - - date_time_tmp = date_time - + + date_time_tmp = date_time + call augment_date_time( -(dtstep_assim/2 + ae_time_offset), date_time_tmp ) - + ! get tmp file name and remove file if it exists - + call date_and_time(date_string, time_string) ! f90 intrinsic function - + tmpfname = trim(work_path) // '/' // 'tmp.' // trim(exp_id) & // '.' // date_string // time_string - cmd = '/bin/rm -f ' // tmpfname - + cmd = '/bin/rm -f ' // tmpfname + call Execute_command_line(trim(cmd)) - + ! identify all files within current assimilation interval ! (list all files within hourly intervals) - + do i=1,(dtstep_assim/3600) - + write (YYYYMM,'(i6.6)') date_time_tmp%year*100 + date_time_tmp%month write (DDHH, '(i4.4)') date_time_tmp%day *100 + date_time_tmp%hour - + ! TUW files time stamped with hour only. Update for EUMETSAT BUFR cmd = 'ls ' // trim(this_obs_param%path) // '/' // YYYYMM(1:4) // & '.' // YYYYMM(5:6) // '/' // trim(this_obs_param%name) & - // YYYYMM // DDHH(1:2) // '.' // DDHH(3:4) !// '??' - + // YYYYMM // DDHH(1:2) // '.' // DDHH(3:4) !// '??' + if (this_obs_param%descr(1:10)=='ASCAT_SM_A') then - + cmd = trim(cmd) // '_A.bin' - + elseif (this_obs_param%descr(1:10)=='ASCAT_SM_D') then - + cmd = trim(cmd) // '_D.bin' - + else - + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown descr') - + end if cmd = trim(cmd) // ' >> ' // trim(tmpfname) - + call Execute_command_line(trim(cmd)) - - + + call augment_date_time( 3600, date_time_tmp ) - + end do - + ! find out how many need to be read - + tmpfname2 = trim(tmpfname) // '.wc' - + cmd = 'wc -w ' // trim(tmpfname) // ' > ' // trim(tmpfname2) - + call Execute_command_line(trim(cmd)) - + open(10, file=tmpfname2, form='formatted', action='read') - + read(10,*) N_files - + close(10,status='delete') - + ! load file names into "fnames" - + open(10, file=tmpfname, form='formatted', action='read') - + if (N_files>0) then - + allocate(fnames(N_files)) - + do i=1,N_files read(10,'(a)') fnames(i) end do - + end if - + close(10,status='delete') - + ! read observations: ! ! 1.) read N_tmp observations and their lat/lon info from file @@ -1430,13 +1437,13 @@ subroutine read_obs_sm_ASCAT( & ! 1.) read N_tmp observations and their lat/lon info from file if (N_files>0) then - + call read_sm_ASCAT_bin( & N_files, fnames, & N_tmp, tmp_lon, tmp_lat, tmp_obs ) - + if (logit) then - + write (logunit,*) 'read_obs_sm_ASCAT: read ', N_tmp, & ' at date_time = ', date_time, ' from ' do i=1,N_files @@ -1447,13 +1454,13 @@ subroutine read_obs_sm_ASCAT( & end if deallocate(fnames) - + else - + N_tmp = 0 - + end if - + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! SOME QC SHOULD BE DONE HERE!!! @@ -1461,7 +1468,7 @@ subroutine read_obs_sm_ASCAT( & ! MAKE SURE no-data-values ARE DEALT WITH ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - + ! ---------------------------------------------------------------- ! ! 2.) for each observation @@ -1469,87 +1476,87 @@ subroutine read_obs_sm_ASCAT( & ! b) determine tile within grid cell that contains lat/lon if (N_tmp>0) then - + allocate(tmp_tile_num(N_tmp)) - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_tmp, tmp_lat, tmp_lon, & this_obs_param, & tmp_tile_num ) - - + + ! ---------------------------------------------------------------- ! ! 3.) compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + sm_ASCAT = 0. N_obs_in_tile = 0 - + do i=1,N_tmp - + ind = tmp_tile_num(i) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this step eliminates obs outside domain - + sm_ASCAT(ind) = sm_ASCAT(ind) + tmp_obs(i) - + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 - + end if - + end do - + ! normalize - + do i=1,N_catd - + if (N_obs_in_tile(i)>1) then - + sm_ASCAT(i) = sm_ASCAT(i)/real(N_obs_in_tile(i)) - + elseif (N_obs_in_tile(i)==0) then - + sm_ASCAT(i) = this_obs_param%nodata - + end if - + end do - + ! clean up - + if (associated(tmp_tile_num)) deallocate(tmp_tile_num) - + ! -------------------------------- - + ! set observation error standard deviation do i=1,N_catd std_sm_ASCAT(i) = this_obs_param%errstd enddo ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if - + ! clean up - + if (associated(tmp_obs)) deallocate(tmp_obs) if (associated(tmp_lon)) deallocate(tmp_lon) if (associated(tmp_lat)) deallocate(tmp_lat) end subroutine read_obs_sm_ASCAT - + ! **************************************************************************** subroutine read_obs_sm_ASCAT_EUMET( & @@ -1557,11 +1564,11 @@ subroutine read_obs_sm_ASCAT_EUMET( & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, ASCAT_sm, ASCAT_sm_std, ASCAT_lon, ASCAT_lat, ASCAT_time ) - + !--------------------------------------------------------------------- - ! - ! Routine to read in ASCAT surface degree of saturation (sfds) obs from - ! BUFR files for both ascending and descending passes. + ! + ! Routine to read in ASCAT surface degree of saturation (sfds) obs from + ! BUFR files for both ascending and descending passes. ! ! ASCAT_sm and ASCAT_sm_std outputs from this subroutine are in wetness fraction (i.e., 0-1) units! ! @@ -1569,11 +1576,11 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! Soil moisture derived from re-sampled (spatially averaged) backscatter (sigma0) values ! on a 25-km orbit swath grid. Input data files are in BUFR file format. ! - ! EUMETSAT BUFR files contain data for both ascending and descending half-orbits. - ! The BUFR field DOMO ("Direction of motion of moving observing platform") could be used to + ! EUMETSAT BUFR files contain data for both ascending and descending half-orbits. + ! The BUFR field DOMO ("Direction of motion of moving observing platform") could be used to ! separate Asc and Desc. (The BUFR files do not contain any explicit orbit indicator variable.) ! According to Pamela Schoebel-Pattiselanno, EUMETSAT User Services Helpdesk: - ! "When the value (of DOMO) is between 180 and 270 degrees, it is the descending part + ! "When the value (of DOMO) is between 180 and 270 degrees, it is the descending part ! of the orbit. When it is between 270 and 360 degrees, it is the ascending part." ! ! Q. Liu, Nov 2019 - based on read_obs_sm_ASCAT @@ -1581,47 +1588,47 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! A. Fox, reichle, Feb 2024 - added ASCAT obs mask ! ! -------------------------------------------------------------------- - - use netcdf + + use netcdf implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim, N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + logical, intent(out) :: found_obs real, intent(out), dimension(N_catd) :: ASCAT_sm ! sfds obs [fraction] (i.e., 0-1) real, intent(out), dimension(N_catd) :: ASCAT_sm_std ! sfds obs err std [fraction] (i.e., 0-1) real, intent(out), dimension(N_catd) :: ASCAT_lon, ASCAT_lat real*8, intent(out), dimension(N_catd) :: ASCAT_time ! J2000 seconds - + ! --------------- - + ! Each obs file contains ~100-110 minutes (1 full orbit?) of observations (dt_ASCAT_obsfile). - ! File name indicates start time of swath. - + ! File name indicates start time of swath. + integer, parameter :: dt_ASCAT_obsfile = 110*60 ! seconds - + integer, parameter :: N_fnames_max = 15 ! max number of files per day character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 - + character( 15) :: str_date_time character( 80) :: fname_of_fname_list character(300) :: tmpfname @@ -1629,11 +1636,11 @@ subroutine read_obs_sm_ASCAT_EUMET( & type(date_time_type) :: date_time_tmp type(date_time_type) :: date_time_low, date_time_low_fname type(date_time_type) :: date_time_up - + integer :: ii, ind, N_tmp, N_files, kk, N_obs, N_fnames, N_fnames_tmp, obs_dir_hier - + character(300), dimension(:), allocatable :: fnames, tmpfnames - + ! -------------------- ! ! variables for BUFR read @@ -1653,14 +1660,14 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! variables for obs mask read integer(kind=1), dimension(:,:), allocatable :: mask_data - + real :: mask_ll_lon, mask_ll_lat, mask_dlon, mask_dlat - - integer :: ncid, ierr + + integer :: ncid, ierr integer :: mask_N_lon, mask_N_lat, mask_lon_ind, mask_lat_ind integer :: lon_dimid, lat_dimid integer :: mask_varid, ll_lon_varid, ll_lat_varid, dlon_varid, dlat_varid - + character(300) :: mask_filename logical :: file_exists @@ -1678,68 +1685,68 @@ subroutine read_obs_sm_ASCAT_EUMET( & integer, dimension(:), pointer :: tmp_tile_num - integer, dimension(N_catd) :: N_obs_in_tile + integer, dimension(N_catd) :: N_obs_in_tile character(len=*), parameter :: Iam = 'read_obs_sm_ASCAT_EUMET' character(len=400) :: err_msg ! ------------------------------------------------------------------- - - nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num, tmp_jtime ) - + + nullify( tmp_obs, tmp_lat, tmp_lon, tmp_tile_num, tmp_jtime ) + ! --------------- - + ! initialize - + found_obs = .false. - - ! find files that are within half-open interval + + ! find files that are within half-open interval ! (date_time-dtstep_assim/2,date_time+dtstep_assim/2] - date_time_low = date_time + date_time_low = date_time call augment_date_time( -(dtstep_assim/2), date_time_low) - date_time_up = date_time + date_time_up = date_time call augment_date_time( (dtstep_assim/2), date_time_up) - - ! calculate "extra" date_time_low to catch files w/ swath start time stamps before window + + ! calculate "extra" date_time_low to catch files w/ swath start time stamps before window ! but containing relevant obs - + date_time_low_fname = date_time_low call augment_date_time( -dt_ASCAT_obsfile, date_time_low_fname) - + ! ---------------------------------------------------------------- ! ! read file with list of ASCAT file names for first day - + fname_of_fname_list = 'dummy' ! make sure it is properly defined in obs_param nml obs_dir_hier = 1 - + call read_obs_fnames( date_time_low_fname, this_obs_param, & fname_of_fname_list, N_fnames_max, & N_fnames, fname_list(1:N_fnames_max), obs_dir_hier ) - + ! if needed, read file with list of ASCAT file names for second day and add ! file names into "fname_list" - + if (date_time_low_fname%day /= date_time_up%day) then - + call read_obs_fnames( date_time_up, this_obs_param, & fname_of_fname_list, N_fnames_max, & N_fnames_tmp, fname_list((N_fnames+1):(N_fnames+N_fnames_max)), obs_dir_hier ) - + N_fnames = N_fnames + N_fnames_tmp - + end if - + tmpfnames = fname_list(1:N_fnames) - + ! ---------------------------------------------------------------- ! ! find files that have obs within assimilation window - - N_tmp = 0 - + + N_tmp = 0 + do kk = 1,N_fnames tmpfname = fname_list(kk) @@ -1748,36 +1755,36 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! ! e.g. Y2019/M07/M01-ASCA-ASCSMO02-NA-5.0-20190702075700.000000000Z-20190702084627-1350204.bfr ! - ! 12345678901234567890123456789012345678901234567890123456789012345678901234567890 + ! 12345678901234567890123456789012345678901234567890123456789012345678901234567890 ! 1 2 3 4 5 6 7 str_date_time = tmpfname(36:49) - + read(str_date_time( 1: 4), *) date_time_tmp%year read(str_date_time( 5: 6), *) date_time_tmp%month read(str_date_time( 7: 8), *) date_time_tmp%day read(str_date_time( 9:10), *) date_time_tmp%hour read(str_date_time(11:12), *) date_time_tmp%min read(str_date_time(13:14), *) date_time_tmp%sec - + if ( datetime_lt_refdatetime( date_time_low_fname, date_time_tmp ) .and. & - datetime_le_refdatetime( date_time_tmp, date_time_up ) ) then + datetime_le_refdatetime( date_time_tmp, date_time_up ) ) then N_tmp = N_tmp + 1 tmpfnames(N_tmp) = trim(this_obs_param%path) // '/' // trim(tmpfname) - + end if - + end do - + fnames = tmpfnames(1:N_tmp) N_files = N_tmp - + ! ---------------------------------------------------------------- ! ! loop through files and read obs + metadata into tmp_data - + if (N_files>0) then ! read and process data if files are found @@ -1786,25 +1793,25 @@ subroutine read_obs_sm_ASCAT_EUMET( & allocate(tmp1_lat( max_rec )) allocate(tmp1_obs( max_rec )) allocate(tmp1_jtime(max_rec )) - + allocate(tmp_data( max_obs,14)) - + N_obs = 0 - + do kk = 1,N_files - + ! open bufr file - + call closbf(lnbufr) ! if a file with unit number lnbufr is open in (or "linked" with) BUFR, close it open(lnbufr, file=trim(fnames(kk)), action='read', form='unformatted') - call openbf(lnbufr, 'SEC3', lnbufr) + call openbf(lnbufr, 'SEC3', lnbufr) call mtinfo( trim(this_obs_param%path) // '/BUFR_mastertable/', lnbufr+1, lnbufr+2) - call datelen(10) ! select date/time format with 4-digit year (YYYYMMDDHH) - + call datelen(10) ! select date/time format with 4-digit year (YYYYMMDDHH) + msg_report: do while( ireadmg(lnbufr,subset,idate) == 0 ) - + loop_report: do while( ireadsb(lnbufr) == 0 ) - + ! columns of tmp_data: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 call ufbint(lnbufr,tmp_vdata,14,1,iret,'YEAR MNTH DAYS HOUR MINU SECO SSOM SMPF SMCF ALFR TPCX IWFR CLATH CLONH') @@ -1817,13 +1824,13 @@ subroutine read_obs_sm_ASCAT_EUMET( & end if tmp_data(N_obs,:) = tmp_vdata - + end do loop_report end do msg_report - + call closbf(lnbufr) close(lnbufr) - + end do ! end file loop ! ---------------------------------------------------------------- @@ -1831,13 +1838,13 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! read mask file for ASCAT obs (netcdf format, regular lat/lon grid) mask_filename = trim(this_obs_param%maskpath) // '/' // trim(this_obs_param%maskname) - + if (logit) write (logunit,'(400A)') ' reading mask for ASCAT obs from ', trim(mask_filename) - + ! check if file exists - + inquire(file=mask_filename, exist=file_exists) - + if (.not. file_exists) then err_msg = 'Mask file for ASCAT obs not found!' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) @@ -1852,11 +1859,11 @@ subroutine read_obs_sm_ASCAT_EUMET( & ierr = nf90_inq_varid(ncid, 'll_lat', ll_lat_varid) ierr = nf90_inq_varid(ncid, 'd_lon', dlon_varid) ierr = nf90_inq_varid(ncid, 'd_lat', dlat_varid) - + ! get variable dimension IDs ierr = nf90_inq_dimid(ncid, 'lon', lon_dimid) ierr = nf90_inq_dimid(ncid, 'lat', lat_dimid) - + ! get dimension size ierr = nf90_inquire_dimension(ncid, lon_dimid, len=mask_N_lon) ierr = nf90_inquire_dimension(ncid, lat_dimid, len=mask_N_lat) @@ -1866,38 +1873,38 @@ subroutine read_obs_sm_ASCAT_EUMET( & ierr = nf90_get_var(ncid, ll_lat_varid, mask_ll_lat) ierr = nf90_get_var(ncid, dlon_varid, mask_dlon) ierr = nf90_get_var(ncid, dlat_varid, mask_dlat) - + ! allocate memory for mask allocate(mask_data(mask_N_lon, mask_N_lat)) ! note: lon-by-lat - + ! read mask ierr = nf90_get_var(ncid, mask_varid, mask_data) - + ! close netCDF mask file ierr = nf90_close(ncid) ! ---------------------------------------------------------------- ! ! select obs within assimilation window and from desired orbit direction; apply basic QC based on obs info - - N_tmp = 0 - + + N_tmp = 0 + do kk = 1,N_obs - + date_time_tmp%year = int(tmp_data(kk, 1)) date_time_tmp%month = int(tmp_data(kk, 2)) date_time_tmp%day = int(tmp_data(kk, 3)) date_time_tmp%hour = int(tmp_data(kk, 4)) date_time_tmp%min = int(tmp_data(kk, 5)) date_time_tmp%sec = int(tmp_data(kk, 6)) - + ! skip if record outside of current assim window if ( datetime_lt_refdatetime( date_time_tmp, date_time_low ) .and. & datetime_le_refdatetime( date_time_up, date_time_tmp ) ) cycle - + ! skip if record contains invalid soil moisture value if ( tmp_data(kk, 7) > 100. .or. tmp_data(kk, 7) < 0. ) cycle - + ! to distinguish orbit directions, must read "DOMO" from BUFR file ! ! 180 <= DOMO < 270 : descending @@ -1905,47 +1912,47 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! ! if (index(this_obs_param%descr,'_A') /=0 .and. (tmp_data(kk, 8) < 270 .or. tmp_data(kk, 8) > 360)) cycle ! if (index(this_obs_param%descr,'_D') /=0 .and. (tmp_data(kk, 8) < 180 .or. tmp_data(kk, 8) >= 270)) cycle - + ! skip if processing flag is set if(int(tmp_data(kk, 8)) /= 0) cycle - + ! skip if correction flag is set ("good" values are 0 and 4) if (.not. ( (int(tmp_data(kk, 9)) == 0) .or. (int(tmp_data(kk, 9)) == 4)) ) cycle - + ! skip if land fraction is missing or < 0.9 if(tmp_data(kk, 10) > 1. .or. tmp_data(kk, 10) < 0.9 ) cycle - + ! skip if topographic complexity > 10% if(tmp_data(kk, 11) > 10.) cycle - + ! skip if inundation and wetland fraction > 10% if(tmp_data(kk, 12) > 10.) cycle - + ! find lat/lon indices of ASCAT mask for this observation lat/lon mask_lat_ind = max( min( ceiling((tmp_data(kk, 13) - mask_ll_lat)/mask_dlat), mask_N_lat ), 1) mask_lon_ind = max( min( ceiling((tmp_data(kk, 14) - mask_ll_lon)/mask_dlon), mask_N_lon ), 1) ! skip if masked if (mask_data( mask_lon_ind, mask_lat_ind ) /= 0) cycle ! note: lon-by-lat - + N_tmp = N_tmp + 1 ! passed all QC if (N_tmp > max_rec) then err_msg = 'Too many obs have passed QC - how long is your assimilation window?' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + tmp1_jtime(N_tmp) = datetime_to_J2000seconds( date_time_tmp, J2000_epoch_id ) - - tmp1_lat( N_tmp) = tmp_data(kk, 13) + + tmp1_lat( N_tmp) = tmp_data(kk, 13) tmp1_lon( N_tmp) = tmp_data(kk, 14) tmp1_obs( N_tmp) = tmp_data(kk, 7)/100. ! change units from percent (0-100) to fraction (0-1) - + end do if (logit) then - + write (logunit,*) 'read_obs_sm_ASCAT_EUMET: read ', N_tmp, ' at date_time = ', date_time, ' from:' do ii=1,N_files write (logunit,*) trim(fnames(ii)) @@ -1953,14 +1960,14 @@ subroutine read_obs_sm_ASCAT_EUMET( & write (logunit,*) '----------' write (logunit,*) 'max(obs)=',maxval(tmp1_obs(1:N_tmp)), ', min(obs)=',minval(tmp1_obs(1:N_tmp)), & ', avg(obs)=',sum(tmp1_obs(1:N_tmp))/N_tmp - + end if - + deallocate(fnames) ! copy "good" obs with lat/lon/time into tmp_* (pointers) - allocate(tmp_jtime(N_tmp)) + allocate(tmp_jtime(N_tmp)) allocate(tmp_lon( N_tmp)) allocate(tmp_lat( N_tmp)) allocate(tmp_obs( N_tmp)) @@ -1969,20 +1976,20 @@ subroutine read_obs_sm_ASCAT_EUMET( & tmp_lon = tmp1_lon( 1:N_tmp) tmp_lat = tmp1_lat( 1:N_tmp) tmp_obs = tmp1_obs( 1:N_tmp) - + deallocate(tmp1_jtime) deallocate(tmp1_lon) deallocate(tmp1_lat) deallocate(tmp1_obs) deallocate(tmp_data) deallocate(mask_data) - + else - + N_tmp = 0 - + end if - + ! ---------------------------------------------------------------- ! ! 2.) for each observation @@ -1990,291 +1997,291 @@ subroutine read_obs_sm_ASCAT_EUMET( & ! b) determine tile within grid cell that contains lat/lon if (N_tmp>0) then - + allocate(tmp_tile_num(N_tmp)) - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_tmp, tmp_lat, tmp_lon, & this_obs_param, & tmp_tile_num ) - - + + ! ---------------------------------------------------------------- ! ! 3.) compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + ASCAT_sm = 0. ASCAT_lon = 0. ASCAT_lat = 0. ASCAT_time = 0.0D0 - + N_obs_in_tile = 0 - + do ii=1,N_tmp - + ind = tmp_tile_num(ii) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this step eliminates obs outside domain - + ASCAT_sm( ind) = ASCAT_sm( ind) + tmp_obs( ii) ASCAT_lon( ind) = ASCAT_lon( ind) + tmp_lon( ii) ASCAT_lat( ind) = ASCAT_lat( ind) + tmp_lat( ii) ASCAT_time(ind) = ASCAT_time(ind) + tmp_jtime(ii) - + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 - + end if - + end do ! normalize and set obs error std-dev do ii=1,N_catd - - ! set observation error standard deviation - + + ! set observation error standard deviation + ASCAT_sm_std(ii) = this_obs_param%errstd/100. ! change units from percent (0-100) to fraction (0-1) ! normalize if (N_obs_in_tile(ii)>1) then - + ASCAT_sm( ii) = ASCAT_sm( ii)/real(N_obs_in_tile(ii)) ASCAT_lon( ii) = ASCAT_lon( ii)/real(N_obs_in_tile(ii)) ASCAT_lat( ii) = ASCAT_lat( ii)/real(N_obs_in_tile(ii)) ASCAT_time( ii) = ASCAT_time(ii)/real(N_obs_in_tile(ii),kind(0.0D0)) - + elseif (N_obs_in_tile(ii)==0) then - + ASCAT_sm( ii) = this_obs_param%nodata ASCAT_lon( ii) = this_obs_param%nodata ASCAT_lat( ii) = this_obs_param%nodata ASCAT_time( ii) = real(this_obs_param%nodata,kind(0.0D0)) ASCAT_sm_std(ii) = this_obs_param%nodata - + else - + ! nothing to do if N_obs_in_tile(ii)==1 (and assuming N_obs_in_tile is never negative) - + end if - + end do - + ! clean up - + if (associated(tmp_tile_num)) deallocate(tmp_tile_num) - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if - + ! clean up - + if (associated(tmp_obs)) deallocate(tmp_obs) if (associated(tmp_lon)) deallocate(tmp_lon) if (associated(tmp_lat)) deallocate(tmp_lat) if (associated(tmp_jtime)) deallocate(tmp_jtime) - + end subroutine read_obs_sm_ASCAT_EUMET - + ! *************************************************************************** subroutine read_sm_ASCAT_bin( & N_files, fnames, N_data, lon, lat, sm_ASCAT, ease_col, ease_row ) - + ! read soil moisture data from one or more ASCAT bin files ! ! return ONLY valid data points (ie. excluding no-data-values) - ! + ! ! no QC in addition to what was done in matlab-preprocessing ! ! DRAPER, May 2011 ! based on read_ae_sm_LPRM_bin ! ! DRAPER, July 2012 - ! updated, for inclusion of SDS error in ASCAT file + ! updated, for inclusion of SDS error in ASCAT file ! (error info currently not saved) implicit none - + integer, intent(in) :: N_files - + character(*), dimension(N_files), intent(in) :: fnames - + integer, intent(out) :: N_data - + real, dimension(:), pointer :: lon, lat, sm_ASCAT ! output - + integer, dimension(:), pointer, optional :: ease_col, ease_row ! output - + ! local variables - + logical :: must_stop - + integer, dimension(N_files) :: N_data_tmp - + integer :: i, j, k_off - + integer, dimension(:), allocatable :: tmpintvec real, dimension(:), allocatable :: tmprealvec - + character(len=*), parameter :: Iam = 'read_sm_ASCAT_bin' character(len=400) :: err_msg ! ------------------------------------------------------------- - + ! make sure pointers are not allocated or associated - + must_stop = .false. - + if ( associated(lon) .or. associated(lat) .or. associated(sm_ASCAT) ) then must_stop = .true. end if - - if ( present(ease_col) ) then + + if ( present(ease_col) ) then if (associated(ease_col)) must_stop = .true. end if - + if ( present(ease_row) ) then if (associated(ease_row)) must_stop = .true. end if - + if (must_stop) then err_msg = 'output pointers must not be ' // & 'associated or allocated on input.' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - - + + !write (logunit,*) 'read_ae_sm_LPRM_bin(): output pointers must not be ' // & ! 'associated or allocated on input. STOPPING.' !stop ! !end if - + ! determine number of data to be read from each file - + N_data = 0 - + do j=1,N_files - - ! open file - + + ! open file + open( 10, file=trim(fnames(j)), form='unformatted',convert='big_endian', action='read' ) - + read( 10) N_data_tmp(j) - + close(10, status='keep') - + end do - + ! allocate pointers (must be deallocated outside this subroutine!) - + N_data = sum(N_data_tmp) - + allocate(lon( N_data)) allocate(lat( N_data)) allocate(sm_ASCAT(N_data)) - + if (present(ease_col)) allocate(ease_col(N_data)) if (present(ease_row)) allocate(ease_row(N_data)) - + ! read data into arrays, concatenate data from N_files files - + ! format of AMSR_sm_LPRM_EASE_bin files: ! ! record 1 -- N_data int*4 - ! record 2 -- lon( 1:N_data) real*4 - ! record 3 -- lat( 1:N_data) real*4 - ! record 4 -- sds( 1:N_data) real*4 - ! record 5 -- sds_noise (1:N_data) real*4 + ! record 2 -- lon( 1:N_data) real*4 + ! record 3 -- lat( 1:N_data) real*4 + ! record 4 -- sds( 1:N_data) real*4 + ! record 5 -- sds_noise (1:N_data) real*4 ! record 6 -- ind_i(1:N_data) int*4 zero-based (!) EASE row index ! record 7 -- ind_j(1:N_data) int*4 zero-based (!) EASE col index - + k_off = 0 - + do j=1,N_files - + allocate(tmprealvec(N_data_tmp(j))) - + if (present(ease_col)) allocate(tmpintvec(N_data_tmp(j))) - + open (10, file=trim(fnames(j)), form='unformatted', convert='big_endian', action='read' ) - + ! re-read N_data - - read (10) N_data_tmp(j) - - ! read data as needed - - read (10) tmprealvec; lon(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - read (10) tmprealvec; lat(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - + + read (10) N_data_tmp(j) + + ! read data as needed + + read (10) tmprealvec; lon(k_off+1:k_off+N_data_tmp(j)) = tmprealvec + read (10) tmprealvec; lat(k_off+1:k_off+N_data_tmp(j)) = tmprealvec + read (10) tmprealvec; sm_ASCAT(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - - ! skip record with error - + + ! skip record with error + read (10) ! tmprealvec; er_ASCAT(k_off+1:k_off+N_data_tmp(j)) = tmprealvec - + if (present(ease_col) .and. present(ease_row)) then - + read (10) tmpintvec; ease_col(k_off+1:k_off+N_data_tmp(j)) = tmpintvec read (10) tmpintvec; ease_row(k_off+1:k_off+N_data_tmp(j)) = tmpintvec - + end if - + ! clean up - + close(10, status='keep') - + deallocate(tmprealvec) if (allocated(tmpintvec)) deallocate(tmpintvec) - + ! prepare next j - + k_off = k_off + N_data_tmp(j) - + end do - + ! ------------------------------------- ! - ! eliminate no-data-values - + ! eliminate no-data-values + j = 0 - + do i=1,N_data - + if (sm_ASCAT(i)>0.) then ! any neg is nodata - + j=j+1 - + sm_ASCAT(j) = sm_ASCAT(i) lon(j) = lon(i) lat(j) = lat(i) if (present(ease_col)) ease_col(j) = ease_col(i) if (present(ease_row)) ease_row(j) = ease_row(i) - + end if - + end do - + N_data = j - + end subroutine read_sm_ASCAT_bin ! ***************************************************************** - + subroutine read_obs_LaRC_Tskin( & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & @@ -2282,158 +2289,158 @@ subroutine read_obs_LaRC_Tskin( & found_obs, ts_LARC, std_ts_LARC ) !--------------------------------------------------------------------- - ! - ! Subroutine to read in Tskin from LaRC. - ! Draper, June 2012. + ! + ! Subroutine to read in Tskin from LaRC. + ! Draper, June 2012. ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim, N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: ts_LARC real, intent(out), dimension(N_catd) :: std_ts_LARC logical, intent(out) :: found_obs - + ! --------------- - + real, parameter :: min_Tskin = 200. real, parameter :: max_Tskin = 370. - + real, parameter :: tol = 1.e-2 - - ! Offsets taken - + + ! Offsets taken + integer :: ts_time_offset ! in seconds - + character(6) :: DDHHMM character(6) :: YYYYMM - + type(date_time_type) :: date_time_tmp - + integer :: i, ind, N_tmp, N_files - + ! 24 = max number of files in one cycle - ! (assumes daily assim cycle, and hourly files) - + ! (assumes daily assim cycle, and hourly files) + character(300), dimension(24) :: fnames - + real, dimension(:), pointer :: tmp_obs, tmp_lat, tmp_lon - + integer, dimension(N_catd) :: N_obs_in_tile - + integer, dimension(:), pointer :: tmp_tile_num - + character(300) :: tmp_fname - + logical :: ex - - integer :: MM - + + integer :: MM + character(len=*), parameter :: Iam = 'read_obs_LaRC_Tskin' ! ------------------------------------------------------------------- - - nullify( tmp_obs, tmp_lat, tmp_lon , tmp_tile_num) - + + nullify( tmp_obs, tmp_lat, tmp_lon , tmp_tile_num) + ! --------------- - + ! initialize - + ! time stampes are at start of scan. MET-09 scans much faster than others - + select case (trim(this_obs_param%descr)) - - case ('LaRC_tskin-GOESE','LaRC_tskin-GOESW', 'LaRC_tskin-FY2E-') - - ts_time_offset=26*60/2 + + case ('LaRC_tskin-GOESE','LaRC_tskin-GOESW', 'LaRC_tskin-FY2E-') + + ts_time_offset=26*60/2 case ('LaRC_tskin-MTST2') - - ts_time_offset=28*60/2 - - case ('LaRC_tskin-MET09') - + + ts_time_offset=28*60/2 + + case ('LaRC_tskin-MET09') + ts_time_offset=12*60/2 - - case default - + + case default + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown obs_param%descr') - + end select - + found_obs = .false. - - ! find files that are within half-open interval + + ! find files that are within half-open interval ! [date_time-dtstep_assim/2,date_time+dtstep_assim/2) select case (trim(this_obs_param%descr)) case ('LaRC_tskin-GOESW') MM=00 - case ('LaRC_tskin-GOESE') + case ('LaRC_tskin-GOESE') MM=45 - case ('LaRC_tskin-MET09') + case ('LaRC_tskin-MET09') MM=00 case ('LaRC_tskin-FY2E-') MM=00 case ('LaRC_tskin-MTST2') MM=30 - case default + case default call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown obs_param%descr') end select - - date_time_tmp = date_time - + + date_time_tmp = date_time + call augment_date_time( -(dtstep_assim/2 + ts_time_offset), date_time_tmp ) - + ! identify all files within current assimilation interval ! inquire for the file once per hour - will always be at same minutes past hour - - N_files=0 - + + N_files=0 + do i=1,(dtstep_assim/3600) - + write (YYYYMM, '(i6.6)') date_time_tmp%year*100 + date_time_tmp%month write (DDHHMM, '(i6.6)') date_time_tmp%day*10000 + date_time_tmp%hour*100 + MM - + tmp_fname=trim(this_obs_param%path) // '/' // YYYYMM(1:4) // & '/' // YYYYMM(5:6) // '/' // DDHHMM(1:2) // '/' // & trim(this_obs_param%name) & // YYYYMM // DDHHMM(1:2) // '_' // DDHHMM(3:6) // 'z.nc4' - - inquire(file=trim(tmp_fname),exist=ex) - - if (ex) then - + + inquire(file=trim(tmp_fname),exist=ex) + + if (ex) then + N_files = N_files + 1 - + fnames(N_files) = tmp_fname - + end if - + call augment_date_time( 3600, date_time_tmp ) - + end do - + ! read observations: ! ! 1.) read N_tmp observations and their lat/lon info from file @@ -2443,360 +2450,360 @@ subroutine read_obs_LaRC_Tskin( & ! ---------------------------------------------------------------- ! ! 1.) read N_tmp observations and their lat/lon info from file - + if (N_files>0) then - + call read_LaRC_Tskin_nc4( & this_obs_param, N_files, fnames(1:N_files), & N_tmp, tmp_lon, tmp_lat, tmp_obs ) - + if (logit) then - + write (logunit,*) 'read_obs_LaRC_Tskin: read ', N_tmp, & ' at date_time = ', date_time do i=1,N_files write (logunit,*) trim(fnames(i)) end do write (logunit,*) '----------' - + end if else - + N_tmp = 0 - + end if ! 2.) get tile number for each obs from lat/lon if (N_tmp>0) then - + allocate(tmp_tile_num(N_tmp)) - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_tmp, tmp_lat, tmp_lon, & this_obs_param, & tmp_tile_num ) - + ! ---------------------------------------------------------------- ! ! 3.) compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + ts_LARC = 0. N_obs_in_tile = 0 - + do i=1,N_tmp - + ind = tmp_tile_num(i) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this step eliminates obs outside domain - + ! make sure obs is within allowed range - + if ( (min_Tskin1) then - + ts_LARC(i) = ts_LARC(i)/real(N_obs_in_tile(i)) - + elseif (N_obs_in_tile(i)==0) then - + ts_LARC(i) = this_obs_param%nodata - + end if - + end do - + if (associated(tmp_tile_num)) deallocate(tmp_tile_num) - + do i=1,N_catd - + std_ts_LARC(i) = this_obs_param%errstd - + ! CSD - temporary - only works over Americas ! if keep, base on SZA - - if ( (date_time_tmp%hour .GT. 2 ) .AND. (date_time_tmp%hour .LT. 14.5 )) then + + if ( (date_time_tmp%hour .GT. 2 ) .AND. (date_time_tmp%hour .LT. 14.5 )) then std_ts_LARC(i) = 1.3 ! 1.5 - - else - + + else + std_ts_LARC(i) = 2.1 ! 2.2 - + end if - - ! CSD - end temporary. + + ! CSD - end temporary. ! -------------------------------- - + end do - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if - + ! clean up - + if (associated(tmp_obs)) deallocate(tmp_obs) if (associated(tmp_lon)) deallocate(tmp_lon) if (associated(tmp_lat)) deallocate(tmp_lat) - + end subroutine read_obs_LaRC_Tskin ! *************************************************************************** subroutine read_LaRC_Tskin_nc4( & this_obs_param, N_files, fnames, N_data, lon, lat, ts_LaRC ) - - ! read Tskin from LaRC nc4 files + + ! read Tskin from LaRC nc4 files ! ! Apply QC: ! Screen longitude to retain obs only from closest GEOsat - ! Viewing zenith angle (VZA) + ! Viewing zenith angle (VZA) ! Cloud fraction - ! Sun zenith angle (SZA) - ! + ! Sun zenith angle (SZA) + ! ! returns number of data, and pointers to the lon, lat, and data ! return ONLY valid data points (ie. excluding no-data-values) ! ! Currently tailored to read in LaRC files that have been reprocessed to replace ! integer data with float data (to enable GFIO to be used) - ! - ! If LaRC fixes the files, will need to change: + ! + ! If LaRC fixes the files, will need to change: ! if ( nvars .LT 5 ) - this is included as some LaRC files are missing HIRESTSKIN field ! Replace lower case variable names with appropriate case - ! Replace nodata value (hardwired, as cannot read missing_value from file) + ! Replace nodata value (hardwired, as cannot read missing_value from file) ! Replace time setting with actual minutes (currently rounded down to nearest hour, - ! due to bug/assumption in GFIO read var routines). + ! due to bug/assumption in GFIO read var routines). ! DRAPER, June 2011 - + implicit none - + type(obs_param_type), intent(in) :: this_obs_param integer, intent(in) :: N_files - + character(*), dimension(N_files), intent(in) :: fnames - + integer, intent(out) :: N_data - + real, dimension(:), pointer :: lon, lat, ts_LaRC ! output - + ! local variables - + integer :: i, j, fid, rc integer :: x, x_min, x_max, y, y_min, y_max integer :: YYYYMMDD, HHMMSS integer :: g_nlon, g_nlat, km, lm, nvars, ngatts - + real :: nodata_LARC - + real, parameter :: tol=0.01 ! QC - + real, parameter :: max_fcld = 20. ! max total cloud fraction (%) real, parameter :: max_vza = 60. ! max viewing zenith angle real, parameter :: min_excl_sza = 82. ! min of sza exclusion interval real, parameter :: max_excl_sza = 90. ! max of sza exclusion interval logical :: tskin_ok, fcld_ok, sza_ok, vza_ok - + real, dimension(2) :: lon_range, lat_range real :: dlat, dlon - logical :: first - + logical :: first + real, dimension(:,:), allocatable :: tmp_data, ave_data, sum_data - + integer, dimension(:,:), allocatable :: cnt_data - + real, dimension(:,:), allocatable :: tmp_vza, tmp_sza, tmp_fcld ! ---------------------- ! - ! variables to fix viewing zenith angle bug in some of the 2012 GOES-East files + ! variables to fix viewing zenith angle bug in some of the 2012 GOES-East files ! reprocessed by Ben Scarino ! - ! the "tmp_goesEast_vza_*" parameters point to a dummy "vza.nc4" file that contains + ! the "tmp_goesEast_vza_*" parameters point to a dummy "vza.nc4" file that contains ! the correct vza values (located within the path of the reprocessed GOES-East files) - + character( 40), parameter :: tmp_goesEast_vza_dirname = 'goesEast201205_201304_RM2' character( 40), parameter :: tmp_goesEast_vza_fname = 'int2float/vza.nc4' - - integer, parameter :: tmp_goesEast_vza_YYYYMMDD = 20120629 - integer, parameter :: tmp_goesEast_vza_HHMMSS = 230000 - + + integer, parameter :: tmp_goesEast_vza_YYYYMMDD = 20120629 + integer, parameter :: tmp_goesEast_vza_HHMMSS = 230000 + character(300) :: tmp_fname - + integer :: ind, tmp_fid character(len=*), parameter :: Iam = 'read_LaRC_Tskin_nc4' character(len=400) :: err_msg - + ! ------------------------------------------------------------- - + ! make sure pointers are not allocated or associated - + if ( associated(lon) .or. associated(lat) .or. associated(ts_LaRC) ) then err_msg = 'output pointers must not be ' // & 'associated or allocated on input.' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + !----------------------------------------- ! read data into global arrays, perform QC !----------------------------------------- first = .true. - + do j=1,N_files - + if (logit) write(logunit,'(400A)') 'reading LaRC Tskin obs from ', trim(fnames(j)) call Gfio_Open ( fnames(j), 1, fid, rc ) - + if (rc<0) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'Error opening gfio file') end if - - ! temporary catch for files missing HIRESTSKIN (LaRC will fix this) - + + ! temporary catch for files missing HIRESTSKIN (LaRC will fix this) + call GFIO_DimInquire (fid,g_nlon,g_nlat,km,lm,nvars,ngatts,rc) - + dlon=360./real(g_nlon ) - dlat=180./real(g_nlat-1) - + dlat=180./real(g_nlat-1) + if (rc<0) then err_msg = 'DimInquire error, reading gfio file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - - !if (nvars .LT. 24) then (for original LaRC files) - + + !if (nvars .LT. 24) then (for original LaRC files) + if (nvars .LT. 5) then ! reprocessed files have only 5 variables - + write (logunit,*) 'CSD_LaRC file missing a variable ', fnames(j) - - else - + + else + if (first) then - + ! Attribute cannot be read from LaRC (or reprocessed) files ! no idea why this is not working. ncdump shows attribute is in file ! call GFIO_GetRealAtt ( fid, 'missing_value', 1, nodata_LARC, rc ) !if (rc<0) call stop_it('Get attribute error, reading nodata from gfio file') - + nodata_LARC=-9.99e33 ! For reprocessed files only!!! LaRC use different value - + write(logunit,*) & 'No-data-value manually set for reprocessed files: ', nodata_LARC - + ! get dimensions of grid and allocate temporary arrays - - allocate(tmp_data(g_nlon, g_nlat)) - allocate(ave_data(g_nlon, g_nlat)) - allocate(sum_data(g_nlon, g_nlat)) - allocate(cnt_data(g_nlon, g_nlat)) - allocate(tmp_vza( g_nlon, g_nlat)) - allocate(tmp_sza( g_nlon, g_nlat)) - allocate(tmp_fcld(g_nlon, g_nlat)) - + + allocate(tmp_data(g_nlon, g_nlat)) + allocate(ave_data(g_nlon, g_nlat)) + allocate(sum_data(g_nlon, g_nlat)) + allocate(cnt_data(g_nlon, g_nlat)) + allocate(tmp_vza( g_nlon, g_nlat)) + allocate(tmp_sza( g_nlon, g_nlat)) + allocate(tmp_fcld(g_nlon, g_nlat)) + sum_data=0 cnt_data=0 - + ! get dimensions of subgrid containing data for this disk lat_range=(/-52.0,52.0/) ! maximum range is +/-52.0 for VZA<60 - - ! specify boundary for min and max lon (select disk with lowest VZA) - + + ! specify boundary for min and max lon (select disk with lowest VZA) + select case (trim(this_obs_param%descr)) - + case ('LaRC_tskin-GOESW') - lon_range=(/-175. ,-105. /) - case ('LaRC_tskin-GOESE') - lon_range=(/-105. , -36.875/) - case ('LaRC_tskin-MET09') - lon_range=(/ -36.875, 54. /) + lon_range=(/-175. ,-105. /) + case ('LaRC_tskin-GOESE') + lon_range=(/-105. , -36.875/) + case ('LaRC_tskin-MET09') + lon_range=(/ -36.875, 54. /) case ('LaRC_tskin-FY2E-') ! bounds will not change if not using FY2, as FY2 eastern hemisphere is missing lon_range=(/ 54. , 90. /) case ('LaRC_tskin-MTST2') lon_range=(/ 90. , 180. /) ! avoid crossing the dateline - case default + case default err_msg = 'unknown obs_param%descr' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end select - + ! find index of smallest lon > min_lon and of largest lon <= max_lon - + x_min=floor( (minval(lon_range)-(-180.))/dlon ) + 2 x_max=floor( (maxval(lon_range)-(-180.))/dlon ) + 1 - + ! find index of smallest lat >= min_lat and of largest lat <= max_lat - + y_min=floor( (lat_range(1)-(-90.))/dlat ) + 1 y_max=floor( (lat_range(2)-(-90.))/dlat ) + 1 - - first = .false. - + + first = .false. + endif ! first - + ! cannot use actual time, as obs files are rounded to nearest hour - - ! ORIGINAL LaRC FILES + + ! ORIGINAL LaRC FILES ! -define time as "minutes since YYYYMMDD, HHMM", and have time=0 ! read(fnames(j)(len_trim(fnames(j))-8:len_trim(fnames(j))-5),'(i4)') HHMMSS ! HHMMSS=HHMMSS*100 ! add seconds - - ! lat4d.sh REPROCESSED FILES + + ! lat4d.sh REPROCESSED FILES ! -define time as "hours since YYYYMMDD, HH", and have time=MM/60 - ! The GFIO routine getbegdatetime calculates the time increment in a file by + ! The GFIO routine getbegdatetime calculates the time increment in a file by ! reading in first two values - ! If there is only one time value in the file, the second read (line 240) fails, + ! If there is only one time value in the file, the second read (line 240) fails, ! and incSecs=1 - ! This results in TimeIndex in GFIO_GetVar= (MM in seconds) (rather than 1) - ! (in summary, getbegdatetime assumes that if there are not multiple times in + ! This results in TimeIndex in GFIO_GetVar= (MM in seconds) (rather than 1) + ! (in summary, getbegdatetime assumes that if there are not multiple times in ! a file, time must equal 0 - ! which it does not for reprocessed LaRC files) - ! Get around this by specifying MM=0 regardless of actual time - + ! which it does not for reprocessed LaRC files) + ! Get around this by specifying MM=0 regardless of actual time + read(fnames(j)(len_trim(fnames(j)) -8:len_trim(fnames(j))-7),'(i2)') HHMMSS - + HHMMSS=HHMMSS*10000 ! add seconds, assume minutes are zero - + read(fnames(j)(len_trim(fnames(j))-17:len_trim(fnames(j))-9),'(i8)') YYYYMMDD - + ! read HIRESTSKIN @@ -2807,8 +2814,8 @@ subroutine read_LaRC_Tskin_nc4( & err_msg = 'GetVar error, reading hirestskin from gfio file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - - ! Viewing zenith angle (VZA) + + ! Viewing zenith angle (VZA) call GFIO_GetVar( fid,'vza', & YYYYMMDD, HHMMSS, g_nlon, g_nlat, & @@ -2817,30 +2824,30 @@ subroutine read_LaRC_Tskin_nc4( & err_msg = 'GetVar error, reading vza from gfio file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - - ! GOES-East VZA for end Sep, start Oct 2012 are incorrect in the files + + ! GOES-East VZA for end Sep, start Oct 2012 are incorrect in the files ! reprocessed by Ben Scarino, replace VZA with that from a hard-coded file - + ind = index(fnames(j), trim(tmp_goesEast_vza_dirname)) - + if (ind/=0) then ! extract path to replacement vza file from "fnames(j)" ! (=fnames(j) up to and including "tmp_goesEast_vza_dirname") - + tmp_fname = fnames(j)(1:ind+len_trim(tmp_goesEast_vza_dirname)-1) - + ! append "tmp_goesEast_vza_fname" - + tmp_fname = trim(tmp_fname) // '/' // trim(tmp_goesEast_vza_fname) - + call Gfio_Open (tmp_fname, 1, tmp_fid, rc ) - + if (rc<0) then err_msg = 'Error opening gfio file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + ! read replacement vza call GFIO_GetVar( tmp_fid,'vza', & @@ -2850,13 +2857,13 @@ subroutine read_LaRC_Tskin_nc4( & err_msg = 'GetVar error, reading vza from gfio file (2)' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + call GFIO_Close ( tmp_fid, rc ) - + end if - - ! Sun zenith angle (SZA) + + ! Sun zenith angle (SZA) call GFIO_GetVar( fid,'sza', & YYYYMMDD, HHMMSS, g_nlon, g_nlat, & @@ -2865,7 +2872,7 @@ subroutine read_LaRC_Tskin_nc4( & err_msg = 'GetVar error, reading sza from gfio file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + ! Cloud fraction (FCLD), level 1 is total cloud @@ -2878,15 +2885,15 @@ subroutine read_LaRC_Tskin_nc4( & end if ! calculate mean of all data that passes QC - + do x=x_min, x_max - + do y=y_min, y_max - + ! Tskin must not be no-data-value - + tskin_ok = (abs(tmp_data(x,y) - nodata_LARC) .GE. tol) - + ! fcld must not be no-data-value; fcld<=max_cld fcld_ok = & @@ -2894,191 +2901,191 @@ subroutine read_LaRC_Tskin_nc4( & (abs(tmp_fcld(x,y)-nodata_LARC) .GT. tol) ! sza must not be no-data-value; sza<=min_excl_sza; sza>=max_excl_sza - + sza_ok = & ( & (tmp_sza(x,y) .LE. min_excl_sza) .OR. & (tmp_sza(x,y) .GE. max_excl_sza) & ) .AND. & - (abs(tmp_sza(x,y)-nodata_LARC) .GT. tol) + (abs(tmp_sza(x,y)-nodata_LARC) .GT. tol) ! vza must not be no-data-value; vza<=vza_max - + vza_ok = & (tmp_vza( x,y) .LE. max_vza) .AND. & (abs(tmp_vza( x,y)-nodata_LARC) .GT. tol) - - + + if (tskin_ok .and. fcld_ok .and. sza_ok .and. vza_ok ) then - + sum_data(x,y) =sum_data(x,y) + tmp_data(x,y) - + cnt_data(x,y) = cnt_data(x,y)+1 - + end if - + end do end do - - end if ! nvars - + + end if ! nvars + call GFIO_Close ( fid, rc ) - + end do ! N_files - - - if ( .not. first) then - + + + if ( .not. first) then + ! only calc averages if found file with HIRESTSKIN - temporary for incomplete files - + ! calculate average over appropriate lat/lon range, and count locations with data - - N_data=0 - + + N_data=0 + do x=x_min, x_max do y=y_min, y_max - if(cnt_data(x,y)>0) then - ave_data(x,y)=sum_data(x,y)/cnt_data(x,y) + if(cnt_data(x,y)>0) then + ave_data(x,y)=sum_data(x,y)/cnt_data(x,y) N_data=N_data+1 else ave_data(x,y)=nodata_LARC endif enddo enddo - + ! allocate pointers for return vectors (must be deallocated outside this subroutine!) - + allocate(lon( N_data)) allocate(lat( N_data)) allocate(ts_LaRC(N_data)) - - ! pass return data into vectors - + + ! pass return data into vectors + i=1 do x=x_min, x_max do y=y_min, y_max - if ( abs(ave_data(x,y)-nodata_LARC) .GT. tol ) then - ts_LaRC(i)=ave_data(x,y) + if ( abs(ave_data(x,y)-nodata_LARC) .GT. tol ) then + ts_LaRC(i)=ave_data(x,y) lon(i)=(x-1)*dlon -180. - lat(i)=(y-1)*dlat - 90. + lat(i)=(y-1)*dlat - 90. i=i+1 endif enddo enddo - - ! clean up - - deallocate(tmp_data) - deallocate(sum_data) - deallocate(cnt_data) - deallocate(ave_data) - deallocate(tmp_vza) - deallocate(tmp_sza) - deallocate(tmp_fcld) - + + ! clean up + + deallocate(tmp_data) + deallocate(sum_data) + deallocate(cnt_data) + deallocate(ave_data) + deallocate(tmp_vza) + deallocate(tmp_sza) + deallocate(tmp_fcld) + else - + N_data=0 - + ! OK not to allocate pointers for data? endif - + end subroutine read_LaRC_Tskin_nc4 ! ***************************************************************** - + subroutine read_obs_RedArkOSSE_sm( & date_time, N_catd, tile_coord, this_obs_param, & found_obs, RedArkOSSE_sm, std_RedArkOSSE_sm ) - + ! Read observations of surface soil moisture from Wade Crow's ! synthetic RedArk OSSE 36km soil moisture files. - ! Set flag "found_obs" to true if observations are available + ! Set flag "found_obs" to true if observations are available ! for assimilation. ! ! If there are N > 1 observations in a given tile, ! a "super-observation" is computed by averaging the N observations. ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! reichle, 25 May 2006 ! reichle, 26 Sep 2006 - added iostat to open statement ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: RedArkOSSE_sm real, intent(out), dimension(N_catd) :: std_RedArkOSSE_sm logical, intent(out) :: found_obs - + ! --------------- - + ! locals ! RedArkOSSE synthetic soil moisture obs at 36km are available once a day. ! ! Mapping of the 36km retrievals to catchment/tile space is precomputed ! and stored in p2t_nearest.dat and t2p_nearest.dat - + integer, parameter :: N_p2t = 1120 ! # of 36 km pixels integer, parameter :: N_t2p = 69 ! # of tiles that need "duplicated" obs - + integer, parameter :: HH_obs = 21 ! obs hour of day (UTC) = 3pm CST - + ! initial QC parameters real, parameter :: obs_min = 0.0 ! min allowed obs real, parameter :: obs_max = 0.45 ! max allowed obs - real, parameter :: opac_max = 0.3 ! max allowed vegetation opacity - + real, parameter :: opac_max = 0.3 ! max allowed vegetation opacity + integer, parameter :: qc_failed_obs = -888. - + integer, dimension(N_p2t) :: p2t integer, dimension(N_t2p,2) :: t2p - + character(3) :: DDD character(4) :: YYYY character(300) :: tmpfname - + integer :: i, j, ind, N_tmp, tmp_tile_id, istat - + real :: tmp_real, tmp_opac - + real, dimension(:), allocatable :: tmp_obs integer, dimension(:), allocatable :: tmp_tile_num - - integer, dimension(N_catd) :: N_obs_in_tile + + integer, dimension(N_catd) :: N_obs_in_tile character(len=*), parameter :: Iam = 'read_obs_RedArkOSSE_sm' character(len=400) :: err_msg ! ------------------------------------------------------------------- - + ! initialize - + found_obs = .false. - + ! obs are available only once per day ! (hard-coded b/c time-of-day is not clear from filename) - + if (date_time%hour == HH_obs) then - + ! read observations: ! ! 1.) read observations, p2t, and t2p from files @@ -3087,330 +3094,330 @@ subroutine read_obs_RedArkOSSE_sm( & ! 4.) compute super-obs for each tile from all obs w/in that tile ! ! ---------------------------------------------------------------- - + ! 1.) a) read obs - + write (YYYY,'(i4.4)') date_time%year write (DDD, '(i3.3)') date_time%dofyr - + tmpfname = trim(this_obs_param%path) // '/' // YYYY // & '/' // trim(this_obs_param%name) // YYYY // '.' // DDD - + open(10, file=tmpfname, form='formatted', action='read', iostat=istat) - + if (istat==0) then - + N_tmp = N_p2t + N_t2p - + allocate(tmp_obs(N_tmp)) allocate(tmp_tile_num(N_tmp)) - + do i=1,N_p2t - + read(10,*) tmp_real, tmp_opac - - tmp_obs(i) = tmp_real/100. ! unit conversion - + + tmp_obs(i) = tmp_real/100. ! unit conversion + ! initial QC - + if ( (tmp_opac > opac_max) .or. & (tmp_obs(i) > obs_max) .or. & (tmp_obs(i) < obs_min) ) then - + tmp_obs(i) = qc_failed_obs - + end if - + end do - + close(10,status='keep') - + ! 1.) b) read p2t - + tmpfname = trim(this_obs_param%path) // '/p2t_nearest.dat' - + open(10, file=tmpfname, form='formatted', action='read') - + do i=1,N_p2t - + read(10,*) p2t(i), tmp_tile_id ! col 1: tile_num, col 2: tile_id - + ! check tile_id for consistency - + if (p2t(i)>0) then ! if statement added 1 May 2007, reichle if (tmp_tile_id/=tile_coord(p2t(i))%tile_id) then - + !write (logunit,*) 'read_obs_RedArkOSSE(): something is wrong (p2t)' !!!write (logunit,*) i, p2t(i), tmp_tile_id, tile_coord(p2t(i))%tile_id !stop err_msg = 'something is wrong (p2t)' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if end if - + end do - + close(10,status='keep') - + ! 1.) c) read t2p - + tmpfname = trim(this_obs_param%path) // '/t2p_nearest.dat' - + open(10, file=tmpfname, form='formatted', action='read') - + do i=1,N_t2p - + read(10,*) t2p(i,1:2), tmp_tile_id ! pixel_num, tile_num, tile_id - + if (tmp_tile_id/=tile_coord(t2p(i,2))%tile_id) then err_msg = 'something is wrong (t2p)' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + !write (logunit,*) 'read_obs_RedArkOSSE(): something is wrong (t2p)' !stop ! !end if - + end do - + close(10,status='keep') - + ! ----------------------- ! ! 2.) for each observation determine tile_num ("p2t") - + do i=1,N_p2t - + tmp_tile_num(i) = p2t(i) - + end do - + ! ----------------------- ! ! 3.) duplicate observations for tiles not covered yet ("t2p") - + do i=1,N_t2p - + j = i + N_p2t - + tmp_obs( j) = tmp_obs(t2p(i,1)) - + tmp_tile_num(j) = t2p(i,2) - + end do - - + + ! ------------------------ ! ! 4.) compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + RedArkOSSE_sm = 0. N_obs_in_tile = 0 - + do i=1,N_tmp - + ind = tmp_tile_num(i) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this step eliminates obs outside domain - + if (tmp_obs(i)>0.) then ! this step eliminates no-data - + RedArkOSSE_sm(ind) = RedArkOSSE_sm(ind) + tmp_obs(i) - + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 - + end if - + end if - + end do - + ! normalize - + do i=1,N_catd - + if (N_obs_in_tile(i)>1) then - + RedArkOSSE_sm(i) = RedArkOSSE_sm(i)/real(N_obs_in_tile(i)) - + elseif (N_obs_in_tile(i)==0) then - + RedArkOSSE_sm(i) = this_obs_param%nodata - + end if - + end do - + ! -------------------------------- - + ! set observation error standard deviation - + do i=1,N_catd std_RedArkOSSE_sm(i) = this_obs_param%errstd end do - + ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if ! istat==0 - + ! clean up - + if (allocated(tmp_tile_num)) deallocate(tmp_tile_num) if (allocated(tmp_obs)) deallocate(tmp_obs) - + end if - + end subroutine read_obs_RedArkOSSE_sm - + ! ******************************************************************* - + subroutine read_obs_RedArkOSSE_CLSMsynthSM( date_time, N_catd, & this_obs_param, & found_obs, RedArkOSSE_CLSMsynthSM, std_RedArkOSSE_CLSMsynthSM ) - + ! Read synthetic observations of surface soil moisture from CLSM - ! RedArkOSSE integration (generated in matlab from innov output with + ! RedArkOSSE integration (generated in matlab from innov output with ! get_RedArk_CLSM_synth_retrievals.m) ! ! synthetic RedArk OSSE 36km soil moisture files. - ! Set flag "found_obs" to true if observations are available + ! Set flag "found_obs" to true if observations are available ! for assimilation. ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! reichle, 16 Feb 2006 - ! reichle, 5 Apr 2007 - use obs std from default nml input + ! reichle, 5 Apr 2007 - use obs std from default nml input ! (instead of reading from file) ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: N_catd - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: RedArkOSSE_CLSMsynthSM real, intent(out), dimension(N_catd) :: std_RedArkOSSE_CLSMsynthSM logical, intent(out) :: found_obs - + ! --------------- - + ! locals ! RedArkOSSE CLSM synthetic soil moisture obs are available once a day. - + integer, parameter :: HH_obs = 21 ! obs hour of day (UTC) = 3pm CST - + real, parameter :: nodata = -9999. - + ! initial QC parameters - + real, parameter :: obs_min = 0.0 ! min allowed obs real, parameter :: obs_max = 0.5 ! max allowed obs - + character(2) :: MM, DD character(4) :: YYYY, HHMM character(300) :: tmpfname - + integer :: i, N_tmp, istat, tmp_tilenum - + real :: tmp_obs, tmp_obs_std - + ! ------------------------------------------------------------------- - + ! initialize - + found_obs = .false. - + ! obs are available only once per day ! (hard-coded b/c time-of-day is not very clear in RedArkOSSE - + if (date_time%hour == HH_obs) then - + write (YYYY,'(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month write (DD, '(i2.2)') date_time%day write (HHMM,'(i4.4)') 100*HH_obs - + tmpfname = trim(this_obs_param%path) // '/Y' // YYYY // '/M' // MM // & '/' // trim(this_obs_param%name) // YYYY // MM // DD // '_' // HHMM - + open(10, file=tmpfname, form='formatted', action='read', iostat=istat) - + if (istat==0) then - + read(10,*) N_tmp read(10,*) ! tmp_obs_std - ! do NOT use obs std from file + ! do NOT use obs std from file ! (s.t. "wrong" obs std can be specified conveniently in nml file) ! reichle, 5 Apr 2007 - + tmp_obs_std = this_obs_param%errstd - + RedArkOSSE_CLSMsynthSM( 1:N_catd) = nodata std_RedArkOSSE_CLSMsynthSM(1:N_catd) = tmp_obs_std - + do i=1,N_tmp - + read(10,*) tmp_tilenum, tmp_obs - + ! initial QC - + tmp_obs = min( obs_max, tmp_obs ) tmp_obs = max( obs_min, tmp_obs ) - + RedArkOSSE_CLSMsynthSM( tmp_tilenum ) = tmp_obs - + end do - + close(10,status='keep') - + found_obs = .true. - + end if ! istat==0 - + end if - + end subroutine read_obs_RedArkOSSE_CLSMsynthSM ! ***************************************************************** - + subroutine read_obs_VivianaOK_CLSMsynthSM( date_time, N_catd, & this_obs_param, & found_obs, VivianaOK_CLSMsynthSM, std_VivianaOK_CLSMsynthSM ) - + ! Read synthetic observations of surface soil moisture from CLSM - ! integration over Viviana's OK domain + ! integration over Viviana's OK domain ! - ! Set flag "found_obs" to true if observations are available + ! Set flag "found_obs" to true if observations are available ! for assimilation. ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! vmaggion + reichle, 4 Aug 2008: @@ -3418,313 +3425,313 @@ subroutine read_obs_VivianaOK_CLSMsynthSM( date_time, N_catd, & ! use obs std from default nml input (instead of reading from file) ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: N_catd - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: VivianaOK_CLSMsynthSM real, intent(out), dimension(N_catd) :: std_VivianaOK_CLSMsynthSM logical, intent(out) :: found_obs - + ! --------------- - + ! locals ! VivianaOK CLSM synthetic soil moisture obs are available once a day. - + integer, parameter :: HH_obs = 12 ! obs hour of day (UTC) = 6am CST - + real, parameter :: nodata = -9999. - + ! initial QC parameters - + real, parameter :: obs_min = 0.0 ! min allowed obs real, parameter :: obs_max = 0.5 ! max allowed obs - + character(2) :: MM, DD, HH character(4) :: YYYY character(300) :: tmpfname - + integer :: i, istat - + real :: tmp_obs, tmp_obs_std - + ! ------------------------------------------------------------------- - + ! initialize - + found_obs = .false. - + ! obs are available only once per day ! (hard-coded b/c time-of-day is not very clear in VivianaOK - + if (date_time%hour == HH_obs) then - + write (YYYY,'(i4.4)') date_time%year write (MM, '(i2.2)') date_time%month write (DD, '(i2.2)') date_time%day write (HH, '(i2.2)') HH_obs - + tmpfname = trim(this_obs_param%path) // '/' // & trim(this_obs_param%name) // YYYY // MM // DD // '_' // HH // 'z.txt' - + open(10, file=tmpfname, form='formatted', action='read', iostat=istat) - + if (istat==0) then - - ! do NOT use obs std from file + + ! do NOT use obs std from file ! (s.t. "wrong" obs std can be specified conveniently in nml file) - + tmp_obs_std = this_obs_param%errstd - + VivianaOK_CLSMsynthSM( 1:N_catd) = nodata std_VivianaOK_CLSMsynthSM(1:N_catd) = tmp_obs_std - + do i=1,N_catd - + read(10,*) tmp_obs - + ! initial QC - + tmp_obs = min( obs_max, tmp_obs ) tmp_obs = max( obs_min, tmp_obs ) - + VivianaOK_CLSMsynthSM(i) = tmp_obs - + if (abs(tmp_obs-nodata)>1e-4) found_obs = .true. - + end do - + close(10,status='keep') - + end if ! istat==0 - + end if - + end subroutine read_obs_VivianaOK_CLSMsynthSM ! ***************************************************************** - + subroutine read_obs_RedArkOSSE_truth( & date_time, N_catd, this_obs_param, & found_obs, RedArkOSSE_truth, std_RedArkOSSE_truth ) - + ! Read "observations" of true surface soil moisture from Wade Crow's - ! truth soil moisture in CLSM catchment space. - ! Set flag "found_obs" to true if observations are available + ! truth soil moisture in CLSM catchment space. + ! Set flag "found_obs" to true if observations are available ! for assimilation. ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! reichle, 18 Sep 2006 ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: N_catd - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: RedArkOSSE_truth real, intent(out), dimension(N_catd) :: std_RedArkOSSE_truth logical, intent(out) :: found_obs - + ! --------------- - + ! locals ! RedArkOSSE soil moisture truth available once a day. ! ! Truth is stored directly in Catchment space - + integer, parameter :: HH_obs = 21 ! obs hour of day (UTC) = 3pm CST - + character(2), parameter :: HH = '15' ! might be CST??? - + character(3) :: DDD character(4) :: YYYY character(300) :: tmpfname - + integer :: i, istat - + real :: tmp_real - + ! ------------------------------------------------------------------- - + ! initialize - + found_obs = .false. - + ! obs are available only once per day ! (hard-coded b/c time-of-day is not clear from filename) - + if (date_time%hour == HH_obs) then - + ! read observations from file - + write (YYYY,'(i4.4)') date_time%year write (DDD, '(i3.3)') date_time%dofyr - + ! file name - + tmpfname = trim(this_obs_param%path) // '/' // YYYY // & '/' // trim(this_obs_param%name) // YYYY // '.' // DDD // '.' // HH - + ! open file and read obs if available - + open(10, file=tmpfname, form='formatted', action='read', iostat=istat) - + if (istat==0) then - + do i=1,N_catd - + read(10,*) tmp_real - - RedArkOSSE_truth(i) = tmp_real/100. ! unit conversion - + + RedArkOSSE_truth(i) = tmp_real/100. ! unit conversion + end do - + close(10,status='keep') - + ! set observation error standard deviation - + do i=1,N_catd std_RedArkOSSE_truth(i) = this_obs_param%errstd end do - + found_obs = .true. - + end if end if end subroutine read_obs_RedArkOSSE_truth ! ***************************************************************** - + subroutine read_obs_isccp_tskin_gswp2_v1( & date_time, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, isccp_tskin_gswp2_v1, std_isccp_tskin_gswp2_v1 ) - + ! Read observations of land skin temperature from ISCCP data ! produced by Sarith on GSWP-2 grid - ! Set flag "found_obs" to true if observations are available + ! Set flag "found_obs" to true if observations are available ! for assimilation. ! ! If there are N > 1 observations in a given tile, ! a "super-observation" is computed by averaging the N observations ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! reichle, 26 Sep 2005 ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: isccp_tskin_gswp2_v1 real, intent(out), dimension(N_catd) :: std_isccp_tskin_gswp2_v1 logical, intent(out) :: found_obs - + ! --------------- - + ! locals - + integer, parameter :: N_gswp2_compressed = 15238 - - ! land_i_gswp2 and land_j_gswp2 as stored in + + ! land_i_gswp2 and land_j_gswp2 as stored in ! ISCCP_Tskin_GSWP2_grid_V1 files (by Sarith) follow the GSWP2 convention ! for grid orientation, that is counting from north-to-south ! and from west-to-east real, parameter :: minlon_gswp2 = -180.5 real, parameter :: maxlat_gswp2 = 90.5 - + real, parameter :: dx_gswp2 = 1. real, parameter :: dy_gswp2 = 1. ! parameters for initial quality control and no-data-value treatment - + real, parameter :: tskin_min = 200. real, parameter :: tskin_max = 400. - + ! ISCCP_Tskin_GSWP2_grid_V1 files are available every 3h - + character(2) :: HH character(4) :: YYYY, MMDD character(300) :: fname - + integer :: i, j, ind, istat, N_tmp - + real :: tsclr - + integer :: land_i_gswp2, land_j_gswp2 - + integer, dimension(N_gswp2_compressed) :: tmp_tile_num - + real, dimension(N_gswp2_compressed) :: tmp_obs, tmp_lat, tmp_lon - - integer, dimension(N_catd) :: N_obs_in_tile - + + integer, dimension(N_catd) :: N_obs_in_tile + ! ------------------------------------------------------------------- - + ! initialize - + found_obs = .false. ! assemble file name - + write (YYYY,'(i4.4)') date_time%year write (MMDD,'(i4.4)') date_time%month*100 + date_time%day write (HH, '(i2.2)') date_time%hour - + fname = trim(this_obs_param%path) // '/' // '/Y' // YYYY // & '/M' // MMDD(1:2) // '/' // trim(this_obs_param%name) & // YYYY // MMDD // '_' // HH // 'z.bin' - + open(10, file=fname, form='unformatted', convert='big_endian', & action='read', status='old', iostat=istat) - + ! read observations: ! ! 1.) read N_tmp observations and their lat/lon info from file @@ -3736,209 +3743,209 @@ subroutine read_obs_isccp_tskin_gswp2_v1( & ! ---------------------------------------------------------------- ! ! 1.) read N_tmp observations and their lat/lon info from file - + if (istat==0) then - + j = 0 - + do i=1,N_gswp2_compressed - + read (10) land_i_gswp2, land_j_gswp2, tsclr - + ! eliminate no-data-values (specify range of acceptable tskin) - - if ( (tsclr>tskin_min) .and. & - (tsclrtskin_min) .and. & + (tsclr0) then - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_tmp, tmp_lat(1:N_tmp), tmp_lon(1:N_tmp), & this_obs_param, & tmp_tile_num(1:N_tmp) ) - + ! ---------------------------------------------------------------- ! ! 3.) compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + isccp_tskin_gswp2_v1 = 0. N_obs_in_tile = 0 - + do i=1,N_tmp - + ind = tmp_tile_num(i) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this step eliminates obs outside domain - + isccp_tskin_gswp2_v1(ind) = isccp_tskin_gswp2_v1(ind) + tmp_obs(i) - + N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 - + end if - + end do - + ! normalize - + do i=1,N_catd - + if (N_obs_in_tile(i)>1) then - + isccp_tskin_gswp2_v1(i) = & isccp_tskin_gswp2_v1(i)/real(N_obs_in_tile(i)) - + elseif (N_obs_in_tile(i)==0) then - + isccp_tskin_gswp2_v1(i) = this_obs_param%nodata - + end if - + end do - - + + ! -------------------------------- - + ! set observation error standard deviation - + do i=1,N_catd std_isccp_tskin_gswp2_v1(i) = this_obs_param%errstd end do - + ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if - + end subroutine read_obs_isccp_tskin_gswp2_v1 - + ! ***************************************************************** - + subroutine read_obs_isccp_tskin_ceop3n4( & date_time, N_catd, tile_coord, & this_obs_param, & found_obs, isccp_tskin_gswp2_v1, std_isccp_tskin_gswp2_v1 ) - + ! *** ONLY for "CEOP3n4_by_tile_FV_144x91" domain *** ! ! Read observations of land skin temperature from ISCCP data - ! produced by Sarith on GSWP-2 grid - ! Set flag "found_obs" to true if observations are available + ! produced by Sarith on GSWP-2 grid + ! Set flag "found_obs" to true if observations are available ! for assimilation (even if only "nodata" values in file...). ! ! inputs to this subroutine: - ! date_time = current model date and time + ! date_time = current model date and time ! N_catd = number of catchments in domain ! ! reichle, 27 Jan 2009 ! ! -------------------------------------------------------------------- - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! outputs: - + real, intent(out), dimension(N_catd) :: isccp_tskin_gswp2_v1 real, intent(out), dimension(N_catd) :: std_isccp_tskin_gswp2_v1 logical, intent(out) :: found_obs - + ! --------------- - + ! locals - + integer, parameter :: N_gswp2_compressed = 15238 - + ! parameters for initial quality control and no-data-value treatment - + real, parameter :: tskin_min = 200. real, parameter :: tskin_max = 400. integer, parameter :: N_tiles = 41 - ! the following mapping is from + ! the following mapping is from ! /land/l_data/CEOP/EOP3n4/coord/map_CEOP3n4_to_ISCCP_GSWP2.txt - ! produced with + ! produced with ! land01:/home/reichle/GMAO/station_data/CEOP/EOP3n4/matlab/map_GEOS5_to_ISCCP_GSWP2.m ! ISCCP_Tskin_GSWP2_grid_V1 files are available every 3h - + character(2) :: HH character(4) :: YYYY, MMDD character(300) :: fname - + integer :: i, istat - + real, dimension(N_gswp2_compressed) :: tsclr - + real :: tmp_obs - + integer :: land_i_gswp2, land_j_gswp2 character(len=*), parameter :: Iam = 'read_obs_isccp_tskin_ceop3n4' character(len=400) :: err_msg - + integer, dimension(2,N_tiles) :: GEOS5_to_ISCCP ! ------------------------------------------ - - + + GEOS5_to_ISCCP = reshape( & - (/ & - 64402, 96, & + (/ & + 64402, 96, & 68663, 1687, & 68677, 1686, & 68771, 1792, & @@ -3982,61 +3989,61 @@ subroutine read_obs_isccp_tskin_ceop3n4( & /), & shape(GEOS5_to_ISCCP) & ) - - + + ! ------------------------------------------------------------------- if (N_catd/=N_tiles) then err_msg = 'error 1 -- use only for CEOP3n4_by_tile_144x91 domain' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + !write (logunit,*) 'error 1 -- use only for CEOP3n4_by_tile_144x91 domain' !stop !end if - + ! initialize - + found_obs = .false. isccp_tskin_gswp2_v1 = this_obs_param%nodata std_isccp_tskin_gswp2_v1 = this_obs_param%nodata - + ! assemble file name - + write (YYYY,'(i4.4)') date_time%year write (MMDD,'(i4.4)') date_time%month*100 + date_time%day write (HH, '(i2.2)') date_time%hour - + fname = trim(this_obs_param%path) // '/' // '/Y' // YYYY // & '/M' // MMDD(1:2) // '/' // trim(this_obs_param%name) & // YYYY // MMDD // '_' // HH // 'z.bin' - + open(10, file=fname, form='unformatted', convert='big_endian', & action='read', status='old', iostat=istat) - + if (istat==0) then - + ! set found_obs=true (even if only nodata values in file) - + found_obs = .true. - + ! read observations from ISCCP file - + do i=1,N_gswp2_compressed - + read (10) land_i_gswp2, land_j_gswp2, tsclr(i) - + end do - + ! extract obs for CEOP3n4_by_tile_FV_144x91 domain - + do i=1,N_tiles - + ! double-check tile ids - - if (GEOS5_to_ISCCP(1,i)/=tile_coord(i)%tile_id) then - + + if (GEOS5_to_ISCCP(1,i)/=tile_coord(i)%tile_id) then + !write (logunit,*) 'error 2 -- only for CEOP3n4_by_tile_144x91 domain' !stop @@ -4044,53 +4051,53 @@ subroutine read_obs_isccp_tskin_ceop3n4( & call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) else - + tmp_obs = tsclr(GEOS5_to_ISCCP(2,i)) - + ! basic QC - - if ( (tmp_obs>tskin_min) .and. (tmp_obstskin_min) .and. (tmp_obs0) then - + tmp_obs = tsclr(GEOS5_to_ISCCP(2,i)) - + else - + tmp_obs = this_obs_param%nodata - + end if - + ! basic QC - - if ( (tmp_obs>tskin_min) .and. (tmp_obstskin_min) .and. (tmp_obsN_files_max) then err_msg = 'too many files found' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + fnames(k) = tmpfname - + end if - + call augment_date_time( dtstep_file, date_time_low ) - + end do ! end do while time step loop - + N_files = k ! --------------------------------------------------------------- @@ -4633,104 +4640,104 @@ subroutine read_obs_SMOS( date_time, N_catd, this_obs_param, & allocate(start_time( N_files, 5)) allocate(end_time( N_files, 5)) - + allocate(Asc_flag( N_files )) allocate(N_obs_tmp( N_files )) allocate(N_ang_tmp( N_files )) - + ! open files, read and interpret headers - + do k=1,N_files - + unitnum(k) = unitnum_off + k - + open(unitnum(k), file=trim(fnames(k)),form='unformatted', convert='big_endian',status='old', & access='SEQUENTIAL', iostat=istat) - + if (istat/=0) then - + err_msg = 'could not open file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + else - + if (logit) write (logunit,'(400A)') 'reading file ' // trim(fnames(k)) - + end if - + read(unitnum(k)) Asc_flag( k ) ! could also read version no. read(unitnum(k)) start_time(k,:) read(unitnum(k)) end_time( k,:) read(unitnum(k)) N_obs_tmp( k ), N_ang_tmp(k) - + if (logit) write (logunit,*) ' Asc_flag = ', Asc_flag( k ) if (logit) write (logunit,*) ' start_time = ', start_time( k,:) if (logit) write (logunit,*) ' end_time = ', end_time( k,:) if (logit) write (logunit,*) ' N_obs_tmp = ', N_obs_tmp( k ) if (logit) write (logunit,*) ' N_ang_tmp = ', N_ang_tmp( k ) - + end do - + ! make sure N_ang is same in all files - + N_ang = N_ang_tmp(1) - + do k=2,N_files - + if ( N_obs_tmp(k)>0 .and. (N_ang/=N_ang_tmp(k)) ) then err_msg = 'angles differ between files' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + end do - + ! allocate data variables - + N_obs = sum(N_obs_tmp(1:N_files)) - + if (N_obs>0) then - + allocate(tmp_ang(N_ang)) - + allocate(tmp_lon(N_obs)) allocate(tmp_lat(N_obs)) - + allocate(tmp_obs(N_obs)) allocate(tmp_std(N_obs)) allocate(tmp_cnt(N_obs)) allocate(tmp_file_ind(N_obs)) - + if (SM_files) allocate(tmp_DQX(N_obs)) - + ! loop through files and read data - + ind_end = 0 - + do k=1,N_files - + if (N_obs_tmp(k)>0) then - + ind_start = ind_end + 1 - + ind_end = ind_start + N_obs_tmp(k) - 1 - + ! record to which file each obs belongs - + tmp_file_ind(ind_start:ind_end) = k ! continue with reading file (scalars in header were read above) - + read(unitnum(k)) tmp_ang ! assume same angles in all files read(unitnum(k)) tmp_lon(ind_start:ind_end) - read(unitnum(k)) tmp_lat(ind_start:ind_end) - - + read(unitnum(k)) tmp_lat(ind_start:ind_end) + + if (SM_files) then - + ! read SMOS SM file - + read(unitnum(k)) tmp_obs(ind_start:ind_end) ! 1 SM read(unitnum(k)) ! 2 ST read(unitnum(k)) ! 3 tau @@ -4741,178 +4748,178 @@ subroutine read_obs_SMOS( date_time, N_catd, this_obs_param, & read(unitnum(k)) ! 8 tau RSTD read(unitnum(k)) tmp_std(ind_start:ind_end) ! 9 std-dev SM read(unitnum(k)) tmp_cnt(ind_start:ind_end) ! 10 count SM - + else - + ! read SMOS Tb file - + ! first time obs are read: figure out angle and polarization of interest - + if (ind_start==1) then - + ! find the index for the angle of interest - ! NOTE: after processing of namelist inputs, each species + ! NOTE: after processing of namelist inputs, each species ! has a unique angle (see subroutine read_ens_upd_inputs()) - + ind_angle = -9999 - + do i=1,N_ang - + if (abs(tmp_ang(i)-this_obs_param%ang(1))<0.01) ind_angle = i - + end do - + if (ind_angle<0) then err_msg = 'Problem with incidence angle' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + ! need h-pol or v-pol? - + if (this_obs_param%pol==1) then - + hpol = .true. - + elseif (this_obs_param%pol==2) then - + hpol = .false. - + else - + err_msg = 'Problem with polarization' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + end if - + ! for each field, loop over all angles, keep data at angle of interest - + do i=1,N_ang ! (1) Tbh - + if ( (i==ind_angle) .and. ( hpol) ) then read(unitnum(k)) tmp_obs(ind_start:ind_end) else read(unitnum(k)) end if - + end do - + do i=1,N_ang ! (2) Tbv - + if ( (i==ind_angle) .and. (.not. hpol) ) then read(unitnum(k)) tmp_obs(ind_start:ind_end) else read(unitnum(k)) end if - + end do - + do i=1,N_ang ! (3) std-dev Tbh - + if ( (i==ind_angle) .and. ( hpol) ) then read(unitnum(k)) tmp_std(ind_start:ind_end) else read(unitnum(k)) end if - + end do - + do i=1,N_ang ! (4) std-dev Tbv - + if ( (i==ind_angle) .and. (.not. hpol) ) then read(unitnum(k)) tmp_std(ind_start:ind_end) else read(unitnum(k)) end if - + end do - + do i=1,N_ang ! (5) count Tbh - + if ( (i==ind_angle) .and. ( hpol) ) then read(unitnum(k)) tmp_cnt(ind_start:ind_end) else read(unitnum(k)) end if - + end do - + do i=1,N_ang ! (6) count Tbv - + if ( (i==ind_angle) .and. (.not. hpol) ) then read(unitnum(k)) tmp_cnt(ind_start:ind_end) else read(unitnum(k)) end if - + end do - + ! additional fields in file that are not currently read: ! ! (7-8) RA Tbh-Tbv ! (9-16) repeat the above for T3 and T4 - + end if ! if SM_files end if ! if N_obs_tmp(k)>0 - + end do ! loop through files - + ! ------------------------------------------------- ! ! eliminate no-data-values and data that fail initial QC ! and keep track how many obs survived from each file - N_obs_tmp = 0 ! re-use N_obs_tmp - + N_obs_tmp = 0 ! re-use N_obs_tmp + j=0 - + do n=1,N_obs - + if (SM_files) then - + keep_data = & (tmp_obs( n) > SM_min) .and. & ! incl: any neg is nodata (tmp_obs( n) < SM_max) .and. & (tmp_DQX( n) > SM_DQX_min) .and. & (tmp_DQX( n) < SM_DQX_max) .and. & (tmp_std( n) < SM_std_max) .and. & - (tmp_cnt( n) >= SM_cnt_min) - + (tmp_cnt( n) >= SM_cnt_min) + else - + keep_data = & (tmp_obs( n) > Tb_min) .and. & ! incl: any neg is nodata (tmp_obs( n) < Tb_max) .and. & (tmp_std( n) < Tb_std_max) .and. & - (tmp_cnt( n) >= Tb_cnt_min) - + (tmp_cnt( n) >= Tb_cnt_min) + end if - + if (keep_data) then - + j=j+1 - + tmp_obs(j) = tmp_obs(n) tmp_lon(j) = tmp_lon(n) tmp_lat(j) = tmp_lat(n) N_obs_tmp(tmp_file_ind(n)) = N_obs_tmp(tmp_file_ind(n)) + 1 - + end if - + end do - + N_obs = j ! Note: This is NOT the final number of valid obs in "SMOS_data"! - + if (SM_files) deallocate(tmp_DQX) - + deallocate(tmp_std) deallocate(tmp_cnt) - - deallocate(tmp_file_ind) - + + deallocate(tmp_file_ind) + ! ------------------------------------------------- ! ! map obs to tiles @@ -4920,24 +4927,24 @@ subroutine read_obs_SMOS( date_time, N_catd, this_obs_param, & ! for each observation ! a) determine grid cell that contains lat/lon ! b) determine tile within grid cell that contains lat/lon - + if (N_obs>0) then - + allocate(tmp_tile_num(N_obs)) - + ! temporarily shift lat/lon of obs for computation of nearest tile to ! avoid ambiguous assignment of M09 model tile within M36 obs grid cell - ! (center of M36 grid cell is equidistant from at least two M09 model + ! (center of M36 grid cell is equidistant from at least two M09 model ! tiles) -- reichle, 23 Aug 2013 - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_obs, tmp_lat, tmp_lon, & this_obs_param, & tmp_tile_num, & tmp_shift_lat, tmp_shift_lon ) - - ! make sure center-of-mass of tile that administers obs + + ! make sure center-of-mass of tile that administers obs ! is within EASEv2 M36 obs grid cell, discard obs otherwise ! (by setting tmp_tile_num to negative value) ! - reichle, 31 Jan 2014 @@ -4948,171 +4955,171 @@ subroutine read_obs_SMOS( date_time, N_catd, this_obs_param, & ! distortion of EASE grid cells at high latitudes, and/or Gabrielle's ! use of an M36 innovations integration to derive Tb scaling files ! for the M09 (SMAP) system. - ! The problem with the piece of code is that it throws out far too many - ! obs (at all latitudes) if the tile space is coarse, e.g., that of the + ! The problem with the piece of code is that it throws out far too many + ! obs (at all latitudes) if the tile space is coarse, e.g., that of the ! 1/2 deg Lat/Lon grids of MERRA or MERRA-2. ! The piece of code may no longer be needed because of improvements in ! the obs readers and in get_obs_pred(), but the impact of removing the ! code on the SMAP L4_SM system is not clear. At this time, just prior - ! to finalizing the L4_SM "validated release", keep the code for the - ! EASEv2 M09 and M36 tile spaces that are relevant for the L4_SM + ! to finalizing the L4_SM "validated release", keep the code for the + ! EASEv2 M09 and M36 tile spaces that are relevant for the L4_SM ! system, but drop it for all other tile spaces. ! - reichle, 3 Feb 2016 - + if ( & (index(tile_grid_d%gridtype, 'EASEv2_M09') /=0) .or. & (index(tile_grid_d%gridtype, 'EASEv2_M36') /=0) ) then - + do ii=1,N_obs - + if (tmp_tile_num(ii)>0) then - + call ease_convert('EASEv2_M36', & tile_coord(tmp_tile_num(ii))%com_lat, & tile_coord(tmp_tile_num(ii))%com_lon, & M36_col_ind_tile, M36_row_ind_tile ) - + call ease_convert('EASEv2_M36', & tmp_lat(ii), & tmp_lon(ii), & M36_col_ind_obs, M36_row_ind_obs ) - + if ( (nint(M36_col_ind_tile)/=nint(M36_col_ind_obs)) .or. & (nint(M36_row_ind_tile)/=nint(M36_row_ind_obs)) ) & tmp_tile_num(ii) = -9999 - + end if - + end do - + end if ! compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + SMOS_data = 0. SMOS_lon = 0. SMOS_lat = 0. N_obs_in_tile = 0 - + do i=1,N_obs - + ind_tile = tmp_tile_num(i) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind_tile>0) then ! this step eliminates obs outside domain - + SMOS_data(ind_tile) = SMOS_data(ind_tile) + tmp_obs(i) SMOS_lon( ind_tile) = SMOS_lon( ind_tile) + tmp_lon(i) SMOS_lat( ind_tile) = SMOS_lat( ind_tile) + tmp_lat(i) - + N_obs_in_tile(ind_tile) = N_obs_in_tile(ind_tile) + 1 - + end if - + end do - + ! normalize - + do i=1,N_catd - + if (N_obs_in_tile(i)>1) then - + tmpreal = real(N_obs_in_tile(i)) SMOS_data(i) = SMOS_data(i)/tmpreal SMOS_lon( i) = SMOS_lon( i)/tmpreal SMOS_lat( i) = SMOS_lat( i)/tmpreal - + elseif (N_obs_in_tile(i)==0) then - + SMOS_data(i) = this_obs_param%nodata SMOS_lon( i) = this_obs_param%nodata SMOS_lat( i) = this_obs_param%nodata - + end if - + end do - + ! clean up - + deallocate(tmp_tile_num) - + ! -------------------------------- - + ! set observation error standard deviation - + do i=1,N_catd std_SMOS_data(i) = this_obs_param%errstd end do - + ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if - + deallocate(tmp_ang) - + deallocate(tmp_lon) deallocate(tmp_lat) - + deallocate(tmp_obs) - + end if - + do k=1,N_files - + close(unitnum(k),status='keep') - + end do - + deallocate(unitnum) deallocate(start_time) - deallocate(end_time) - deallocate(Asc_flag) - deallocate(N_ang_tmp) - + deallocate(end_time) + deallocate(Asc_flag) + deallocate(N_ang_tmp) + end if ! if N_files>0 ! ------------------------------------------------- ! ! write "obslog" file - + if (write_obslog) then - + YYYYMMDD_HHMMSSz = date_time2string(date_time) - + tmpstr80 = 'read_obs_SMOS()' ! name of this subroutine - + do k=1,N_files - + write (tmpstr12,'(i12)') N_obs_tmp(k) ! convert integer to string - + call add_to_obslog( YYYYMMDD_HHMMSSz, this_obs_param%descr, tmpstr80, & tmpstr12, fnames(k) ) - + end do - + end if ! clean up - if (N_files>0) deallocate(N_obs_tmp) - + if (N_files>0) deallocate(N_obs_tmp) + deallocate(fnames) - + if (logit) write (logunit,*) 'read_obs_SMOS(): done.' - + end subroutine read_obs_SMOS @@ -5123,16 +5130,16 @@ subroutine read_obs_MODIS_SCF( & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, MODIS_obs, std_MODIS_obs, MODIS_lon, MODIS_lat ) - + ! read MODIS snow cover fraction observations on MODIS 0.05-degree climate ! modeling grid (CMG) ! ! Terra: MOD10C1 ! Aqua: MYD10C1 ! - ! For now, assume that MODIS resolution is finer than Catchment tile space + ! For now, assume that MODIS resolution is finer than Catchment tile space ! and super-ob data to tile space, with lat/lon coords of obs matching - ! lat/lon coords of tiles + ! lat/lon coords of tiles ! ! reichle, 18 Oct 2023 ! @@ -5143,20 +5150,20 @@ subroutine read_obs_MODIS_SCF( & ! inputs type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim ! [seconds] integer, intent(in) :: N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! output logical, intent(out) :: found_obs @@ -5164,27 +5171,27 @@ subroutine read_obs_MODIS_SCF( & real, dimension(N_catd), intent(out) :: MODIS_obs, std_MODIS_obs, MODIS_lon, MODIS_lat ! ------------------------------------------------------------------------ - + ! locals integer, parameter :: dtstep_assim_max = 21600 ! [seconds] avoid assim window spanning >=180 deg lon real, parameter :: CMG_dlon = 0.05 ! [degrees] longitude spacing of MODIS CMG grid real, parameter :: CMG_dlat = 0.05 ! [degrees] latitude spacing of MODIS CMG grid - + real, parameter :: CMG_ll_lon = -180. ! [degrees] lower-left longitude of MODIS CMG grid real, parameter :: CMG_ur_lat = 90. ! [degrees] upper-right latitude of MODIS CMG grid - + character(7) :: MODIS_product_ID - + character(4) :: YYYY - character(3) :: DDD ! day of year - + character(3) :: DDD ! day of year + character(400) :: fname real :: overpass_hour, tmp_delta, tmp_real, max_delta_lon - - integer :: N_files, N_lon, N_lat, N_tmp, nn, kk, ind + + integer :: N_files, N_lon, N_lat, N_tmp, nn, kk, ind integer :: N_CMG_obs, N_good_data, tmp_ind_start, tmp_ind_last type(date_time_type) :: date_time_beg, date_time_end @@ -5192,7 +5199,7 @@ subroutine read_obs_MODIS_SCF( & real :: lon_beg, lon_end real :: lon_beg_MODIS, lon_end_MODIS - + integer :: delta_day_beg, delta_day_end integer :: delta_day_beg_MODIS, delta_day_end_MODIS @@ -5210,80 +5217,80 @@ subroutine read_obs_MODIS_SCF( & character(len=*), parameter :: Iam = 'read_obs_MODIS_SCF' character(len=400) :: err_msg - + ! ---------------------------------------------------------------------------------- ! - ! restrict assimilation time step to max allowed + ! restrict assimilation time step to max allowed + + if (dtstep_assim > dtstep_assim_max) then - if (dtstep_assim > dtstep_assim_max) then - err_msg = 'dtstep_assim exceeds max allowed' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - - end if - + end if + + ! initialize - + found_obs = .false. - + ! identify MODIS product and overpass hour MODIS_product_ID = this_obs_param%name(1:7) - + select case (MODIS_product_ID) - + case('MOD10C1'); overpass_hour = 10.5 ! [hours] Terra: 10:30am local time - + case('MYD10C1'); overpass_hour = 13.5 ! [hours] Aqua: 1:30pm local time - + end select - - + + ! determine beginning and end of assimilation window - + date_time_beg = date_time date_time_end = date_time - + call augment_date_time( -dtstep_assim/2, date_time_beg ) call augment_date_time( dtstep_assim/2, date_time_end ) - - + + ! determine longitude band associated with assimilation window and local overpass hour ! ! observations will be returned only for tiles with lon_end < tile_coord%com_lon <= lon_beg - + call localtime2longitude( date_time_beg, overpass_hour, lon_beg, delta_day_beg ) call localtime2longitude( date_time_end, overpass_hour, lon_end, delta_day_end ) - - - ! determine (rough) longitude band for which MODIS obs need to be read + + + ! determine (rough) longitude band for which MODIS obs need to be read ! - ! --> because tiles have a non-zero extent, need to read MODIS obs in CMG grid cells located + ! --> because tiles have a non-zero extent, need to read MODIS obs in CMG grid cells located ! in a wider band (lon_min-delta:lon_max+delta), ! where delta should be somewhat larger than max( tile_coord%max_lon - tile_coord%min_lon ) - + tmp_delta = 3.*maxval( tile_coord(1:N_catd)%max_lon - tile_coord(1:N_catd)%min_lon ) ! [degrees] - + tmp_delta = tmp_delta/360.*86400. ! [seconds] - + date_time_beg_MODIS = date_time_beg date_time_end_MODIS = date_time_end call augment_date_time( -nint(tmp_delta), date_time_beg_MODIS ) call augment_date_time( nint(tmp_delta), date_time_end_MODIS ) - + call localtime2longitude( date_time_beg_MODIS, overpass_hour, lon_beg_MODIS, delta_day_beg_MODIS ) call localtime2longitude( date_time_end_MODIS, overpass_hour, lon_end_MODIS, delta_day_end_MODIS ) - + ! adjust date_time_*_MODIS to reflect calendar date at lon_*_MODIS - + call augment_date_time( delta_day_beg_MODIS*86400, date_time_beg_MODIS ) call augment_date_time( delta_day_end_MODIS*86400, date_time_end_MODIS ) - + ! put together arguments for call(s) to read_MODIS_SCF_hdf() - + lon_min_vec = MAPL_UNDEF lon_max_vec = MAPL_UNDEF @@ -5293,23 +5300,23 @@ subroutine read_obs_MODIS_SCF( & N_lon_vec = 0 if (lon_end_MODIS < lon_beg_MODIS) then - + if ( date_time_end_MODIS%dofyr == date_time_beg_MODIS%dofyr ) then - - ! need only one daily MODIS file and longitude band - - N_files = 1 - - lon_min_vec(1) = lon_end_MODIS + + ! need only one daily MODIS file and longitude band + + N_files = 1 + + lon_min_vec(1) = lon_end_MODIS lon_max_vec(1) = lon_beg_MODIS - + year_vec( 1) = date_time_beg_MODIS%year dofyr_vec( 1) = date_time_beg_MODIS%dofyr - + else - + ! this should never happen for dtstep_assim_max=21600 and overpass_hour=10:30am or 1:30pm - + write (logunit,*) 'overpass_hour = ', overpass_hour write (logunit,*) 'date_time = ', date_time write (logunit,*) 'date_time_beg = ', date_time_beg @@ -5320,47 +5327,47 @@ subroutine read_obs_MODIS_SCF( & write (logunit,*) 'lon_end = ', lon_end write (logunit,*) 'lon_beg_MODIS = ', lon_beg_MODIS write (logunit,*) 'lon_end_MODIS = ', lon_end_MODIS - + err_msg = 'encountered unexpected condition for longitude band!!!' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + else - - ! longitude band wraps around dateline, two daily MODIS files are needed + + ! longitude band wraps around dateline, two daily MODIS files are needed ! (this could also occur if lon_*_MODIS=180., which would result in an ! empty first longitude band, but because of tmp_delta>0, this should ! never happen) - + N_files = 2 - - lon_min_vec(1) = lon_end_MODIS + + lon_min_vec(1) = lon_end_MODIS lon_max_vec(1) = 179.999 ! use 179.999 such that zero-based last_ind<=N_lon-1 (see below) - + year_vec( 1) = date_time_end_MODIS%year dofyr_vec( 1) = date_time_end_MODIS%dofyr - - lon_min_vec(2) = -180. + + lon_min_vec(2) = -180. lon_max_vec(2) = lon_beg_MODIS - + year_vec( 2) = date_time_beg_MODIS%year dofyr_vec( 2) = date_time_beg_MODIS%dofyr - + end if - - ! verify that longitude bands do not exceed max expected expected width + + ! verify that longitude bands do not exceed max expected expected width ! (add 0.1 degree of wiggle room) - + max_delta_lon = real(dtstep_assim_max + 2*nint(tmp_delta))/86400.*360. + 0.1 ! [degree] - + do nn=1,N_files - + if ( lon_max_vec(nn) - lon_min_vec(nn) > max_delta_lon ) then - + err_msg = 'longitude band too wide' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if end do @@ -5371,23 +5378,23 @@ subroutine read_obs_MODIS_SCF( & lat_min = minval( tile_coord(1:N_catd)%min_lat ) lat_max = maxval( tile_coord(1:N_catd)%max_lat ) - + ! determine N_lat and N_lon[_vec] (# CMG grid cells in lat/lon bands ) start_ind(1) = (CMG_ur_lat - lat_max )/CMG_dlat last_ind(1) = (CMG_ur_lat - lat_min )/CMG_dlat - + N_lat = last_ind(1) - start_ind(1) + 1 - + start_ind = (lon_min_vec - CMG_ll_lon)/CMG_dlon last_ind = (lon_max_vec - CMG_ll_lon)/CMG_dlon - - N_lon_vec = last_ind - start_ind + 1 - + + N_lon_vec = last_ind - start_ind + 1 + N_lon = sum( N_lon_vec(1:N_files) ) - - + + ! ! dbg ! ! write (*,*) '###############################################################################' ! ! dbg ! ! write (*,*) Iam // '():' ! ! dbg ! ! write (*,*) 'lon_min_vec = ', lon_min_vec @@ -5400,57 +5407,57 @@ subroutine read_obs_MODIS_SCF( & ! ! dbg ! ! write (*,*) 'dofyr_vec = ', dofyr_vec ! ! dbg ! ! write (*,*) 'N_lon_vec = ', N_lon_vec ! ! dbg ! ! write (*,*) 'N_lat = ', N_lat - ! ! dbg ! ! write (*,*) '###############################################################################' + ! ! dbg ! ! write (*,*) '###############################################################################' ! allocate arrays for MODIS CMG data (max size that could possibly be needed for obs from both files) - + N_tmp = N_lon*N_lat allocate( CMG_obs(N_tmp) ) allocate( CMG_lon(N_tmp) ) allocate( CMG_lat(N_tmp) ) - - - ! read MODIS SCF obs + + + ! read MODIS SCF obs ! ! - (renamed) files currently located at /discover/nobackup/projects/S2SHMA/MODIS/MOD10C1_V61/ (2010-2022) ! ! - in ensupd nml file, specify the file "name" according to the following template: - ! + ! ! %name = 'MOD10C1.Ayyyyddd.061.hdf' ! ! 1 2 ! 123456789012345678901234 - ! + ! ! MOD10C1 = MODIS product name ! .A = "acquisition time" indicator - ! yyyyddd = placeholder for year/day-of-year + ! yyyyddd = placeholder for year/day-of-year ! .061 = version (Collection) indicator - ! .hdf = file name extension + ! .hdf = file name extension ! ! Assuming the MODIS file naming convention remains unchanged, the version can then ! be specified in the nml file. - + N_CMG_obs = 0 ! initialize counter for "good" obs returned by calls to read_MODIS_SCF_hdf() - + do nn=1,N_files ! loop through longitude bands - + ! determine MODIS file name(s) - + write (YYYY,'(i4.4)') year_vec( nn) write (DDD, '(i3.3)') dofyr_vec(nn) - - fname = & + + fname = & trim(this_obs_param%path) // '/' // YYYY // '/' // & this_obs_param%name(1:9) // YYYY // DDD // this_obs_param%name(17:24) - - ! determine sub-array of CMG_* - + + ! determine sub-array of CMG_* + tmp_ind_start = N_CMG_obs + 1 tmp_ind_last = N_CMG_obs + N_lon_vec(nn)*N_lat - + call read_MODIS_SCF_hdf( fname, & lon_min_vec(nn), lon_max_vec(nn), lat_min, lat_max, & N_good_data, & @@ -5459,114 +5466,114 @@ subroutine read_obs_MODIS_SCF( & CMG_obs(tmp_ind_start:tmp_ind_last) ) N_CMG_obs = N_CMG_obs + N_good_data - + end do - - ! return if no MODIS obs were found (found_obs=.false. per initialization above) + + ! return if no MODIS obs were found (found_obs=.false. per initialization above) if (N_CMG_obs==0) return - ! map to tile space - + ! map to tile space + allocate(tmp_tile_num(N_CMG_obs)) - + call get_tile_num_for_obs( N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, & tile_num_in_cell_ij, & N_CMG_obs, CMG_lat(1:N_CMG_obs), CMG_lon(1:N_CMG_obs), & this_obs_param, & tmp_tile_num ) - - + + std_MODIS_obs = this_obs_param%errstd - + MODIS_obs = 0. MODIS_lon = 0. MODIS_lat = 0. N_obs_in_tile = 0 - + do kk=1,N_CMG_obs - + ind = tmp_tile_num(kk) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind>0) then ! this condition eliminates obs outside domain - + MODIS_obs( ind) = MODIS_obs( ind) + CMG_obs(kk) MODIS_lon( ind) = MODIS_lon( ind) + CMG_lon(kk) MODIS_lat( ind) = MODIS_lat( ind) + CMG_lat(kk) N_obs_in_tile(ind) = N_obs_in_tile(ind) + 1 - + end if - + end do - + ! normalize - + do kk=1,N_catd if (N_obs_in_tile(kk)>0) then - + tmp_real = real(N_obs_in_tile(kk)) MODIS_obs(kk) = MODIS_obs(kk)/tmp_real MODIS_lon(kk) = MODIS_lon(kk)/tmp_real MODIS_lat(kk) = MODIS_lat(kk)/tmp_real - + else if (N_obs_in_tile(kk)==0) then - + MODIS_obs(kk) = this_obs_param%nodata MODIS_lon(kk) = this_obs_param%nodata MODIS_lat(kk) = this_obs_param%nodata - + end if end do - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + deallocate(tmp_tile_num) - - deallocate(CMG_obs) + + deallocate(CMG_obs) deallocate(CMG_lon) deallocate(CMG_lat) - ! to avoid double-counting of MODIS CMG obs, remove obs for tiles with center-of-mass longitude + ! to avoid double-counting of MODIS CMG obs, remove obs for tiles with center-of-mass longitude ! falling outside the longitude band associated with assimilation window - + if (lon_end < lon_beg) then - + ! need only one longitude band (keep obs when lon_end<=com_lon<=lon_beg) - + do kk=1,N_catd - + if ( (tile_coord(kk)%com_lonlon_beg) ) & MODIS_obs(kk) = this_obs_param%nodata - + end do - + else - + ! longitude band wraps around dateline (keep obs when -180<=com_lon<=lon_beg or lon_end<=com_lon<=180) - + do kk=1,N_catd - + if ( (tile_coord(kk)%com_lon>lon_beg) .and. (tile_coord(kk)%com_lon= 24.) ) then - + err_msg = 'input local_hour falls outside allowed range of 0:24' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + ! initialize - + delta_day = 0 - + ! determine fractional UTC hour and time difference with local_hour - + UTC_hour = ( date_time%hour*3600 + date_time%min*60 + date_time%sec )/3600. ! 0 <= UTC_hour < 24 - + time_diff = local_hour - UTC_hour - + ! enforce -12. < time_diff <= 12. and determine associated date difference, if any - + if (time_diff <= -12.) then - + delta_day = 1 - time_diff = time_diff + 24. - + time_diff = time_diff + 24. + elseif (time_diff > 12.) then - + delta_day = -1 - time_diff = time_diff - 24. - + time_diff = time_diff - 24. + end if - + ! determine longitude - - longitude = time_diff/24.*360. - + + longitude = time_diff/24.*360. + end subroutine localtime2longitude - + ! ***************************************************************** - + subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & N_good_data, CMG_lon, CMG_lat, CMG_SCF ) - - ! read snow cover area fraction (SCF) obs from daily MODIS Terra or Aqua M?D10C1, version 6.1 + + ! read snow cover area fraction (SCF) obs from daily MODIS Terra or Aqua M?D10C1, version 6.1 ! - Terra: https://nsidc.org/data/mod10c1/versions/61 ! - Aqua: https://nsidc.org/data/myd10c1/versions/61 ! - daily data with spatial resolution of 0.05 deg on MODIS climate modeling grid (CMG) ! - Terra: missing days 2016 d. 50-58 - ! - data are read for the requested lat/lon range + ! - data are read for the requested lat/lon range ! - apply QC ! - ! reichle, 20 Oct 2023 + ! reichle, 20 Oct 2023 ! ! ------------------------------------------------------------------------------------------------- - + implicit none - + character(*), intent(in) :: fname ! MODIS file name with full path - + real, intent(in) :: lon_min, lon_max ! -180 <= lon_* <= 180 real, intent(in) :: lat_min, lat_max ! -90 <= lat_* <= 90 - + integer, intent(out) :: N_good_data - + real, dimension(:), intent(out) :: CMG_lon, CMG_lat, CMG_SCF ! NOTE: 1-dim array on CMG grid - + ! ------------------------------------------------- ! local parameters - + ! ll/ur_lon/lat simply indicate the extent of the MODIS CMG grid ! ! index increases from (-180,90) to (180,-90) (southward and eastward) - + integer, parameter :: CMG_N_lon = 7200 integer, parameter :: CMG_N_lat = 3600 - - real, parameter :: CMG_ll_lon = -180. + + real, parameter :: CMG_ll_lon = -180. real, parameter :: CMG_ll_lat = -90. - + real, parameter :: CMG_ur_lon = 180. real, parameter :: CMG_ur_lat = 90. - + real, parameter :: CMG_dlon = 0.05 real, parameter :: CMG_dlat = 0.05 - integer, parameter :: N_fields = 3 + integer, parameter :: N_fields = 3 character(22), dimension(N_fields), parameter :: field_names = (/ & 'Day_CMG_Snow_Cover ', & ! 1 @@ -5699,20 +5706,20 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & ! 1234567890123456789012 ! 1 2 - integer, parameter :: SCF_nodata = -9999. - + integer, parameter :: SCF_nodata = -9999. + integer(KIND=2), parameter :: qc_snow_cover_max = 100 ! exclude lake ice, night, inland water, ocean, etc - integer(KIND=2), parameter :: qc_clear_index_min = 20 ! ensure sufficiently clear conditions - integer(KIND=2), parameter :: qc_snow_spatial_max = 2 ! data quality (0=best, 1=good, 2=OK, 3=poor, 4=other) - + integer(KIND=2), parameter :: qc_clear_index_min = 20 ! ensure sufficiently clear conditions + integer(KIND=2), parameter :: qc_snow_spatial_max = 2 ! data quality (0=best, 1=good, 2=OK, 3=poor, 4=other) + integer, parameter :: DFACC_READ = 1 ! from hdf.inc - - integer, parameter :: uint8_min = 0 + + integer, parameter :: uint8_min = 0 integer, parameter :: uint8_max = 255 ! local variables - + integer :: N_lon, N_lat, N_tmp, ii, jj, kk, nn real, dimension(:), allocatable :: lon_c @@ -5722,14 +5729,14 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & real, dimension(:), allocatable :: lat_ind integer, dimension(2) :: start, edge, stride, last, dimsizes - + logical :: file_exists, keep_data - + integer :: status, sd_id, sds_id, sds_index - + integer :: sfstart, sfn2index, sfselect, sfginfo integer :: sfrdata, sfendacc, sfend - + character(64) :: sds_name integer :: rank, data_type, num_attrs @@ -5737,20 +5744,20 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & integer(KIND=2), dimension(:,:), allocatable :: Snow_Cover integer(KIND=2), dimension(:,:), allocatable :: Clear_Index integer(KIND=2), dimension(:,:), allocatable :: Snow_Spatial_QA - + character(1), dimension(:,:), allocatable :: tmp_char1 character(len=*), parameter :: Iam = 'read_MODIS_SCF_hdf' character(len=400) :: err_msg - + ! ------------------------------------------------------------------------- ! ! make sure file exists - + inquire( file=trim(fname), exist=file_exists ) - + if (.not. file_exists ) then - + if (logit) then write (logunit,'(400A)') trim(Iam), ': cannot find file ', trim(fname) write (logunit,* ) 'not reading MODIS SCF obs' @@ -5759,97 +5766,97 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & N_good_data = 0 return - + end if - + ! ensure lon_* and lat_* inputs are within range - + if ( (lon_min < CMG_ll_lon) .or. & - (lon_max > CMG_ur_lon) .or. & + (lon_max > CMG_ur_lon) .or. & (lat_min < CMG_ll_lat) .or. & (lat_max > CMG_ur_lat) ) then err_msg = 'lat/lon min/max inputs out of range' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + ! determine MODIS CMG array indices for requested lat/lon_min/max range - - ! MODIS CMG hdf files: + + ! MODIS CMG hdf files: ! - lon-by-lat (NOTE: The metadata "s.Vgroup(1).Vgroup(1).SDS(1).Dims.Size" returned - ! by Matlab's hdfinfo() confusingly suggests that CMG files are + ! by Matlab's hdfinfo() confusingly suggests that CMG files are ! lat-by-lon, which is not the case!) ! - index values increase eastward and southward start(1) = (lon_min - CMG_ll_lon)/CMG_dlon ! 0-based [as required for hdf read] start(2) = (CMG_ur_lat - lat_max )/CMG_dlat ! 0-based [as required for hdf read] - + last(1) = (lon_max - CMG_ll_lon)/CMG_dlon ! 0-based [as required for hdf read] last(2) = (CMG_ur_lat - lat_min )/CMG_dlat ! 0-based [as required for hdf read] - + N_lon = last(1) - start(1) + 1 - N_lat = last(2) - start(2) + 1 - + N_lat = last(2) - start(2) + 1 + edge(1) = N_lon edge(2) = N_lat - + stride(1) = 1 stride(2) = 1 - + ! ! dbg ! ! write (*,*) '###############################################################################' ! ! dbg ! ! write (*,*) Iam // '():' ! ! dbg ! ! write (*,*) 'size(CMG_SCF), N_lon, N_lat = ', size(CMG_SCF), N_lon, N_lat ! ! dbg ! ! write (*,*) 'lon_min, lon_max = ', lon_min, lon_max - ! ! dbg ! ! write (*,*) 'lat_min, lat_max = ', lat_min, lat_max + ! ! dbg ! ! write (*,*) 'lat_min, lat_max = ', lat_min, lat_max ! ! dbg ! ! write (*,*) 'start [lon, lat] = ', start ! ! dbg ! ! write (*,*) 'last [lon, lat] = ', last ! ! dbg ! ! write (*,*) 'edge [lon, lat] = ', edge - ! ! dbg ! ! write (*,*) '###############################################################################' - + ! ! dbg ! ! write (*,*) '###############################################################################' + ! checks array dimensions - + N_tmp = N_lon*N_lat - + if ( (N_tmp /= size(CMG_lon)) .or. & (N_tmp /= size(CMG_lat)) .or. & (N_tmp /= size(CMG_SCF)) ) then - + err_msg = 'inconsistent array dimensions' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + ! check array bounds - + if ( ( start(1) < 0 ) .or. ( start(1) > CMG_N_lon - 1 ) .or. & ( start(2) < 0 ) .or. ( start(2) > CMG_N_lat - 1 ) .or. & ( last( 1) < 0 ) .or. ( last( 1) > CMG_N_lon - 1 ) .or. & ( last( 2) < 0 ) .or. ( last( 2) > CMG_N_lat - 1 ) .or. & - ( start(1) > last(1) ) .or. ( start(2) > last(2) ) & - ) then - + ( start(1) > last(1) ) .or. ( start(2) > last(2) ) & + ) then + err_msg = 'start/edge indices out of bounds' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if if (logit) then write (logunit,'(400A)') trim(Iam), '(): reading MODIS SCF obs from ', trim(fname) write (logunit,* ) ' N_lon, N_lat, lon_min, lon_max, lat_min, lat_max = ' - write (logunit,* ) N_lon, N_lat, lon_min, lon_max, lat_min, lat_max + write (logunit,* ) N_lon, N_lat, lon_min, lon_max, lat_min, lat_max end if - + ! allocate arrays allocate(lon_c( N_lon)) allocate(lat_c( N_lat)) - + allocate(lon_ind(N_lon)) allocate(lat_ind(N_lat)) @@ -5860,47 +5867,47 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & allocate(tmp_char1( N_lon,N_lat)) ! -------------------------- - - ! determine center lat/lon of CMG cells + + ! determine center lat/lon of CMG cells lon_ind = (/(ii, ii=0, N_lon-1, 1)/) ! =0:(N_lon-1) lat_ind = (/(jj, jj=0, N_lat-1, 1)/) ! =0:(N_lat-1) - ! lat_c, lon_c are lat/lon at center of CMG grid cell - + ! lat_c, lon_c are lat/lon at center of CMG grid cell + lon_c = CMG_ll_lon + 0.5*CMG_dlon + (start(1)+lon_ind)*CMG_dlon lat_c = CMG_ur_lat - 0.5*CMG_dlat - (start(2)+lat_ind)*CMG_dlat - + ! -------------------------- - - ! open hdf file (read-only) and initialize SD interface - + + ! open hdf file (read-only) and initialize SD interface + sd_id = sfstart(trim(fname), DFACC_READ) if (sd_id<0) then - + err_msg = 'cannot sfstart (open) file: ' // trim(fname) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + ! read data do nn=1,N_fields - + sds_index = sfn2index( sd_id, trim(field_names(nn)) ) sds_id = sfselect( sd_id, sds_index ) status = sfginfo( sds_id, sds_name, rank, dimsizes, data_type, num_attrs ) - - if ( (dimsizes(1)/=CMG_N_lon) .or. (dimsizes(2)/=CMG_N_lat) ) then + + if ( (dimsizes(1)/=CMG_N_lon) .or. (dimsizes(2)/=CMG_N_lat) ) then err_msg = 'dimensions in hdf file doe not match CMG_N_lon and/or CMG_N_lat' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + ! ######################################################### ! ! dbg ! ! write (*,*) 'sds_name = ', trim(sds_name) @@ -5908,28 +5915,28 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & ! ! dbg ! ! write (*,*) 'dimsizes = ', dimsizes ! ######################################################### - + select case (nn) - + case (1); status = sfrdata( sds_id, start, stride, edge, tmp_char1 ); Snow_Cover = ichar(tmp_char1,2) case (2); status = sfrdata( sds_id, start, stride, edge, tmp_char1 ); Clear_Index = ichar(tmp_char1,2) case (3); status = sfrdata( sds_id, start, stride, edge, tmp_char1 ); Snow_Spatial_QA = ichar(tmp_char1,2) - + case default; call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown field') - + end select - + if (status/=0) then - + err_msg = 'error reading data from hdf file: ' // trim(fname) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if status = sfendacc(sds_id) ! terminate access to SDS (field) - + end do - + status = sfend(sd_id) ! close hdf file and SD interface @@ -5940,10 +5947,10 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & ! ! dbg ! ! write (*,*) Clear_Index(1:10,1:5) ! ! dbg ! ! write (*,*) 'Snow_Spatial_QA(1:10,1:5)=' ! ! dbg ! ! write (*,*) Snow_Spatial_QA(1:10,1:5) - + ! ! dbg ! ! !if ( (lon_min>-138) .and. (lon_min<-134) ) then ! ! dbg ! ! if (.false.) then - ! ! dbg ! ! do jj=1,N_lat + ! ! dbg ! ! do jj=1,N_lat ! ! dbg ! ! write(997) Snow_Cover( jj,:) ! ! dbg ! ! write(998) Clear_Index( jj,:) ! ! dbg ! ! write(999) Snow_Spatial_QA(jj,:) @@ -5963,16 +5970,16 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & any(Clear_Index > uint8_max) .or. & any(Snow_Spatial_QA < uint8_min) .or. & any(Snow_Spatial_QA > uint8_max) ) then - + err_msg = 'unexpected range in data from hdf file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + ! ------------------------------------- ! ! apply QC and put SCF obs into output array - + CMG_lon = SCF_nodata ! initialize CMG_lat = SCF_nodata ! initialize CMG_SCF = SCF_nodata ! initialize @@ -5981,41 +5988,41 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & do ii=1,N_lon do jj=1,N_lat - - ! note: Snow_Cover >= 0 per range check above, no need to check for minimum - - keep_data = & - (Snow_Cover( ii,jj) <= qc_snow_cover_max ) .and. & ! 0<=SCF<=100 (1) + + ! note: Snow_Cover >= 0 per range check above, no need to check for minimum + + keep_data = & + (Snow_Cover( ii,jj) <= qc_snow_cover_max ) .and. & ! 0<=SCF<=100 (1) (Clear_Index( ii,jj) > qc_clear_index_min ) .and. & ! sufficiently clear sky (2) (Snow_Spatial_QA(ii,jj) <= qc_snow_spatial_max) ! keep "best", "good", or "OK" quality (3) - + ! (1) excludes "lake ice", "night", "inland water", "ocean", "cloud obscured water", "data not mapped", "fill" ! (2) clear_index>100 already removed via qc_snow_cover_max ! (3) excludes Antarctica - + if (keep_data) then kk = kk + 1 - + ! raw SCF value is for clear portion of grid cell only, need to normalize with Clear_Index - - CMG_SCF(kk) = real(Snow_Cover(ii,jj))/real(Clear_Index(ii,jj)) - + + CMG_SCF(kk) = real(Snow_Cover(ii,jj))/real(Clear_Index(ii,jj)) + CMG_lon(kk) = lon_c(ii) CMG_lat(kk) = lat_c(jj) - + end if - + end do end do - + N_good_data = kk - + if (logit) write (logunit,*) 'N_good_data = ', N_good_data - + deallocate(lat_c) deallocate(lon_c) - + deallocate(lat_ind) deallocate(lon_ind) @@ -6026,43 +6033,43 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & end subroutine read_MODIS_SCF_hdf ! ***************************************************************** - + subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & dtstep_assim, tile_coord, tile_grid_d, & N_tile_in_cell_ij, tile_num_in_cell_ij, write_obslog, & found_obs, SMAP_data, std_SMAP_data, SMAP_lon, SMAP_lat, SMAP_time ) - - ! read freeze/thaw (FT) data within the assimilation window from one or more + + ! read freeze/thaw (FT) data within the assimilation window from one or more ! SMAP L2_SM_AP half-orbit h5 files (or L3_FT_A files -- TO BE IMPLEMENTED) ! - ! this subroutine reads each species independently of the others; + ! this subroutine reads each species independently of the others; ! TO DO: what if more than one flavor of FT data is read from SMAP? ! ! reichle, 14 Nov 2014 - + implicit none - + ! inputs: - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: N_catd, dtstep_assim type(obs_param_type), intent(in) :: this_obs_param type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input logical, intent(in) :: write_obslog ! outputs: - + logical, intent(out) :: found_obs real, intent(out), dimension(N_catd) :: SMAP_data, std_SMAP_data @@ -6072,43 +6079,43 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & ! -------------------------------------- ! local variables - + character(len=*), parameter :: Iam = 'read_obs_SMAP_FT' character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 ! SMAP data are stored in Yyyyy/Mmm/Ddd directories ! - ! Within each directory, there must be an ASCII text file that lists the h5 file + ! Within each directory, there must be an ASCII text file that lists the h5 file ! names of all files in the directory that are suitable for assimilation. ! These ASCII files must be named as follows: - + character(80), parameter :: fname_of_fname_list_L2_SM_AP_A = 'SMAP_L2_SM_AP_A_list.txt' character(80), parameter :: fname_of_fname_list_L2_SM_AP_D = 'SMAP_L2_SM_AP_D_list.txt' character(80), parameter :: fname_of_fname_list_L3_FT_A = 'SMAP_L3_FT_A_list.txt' - + logical, parameter :: tmp_debug = .false. - + real, parameter :: FT_min = 0. ! min allowed FT real, parameter :: FT_max = 1. ! max allowed FT integer, parameter :: dt_halforbit = 50*60 ! seconds - + integer, parameter :: N_halforbits_max = 15 ! max number of half-orbits per day - + integer, parameter :: dtstep_assim_max = 10800 ! max allowed dtstep_assim [seconds] - + ! get max number of files to be read (use 2*dt_halforbit just to be safe) - + integer, parameter :: N_fnames_max = & ceiling(real(N_halforbits_max)/86400.*real(dtstep_assim_max+2*dt_halforbit)) - + ! -------------------------------------------- type(hdf5read) :: h5r - + logical :: L2AP_files, keep_data, file_exists - + type(date_time_type) :: date_time_low, date_time_upp type(date_time_type) :: date_time_low_fname, date_time_tmp @@ -6117,15 +6124,15 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & integer :: dset_rank integer :: ind_tile, ind_start, ind_end - real :: M09_col_ind_tile, M09_row_ind_tile + real :: M09_col_ind_tile, M09_row_ind_tile real :: M09_col_ind_obs, M09_row_ind_obs real :: tmpreal - + real*8 :: J2000_seconds_low, J2000_seconds_upp character( 12) :: tmpstr12 character( 15) :: SMAP_date_time - character( 16) :: YYYYMMDD_HHMMSSz + character( 16) :: YYYYMMDD_HHMMSSz character( 80) :: fname_of_fname_list, tmpstr80 character(300) :: fname_tmp, tmp_err_msg @@ -6139,7 +6146,7 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & integer, dimension(N_catd) :: N_obs_in_tile real, dimension(:), allocatable :: tmp_lon, tmp_lat - + real, dimension(:), allocatable :: tmp_ft real*8, dimension(:), allocatable :: tmp_time @@ -6151,100 +6158,100 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & character(len=400) :: err_msg ! ------------------------------------------------------------------- - + ! check inputs - + ! the subroutine makes sense only if dtstep_assim <= 3 hours ! ! (this avoids that more than 2 different Yyyyy/Mmm/Ddd directories are needed - ! and that the time mismatch between the observed Tb and the model forecast Tb + ! and that the time mismatch between the observed Tb and the model forecast Tb ! becomes excessive; in future, add time interpolation of forecast Tb) - + if (dtstep_assim > dtstep_assim_max) then err_msg = 'dtstep_assim must not exceed 3 hours' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + ! initialize - + found_obs = .false. - - ! read L2_SM_AP or L3_FT_A files? - + + ! read L2_SM_AP or L3_FT_A files? + if (index(this_obs_param%descr,'L2AP') /= 0) then - + L2AP_files = .true. ! read SMAP L2_SM_AP files - + elseif (index(this_obs_param%descr,'L3FT') /= 0) then - + L2AP_files = .false. ! read SMAP L3_FT_A files else - + err_msg = 'cannot interpret %descr' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - - + + ! determine the name of the file that contain the relevant list of file names - + if (L2AP_files .and. this_obs_param%orbit==1) then - + fname_of_fname_list = trim(fname_of_fname_list_L2_SM_AP_A) - + elseif (L2AP_files .and. this_obs_param%orbit==2) then - + fname_of_fname_list = trim(fname_of_fname_list_L2_SM_AP_D) - + else - + fname_of_fname_list = trim(fname_of_fname_list_L3_FT_A) - + end if - + ! --------------------------- ! ! define h5 data set names - + if (L2AP_files) then - + ! L2_SM_AP - + dset_name_lon = '/Soil_Moisture_Retrieval_Data/longitude' dset_name_lat = '/Soil_Moisture_Retrieval_Data/latitude' - + dset_name_time = '/Soil_Moisture_Retrieval_Data/spacecraft_overpass_time_seconds' dset_name_ft = '/Soil_Moisture_Retrieval_Data/freeze_thaw_fraction' - + dset_name_ft_qual_flag = '/Soil_Moisture_Retrieval_Data/NOT_YET_IMPLEMENTED' else - + ! L3FT err_msg = 'L3FT NOT YET IMPLEMENTED' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + ! --------------------------- ! ! Define search interval for obs ! ! [ date_time - dtstep_assim/2 + one_second, date_time + dtstep_assim/2 ] - + ! lower boundary - + date_time_low = date_time - + call augment_date_time( -dtstep_assim/2+1, date_time_low ) - + J2000_seconds_low = datetime_to_J2000seconds( date_time_low, J2000_epoch_id ) ! upper boundary - + date_time_upp = date_time call augment_date_time( dtstep_assim/2, date_time_upp ) @@ -6257,132 +6264,132 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & ! assimilation window [date_time_low,date_time_upp] date_time_low_fname = date_time_low - + call augment_date_time( -dt_halforbit, date_time_low_fname ) - + ! read file with list of SMAP file names for first day - + call read_obs_fnames( date_time_low_fname, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames, fname_list(1:N_halforbits_max) ) - + ! if needed, read file with list of SMAP file names for second day and add ! file names into "fname_list" - + if (date_time_low_fname%day /= date_time_upp%day) then - + call read_obs_fnames( date_time_upp, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames_tmp, fname_list((N_fnames+1):(N_fnames+N_halforbits_max)) ) - + N_fnames = N_fnames + N_fnames_tmp - + end if ! ------------------------------------------------------------------ - + ! extract names of files that could contain data within the assimilation ! window ! ! sample file names: ! ! Yyyyy/Mmm/Ddd/SMAP_L2_SM_AP_03073_D_20010730T193828_D04003_000.h5 - ! ||||||||||||||| + ! ||||||||||||||| ! counter: 1 2 3 4 5 6 ! 1234567890123456789012345678901234567890123456789012345678901234567890 if (L2AP_files) then - + ind_start = 37 ind_end = 51 - + else - + ind_start = -9999999 ! TO DO: IMPLEMENT F3FT ind_end = -9999999 - + end if kk = 0 - + do ii=1,N_fnames SMAP_date_time = fname_list(ii)(ind_start:ind_end) - + date_time_tmp = SMAPdatetime_to_DateTimeType( SMAP_date_time ) ! check whether: date_time_low_fname < date_time_tmp <= date_time_upp if ( datetime_lt_refdatetime( date_time_low_fname, date_time_tmp ) .and. & - datetime_le_refdatetime( date_time_tmp, date_time_upp ) ) then - + datetime_le_refdatetime( date_time_tmp, date_time_upp ) ) then + kk = kk+1 - ! there can be no more than N_fnames_max files that have data falling into the + ! there can be no more than N_fnames_max files that have data falling into the ! assimilation window (see also dtstep_assim_max) - + if (kk>N_fnames_max) then err_msg = 'too many files match assimilation window' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + fname_list(kk) = fname_list(ii) - + end if - + end do ! ii=1,N_fnames - + N_fnames = kk - + ! the "N_fnames" files of interest are now the first "N_fnames" entries ! in "fname_list" - + ! --------------------------------------------------------- ! ! read and process data if files were found - + if (N_fnames==0) then - + ! no data files found - + SMAP_data = this_obs_param%nodata SMAP_lon = this_obs_param%nodata SMAP_lat = this_obs_param%nodata SMAP_time = real(this_obs_param%nodata,kind(0.0D0)) std_SMAP_data = this_obs_param%nodata - + else - + ! initialize outputs - + SMAP_data = 0. SMAP_lon = 0. SMAP_lat = 0. SMAP_time = 0.0D0 N_obs_in_tile = 0 ! for normalization after mapping to tile and super-obs - + ! loop through files - + do kk=1,N_fnames - + ! open file - + fname_tmp = trim(this_obs_param%path) // '/' // fname_list(kk) - + if (logit) write(logunit,'(400A)') 'reading file: ' // trim(fname_tmp) - + inquire(file=fname_tmp, exist=file_exists) - + if (.not. file_exists) then - + err_msg = 'file does NOT exist' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if call h5r%openFile(fname_tmp) - + ! ------------------------------ ! read h5 datasets @@ -6390,65 +6397,65 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & tmp_err_msg = trim(Iam) // ': inconsistent dataset lengths' ! LONGITUDE: query dataset, record size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_lon) call h5r%queryDataset(dset_name_lon, dset_rank, dset_size) - + N_obs_tmp = dset_size(1) - + allocate(tmp_lon(N_obs_tmp)) - + call h5r%readDataset(tmp_lon) - + ! LATITUDE: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_lat) call h5r%queryDataset(dset_name_lat, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_lat(N_obs_tmp)) call h5r%readDataset(tmp_lat) - + ! TIME: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_time) - + call h5r%queryDataset(dset_name_time, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_time(N_obs_tmp)) call h5r%readDataset(tmp_time) ! FT: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_ft) - + call h5r%queryDataset(dset_name_ft, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_ft(N_obs_tmp)) - + call h5r%readDataset(tmp_ft) - + !! FT_QUAL_FLAG: query dataset, check size, allocate space, read data - ! + ! !if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_ft_qual_flag) ! !call h5r%queryDataset(dset_name_ft_qual_flag, dset_rank, dset_size) - ! + ! !if (N_obs_tmp/=dset_size(1)) then ! call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) !end if @@ -6457,62 +6464,62 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & ! !call h5r%readDataset(tmp_ft_qual_flag) - + ! close file - + call h5r%closeFile if (logit) write(logunit,*) 'done reading file' ! -------------------------------------------------------- ! - ! eliminate obs outside desired time window, no-data-values and data + ! eliminate obs outside desired time window, no-data-values and data ! that fail initial QC ! keep track of how many obs survived from current file ! ##################################### - ! TO DO: + ! TO DO: ! - use quality flag once available ! ##################################### - + jj = 0 - + if (L2AP_files) then - + do nn=1,N_obs_tmp - + !(mod(tmp_tb_qual_flag_1(nn),2)==0) .and. & ! lowest bit must be 0 - + keep_data = & (tmp_time(nn) > J2000_seconds_low) .and. & (tmp_time(nn) <= J2000_seconds_upp) .and. & (tmp_ft(nn) > FT_min) .and. & ! elim neg nodata (tmp_ft(nn) < FT_max) ! elim unphysically large value - + if (keep_data) then - + jj=jj+1 - + tmp_lon( jj) = tmp_lon( nn) tmp_lat( jj) = tmp_lat( nn) tmp_ft( jj) = tmp_ft( nn) tmp_time(jj) = tmp_time(nn) - + end if - + end do else - + ! L3FT - + err_msg = 'L3FT NOT YET IMPLEMENTED' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + end if ! (L2AP_files) - + N_obs_kept(kk) = jj - + ! ------------------------------------------------- ! ! map obs to tiles @@ -6520,11 +6527,11 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & ! for each observation ! a) determine grid cell that contains lat/lon ! b) determine tile within grid cell that contains lat/lon - + if (N_obs_kept(kk)>0) then - + allocate(tmp_tile_num(N_obs_kept(kk))) - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_obs_kept(kk), & @@ -6532,15 +6539,15 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & tmp_lon(1:N_obs_kept(kk)), & this_obs_param, & tmp_tile_num ) - - ! make sure center-of-mass of tile that administers obs + + ! make sure center-of-mass of tile that administers obs ! is within EASEv2 M09 obs grid cell, discard obs otherwise ! (by setting tmp_tile_num to negative value) ! - ! this step eliminates SMAP obs that fall outside of the GEOS-5 + ! this step eliminates SMAP obs that fall outside of the GEOS-5 ! land mask at M09 resolution ! - ! not sure what this does if model is run on tile space other than + ! not sure what this does if model is run on tile space other than ! EASEv2_M09 - reichle, 11 Nov 2014 ! ! It is not 100 percent clear why this piece code had been added @@ -6549,125 +6556,125 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & ! distortion of EASE grid cells at high latitudes, and/or Gabrielle's ! use of an M36 innovations integration to derive Tb scaling files ! for the M09 (SMAP) system. - ! The problem with the piece of code is that it throws out far too many - ! obs (at all latitudes) if the tile space is coarse, e.g., that of the + ! The problem with the piece of code is that it throws out far too many + ! obs (at all latitudes) if the tile space is coarse, e.g., that of the ! 1/2 deg Lat/Lon grids of MERRA or MERRA-2. ! The piece of code may no longer be needed because of improvements in ! the obs readers and in get_obs_pred(), but the impact of removing the ! code on the SMAP L4_SM system is not clear. At this time, just prior - ! to finalizing the L4_SM "validated release", keep the code for the - ! EASEv2 M09 and M36 tile spaces that are relevant for the L4_SM + ! to finalizing the L4_SM "validated release", keep the code for the + ! EASEv2 M09 and M36 tile spaces that are relevant for the L4_SM ! system, but drop it for all other tile spaces. ! - reichle, 3 Feb 2016 - + if ( & (index(tile_grid_d%gridtype, 'EASEv2_M09') /=0) .or. & (index(tile_grid_d%gridtype, 'EASEv2_M36') /=0) ) then - + do ii=1,N_obs_kept(kk) - + if (tmp_tile_num(ii)>0) then - + call ease_convert('EASEv2_M09', & tile_coord(tmp_tile_num(ii))%com_lat, & tile_coord(tmp_tile_num(ii))%com_lon, & M09_col_ind_tile, M09_row_ind_tile ) - + call ease_convert('EASEv2_M09', & tmp_lat(ii), & tmp_lon(ii), & M09_col_ind_obs, M09_row_ind_obs ) - + if ( (nint(M09_col_ind_tile)/=nint(M09_col_ind_obs)) .or. & (nint(M09_row_ind_tile)/=nint(M09_row_ind_obs)) ) & tmp_tile_num(ii) = -9999 - + end if - + end do end if ! compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + do ii=1,N_obs_kept(kk) - + ind_tile = tmp_tile_num(ii) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind_tile>0) then ! this step eliminates obs outside domain - + SMAP_data(ind_tile) = SMAP_data(ind_tile) + tmp_ft( ii) SMAP_lon( ind_tile) = SMAP_lon( ind_tile) + tmp_lon( ii) SMAP_lat( ind_tile) = SMAP_lat( ind_tile) + tmp_lat( ii) SMAP_time(ind_tile) = SMAP_time(ind_tile) + tmp_time(ii) - + N_obs_in_tile(ind_tile) = N_obs_in_tile(ind_tile) + 1 - + end if - + end do - + deallocate(tmp_tile_num) - + end if - + ! clean up - + if (allocated(tmp_lon )) deallocate(tmp_lon ) if (allocated(tmp_lat )) deallocate(tmp_lat ) if (allocated(tmp_time )) deallocate(tmp_time ) if (allocated(tmp_ft )) deallocate(tmp_ft ) if (allocated(tmp_ft_qual_flag)) deallocate(tmp_ft_qual_flag) - + end do ! kk=1,N_fnames - + ! normalize - + do ii=1,N_catd - + if (N_obs_in_tile(ii)>1) then - + tmpreal = real(N_obs_in_tile(ii)) - + SMAP_data(ii) = SMAP_data(ii)/ tmpreal SMAP_lon( ii) = SMAP_lon( ii)/ tmpreal SMAP_lat( ii) = SMAP_lat( ii)/ tmpreal SMAP_time(ii) = SMAP_time(ii)/real(tmpreal,kind(0.0D0)) - + elseif (N_obs_in_tile(ii)==0) then - + SMAP_data(ii) = this_obs_param%nodata SMAP_lon( ii) = this_obs_param%nodata SMAP_lat( ii) = this_obs_param%nodata SMAP_time(ii) = real(this_obs_param%nodata,kind(0.0D0)) - + end if - + end do - + ! -------------------------------- - + ! set observation error standard deviation - + do ii=1,N_catd std_SMAP_data(ii) = this_obs_param%errstd end do - + ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if ! N_fnames==0 - + ! ------------------------------------------------- ! ! write "obslog" file @@ -6677,39 +6684,39 @@ subroutine read_obs_SMAP_FT( date_time, N_catd, this_obs_param, & YYYYMMDD_HHMMSSz = date_time2string(date_time) tmpstr80 = 'read_obs_SMAP_FT()' ! name of this subroutine - + do kk=1,N_fnames - + fname_tmp = trim(this_obs_param%path) // fname_list(kk) - + write (tmpstr12,'(i12)') N_obs_kept(kk) ! convert integer to string - + call add_to_obslog( YYYYMMDD_HHMMSSz, this_obs_param%descr, tmpstr80, & tmpstr12, fname_tmp ) - + end do - + end if - + ! clean up - + if (logit) write (logunit,*) 'read_obs_SMAP_FT(): done.' - + end subroutine read_obs_SMAP_FT - + ! ***************************************************************** ! ***************************************************************** - + subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & dtstep_assim, tile_coord, tile_grid_d, & N_tile_in_cell_ij, tile_num_in_cell_ij, write_obslog, & found_obs, SMAP_data, std_SMAP_data, SMAP_lon, SMAP_lat, SMAP_time ) - - ! read brightness temperature data within the assimilation window from one or more + + ! read brightness temperature data within the assimilation window from one or more ! SMAP half-orbit h5 files (L1C, L1CE (Enhanced), or L2AP) ! - ! this subroutine reads each species independently of the others; + ! this subroutine reads each species independently of the others; ! see subroutine turn_off_assim_SMAP_L1CTb() for avoiding the assimilation ! of redundant L1C_Tb observations when corresponding L2AP_Tb obs are assimilated ! @@ -6721,14 +6728,14 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! reichle, 31 May 2016: added stats check for L1C fore-minus-aft Tb differences ! reichle, 26 Dec 2017: added functionality for L1CE (Enhanced) files, incl. thinning ! reichle, 23 Jan 2018: removed stats check for L1C fore-minus-aft Tb differences; - ! use avg fore/aft timestamp so that fore and aft Tbs for same + ! use avg fore/aft timestamp so that fore and aft Tbs for same ! location are never used in different assimilation windows ! reichle, 22 Apr 2020: resurrected check for L1C fore-minus-aft Tb differences ! after antenna-scan-angle (ASA) issues continued and the ! SMAP Project declined to address these issues in L1 ops implicit none - + ! inputs: type(date_time_type), intent(in) :: date_time @@ -6738,18 +6745,18 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & type(obs_param_type), intent(in) :: this_obs_param type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input logical, intent(in) :: write_obslog ! outputs: - + logical, intent(out) :: found_obs real, intent(out), dimension(N_catd) :: SMAP_data, std_SMAP_data @@ -6759,49 +6766,49 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! -------------------------------------- ! local variables - + character(4), parameter :: J2000_epoch_id = 'TT12' ! see date_time_util.F90 - + ! SMAP data are stored in Yyyyy/Mmm/Ddd directories ! - ! Within each directory, there must be an ASCII text file that lists the h5 file + ! Within each directory, there must be an ASCII text file that lists the h5 file ! names of all files in the directory that are suitable for assimilation. ! These ASCII files must be named as follows: - character(80), parameter :: fname_of_fname_list_L1C_TB_A = 'SMAP_L1C_TB_A_list.txt' - character(80), parameter :: fname_of_fname_list_L1C_TB_D = 'SMAP_L1C_TB_D_list.txt' + character(80), parameter :: fname_of_fname_list_L1C_TB_A = 'SMAP_L1C_TB_A_list.txt' + character(80), parameter :: fname_of_fname_list_L1C_TB_D = 'SMAP_L1C_TB_D_list.txt' + + character(80), parameter :: fname_of_fname_list_L1C_TB_E_A = 'SMAP_L1C_TB_E_A_list.txt' + character(80), parameter :: fname_of_fname_list_L1C_TB_E_D = 'SMAP_L1C_TB_E_D_list.txt' - character(80), parameter :: fname_of_fname_list_L1C_TB_E_A = 'SMAP_L1C_TB_E_A_list.txt' - character(80), parameter :: fname_of_fname_list_L1C_TB_E_D = 'SMAP_L1C_TB_E_D_list.txt' - character(80), parameter :: fname_of_fname_list_L2_SM_AP_A = 'SMAP_L2_SM_AP_A_list.txt' character(80), parameter :: fname_of_fname_list_L2_SM_AP_D = 'SMAP_L2_SM_AP_D_list.txt' - + logical, parameter :: tmp_debug = .false. - + real, parameter :: Tb_min = 100.0 ! min allowed Tb real, parameter :: Tb_max = 320.0 ! max allowed Tb real, parameter :: max_std_tb_fore_minus_aft = 20. ! max std-dev L1C[E] fore-minus-aft Tb diffs integer, parameter :: L1CE_spacing = 3 ! thinning of L1C_TB_E in units of 9-km indices ("3" => 27 km) - + ! temporarily shift lat/lon of obs for computation of nearest tile to ! avoid ambiguous assignment of M09 model tile within M36 obs grid cell - ! (center of M36 grid cell is equidistant from at least two M09 model + ! (center of M36 grid cell is equidistant from at least two M09 model ! tiles) -- reichle, 23 Aug 2013 - + real, parameter :: tmp_shift_lon = 0.01 real, parameter :: tmp_shift_lat = 0.005 integer, parameter :: dt_halforbit = 50*60 ! seconds - + integer, parameter :: N_halforbits_max = 15 ! max number of half-orbits per day integer, parameter :: dtstep_assim_max = 10800 ! max allowed dtstep_assim [seconds] - + ! get max number of files to be read (use 2*dt_halforbit just to be safe) - + integer, parameter :: N_fnames_max = & ceiling(real(N_halforbits_max)/86400.*real(dtstep_assim_max+2*dt_halforbit)) @@ -6820,7 +6827,7 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & integer :: dset_rank integer :: ind_tile, ind_start, ind_end - real :: M36_col_ind_tile, M36_row_ind_tile + real :: M36_col_ind_tile, M36_row_ind_tile real :: M36_col_ind_obs, M36_row_ind_obs real :: tmpreal real :: tmp_tb_diff, tmp_tb_diff_Sum, tmp_tb_diff_SumOfSq, tmp_var @@ -6830,7 +6837,7 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & character( 2) :: orbit_tag character( 12) :: tmpstr12 character( 15) :: SMAP_date_time - character( 16) :: YYYYMMDD_HHMMSSz + character( 16) :: YYYYMMDD_HHMMSSz character( 80) :: fname_of_fname_list, tmpstr80 character(300) :: fname_tmp, tmp_err_msg @@ -6848,62 +6855,62 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & real, dimension(:), allocatable :: tmp_lon, tmp_lat integer, dimension(:), allocatable :: tmp_col, tmp_row - + real, dimension(:), allocatable :: tmp_tb_1 real, dimension(:), allocatable :: tmp_tb_2 - + real*8, dimension(:), allocatable :: tmp_time_1 real*8, dimension(:), allocatable :: tmp_time_2 integer, dimension(:), allocatable :: tmp_tb_qual_flag_1 integer, dimension(:), allocatable :: tmp_tb_qual_flag_2 - + integer, dimension(:), allocatable :: tmp_tile_num character(len=*), parameter :: Iam = 'read_obs_SMAP_halforbit_Tb' character(len=400) :: err_msg ! ------------------------------------------------------------------- - + ! check inputs - + ! the subroutine makes sense only if dtstep_assim <= 3 hours ! ! (this avoids that more than 2 different Yyyyy/Mmm/Ddd directories are needed - ! and that the time mismatch between the observed Tb and the model forecast Tb + ! and that the time mismatch between the observed Tb and the model forecast Tb ! becomes excessive; in future, add time interpolation of forecast Tb) - + if (dtstep_assim > dtstep_assim_max) then err_msg = 'dtstep_assim must not exceed 3 hours' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + ! initialize - + found_obs = .false. - - ! read Tbs from L1C_TB, L1C_TB_E, or L2_SM_AP files? + + ! read Tbs from L1C_TB, L1C_TB_E, or L2_SM_AP files? L1CE_thinning = .false. ! initialize if (index(this_obs_param%descr,'L1C') /= 0) then - - if (index(this_obs_param%descr,'_E') /= 0) then - + + if (index(this_obs_param%descr,'_E') /= 0) then + L1CE_files = .true. ! read SMAP L1C_TB_E (Enhanced) files L1C_files = .false. - - if (index(this_obs_param%descr,'_E27') /= 0) L1CE_thinning = .true. - + + if (index(this_obs_param%descr,'_E27') /= 0) L1CE_thinning = .true. + else - + L1C_files = .true. ! read SMAP L1C_TB (standard) files - L1CE_files = .false. - + L1CE_files = .false. + end if - + elseif (index(this_obs_param%descr,'L2AP') /= 0) then - + L1C_files = .false. ! read SMAP L2_SM_AP files L1CE_files = .false. @@ -6913,101 +6920,101 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - - - ! read ascending or descending files? - + + + ! read ascending or descending files? + tmp_err_msg = 'inconsistent %descr and %orbit' - + if (index(this_obs_param%descr,'_A') /=0 ) then - + orbit_tag = '_A' - + if (this_obs_param%orbit/=1) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + elseif (index(this_obs_param%descr,'_D') /=0 ) then - + orbit_tag = '_D' - + if (this_obs_param%orbit/=2) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + else - + err_msg = 'unknown %descr or %orbit' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if ! determine the name of the file that contains the relevant list of file names - + if (this_obs_param%orbit==1) then ! ascending - + if (L1C_files) then - + fname_of_fname_list = trim(fname_of_fname_list_L1C_TB_A) - + elseif (L1CE_files) then - + fname_of_fname_list = trim(fname_of_fname_list_L1C_TB_E_A) - + else - + fname_of_fname_list = trim(fname_of_fname_list_L2_SM_AP_A) - + end if - + elseif (this_obs_param%orbit==2) then ! descending - + if (L1C_files) then - + fname_of_fname_list = trim(fname_of_fname_list_L1C_TB_D) - + elseif (L1CE_files) then - + fname_of_fname_list = trim(fname_of_fname_list_L1C_TB_E_D) - + else - + fname_of_fname_list = trim(fname_of_fname_list_L2_SM_AP_D) - + end if - + else - + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown obs type') end if - + ! need h-pol or v-pol? - + if (this_obs_param%pol==1) then - + hpol = .true. - + elseif (this_obs_param%pol==2) then - + hpol = .false. - + else - + err_msg = 'Problem with polarization' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if ! --------------------------- ! ! define h5 data set names - + if (L1C_files .or. L1CE_files) then - + ! L1C_TB or L1_C_TB_E - + dset_name_lon = '/Global_Projection/cell_lon' dset_name_lat = '/Global_Projection/cell_lat' @@ -7018,53 +7025,53 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & dset_name_time_2 = '/Global_Projection/cell_tb_time_seconds_aft' if (hpol) then - + dset_name_tb_1 = '/Global_Projection/cell_tb_h_fore' dset_name_tb_2 = '/Global_Projection/cell_tb_h_aft' - + dset_name_tb_qual_flag_1 = '/Global_Projection/cell_tb_qual_flag_h_fore' dset_name_tb_qual_flag_2 = '/Global_Projection/cell_tb_qual_flag_h_aft' - + else - + dset_name_tb_1 = '/Global_Projection/cell_tb_v_fore' dset_name_tb_2 = '/Global_Projection/cell_tb_v_aft' - + dset_name_tb_qual_flag_1 = '/Global_Projection/cell_tb_qual_flag_v_fore' dset_name_tb_qual_flag_2 = '/Global_Projection/cell_tb_qual_flag_v_aft' - + end if - - else - + + else + ! L2_SM_AP - + dset_name_lon = '/Soil_Moisture_Retrieval_Data/longitude' dset_name_lat = '/Soil_Moisture_Retrieval_Data/latitude' - + dset_name_time_1 = '/Soil_Moisture_Retrieval_Data/spacecraft_overpass_time_seconds' dset_name_time_2 = '' ! *not* used if (hpol) then - + dset_name_tb_1 = '/Soil_Moisture_Retrieval_Data/tb_h_disaggregated' dset_name_tb_2 = '' ! *not* used - + dset_name_tb_qual_flag_1 = '/Soil_Moisture_Retrieval_Data/tb_h_disaggregated_qual_flag' dset_name_tb_qual_flag_2 = '' ! *not* used - + else - + dset_name_tb_1 = '/Soil_Moisture_Retrieval_Data/tb_v_disaggregated' dset_name_tb_2 = '' ! *not* used - + dset_name_tb_qual_flag_1 = '/Soil_Moisture_Retrieval_Data/tb_v_disaggregated_qual_flag' dset_name_tb_qual_flag_2 = '' ! *not* used - + end if - + end if - + ! --------------------------- ! ! Define search interval for obs @@ -7072,15 +7079,15 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! [ date_time - dtstep_assim/2 + one_second, date_time + dtstep_assim/2 ] ! lower boundary - + date_time_low = date_time - + call augment_date_time( -dtstep_assim/2+1, date_time_low ) J2000_seconds_low = datetime_to_J2000seconds( date_time_low, J2000_epoch_id ) ! upper boundary - + date_time_upp = date_time call augment_date_time( dtstep_assim/2, date_time_upp ) @@ -7093,72 +7100,72 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! assimilation window [date_time_low,date_time_upp] date_time_low_fname = date_time_low - + call augment_date_time( -dt_halforbit, date_time_low_fname ) - + ! read file with list of SMAP file names for first day - + call read_obs_fnames( date_time_low_fname, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames, fname_list(1:N_halforbits_max) ) - + ! if needed, read file with list of SMAP file names for second day and add ! file names into "fname_list" - + if (date_time_low_fname%day /= date_time_upp%day) then - + call read_obs_fnames( date_time_upp, this_obs_param, & fname_of_fname_list, N_halforbits_max, & N_fnames_tmp, fname_list((N_fnames+1):(N_fnames+N_halforbits_max)) ) - + N_fnames = N_fnames + N_fnames_tmp - + end if ! ------------------------------------------------------------------ - + ! extract names of files that could contain data within the assimilation ! window ! ! sample file names: - ! ||||||||||||||| + ! ||||||||||||||| ! Yyyyy/Mmm/Ddd/SMAP_L1C_TB_03027_D_20010727T160914_D04003_000.h5 ! - ! ||||||||||||||| + ! ||||||||||||||| ! Yyyyy/Mmm/Ddd/SMAP_L1C_TB_E_03027_D_20010727T160914_D04003_000.h5 ! Yyyyy/Mmm/Ddd/SMAP_L2_SM_AP_03073_D_20010730T193828_D04003_000.h5 - ! ||||||||||||||| + ! ||||||||||||||| ! counter: 1 2 3 4 5 6 ! 1234567890123456789012345678901234567890123456789012345678901234567890 if (L1C_files) then - + ind_start = 35 ind_end = 49 else ! L1CE or L2AP files - + ind_start = 37 ind_end = 51 - + end if kk = 0 - + do ii=1,N_fnames SMAP_date_time = fname_list(ii)(ind_start:ind_end) - + date_time_tmp = SMAPdatetime_to_DateTimeType( SMAP_date_time ) ! check whether: date_time_low_fname < date_time_tmp <= date_time_upp if ( datetime_lt_refdatetime( date_time_low_fname, date_time_tmp ) .and. & - datetime_le_refdatetime( date_time_tmp, date_time_upp ) ) then - + datetime_le_refdatetime( date_time_tmp, date_time_upp ) ) then + kk = kk+1 - ! there can be no more than N_fnames_max files that have data falling into the + ! there can be no more than N_fnames_max files that have data falling into the ! assimilation window (see also dtstep_assim_max) if (kk>N_fnames_max) then @@ -7167,62 +7174,62 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & end if fname_list(kk) = fname_list(ii) - + end if - + end do ! ii=1,N_fnames - + N_fnames = kk ! the "N_fnames" files of interest are now the first "N_fnames" entries ! in "fname_list" - + ! --------------------------------------------------------- ! ! read and process data if files were found if (N_fnames==0) then - + ! no data files found - + SMAP_data = this_obs_param%nodata SMAP_lon = this_obs_param%nodata SMAP_lat = this_obs_param%nodata SMAP_time = real(this_obs_param%nodata,kind(0.0D0)) std_SMAP_data = this_obs_param%nodata - + else - + ! initialize outputs - + SMAP_data = 0. SMAP_lon = 0. SMAP_lat = 0. SMAP_time = 0.0D0 N_obs_in_tile = 0 ! for normalization after mapping to tile and super-obs - + ! loop through files - + do kk=1,N_fnames - + ! open file - + fname_tmp = trim(this_obs_param%path) // '/' // fname_list(kk) - + if (logit) write(logunit,'(400A)') 'reading file: ' // trim(fname_tmp) - + inquire(file=fname_tmp, exist=file_exists) - + if (.not. file_exists) then - + err_msg = 'file does NOT exist' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if call h5r%openFile(fname_tmp) - + ! ------------------------------ ! read h5 datasets @@ -7230,27 +7237,27 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & tmp_err_msg = 'read_obs_SMAP_halforbit_Tb(): inconsistent dataset lengths' ! LONGITUDE: query dataset, record size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_lon) call h5r%queryDataset(dset_name_lon, dset_rank, dset_size) - + N_obs_tmp = dset_size(1) - + allocate(tmp_lon(N_obs_tmp)) - + call h5r%readDataset(tmp_lon) - + ! LATITUDE: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_lat) call h5r%queryDataset(dset_name_lat, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_lat(N_obs_tmp)) call h5r%readDataset(tmp_lat) @@ -7258,73 +7265,73 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & if (L1CE_thinning) then ! need to read column and row indices ! COLUMN: query dataset, record size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_col) - + call h5r%queryDataset(dset_name_col, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_col(N_obs_tmp)) - + call h5r%readDataset(tmp_col) - + ! ROW: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_row) - + call h5r%queryDataset(dset_name_row, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_row(N_obs_tmp)) - + call h5r%readDataset(tmp_row) end if - + ! TIME_1: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_time_1) - + call h5r%queryDataset(dset_name_time_1, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_time_1(N_obs_tmp)) call h5r%readDataset(tmp_time_1) ! TB_1: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_tb_1) - + call h5r%queryDataset(dset_name_tb_1, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_tb_1(N_obs_tmp)) call h5r%readDataset(tmp_tb_1) ! TB_QUAL_FLAG_1: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_tb_qual_flag_1) call h5r%queryDataset(dset_name_tb_qual_flag_1, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_tb_qual_flag_1(N_obs_tmp)) call h5r%readDataset(tmp_tb_qual_flag_1) @@ -7332,30 +7339,30 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! for L1C_TB or L1C_TB_E files also read "aft" if (L1C_files .or. L1CE_files) then - + ! *_1 = "fore" ! *_2 = "aft" ! TIME_2: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_time_2) call h5r%queryDataset(dset_name_time_2, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_time_2(N_obs_tmp)) call h5r%readDataset(tmp_time_2) ! TB_2: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_tb_2) - + call h5r%queryDataset(dset_name_tb_2, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if @@ -7363,43 +7370,43 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & allocate(tmp_tb_2(N_obs_tmp)) call h5r%readDataset(tmp_tb_2) - + ! TB_QUAL_FLAG_2: query dataset, check size, allocate space, read data - + if (tmp_debug .and. logit) write(logunit,*) trim(dset_name_tb_qual_flag_2) - + call h5r%queryDataset(dset_name_tb_qual_flag_2, dset_rank, dset_size) - + if (N_obs_tmp/=dset_size(1)) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, tmp_err_msg) end if - + allocate(tmp_tb_qual_flag_2(N_obs_tmp)) - + call h5r%readDataset(tmp_tb_qual_flag_2) end if - + ! close file - + call h5r%closeFile if (logit) write(logunit,*) 'done reading file' ! -------------------------------------------------------- ! - ! eliminate obs outside desired time window, no-data-values and data + ! eliminate obs outside desired time window, no-data-values and data ! that fail initial QC ! keep track of how many obs survived from current file ! ##################################### - ! TO DO: + ! TO DO: ! - refine use of quality flag? - ! - how to QC for proximity to water body? - ! - anything else missing that is done for SMOS + ! - how to QC for proximity to water body? + ! - anything else missing that is done for SMOS ! (either in preproc or in obs reader)?? ! ##################################### - + jj = 0 if (L1C_files .or. L1CE_files) then @@ -7407,7 +7414,7 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! initialize stats check for fore-minus-aft Tb diffs mm = 0 - + tmp_tb_diff_Sum = 0. tmp_tb_diff_SumOfSq = 0. @@ -7419,52 +7426,52 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & (mod(tmp_tb_qual_flag_1(nn),2)==0) .and. & ! lowest bit must be 0 (tmp_tb_1(nn) > Tb_min) .and. & ! elim neg nodata (tmp_tb_1(nn) < Tb_max) ! elim huge pos nodata - + keep_data_2 = & (mod(tmp_tb_qual_flag_2(nn),2)==0) .and. & ! lowest bit must be 0 (tmp_tb_2(nn) > Tb_min) .and. & ! elim neg nodata (tmp_tb_2(nn) < Tb_max) ! elim huge pos nodata ! thinning of L1C_TB_E obs - + if (L1CE_thinning) then - + tmp_keep = & ( mod(tmp_col(nn),L1CE_spacing) == 1 ) .and. & - ( mod(tmp_row(nn),L1CE_spacing) == 1 ) - + ( mod(tmp_row(nn),L1CE_spacing) == 1 ) + keep_data_1 = keep_data_1 .and. tmp_keep keep_data_2 = keep_data_2 .and. tmp_keep - + end if ! compute fore and aft average, put into "tb_1", "tmp_time_1" - + if (keep_data_1 .and. keep_data_2) then - + ! Compute stats for fore-minus-aft Tb differences. - ! Excessive diffs are found in bad L1C_TB files, which occur + ! Excessive diffs are found in bad L1C_TB files, which occur ! occasionally due to bad ANT_AZ files in L1B processing. ! Includes ALL SURFACES!!! ! - reichle, 22 Apr 2020 (resurrected) ! - reichle, 16 Oct 2020 (bug fix: do stats first, then avg) - + mm=mm+1 - - tmp_tb_diff = tmp_tb_1(nn) - tmp_tb_2(nn) - + + tmp_tb_diff = tmp_tb_1(nn) - tmp_tb_2(nn) + tmp_tb_diff_Sum = tmp_tb_diff_Sum + tmp_tb_diff tmp_tb_diff_SumOfSq = tmp_tb_diff_SumOfSq + tmp_tb_diff**2 - + ! put average of "fore" and "aft" into "tb_1", "tmp_time_1" - + tmp_tb_1( nn) = 0.5 *( tmp_tb_1( nn) + tmp_tb_2( nn) ) - tmp_time_1(nn) = 0.5D0*( tmp_time_1(nn) + tmp_time_2(nn) ) - + tmp_time_1(nn) = 0.5D0*( tmp_time_1(nn) + tmp_time_2(nn) ) + elseif (keep_data_2) then ! put "aft" data into "tb_1", "tmp_time_1" - + tmp_tb_1( nn) = tmp_tb_2( nn) tmp_time_1(nn) = tmp_time_2(nn) @@ -7473,78 +7480,78 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! nothing to do here ! - if only keep_data_1 is true (tmp_tb_1 and tmp_time_1 are already as needed) ! - if both keep_data_1 and keep_data_2 are false (next if block ignores data) - + end if ! apply QC and thinning, ensure that time stamp is within assimilation window - - if ( (keep_data_1 .or. keep_data_2) .and. & + + if ( (keep_data_1 .or. keep_data_2) .and. & (tmp_time_1(nn) > J2000_seconds_low) .and. & (tmp_time_1(nn) <= J2000_seconds_upp) & ) then - + jj=jj+1 - + tmp_lon( jj) = tmp_lon(nn) tmp_lat( jj) = tmp_lat(nn) tmp_tb_1( jj) = tmp_tb_1( nn) tmp_time_1(jj) = tmp_time_1(nn) - + end if - + end do - + ! finalize stats check for fore-minus-aft differences (ALL SURFACES!!!) ! - reichle, 22 Apr 2020 (resurrected) - + if (mm>1) then - + tmp_var = ( tmp_tb_diff_SumOfSq - (tmp_tb_diff_Sum**2)/real(mm) )/(real(mm-1)) - + if ( tmp_var > max_std_tb_fore_minus_aft**2 ) then - - write(err_msg, '(e12.5)') sqrt(tmp_var) - + + write(err_msg, '(e12.5)') sqrt(tmp_var) + err_msg = & 'Ignoring ALL obs in halforbit file b/c of excessive std-dev in fore-minus-aft Tbs. ' // & 'std-dev( tb_fore - tb_aft ) = ' // trim(err_msg) - + call ldas_warn(LDAS_GENERIC_WARNING, Iam, err_msg) - + jj = 0 ! results in N_obs_kept=0 below - + end if end if else ! L2_SM_AP - + do nn=1,N_obs_tmp - + keep_data_1 = & (mod(tmp_tb_qual_flag_1(nn),2)==0) .and. & ! lowest bit must be 0 (tmp_time_1(nn) > J2000_seconds_low) .and. & (tmp_time_1(nn) <= J2000_seconds_upp) .and. & (tmp_tb_1(nn) > Tb_min) .and. & ! elim neg nodata (tmp_tb_1(nn) < Tb_max) ! elim huge pos nodata - + if (keep_data_1) then - + jj=jj+1 - + tmp_lon( jj) = tmp_lon( nn) tmp_lat( jj) = tmp_lat( nn) tmp_tb_1( jj) = tmp_tb_1( nn) tmp_time_1(jj) = tmp_time_1(nn) - + end if - + end do - + end if ! (L1C_files .or. L1CE_files) - + N_obs_kept(kk) = jj - + ! ------------------------------------------------- ! ! map obs to tiles @@ -7552,20 +7559,20 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! for each observation ! a) determine grid cell that contains lat/lon ! b) determine tile within grid cell that contains lat/lon - + if (N_obs_kept(kk)>0) then - + allocate(tmp_tile_num(N_obs_kept(kk))) - + ! shift M36 obs lat/lon for proper assignment of M09 tile? - + if ( L1C_files .and. (index(tile_grid_d%gridtype, 'EASEv2_M09') /=0 .or. index(tile_grid_d%gridtype, 'EASEv2-M09') /=0 )) then - + ! temporarily shift lat/lon of obs for computation of nearest tile to ! avoid ambiguous assignment of M09 model tile within M36 obs grid cell - ! (center of M36 grid cell is equidistant from at least two M09 model + ! (center of M36 grid cell is equidistant from at least two M09 model ! tiles) -- reichle, 23 Aug 2013 - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_obs_kept(kk), & @@ -7574,9 +7581,9 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & this_obs_param, & tmp_tile_num, & tmp_shift_lat, tmp_shift_lon ) - + else - + call get_tile_num_for_obs(N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & N_obs_kept(kk), & @@ -7584,17 +7591,17 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & tmp_lon(1:N_obs_kept(kk)), & this_obs_param, & tmp_tile_num) - + end if - - ! make sure center-of-mass of tile that administers obs + + ! make sure center-of-mass of tile that administers obs ! is within EASEv2 M36 obs grid cell, discard obs otherwise ! (by setting tmp_tile_num to negative value) ! - ! this step eliminates SMAP obs that fall outside of the GEOS-5 + ! this step eliminates SMAP obs that fall outside of the GEOS-5 ! land mask at M09 resolution ! - ! not sure what this does if model is run on tile space other than + ! not sure what this does if model is run on tile space other than ! EASEv2_M09 - reichle, 11 Nov 2014 ! ! Bug fix - index used within loop should be "ii" (not "kk"). @@ -7607,71 +7614,71 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & ! distortion of EASE grid cells at high latitudes, and/or Gabrielle's ! use of an M36 innovations integration to derive Tb scaling files ! for the M09 (SMAP) system. - ! The problem with the piece of code is that it throws out far too many - ! obs (at all latitudes) if the tile space is coarse, e.g., that of the + ! The problem with the piece of code is that it throws out far too many + ! obs (at all latitudes) if the tile space is coarse, e.g., that of the ! 1/2 deg Lat/Lon grids of MERRA or MERRA-2. ! The piece of code may no longer be needed because of improvements in ! the obs readers and in get_obs_pred(), but the impact of removing the ! code on the SMAP L4_SM system is not clear. At this time, just prior - ! to finalizing the L4_SM "validated release", keep the code for the - ! EASEv2 M09 and M36 tile spaces that are relevant for the L4_SM + ! to finalizing the L4_SM "validated release", keep the code for the + ! EASEv2 M09 and M36 tile spaces that are relevant for the L4_SM ! system, but drop it for all other tile spaces. ! - reichle, 3 Feb 2016 - + if ( & (index(tile_grid_d%gridtype, 'EASEv2_M09') /=0) .or. & (index(tile_grid_d%gridtype, 'EASEv2_M36') /=0) ) then - + do ii=1,N_obs_kept(kk) - + if (tmp_tile_num(ii)>0) then - + call ease_convert('EASEv2_M36', & tile_coord(tmp_tile_num(ii))%com_lat, & tile_coord(tmp_tile_num(ii))%com_lon, & M36_col_ind_tile, M36_row_ind_tile ) - + call ease_convert('EASEv2_M36', & tmp_lat(ii), & tmp_lon(ii), & M36_col_ind_obs, M36_row_ind_obs ) - + if ( (nint(M36_col_ind_tile)/=nint(M36_col_ind_obs)) .or. & (nint(M36_row_ind_tile)/=nint(M36_row_ind_obs)) ) & tmp_tile_num(ii) = -9999 - + end if - + end do - + end if ! compute super-obs for each tile from all obs w/in that tile ! (also eliminate observations that are not in domain) - + do ii=1,N_obs_kept(kk) - + ind_tile = tmp_tile_num(ii) ! 1<=tmp_tile_num<=N_catd (unless nodata) - + if (ind_tile>0) then ! this step eliminates obs outside domain - + SMAP_data(ind_tile) = SMAP_data(ind_tile) + tmp_tb_1( ii) SMAP_lon( ind_tile) = SMAP_lon( ind_tile) + tmp_lon( ii) SMAP_lat( ind_tile) = SMAP_lat( ind_tile) + tmp_lat( ii) SMAP_time(ind_tile) = SMAP_time(ind_tile) + tmp_time_1(ii) - + N_obs_in_tile(ind_tile) = N_obs_in_tile(ind_tile) + 1 - + end if - + end do deallocate(tmp_tile_num) - + end if - + ! clean up - + if (allocated(tmp_lon )) deallocate(tmp_lon ) if (allocated(tmp_lat )) deallocate(tmp_lat ) if (allocated(tmp_col )) deallocate(tmp_col ) @@ -7682,55 +7689,55 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & if (allocated(tmp_tb_2 )) deallocate(tmp_tb_2 ) if (allocated(tmp_tb_qual_flag_1)) deallocate(tmp_tb_qual_flag_1) if (allocated(tmp_tb_qual_flag_2)) deallocate(tmp_tb_qual_flag_2) - + end do ! kk=1,N_fnames - + ! normalize - + do ii=1,N_catd - + if (N_obs_in_tile(ii)>1) then - + tmpreal = real(N_obs_in_tile(ii)) - + SMAP_data(ii) = SMAP_data(ii)/ tmpreal SMAP_lon( ii) = SMAP_lon( ii)/ tmpreal SMAP_lat( ii) = SMAP_lat( ii)/ tmpreal SMAP_time(ii) = SMAP_time(ii)/real(tmpreal,kind(0.0D0)) - + elseif (N_obs_in_tile(ii)==0) then - + SMAP_data(ii) = this_obs_param%nodata SMAP_lon( ii) = this_obs_param%nodata SMAP_lat( ii) = this_obs_param%nodata SMAP_time(ii) = real(this_obs_param%nodata,kind(0.0D0)) - + end if - + end do - + ! -------------------------------- - + ! set observation error standard deviation - + do ii=1,N_catd std_SMAP_data(ii) = this_obs_param%errstd end do - + ! -------------------------------- - + if (any(N_obs_in_tile>0)) then - + found_obs = .true. - - else - + + else + found_obs = .false. - + end if - + end if ! N_fnames==0 - + ! ------------------------------------------------- ! ! write "obslog" file @@ -7740,54 +7747,54 @@ subroutine read_obs_SMAP_halforbit_Tb( date_time, N_catd, this_obs_param, & YYYYMMDD_HHMMSSz = date_time2string(date_time) tmpstr80 = 'read_obs_SMAP_halforbit_Tb()' ! name of this subroutine - + do kk=1,N_fnames - + fname_tmp = trim(this_obs_param%path) // fname_list(kk) - + write (tmpstr12,'(i12)') N_obs_kept(kk) ! convert integer to string - + call add_to_obslog( YYYYMMDD_HHMMSSz, this_obs_param%descr, tmpstr80, & tmpstr12, fname_tmp ) - + end do end if ! clean up - + if (logit) write (logunit,*) 'read_obs_SMAP_halforbit_Tb(): done.' - + end subroutine read_obs_SMAP_halforbit_Tb - + ! ***************************************************************** subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observations_l) - - ! this subroutine turns off the assimilation of *individual* L1C_Tb obs - ! if corresponding L2AP_Tb obs are assimilated - + + ! this subroutine turns off the assimilation of *individual* L1C_Tb obs + ! if corresponding L2AP_Tb obs are assimilated + ! rationale: SMAP L1C_Tb obs are used to derive the disaggregated L2AP_Tb obs ! and should not be assimilated along with L2AP_Tb - + ! L1C_Tb obs are on the EASEv2 36 km grid ("M36") ! L2AP_Tb obs are on the EASEv2 9 km grid ("M09") ! reichle, 17 Jan 2014 ! reichle, 5 May 2015 - added *ascending* L2AP - + implicit none - + integer, intent(in) :: N_obs_param, N_obsl - - type(obs_param_type), dimension(N_obs_param), intent(in) :: obs_param - - type(obs_type), dimension(N_obsl), intent(inout) :: Observations_l - - ! local variables - logical :: turnoff_L1C_Tbh_A, turnoff_L1C_Tbv_A - logical :: turnoff_L1C_Tbh_D, turnoff_L1C_Tbv_D + type(obs_param_type), dimension(N_obs_param), intent(in) :: obs_param + + type(obs_type), dimension(N_obsl), intent(inout) :: Observations_l + + ! local variables + + logical :: turnoff_L1C_Tbh_A, turnoff_L1C_Tbv_A + logical :: turnoff_L1C_Tbh_D, turnoff_L1C_Tbv_D integer :: n, ii, j, N_obsf, N_cols, N_rows @@ -7802,20 +7809,20 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation integer :: ind_L1C_Tbv_D, ind_L2AP_Tbv_D real :: col, row - + integer, dimension(numprocs) :: N_obsl_vec, tmp_low_ind - + type(obs_type), dimension(:), allocatable :: Observations_f - + logical, dimension(:,:), allocatable :: mask_h_A, mask_v_A logical, dimension(:,:), allocatable :: mask_h_D, mask_v_D character(len=*), parameter :: Iam = 'turn_off_assim_SMAP_L1CTb' - ! ------------------------------------------------------------------------------ - + ! ------------------------------------------------------------------------------ + ! determine species and index numbers of conflicting obs species - + species_L1C_Tbh_A = -9999 species_L1C_Tbh_D = -9999 species_L1C_Tbv_A = -9999 @@ -7833,39 +7840,39 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation ind_L2AP_Tbh_D = -9999 ind_L2AP_Tbv_A = -9999 ind_L2AP_Tbv_D = -9999 - + do j=1,N_obs_param ! reichle, 1 Feb 2018: Only worry about assimilated species. - ! NOTE: At most one of L1C_Tb, L1C_Tb E09, or L1C_Tb E27 can be + ! NOTE: At most one of L1C_Tb, L1C_Tb E09, or L1C_Tb E27 can be ! assimilated at a time, checked in read_ens_upd_inputs(). - if (obs_param(j)%assim) then - + if (obs_param(j)%assim) then + select case (trim(obs_param(j)%descr)) - - case('SMAP_L1C_Tbh_A','SMAP_L1C_Tbh_E09_A','SMAP_L1C_Tbh_E27_A'); species_L1C_Tbh_A =obs_param(j)%species; ind_L1C_Tbh_A =j - case('SMAP_L1C_Tbh_D','SMAP_L1C_Tbh_E09_D','SMAP_L1C_Tbh_E27_D'); species_L1C_Tbh_D =obs_param(j)%species; ind_L1C_Tbh_D =j - case('SMAP_L1C_Tbv_A','SMAP_L1C_Tbv_E09_A','SMAP_L1C_Tbv_E27_A'); species_L1C_Tbv_A =obs_param(j)%species; ind_L1C_Tbv_A =j - case('SMAP_L1C_Tbv_D','SMAP_L1C_Tbv_E09_D','SMAP_L1C_Tbv_E27_D'); species_L1C_Tbv_D =obs_param(j)%species; ind_L1C_Tbv_D =j - case('SMAP_L2AP_Tbh_A'); species_L2AP_Tbh_A=obs_param(j)%species; ind_L2AP_Tbh_A=j - case('SMAP_L2AP_Tbh_D'); species_L2AP_Tbh_D=obs_param(j)%species; ind_L2AP_Tbh_D=j - case('SMAP_L2AP_Tbv_A'); species_L2AP_Tbv_A=obs_param(j)%species; ind_L2AP_Tbv_A=j - case('SMAP_L2AP_Tbv_D'); species_L2AP_Tbv_D=obs_param(j)%species; ind_L2AP_Tbv_D=j + + case('SMAP_L1C_Tbh_A','SMAP_L1C_Tbh_E09_A','SMAP_L1C_Tbh_E27_A'); species_L1C_Tbh_A =obs_param(j)%species; ind_L1C_Tbh_A =j + case('SMAP_L1C_Tbh_D','SMAP_L1C_Tbh_E09_D','SMAP_L1C_Tbh_E27_D'); species_L1C_Tbh_D =obs_param(j)%species; ind_L1C_Tbh_D =j + case('SMAP_L1C_Tbv_A','SMAP_L1C_Tbv_E09_A','SMAP_L1C_Tbv_E27_A'); species_L1C_Tbv_A =obs_param(j)%species; ind_L1C_Tbv_A =j + case('SMAP_L1C_Tbv_D','SMAP_L1C_Tbv_E09_D','SMAP_L1C_Tbv_E27_D'); species_L1C_Tbv_D =obs_param(j)%species; ind_L1C_Tbv_D =j + case('SMAP_L2AP_Tbh_A'); species_L2AP_Tbh_A=obs_param(j)%species; ind_L2AP_Tbh_A=j + case('SMAP_L2AP_Tbh_D'); species_L2AP_Tbh_D=obs_param(j)%species; ind_L2AP_Tbh_D=j + case('SMAP_L2AP_Tbv_A'); species_L2AP_Tbv_A=obs_param(j)%species; ind_L2AP_Tbv_A=j + case('SMAP_L2AP_Tbv_D'); species_L2AP_Tbv_D=obs_param(j)%species; ind_L2AP_Tbv_D=j end select - + end if end do - + ! determine whether there is possibly a redundancy turnoff_L1C_Tbh_A = .false. turnoff_L1C_Tbh_D = .false. turnoff_L1C_Tbv_A = .false. turnoff_L1C_Tbv_D = .false. - + if (ind_L1C_Tbh_A>0 .and. ind_L2AP_Tbh_A>0) turnoff_L1C_Tbh_A = .true. if (ind_L1C_Tbh_D>0 .and. ind_L2AP_Tbh_D>0) turnoff_L1C_Tbh_D = .true. if (ind_L1C_Tbv_A>0 .and. ind_L2AP_Tbv_A>0) turnoff_L1C_Tbv_A = .true. @@ -7874,186 +7881,186 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation ! ----------------------------------------------------------------------- ! proceed only if there is work to be done - + if ( turnoff_L1C_Tbh_A .or. turnoff_L1C_Tbv_A .or. & turnoff_L1C_Tbh_D .or. turnoff_L1C_Tbv_D ) then - - ! gather local obs - + + ! gather local obs + #ifdef LDAS_MPI - + call MPI_Gather( & N_obsl, 1, MPI_integer, & - N_obsl_vec, 1, MPI_integer, & - 0, mpicomm, mpierr ) - + N_obsl_vec, 1, MPI_integer, & + 0, mpicomm, mpierr ) + #else - + N_obsl_vec(1) = N_obsl - + #endif - + if (root_proc) then - + N_obsf = sum(N_obsl_vec) - + allocate(Observations_f(N_obsf)) - + tmp_low_ind(1) = 1 - + do n=1,numprocs-1 - + tmp_low_ind(n+1) = tmp_low_ind(n) + N_obsl_vec(n) - + end do - + end if - + #ifdef LDAS_MPI - + call MPI_GATHERV( & Observations_l, N_obsl, MPI_obs_type, & Observations_f, N_obsl_vec, tmp_low_ind-1, MPI_obs_type, & - 0, mpicomm, mpierr ) - -#else - + 0, mpicomm, mpierr ) + +#else + Observations_f = Observations_l - + #endif - - ! --------------------------------------------------------- + + ! --------------------------------------------------------- ! - ! assemble 36 km EASEv2 mask of L2AP_Tb obs + ! assemble 36 km EASEv2 mask of L2AP_Tb obs call ease_extent( 'EASEv2_M36', N_cols, N_rows ) - + allocate( mask_h_A(N_cols,N_rows) ) allocate( mask_h_D(N_cols,N_rows) ) allocate( mask_v_A(N_cols,N_rows) ) allocate( mask_v_D(N_cols,N_rows) ) - + mask_h_A = .false. ! initialize mask_h_D = .false. ! initialize mask_v_A = .false. ! initialize mask_v_D = .false. ! initialize - + if (root_proc) then - + ! mask for H-pol ascending - + if (turnoff_L1C_Tbh_A) then - + do ii=1,N_obsf - - if (Observations_f(ii)%species==species_L2AP_Tbh_A) then - + + if (Observations_f(ii)%species==species_L2AP_Tbh_A) then + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; ! note conversion to one-based indices - + mask_h_A(nint(col)+1,nint(row)+1) = .true. - + end if - + end do end if ! mask for H-pol descending - + if (turnoff_L1C_Tbh_D) then - + do ii=1,N_obsf - - if (Observations_f(ii)%species==species_L2AP_Tbh_D) then - + + if (Observations_f(ii)%species==species_L2AP_Tbh_D) then + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; ! note conversion to one-based indices - + mask_h_D(nint(col)+1,nint(row)+1) = .true. - + end if - + end do end if - + ! mask for V-pol ascending - + if (turnoff_L1C_Tbv_A) then - + do ii=1,N_obsf - - if (Observations_f(ii)%species==species_L2AP_Tbv_A) then - + + if (Observations_f(ii)%species==species_L2AP_Tbv_A) then + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; ! note conversion to one-based indices - + mask_v_A(nint(col)+1,nint(row)+1) = .true. - + end if - + end do - + end if ! mask for V-pol descending - + if (turnoff_L1C_Tbv_D) then - + do ii=1,N_obsf - - if (Observations_f(ii)%species==species_L2AP_Tbv_D) then - + + if (Observations_f(ii)%species==species_L2AP_Tbv_D) then + call ease_convert('EASEv2_M36', Observations_f(ii)%lat, Observations_f(ii)%lon, & col, row) ! set mask=.true. for the M36 grid cell that contains the L2AP_Tb obs; ! note conversion to one-based indices - + mask_v_D(nint(col)+1,nint(row)+1) = .true. - + end if - + end do - + end if deallocate(Observations_f) - + end if ! (root_proc) - + ! MPI broadcast masks - + call MPI_Bcast(mask_h_A, N_cols*N_rows, MPI_LOGICAL, 0, mpicomm, mpierr) call MPI_Bcast(mask_h_D, N_cols*N_rows, MPI_LOGICAL, 0, mpicomm, mpierr) call MPI_Bcast(mask_v_A, N_cols*N_rows, MPI_LOGICAL, 0, mpicomm, mpierr) call MPI_Bcast(mask_v_D, N_cols*N_rows, MPI_LOGICAL, 0, mpicomm, mpierr) - - ! --------------------------------------------------------- + + ! --------------------------------------------------------- ! ! apply H-pol masks - + if (turnoff_L1C_Tbh_A) then - + do ii=1,N_obsl - - if (Observations_l(ii)%species==species_L1C_Tbh_A) then - + + if (Observations_l(ii)%species==species_L1C_Tbh_A) then + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) - + ! note conversion to one-based indices - + if (mask_h_A(nint(col)+1,nint(row)+1)) Observations_l(ii)%assim = .false. end if @@ -8063,16 +8070,16 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation end if if (turnoff_L1C_Tbh_D) then - + do ii=1,N_obsl - - if (Observations_l(ii)%species==species_L1C_Tbh_D) then - + + if (Observations_l(ii)%species==species_L1C_Tbh_D) then + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) - + ! note conversion to one-based indices - + if (mask_h_D(nint(col)+1,nint(row)+1)) Observations_l(ii)%assim = .false. end if @@ -8080,87 +8087,87 @@ subroutine turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observation end do end if - + ! apply V-pol masks - + if (turnoff_L1C_Tbv_A) then - + do ii=1,N_obsl - - if (Observations_l(ii)%species==species_L1C_Tbv_A) then - + + if (Observations_l(ii)%species==species_L1C_Tbv_A) then + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) - + ! note conversion to one-based indices - + if (mask_v_A(nint(col)+1,nint(row)+1)) Observations_l(ii)%assim = .false. end if - + end do - + end if - + if (turnoff_L1C_Tbv_D) then - + do ii=1,N_obsl - - if (Observations_l(ii)%species==species_L1C_Tbv_D) then - + + if (Observations_l(ii)%species==species_L1C_Tbv_D) then + call ease_convert('EASEv2_M36', & Observations_l(ii)%lat, Observations_l(ii)%lon, col, row) - + ! note conversion to one-based indices - + if (mask_v_D(nint(col)+1,nint(row)+1)) Observations_l(ii)%assim = .false. end if - + end do - + end if - + ! clean up - + deallocate( mask_h_A ) deallocate( mask_h_D ) deallocate( mask_v_A ) deallocate( mask_v_D ) - + end if ! (turnoff_L1C_Tbh_A .or. turnoff_L1C_Tbv_A .or. turnoff_L1C_Tbh_D .or. turnoff_L1C_Tbv_D) - + end subroutine turn_off_assim_SMAP_L1CTb - + ! ***************************************************************** - + subroutine read_obs_fnames( date_time, this_obs_param, & fname_of_fname_list, N_max, N_fnames, fname_list, & obs_dir_hier ) - + ! read the file within an obs Yyyyy/Mmm/Ddd directory that lists ! the obs file names, preface file names with "Yyyyy/Mmm/Ddd", ! and return in "fname_list" ! - ! optional input argument: - ! obs_dir_hier==1 : preface file names with "Yyyyy/Mmm" instead + ! optional input argument: + ! obs_dir_hier==1 : preface file names with "Yyyyy/Mmm" instead ! ! this subroutine is needed when obs file names cannot be predicted ! and must be provided in a short text file that lists the file names ! (e.g., SMAP Tb or soil moisture h5 files, ASCAT soil moisture BUFR files) ! ! reichle, 3 Jan 2014 - ! reichle, 8 Jun 2017: Use "%flistpath" and "%flistname" from "obs_param_type". + ! reichle, 8 Jun 2017: Use "%flistpath" and "%flistname" from "obs_param_type". ! A M Fox, reichle, 22 Sep 2023: added optional argument obs_dir_hier ! ! --------------------------------------------------------------------------------- - + implicit none - + type(date_time_type), intent(in) :: date_time - + type(obs_param_type), intent(in) :: this_obs_param - + character( *), intent(in) :: fname_of_fname_list integer, intent(in) :: N_max @@ -8168,11 +8175,11 @@ subroutine read_obs_fnames( date_time, this_obs_param, & integer, intent(out) :: N_fnames character(100), dimension(N_max), intent(out) :: fname_list - + integer, optional, intent(in) :: obs_dir_hier - + ! local variables - + character(300) :: fname character(200) :: fpath_tmp character( 80) :: fname_tmp @@ -8189,11 +8196,11 @@ subroutine read_obs_fnames( date_time, this_obs_param, & character(len=400) :: err_msg ! --------------------------------------------------------------------- - + write (YYYY,'(i4.4)') date_time%year write (MM ,'(i2.2)') date_time%month write (DD ,'(i2.2)') date_time%day - + YYYYMMDDdir = 'Y' // YYYY // '/M' // MM // '/D' // DD // '/' YYYYMMdir = 'Y' // YYYY // '/M' // MM // '/' @@ -8202,33 +8209,33 @@ subroutine read_obs_fnames( date_time, this_obs_param, & fpath_tmp = this_obs_param%path fname_tmp = fname_of_fname_list - ! use inputs from ensupd nml file if not empty + ! use inputs from ensupd nml file if not empty if (len_trim(this_obs_param%flistpath)>0) fpath_tmp = this_obs_param%flistpath if (len_trim(this_obs_param%flistname)>0) fname_tmp = this_obs_param%flistname - + fname = trim(fpath_tmp) // '/' // YYYYMMDDdir // trim(fname_tmp) - + open( 10, file=fname, form='formatted', action='read', iostat=istat) - + if (istat/=0) then err_msg = 'cannot open file ' // trim(fname) call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + if (logit) write(logunit,'(400A)') 'reading file ' // trim(fname) - + ! read list of file names istat = 0 ii = 0 - + do while (istat==0) - + read(10,*,iostat=istat) tmpstr80 - + if (istat==0) then - + ii = ii+1 if (ii>N_max) then @@ -8237,33 +8244,33 @@ subroutine read_obs_fnames( date_time, this_obs_param, & end if ! preface file names with "Yyyyy/Mmm/Ddd" (default) - + fname_list(ii) = YYYYMMDDdir // trim(tmpstr80) - + if (present(obs_dir_hier)) then - + if (obs_dir_hier == 1) then - + ! preface file names with "Yyyyy/Mmm" - + fname_list(ii) = YYYYMMdir // trim(tmpstr80) - + else - + err_msg = 'Unrecognized value of optional argument obs_dir_hier' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if end if - + else - + exit - + end if - + end do - + close(10, status='keep') N_fnames = ii @@ -8273,21 +8280,21 @@ end subroutine read_obs_fnames ! ***************************************************************** type(date_time_type) function SMAPdatetime_to_DateTimeType( SMAP_date_time ) - + ! convert SMAP date/time strings to LDASsa date_time_type ! reichle, 6 Jan 2014 - + implicit none - + ! input arguments - + character(len=*) :: SMAP_date_time - + ! local variables - + type(date_time_type) :: date_time - + character(4) :: YYYY character(2) :: MM, DD, HH, MI, SS @@ -8297,15 +8304,15 @@ type(date_time_type) function SMAPdatetime_to_DateTimeType( SMAP_date_time ) integer :: HH_is, HH_ie integer :: MI_is, MI_ie integer :: SS_is, SS_ie - + character(len=*), parameter :: Iam = 'SMAPdatetime_to_DateTimeType' character(len=400) :: err_msg ! --------------------------------------------------- - + if ( (len_trim(SMAP_date_time)==15) .and. & (SMAP_date_time( 9: 9)=='T') ) then - + ! format: "20010727T160914" (yyyymmddThhmmss) YYYY_is = 1; YYYY_ie = 4; @@ -8314,7 +8321,7 @@ type(date_time_type) function SMAPdatetime_to_DateTimeType( SMAP_date_time ) HH_is = 10; HH_ie = 11; MI_is = 12; MI_ie = 13; SS_is = 14; SS_ie = 15; - + else if ( & (len_trim(SMAP_date_time)==24) .and. & (SMAP_date_time( 5: 5)=='-') .and. & @@ -8324,45 +8331,45 @@ type(date_time_type) function SMAPdatetime_to_DateTimeType( SMAP_date_time ) (SMAP_date_time(17:17)==':') .and. & (SMAP_date_time(20:20)=='.') .and. & (SMAP_date_time(24:24)=='Z') ) then - + ! format: "2001-07-27T16:09:14.567Z" - + YYYY_is = 1; YYYY_ie = 4; MM_is = 6; MM_ie = 7; DD_is = 9; DD_ie = 10; HH_is = 12; HH_ie = 13; MI_is = 15; MI_ie = 16; SS_is = 18; SS_ie = 19; - - else + + else err_msg = 'invalid SMAPdatetime string' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end if - + YYYY = SMAP_date_time(YYYY_is:YYYY_ie) MM = SMAP_date_time(MM_is :MM_ie ) DD = SMAP_date_time(DD_is :DD_ie ) HH = SMAP_date_time(HH_is :HH_ie ) MI = SMAP_date_time(MI_is :MI_ie ) SS = SMAP_date_time(SS_is :SS_ie ) - + read (YYYY,*) date_time%year read (MM ,*) date_time%month read (DD ,*) date_time%day read (HH ,*) date_time%hour read (MI ,*) date_time%min read (SS ,*) date_time%sec - + call get_dofyr_pentad( date_time ) SMAPdatetime_to_DateTimeType = date_time - + end function SMAPdatetime_to_DateTimeType - + ! ***************************************************************** - + subroutine read_obs( & work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & @@ -8370,7 +8377,7 @@ subroutine read_obs( & this_obs_param, write_obslog, & found_obs, scaled_obs, & tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, tmp_time, tmp_assim ) - + ! read observations and optionally scale observations to model clim ! ! intended to be called by root_proc @@ -8378,10 +8385,10 @@ subroutine read_obs( & ! 10 Jun 2011 - removed model-based QC for MPI re-structuring (now done ! in connection with get_obs_pred()) ! 22 Nov 2011 - minor clean-up, renamed scale_obs_*() subroutines - ! 23 Aug 2013 - added possibility of using lat/lon of obs from reader + ! 23 Aug 2013 - added possibility of using lat/lon of obs from reader ! (rather than using lat/lon from model tile_coord) ! 31 Jan 2014 - added output of time stamp ("tmp_time") - ! 14 Jul 2014 - added summary diagnostic "scaled_obs" + ! 14 Jul 2014 - added summary diagnostic "scaled_obs" ! (indicates whether any of the founds obs were scaled) ! 6 Jun 2016 - added functionality to keep obs that cannot be scaled ! (but will not be assimilated; SMAP/SMOS Tb only for now) @@ -8389,23 +8396,23 @@ subroutine read_obs( & implicit none ! inputs: - + character(*), intent(in) :: work_path character(*), intent(in) :: exp_id type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim, N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(grid_def_type), intent(in) :: tile_grid_d - + integer, dimension(tile_grid_d%N_lon,tile_grid_d%N_lat), intent(in) :: & N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + type(obs_param_type), intent(in) :: this_obs_param logical, intent(in) :: write_obslog @@ -8436,7 +8443,7 @@ subroutine read_obs( & tmp_assim = .true. ! initialize ! initialize lat/lon/time info (may later be overwritten by individual reader) - + tmp_lon = tile_coord%com_lon tmp_lat = tile_coord%com_lat @@ -8445,75 +8452,75 @@ subroutine read_obs( & ! if needed, must be converted by obs reader (e.g., ASCAT/EUMETSAT) tmp_time = datetime_to_J2000seconds(date_time, J2000_epoch_id ) - + ! ----------------------------- - + ! choose appropriate reader - + select case (trim(this_obs_param%descr)) - + case ('ae_l2_sm_a', 'ae_l2_sm_d') - + call read_obs_ae_l2_sm( & work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + ! scale observations to model climatology - + if (this_obs_param%scale .and. found_obs) then scaled_obs = .true. - + call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) - + end if - + case ('ae_sm_LPRM_a_C', 'ae_sm_LPRM_d_C', 'ae_sm_LPRM_a_X', 'ae_sm_LPRM_d_X' ) - + call read_obs_ae_sm_LPRM( & work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + ! scale observations to model climatology - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. - + call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) end if - + case ('ASCAT_SM_A', 'ASCAT_SM_D' ) - + call read_obs_sm_ASCAT( & work_path, exp_id, & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + ! scale observations to model climatology - + if (this_obs_param%scale .and. found_obs) then scaled_obs = .true. call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) - + end if - + case ('ASCAT_META_SM','ASCAT_METB_SM','ASCAT_METC_SM' ) call read_obs_sm_ASCAT_EUMET( & @@ -8522,9 +8529,9 @@ subroutine read_obs( & this_obs_param, & found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, & tmp_time) - + ! scale observations to model climatology - + if (this_obs_param%scale .and. found_obs) then scaled_obs = .true. @@ -8532,181 +8539,181 @@ subroutine read_obs( & call scale_obs_sfmc_zscore( N_catd, tile_coord, & date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & tmp_obs, tmp_std_obs ) - - end if + + end if case ('isccp_tskin_gswp2_v1') - + call read_obs_isccp_tskin_gswp2_v1( & date_time, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. - + call scale_obs_tskin_zscore( N_catd, tile_coord, & date_time, this_obs_param, tmp_obs, tmp_std_obs ) - + end if - + case ('isccp_tskin_ceop3n4') - + call read_obs_isccp_tskin_ceop3n4( & date_time, N_catd, tile_coord, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. call scale_obs_tskin_zscore( N_catd, tile_coord, & date_time, this_obs_param, tmp_obs, tmp_std_obs ) - + end if - + case ('isccp_tskin_ceop3n4_hdASC') - + call read_obs_isccp_ts_ceop3n4_hdASC( & date_time, N_catd, tile_coord, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. - + call scale_obs_tskin_zscore( N_catd, tile_coord, & date_time, this_obs_param, tmp_obs, tmp_std_obs ) end if - + case ('RedArkOSSE_sm') - + call read_obs_RedArkOSSE_sm( & date_time, N_catd, tile_coord, this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + ! scale observations to model climatology (use AMSR-E subroutine) - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. - + call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) end if - + case ('RedArkOSSE_CLSMsynthSM') - + call read_obs_RedArkOSSE_CLSMsynthSM( & date_time, N_catd, this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + ! scale observations to model climatology (use AMSR-E subroutine) - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. - + call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) end if - + case ('RedArkOSSE_truth_50mm','RedArkOSSE_truth_400mm') - + call read_obs_RedArkOSSE_truth( & date_time, N_catd, this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) ! scale observations to model climatology (use AMSR-E subroutine) - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. - + call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) end if ! assimilation NOT implemented for RedArkOSSE_truth obs - + if (this_obs_param%assim) then err_msg = 'assimilation NOT implemented for RedArkOSSE_truth obs' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - - + + case ('VivianaOK_CLSMsynthSM') - + call read_obs_VivianaOK_CLSMsynthSM( & date_time, N_catd, this_obs_param, & found_obs, tmp_obs, tmp_std_obs ) - + ! scale observations to model climatology (use AMSR-E subroutine) - + if (this_obs_param%scale .and. found_obs) then - + scaled_obs = .true. - + call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) - + end if case('SMOS_SM_A','SMOS_SM_D') - + call read_obs_SMOS( & date_time, N_catd, this_obs_param, & dtstep_assim, tile_coord, tile_grid_d, & N_tile_in_cell_ij, tile_num_in_cell_ij, write_obslog, & found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat ) - + ! scale observations to model climatology - + if (this_obs_param%scale .and. found_obs) then scaled_obs = .true. call scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) - + end if case('SMOS_reg_Tbh_A','SMOS_reg_Tbh_D','SMOS_reg_Tbv_A','SMOS_reg_Tbv_D', & 'SMOS_fit_Tbh_A','SMOS_fit_Tbh_D','SMOS_fit_Tbv_A','SMOS_fit_Tbv_D') - + call read_obs_SMOS( & date_time, N_catd, this_obs_param, & dtstep_assim, tile_coord, tile_grid_d, & N_tile_in_cell_ij, tile_num_in_cell_ij, write_obslog, & found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat ) - + ! scale observations to model climatology - + if (this_obs_param%scale .and. found_obs) then scaled_obs = .true. call scale_obs_Tb_zscore( N_catd, tile_coord, date_time, & this_obs_param, tmp_obs, tmp_std_obs, tmp_assim ) - + end if - + case ('MOD10C1','MYD10C1') - + call read_obs_MODIS_SCF( & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat ) - + case('SMAP_L1C_Tbh_A', 'SMAP_L1C_Tbv_A', & 'SMAP_L1C_Tbh_D', 'SMAP_L1C_Tbv_D', & 'SMAP_L1C_Tbh_E09_A', 'SMAP_L1C_Tbv_E09_A', & @@ -8721,21 +8728,21 @@ subroutine read_obs( & dtstep_assim, tile_coord, tile_grid_d, & N_tile_in_cell_ij, tile_num_in_cell_ij, write_obslog, & found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, tmp_time ) - + ! scale observations to model climatology - + if (this_obs_param%scale .and. found_obs) then scaled_obs = .true. - + call scale_obs_Tb_zscore( N_catd, tile_coord, date_time, & this_obs_param, tmp_obs, tmp_std_obs, tmp_assim ) end if - - case('LaRC_tskin-GOESW', 'LaRC_tskin-GOESE', 'LaRC_tskin-MET09', & + + case('LaRC_tskin-GOESW', 'LaRC_tskin-GOESE', 'LaRC_tskin-MET09', & 'LaRC_tskin-FY2E-', 'LaRC_tskin-MTST2') - + call read_obs_LaRC_Tskin( & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & @@ -8745,427 +8752,427 @@ subroutine read_obs( & ! NOT IMPLEMENTED: scale observations to model climatology case default - + err_msg = 'unknown obs_param%descr' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) - + end select - - + + end subroutine read_obs - - + + ! ***************************************************************** - + subroutine scale_obs_sfmc_cdf( N_catd, tile_coord, this_obs_param, & tmp_obs, tmp_std_obs ) - + ! scale sfmc obs to model climatology via cdf matching - ! - ! use matlab functions "get_cdf_match_AMSR.m" and "get_model_and_obs_stats.m" + ! + ! use matlab functions "get_cdf_match_AMSR.m" and "get_model_and_obs_stats.m" ! to create input scaling files ! ! IMPORTANT: Make sure that model and obs data are in the SAME UNITS prior - ! to generating the input scaling files with the matlab routines. - ! Otherwise, the (linear) rescaling outside of the observed + ! to generating the input scaling files with the matlab routines. + ! Otherwise, the (linear) rescaling outside of the observed ! range (between obs_min and obs_max) will fail! For example, - ! ASCAT soil moisture retrievals must first be converted into - ! volumetric units that range roughly from 0 to 0.5 m3/m3 - ! (for some porosity) before they can be rescaled to the + ! ASCAT soil moisture retrievals must first be converted into + ! volumetric units that range roughly from 0 to 0.5 m3/m3 + ! (for some porosity) before they can be rescaled to the ! Catchment model's sfmc (which is in m3/m3). - + ! THIS SUBROUTINE COULD USE WORK... ! - map from scaling parameter domain to model domain ! reichle, 14 Oct 2005 ! reichle, 26 Sep 2006 ! ! bug fix: scaling via linear interpolation outside of [obs_min,obs_max] - ! reichle, 27 Jan 2006 + ! reichle, 27 Jan 2006 ! ! reichle, 22 Nov 2011 - renamed subroutine, minor clean-up, added comments - ! reichle, 9 Nov 2012 - edited to enable use of stats file that do not + ! reichle, 9 Nov 2012 - edited to enable use of stats file that do not ! perfectly match current domain (as in tile_coord) - + implicit none - + integer, intent(in) :: N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(obs_param_type), intent(in) :: this_obs_param - + ! inout - + real, intent(inout), dimension(N_catd) :: tmp_obs real, intent(inout), dimension(N_catd) :: tmp_std_obs - + ! local variables - + real, parameter :: no_data_stats = -9999. - + real, parameter :: tol = 1e-2 - + character(300) :: fname - + integer :: i, j, N_sclprm, N_poly, ind - + real :: edge_min, edge_max, edge_dx - + real :: tmpreal, x, x0, x1, y0, y1 - + integer, dimension(:), allocatable :: tmp_tile_id - + real, dimension(:), allocatable :: std_obs, std_mod, min_obs, max_obs - + real, dimension(:,:), allocatable :: fit_coeff character(len=*), parameter :: Iam = 'scale_obs_sfmc_cdf' character(len=400) :: err_msg - + ! ------------------------------------------------------------------ - + ! read scaling parameters from file - + fname = trim(this_obs_param%scalepath) // '/' // & trim(this_obs_param%scalename) - + if (logit) write (logunit,*) 'scaling obs species ', this_obs_param%species, ':' if (logit) write (logunit,'(400A)') ' reading ', trim(fname) - + open(10, file=fname, form='formatted', action='read') - - read(10,*) N_sclprm, N_poly, edge_min, edge_max, edge_dx - + + read(10,*) N_sclprm, N_poly, edge_min, edge_max, edge_dx + ! minimal consistency check - + if (N_catd>N_sclprm) then err_msg = 'N_sclprm too small' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + allocate( tmp_tile_id( N_sclprm )) allocate( std_obs( N_sclprm )) allocate( std_mod( N_sclprm )) allocate( min_obs( N_sclprm )) allocate( max_obs( N_sclprm )) allocate( fit_coeff( N_sclprm, N_poly+1 )) - + do i=1,N_sclprm - + read (10,*) tmp_tile_id(i), std_obs(i), std_mod(i), & min_obs(i), max_obs(i), fit_coeff(i,:) - + end do - + close(10,status='keep') - + ! -------------------------------------------------------------- - + ! scale observations - + do i=1,N_catd - + ! check for no-data-values in observation ! (any negative number could be no-data-value for observations) - if (tmp_obs(i)>=0.) then - + if (tmp_obs(i)>=0.) then + ! find ind for current tile id in scaling parameters - + if (tmp_tile_id(i)==tile_coord(i)%tile_id) then - + ind = i ! educated guess for global domain - + else - + do ind=1,N_sclprm - + if (tmp_tile_id(ind)==tile_coord(i)%tile_id) exit - + end do - + if (ind>N_sclprm) then err_msg = 'tile_id not found in scaling parameter file' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + end if - - + + ! make sure obs falls within range valid for scaling and ! check for no-data-value in fit parameters - + if ( (tmp_obs(i)>=edge_min) .and. & (tmp_obs(i)<=edge_max) .and. & (abs(fit_coeff(ind,1)-no_data_stats) > tol) .and. & - (abs(fit_coeff(ind,2)-no_data_stats) > tol) ) then - + (abs(fit_coeff(ind,2)-no_data_stats) > tol) ) then + ! evaluate polynomial fit ! ! fit_coeff(ind,1) is coeff for highest order ! fit_coeff(ind,N_poly+1) is constant term - + ! evaluate polynomial at tmp_obs if min_obs<=tmp_obs<=max_obs, ! otherwise at min_obs or max_obs (accordingly) - + x = min( max( tmp_obs(i), min_obs(ind) ), max_obs(ind) ) - - tmpreal = fit_coeff(ind,1) - + + tmpreal = fit_coeff(ind,1) + do j=1,N_poly - + tmpreal = tmpreal*x + fit_coeff(ind,j+1) - + end do - + if (tmp_obs(i)max_obs(ind)) then - + ! linear interpolation between max_obs and max(edges) - ! (NOTE: model and obs data must be in SAME units, + ! (NOTE: model and obs data must be in SAME units, ! reichle, 22 Nov 2011) - + y1 = edge_max y0 = tmpreal x1 = edge_max x0 = max_obs(ind) - + tmp_obs(i) = (y1-y0)/(x1-x0)*( tmp_obs(i) - x0 ) + y0 - + else - + ! accept polynomial fit as is - + tmp_obs(i) = tmpreal - + end if - + ! scale observation error std - + tmp_std_obs(i) = std_mod(ind)/std_obs(ind)*tmp_std_obs(i) - + else - + tmp_obs(i) = this_obs_param%nodata - + end if - + end if ! qc check after scaling - + if ((tmp_obs(i)>edge_max) .or. (tmp_obs(i)N_sclprm) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'something is wrong') end if - + ! -------------------------------------------------------------- - - ! scale observations (at this point all obs are of type + + ! scale observations (at this point all obs are of type ! isccp_tskin_gswp2_v1 because of the way the subroutine is called ! from subroutine read_obs()) - + do i=1,N_catd - + ! check for no-data-values in observation (any neg Tskin is no_data) - + if (tmp_obs(i)>=0.) then - + ! find ind for current tile id in scaling parameters - + do ind=1,N_sclprm - + if (sclprm_tile_id(ind)==tile_coord(i)%tile_id) exit - + end do - + ! sanity check (against accidental use of wrong tile space) - + if ( abs(tile_coord(i)%com_lat-sclprm_lat(ind))>tol .or. & abs(tile_coord(i)%com_lon-sclprm_lon(ind))>tol ) then err_msg = 'something wrong' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + ! check for no-data-values in observation and fit parameters ! (any negative number could be no-data-value for observations) - + if ( sclprm_mean_obs(ind)>0. .and. & sclprm_mean_mod(ind)>0. .and. & sclprm_std_obs(ind)>=0. .and. & sclprm_std_mod(ind)>=0. ) then - + ! scale via standard normal deviates - - tmpreal = sclprm_std_mod(ind)/sclprm_std_obs(ind) - + + tmpreal = sclprm_std_mod(ind)/sclprm_std_obs(ind) + tmp_obs(i) = sclprm_mean_mod(ind) & - + tmpreal*(tmp_obs(i)-sclprm_mean_obs(ind)) - + + tmpreal*(tmp_obs(i)-sclprm_mean_obs(ind)) + ! scale observation error std - + tmp_std_obs(i) = tmpreal*tmp_std_obs(i) - + else - + tmp_obs(i) = this_obs_param%nodata - + end if - + end if - + end do - + deallocate(sclprm_tile_id) - deallocate(sclprm_lon) - deallocate(sclprm_lat) - deallocate(sclprm_mean_obs) - deallocate(sclprm_std_obs) - deallocate(sclprm_mean_mod) - deallocate(sclprm_std_mod) - + deallocate(sclprm_lon) + deallocate(sclprm_lat) + deallocate(sclprm_mean_obs) + deallocate(sclprm_std_obs) + deallocate(sclprm_mean_mod) + deallocate(sclprm_std_mod) + end subroutine scale_obs_tskin_zscore - + ! ***************************************************************** - + subroutine scale_obs_sfmc_zscore( N_catd, tile_coord, & date_time, this_obs_param, tmp_lon, tmp_lat, tmp_time, & tmp_obs, tmp_std_obs ) - + ! scale sfmc or sfds obs to model climatology via standard-normal-deviate (zscore) ! scaling using seasonally varying (pentad) stats - ! Grid information is read from a NetCDF file - ! + ! Grid information is read from a NetCDF file + ! ! Scaling parameters are read from a NetCDF file that contains the following: ! variables: ! int version ; @@ -9199,30 +9206,30 @@ subroutine scale_obs_sfmc_zscore( N_catd, tile_coord, & ! n_data:standard_name = "number of data points" ; ! ! A M Fox, reichle, April 2023 - - use netcdf + + use netcdf implicit none - + integer, intent(in) :: N_catd - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - + type(date_time_type), intent(in) :: date_time - + type(obs_param_type), intent(in) :: this_obs_param real, intent(in), dimension(N_catd) :: tmp_lon, tmp_lat real*8, intent(in), dimension(N_catd) :: tmp_time ! J2000 seconds - + ! inout - + real, intent(inout), dimension(N_catd) :: tmp_obs real, intent(inout), dimension(N_catd) :: tmp_std_obs - + ! ------------------- - + character(300) :: fname - + integer :: i, ind, pp, j_ind, i_ind integer :: ncid, varid, ierr, ierr2 integer :: pentad_dimid, lon_dimid, lat_dimid @@ -9233,44 +9240,44 @@ subroutine scale_obs_sfmc_zscore( N_catd, tile_coord, & integer, dimension(3) :: start, icount logical :: file_exists - + real :: tmpreal, this_lon, this_lat, ll_lon, ll_lat, dlon, dlat - + integer, dimension(:), allocatable :: sclprm_tile_id integer, dimension(:), allocatable :: pentads - + real, dimension(:,:), allocatable :: sclprm_mean_obs, sclprm_std_obs real, dimension(:,:), allocatable :: sclprm_mean_mod, sclprm_std_mod real, dimension(:,:), allocatable :: sclprm_min_mod, sclprm_max_mod - + character(len=*), parameter :: Iam = ' scale_obs_sfmc_zscore' character(len=400) :: err_msg - + ! ------------------------------------------------------------------ - + ! Read scaling parameters from file fname = trim(this_obs_param%scalepath) // '/' // & trim(this_obs_param%scalename) // '.nc4' - + if (logit) write (logunit,*) 'scaling obs species ', this_obs_param%species, ':' if (logit) write (logunit,'(400A)') ' reading ', trim(fname) ! Check if file exists inquire(file=fname, exist=file_exists) - + if (.not. file_exists) then - + err_msg = 'Scaling parameter file not found' - call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) + call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) end if - + ! Determine pentad to use pp = date_time%pentad - + ! Open the NetCDF file ierr = nf90_open(fname, nf90_nowrite, ncid) @@ -9287,11 +9294,11 @@ subroutine scale_obs_sfmc_zscore( N_catd, tile_coord, & ierr = nf90_inq_dimid(ncid, 'pentad', pentad_dimid) ierr = nf90_inq_dimid(ncid, 'lon', lon_dimid) ierr = nf90_inq_dimid(ncid, 'lat', lat_dimid) - + ierr = nf90_inquire_dimension(ncid, pentad_dimid, len = N_pentad) ierr = nf90_inquire_dimension(ncid, lon_dimid, len = N_lon) ierr = nf90_inquire_dimension(ncid, lat_dimid, len = N_lat) - + ! Get the variable IDs ierr = nf90_inq_varid(ncid, 'o_mean', o_mean_varid) @@ -9299,18 +9306,18 @@ subroutine scale_obs_sfmc_zscore( N_catd, tile_coord, & ierr = nf90_inq_varid(ncid, 'm_mean', m_mean_varid) ierr = nf90_inq_varid(ncid, 'm_std', m_std_varid) ierr = nf90_inq_varid(ncid, 'm_min', m_min_varid) - ierr = nf90_inq_varid(ncid, 'm_max', m_max_varid) - + ierr = nf90_inq_varid(ncid, 'm_max', m_max_varid) + ! Read grid variables ierr = nf90_get_var(ncid, ll_lon_varid, ll_lon) ierr = nf90_get_var(ncid, ll_lat_varid, ll_lat) ierr = nf90_get_var(ncid, dlon_varid, dlon) ierr = nf90_get_var(ncid, dlat_varid, dlat) - + start = [1, 1, pp] icount = [N_lat, N_lon, 1 ] - + ! Read mean and std variables allocate(sclprm_mean_obs(N_lat, N_lon), sclprm_std_obs(N_lat, N_lon)) @@ -9327,44 +9334,44 @@ subroutine scale_obs_sfmc_zscore( N_catd, tile_coord, & ! Close the netcdf file ierr = nf90_close(ncid) - + ! -------------------------------------------------------------- - + ! Scale observations (at this point all obs are of same type because ! of the way the subroutine is called from subroutine read_obs() - + do i=1,N_catd - + ! Check for no-data-values in observation (any neg value is no_data) if (tmp_obs(i)>=0.) then - + ! ll_lon and ll_lat refer to lower left corner of grid cell ! (as opposed to the grid point in the center of the grid cell) - + this_lon = tmp_lon(i) this_lat = tmp_lat(i) - + ! Find indices for current tile lat and lon on scaling parameter grid - i_ind = ceiling((this_lon - ll_lon)/dlon) - j_ind = ceiling((this_lat - ll_lat)/dlat) - + i_ind = ceiling((this_lon - ll_lon)/dlon) + j_ind = ceiling((this_lat - ll_lat)/dlat) + ! Check for no-data-values in observation and fit parameters ! (any negative number could be no-data-value for observations) - + if ( sclprm_mean_obs(j_ind, i_ind)>0. .and. & sclprm_mean_mod(j_ind, i_ind)>0. .and. & sclprm_std_obs(j_ind, i_ind)>=0. .and. & sclprm_std_mod(j_ind, i_ind)>=0. ) then - + ! Scale via standard normal deviates - - tmpreal = sclprm_std_mod(j_ind, i_ind)/sclprm_std_obs(j_ind, i_ind) - + + tmpreal = sclprm_std_mod(j_ind, i_ind)/sclprm_std_obs(j_ind, i_ind) + tmp_obs(i) = sclprm_mean_mod(j_ind, i_ind) & - + tmpreal*(tmp_obs(i)-sclprm_mean_obs(j_ind, i_ind)) - + + tmpreal*(tmp_obs(i)-sclprm_mean_obs(j_ind, i_ind)) + ! Check of tmp_obs is within range of model climatology if (tmp_obs(i)=0.) then - + ! find ind for current tile id in scaling parameters ind = -9999 ! initialize to negative integer - + if (N_sclprm==N_catd) then - + ! try an educated guess (to avoid an additional loop) - - if (sclprm_tile_id(i)==tile_coord(i)%tile_id) ind = i - + + if (sclprm_tile_id(i)==tile_coord(i)%tile_id) ind = i + end if - + if (ind<0) then ! the educated guess failed, try again - + do ind=1,N_sclprm - + if (sclprm_tile_id(ind)==tile_coord(i)%tile_id) exit - + end do if (ind>N_sclprm) then @@ -9744,77 +9751,77 @@ subroutine scale_obs_Tb_zscore( N_catd, tile_coord, date_time, this_obs_param, end if end if - + ! check for no-data-values in observation and fit parameters ! (any negative number could be no-data-value for observations) - + if ( sclprm_mean_obs(ind)>0. .and. & sclprm_mean_mod(ind)>0. .and. & sclprm_std_obs( ind)>0. .and. & sclprm_std_mod( ind)>0. ) then - - + + ! sanity check (against accidental use of wrong tile space) - + if ( abs(tile_coord(i)%com_lat-sclprm_lat(ind))>tol .or. & abs(tile_coord(i)%com_lon-sclprm_lon(ind))>tol ) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'something wrong') end if ! scale - + if (scale_mean_only) then - + ! adjust mean only - + tmp_obs(i) = & - tmp_obs(i) - sclprm_mean_obs(ind) + sclprm_mean_mod(ind) - + tmp_obs(i) - sclprm_mean_obs(ind) + sclprm_mean_mod(ind) + else - + ! scale via standard normal deviates - - tmpreal = sclprm_std_mod(ind)/sclprm_std_obs(ind) - + + tmpreal = sclprm_std_mod(ind)/sclprm_std_obs(ind) + tmp_obs(i) = sclprm_mean_mod(ind) & - + tmpreal*(tmp_obs(i)-sclprm_mean_obs(ind)) - + + tmpreal*(tmp_obs(i)-sclprm_mean_obs(ind)) + ! scale observation error std - + tmp_std_obs(i) = tmpreal*tmp_std_obs(i) - + end if - + else - + ! keep unscaled obs, provide only "do not assimilate" info ! reichle, 6 Jun 2016 tmp_assim(i) = .false. - + !tmp_obs(i) = this_obs_param%nodata - + end if - + end if - + end do - + deallocate(sclprm_ang) - + deallocate(sclprm_tile_id) - deallocate(sclprm_lon) - deallocate(sclprm_lat) - - deallocate(sclprm_mean_obs) - deallocate(sclprm_std_obs) - deallocate(sclprm_mean_mod) - deallocate(sclprm_std_mod) - + deallocate(sclprm_lon) + deallocate(sclprm_lat) + + deallocate(sclprm_mean_obs) + deallocate(sclprm_std_obs) + deallocate(sclprm_mean_mod) + deallocate(sclprm_std_mod) + end subroutine scale_obs_Tb_zscore - + ! ***************************************************************** - + subroutine collect_obs( & work_path, exp_id, date_time, dtstep_assim, & N_catl, & @@ -9822,27 +9829,27 @@ subroutine collect_obs( N_catl_vec, low_ind, l2f, & N_obs_param, obs_param, N_obsl_max, write_obslog, & N_obsl, Observations_l, found_obs_f ) - - ! check for observations that must be assimilated, + + ! check for observations that must be assimilated, ! collect into measurement vector ! ! a total of "N_obsl" observations are returned in "Observations_l" ! ! 25 Jul 2005 - rewritten - ! 10 Jun 2011 - re-structured for MPI, incl. removal of model-based QC + ! 10 Jun 2011 - re-structured for MPI, incl. removal of model-based QC ! (now done in connection with get_obs_pred()) ! 24 Dec 2013 - added "obslog" output files (not yet implemented for all readers!) ! 31 Dec 2014 - added time stamp to Observations - ! 21 Mar 2014 - sort Observations_l by tilenum and then species to avoid lay-out + ! 21 Mar 2014 - sort Observations_l by tilenum and then species to avoid lay-out ! dependency for MPI parallel execution implicit none - + character(*), intent(in) :: work_path character(*), intent(in) :: exp_id - + type(date_time_type), intent(in) :: date_time - + integer, intent(in) :: dtstep_assim integer, intent(in) :: N_catl, N_catf @@ -9851,41 +9858,41 @@ subroutine collect_obs( integer, dimension(N_catl), intent(in) :: l2f ! tile_coord_f of catchments in domain (length N_catf) - + type(tile_coord_type), dimension(:), pointer :: tile_coord_f ! input - + type(grid_def_type), intent(in) :: tile_grid_f - + ! N_tile_in_cell_ij and tile_num_in_cell_ij are on the "full" domain ! and guaranteed to be allocated ONLY for the root_proc ! (but may be allocated on all processors depending on obs_param%FOV) - integer, dimension(:,:), pointer :: N_tile_in_cell_ij_f ! input - integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij_f ! input - + integer, dimension(:,:), pointer :: N_tile_in_cell_ij_f ! input + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij_f ! input + integer, intent(in) :: N_obs_param - - type(obs_param_type), dimension(N_obs_param), intent(in) :: obs_param - + + type(obs_param_type), dimension(N_obs_param), intent(in) :: obs_param + integer, intent(in) :: N_obsl_max - + logical, intent(in) :: write_obslog - + integer, intent(out) :: N_obsl type(obs_type), dimension(N_obsl_max), intent(out) :: Observations_l - + logical, intent(out) :: found_obs_f ! ---------------------------- - + ! locals - + logical :: found_obs, scaled_obs, any_scaled_obs integer :: obs_count, species integer :: ii, ind_start, ind_end, N_tmp, this_tilenum, this_tilenum_new - + real, dimension(N_catl) :: tmp_obs, tmp_std_obs, tmp_lon, tmp_lat real, dimension(N_catf) :: tmp_obs_f, tmp_std_obs_f, tmp_lon_f, tmp_lat_f @@ -9896,46 +9903,46 @@ subroutine collect_obs( logical, dimension(N_catf) :: tmp_assim_f integer, dimension(numprocs) :: N_obsl_vec - + integer, dimension(N_obs_param) :: tmp_species - + integer, dimension(:), allocatable :: indx, tilenums - + character( 12) :: tmpstr12 character(300) :: tmpstr300 character(len=*), parameter :: Iam = 'collect_obs' character(len=400) :: err_msg - + ! --------------------------------------------------------------- - if (logit) write (logunit,*) 'collecting observations...' - + if (logit) write (logunit,*) 'collecting observations...' + ! obs_count serves as counter for total number of observations, ! excluding no-data-values - + obs_count = 0 any_scaled_obs = .false. ! initialize - + ! ---------------------------------------------------------- - + do species = 1,N_obs_param - + ! make sure species number here is consistent with ! definitions in nml input file - + if (obs_param(species)%species .ne. species) then call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'something wrong') end if - + if (root_proc) then - - ! subroutine read_obs() reads all observations in obs files - ! (typically global) and returns a vector in (full domain) + + ! subroutine read_obs() reads all observations in obs files + ! (typically global) and returns a vector in (full domain) ! tile space with the values of the observations (at most one ! observation per tile and per species) - + call read_obs( & work_path, exp_id, & date_time, dtstep_assim, N_catf, tile_coord_f, & @@ -9943,52 +9950,52 @@ subroutine collect_obs( obs_param(species), write_obslog, & found_obs, scaled_obs, & tmp_obs_f, tmp_std_obs_f, tmp_lon_f, tmp_lat_f, tmp_time_f, tmp_assim_f) - + if (scaled_obs) any_scaled_obs = .true. end if - + ! put "tmp_obs" (in "tile" space) into "compressed" vector "Observations_l" ! - ! each observation is "managed" ("administered") by the processor that + ! each observation is "managed" ("administered") by the processor that ! contains the tile to which the observation was assigned ! ! the assignment of each observation to tile space does NOT imply that ! only the observations from one (local) processor impact the increments - ! for that processor --> see "halo" below - + ! for that processor --> see "halo" below + #ifdef LDAS_MPI - + call MPI_BCAST(found_obs, 1, MPI_LOGICAL, 0,mpicomm, mpierr) - -#endif - - if (found_obs) then - + +#endif + + if (found_obs) then + ! map from full to local domain - + call f2l_real( N_catf,N_catl,N_catl_vec,low_ind, tmp_obs_f, tmp_obs) call f2l_real( N_catf,N_catl,N_catl_vec,low_ind, tmp_std_obs_f, tmp_std_obs) call f2l_real( N_catf,N_catl,N_catl_vec,low_ind, tmp_lon_f, tmp_lon) call f2l_real( N_catf,N_catl,N_catl_vec,low_ind, tmp_lat_f, tmp_lat) call f2l_real8( N_catf,N_catl,N_catl_vec,low_ind, tmp_time_f, tmp_time) - + call f2l_logical(N_catf,N_catl,N_catl_vec,low_ind, tmp_assim_f, tmp_assim) - + ! NOTE: "Observations" here are l(ocal) obs only - + call put_into_Observations( obs_param(species), N_obsl_max, N_catl, l2f, & tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, tmp_time, tmp_assim, & obs_count, Observations_l ) - + end if - + end do #ifdef LDAS_MPI - + call MPI_BCAST(any_scaled_obs, 1, MPI_LOGICAL, 0,mpicomm, mpierr) #endif @@ -9997,83 +10004,83 @@ subroutine collect_obs( ! ----------------------------------------------------------------- ! - ! sort relevant elements of Observations_l by tilenum and then by species + ! sort relevant elements of Observations_l by tilenum and then by species ! to avoid lay-out dependency for MPI parallel execution ! - reichle, 21 March 2014 - + if (N_obsl>1) then ! sort only if 2 or more obs - + allocate(indx( N_obsl)) - + allocate(tilenums(N_obsl)) - + tilenums = Observations_l(1:N_obsl)%tilenum - - ! get index vector, NOTE: + + ! get index vector, NOTE: ! nr_indexx() does not change input "arr" ! nr_indexx() only works with *real* input "arr" - + call nr_indexx( N_obsl, real(tilenums), indx(1:N_obsl) ) - + ! apply sort by tilenum - + Observations_l(1:N_obsl) = Observations_l(indx(1:N_obsl)) - + ! now make sure that within each tilenum, Observations_l are sorted by species - + ind_start = 1 - + this_tilenum = Observations_l(ind_start)%tilenum - + do ii=2,N_obsl - + this_tilenum_new = Observations_l(ii)%tilenum if ( (this_tilenum_new/=this_tilenum) .or. (ii==N_obsl) ) then - + if ( (this_tilenum_new/=this_tilenum) .and. (ii<=N_obsl) ) then - + ind_end = ii-1 - + else - + ind_end = N_obsl - + end if - - ! Observations_l(ind_start:ind_end) are the complete subset + + ! Observations_l(ind_start:ind_end) are the complete subset ! of (local) obs with the same tilenum - + N_tmp = ind_end - ind_start + 1 if (N_tmp>1) then - + tmp_species(1:N_tmp) = Observations_l(ind_start:ind_end)%species - + ! get index vector for sorting by species (see NOTES above!) - + call nr_indexx( N_tmp, real(tmp_species(1:N_tmp)), indx(1:N_tmp) ) - + ! apply sort by species - + indx(1:N_tmp) = indx(1:N_tmp) + ind_start - 1 ! add offset - + Observations_l(ind_start:ind_end) = Observations_l(indx(1:N_tmp)) end if ! re-initialize - + ind_start = ii - + this_tilenum = this_tilenum_new - + end if - + end do - + ! clean up - + deallocate(indx) deallocate(tilenums) @@ -10082,7 +10089,7 @@ subroutine collect_obs( ! ----------------------------------------------------------------- ! ! check whether number of obs exceeds max number allowed - + if (N_obsl>N_obsl_max) then err_msg = 'N_obsl > N_obsl_max, too many observations' call ldas_abort(LDAS_GENERIC_ERROR, Iam, err_msg) @@ -10090,111 +10097,111 @@ subroutine collect_obs( ! make sure L1C_Tb obs are not assimilated at the same time/place where ! corresponding disaggregated L2AP_Tb observations are assimilated - + call turn_off_assim_SMAP_L1CTb(N_obs_param, obs_param, N_obsl, Observations_l) - + ! gather information about number of observations assigned to each processor - + #ifdef LDAS_MPI - + call MPI_AllGather( & N_obsl, 1, MPI_integer, & - N_obsl_vec, 1, MPI_integer, & - mpicomm, mpierr ) - + N_obsl_vec, 1, MPI_integer, & + mpicomm, mpierr ) + #else - + N_obsl_vec(1) = N_obsl - + #endif - + if (any(N_obsl_vec>0)) then - + found_obs_f = .true. ! found obs somewhere in full domain - + else - - found_obs_f = .false. - + + found_obs_f = .false. + end if ! ----------------------------------------------------------------- - + ! TO DO: scale Observations (note: no scaling for moisture contents) - + ! may be needed for multi-variate assimilation, reichle, 10 Jun 2011 !!do i=1,N_obs - !! + !! !! select case (Observations(i)%species) - !! + !! !! case (5) - !! + !! !! Observations(i)%obs = Observations(i)%obs / scale_temp !! Observations(i)%obsvar = Observations(i)%obsvar /(scale_temp**2) - !! + !! !! end select - !! + !! !!end do - + ! ----------------------------------------------------------------- ! determine total number of observations to be assimilated ! ! NOTE: This number may be different from the total number of obs - ! recorded in the "obslog" file because: + ! recorded in the "obslog" file because: ! (i) the number of obs recorded in the "obslog" file is ! before obs from separate files may have been aggregated ! to "super-obs", and ! (ii) older obs readers may not contribute to the "obslog" file. - + write (tmpstr12,'(i12)') sum(N_obsl_vec) ! convert integer to string - + tmpstr300 = 'collect_obs(): read N_obsf = ' // tmpstr12 // & ' [after obs-based QC, super-obbing' - + if (.not. any_scaled_obs) then - + tmpstr300 = trim(tmpstr300) // ']' - + else - + tmpstr300 = trim(tmpstr300) // ', scaling]' - + end if - + if (logit) write (logunit,'(400A)') trim(tmpstr300) - + ! ------------------------------------- - + end subroutine collect_obs - - + + ! ************************************************************** subroutine put_into_Observations( this_obs_param, N_obs_max, N_catd, l2f, & tmp_obs, tmp_std_obs, tmp_lon, tmp_lat, tmp_time, tmp_assim, & obs_count, Observations ) - + ! Put one type of observations into the general "Observations" vector: - ! throw out no-data-values and keep track of which components have + ! throw out no-data-values and keep track of which components have ! been filled ("obs_count") ! All observations (except no-data-values) are included at this stage ! regardless of whether they will be assimilated or whether only ! innovations will be computed. ! - ! added "l2f" index vector that maps "local" tilenum to tilenum + ! added "l2f" index vector that maps "local" tilenum to tilenum ! within "full" domain ! - reichle, 17 Oct 2011 ! ! reichle, 31 Jan 2014: added "tmp_time" ! reichle, 6 Jun 2016: added flag "tmp_assim" to facilitate retaining unscaled obs - + implicit none - + type(obs_param_type), intent(in) :: this_obs_param - - integer, intent(in) :: N_obs_max, N_catd + + integer, intent(in) :: N_obs_max, N_catd integer, dimension(N_catd), intent(in) :: l2f @@ -10203,184 +10210,184 @@ subroutine put_into_Observations( this_obs_param, N_obs_max, N_catd, l2f, & real*8, dimension(N_catd), intent(in) :: tmp_time logical, dimension(N_catd), intent(in) :: tmp_assim - + integer, intent(inout) :: obs_count - - type(obs_type), dimension(N_obs_max), intent(inout) :: Observations - + + type(obs_type), dimension(N_obs_max), intent(inout) :: Observations + ! ------------------------------------------------------ integer :: i - + real :: nodatavalue, tol - + ! ------------------------------------------------------ - + nodatavalue = this_obs_param%nodata - + tol = abs(nodatavalue*nodata_tolfrac_generic) - + do i=1,N_catd if (abs(tmp_obs(i)-nodatavalue) > tol) then ! check for no-data-value - - obs_count = obs_count+1 ! augment observation counter - + + obs_count = obs_count+1 ! augment observation counter + Observations(obs_count)%obs = tmp_obs(i) - + ! check if std has been set already; if not, use default value ! (no-data-value for std is any negative value) - + if (tmp_std_obs(i) > .0) then - + Observations(obs_count)%obsvar = tmp_std_obs(i)**2 - + else - + Observations(obs_count)%obsvar = this_obs_param%errstd**2 - + end if - + Observations(obs_count)%tilenum = l2f(i) Observations(obs_count)%time = tmp_time(i) Observations(obs_count)%lat = tmp_lat(i) Observations(obs_count)%lon = tmp_lon(i) - - Observations(obs_count)%species = this_obs_param%species - + + Observations(obs_count)%species = this_obs_param%species + Observations(obs_count)%assim = this_obs_param%assim .and. tmp_assim(i) - + Observations(obs_count)%fcst = this_obs_param%nodata Observations(obs_count)%fcstvar = this_obs_param%nodata - + Observations(obs_count)%ana = this_obs_param%nodata Observations(obs_count)%anavar = this_obs_param%nodata - + end if end do - + end subroutine put_into_Observations - - + + ! ***************************************************************** - + subroutine add_to_obslog( & date_time_string, obs_param_descr, subroutine_name, num_obs_string, file_name ) - + implicit none - + character( *), intent(in) :: date_time_string ! format: YYYYMMDD_HHMMSSz character( *), intent(in) :: obs_param_descr character( *), intent(in) :: subroutine_name character( *), intent(in) :: num_obs_string character(*), intent(in) :: file_name - - ! obslog file format (comma-separated values; CSV): + + ! obslog file format (comma-separated values; CSV): ! ! analysis time, obs species descriptor, subroutine name, obs_count, file name - + write (unitnum_obslog,'(500A)') & date_time_string // ', ' // & trim(obs_param_descr) // ', ' // & trim(subroutine_name) // ', ' // & num_obs_string // ', ' // & - trim(file_name) - + trim(file_name) + end subroutine add_to_obslog - + ! ***************************************************************** - + subroutine get_tile_num_for_obs(N_catd, tile_coord, & tile_grid, N_tile_in_cell_ij, tile_num_in_cell_ij, N_latlon, lat, lon, & this_obs_param, & tile_num, & shift_lat, shift_lon ) - - ! find one tile for each obs that "administers" the obs + + ! find one tile for each obs that "administers" the obs ! ! designed to work only for a vector of obs from a *single* species ! (to be called from within obs readers) ! ! - reichle, 2015/02/20 - + implicit none - + integer, intent(in) :: N_catd, N_latlon - + type(tile_coord_type), dimension(:), pointer :: tile_coord ! input - - type(grid_def_type), intent(in) :: tile_grid - + + type(grid_def_type), intent(in) :: tile_grid + integer, dimension(tile_grid%N_lon,tile_grid%N_lat), intent(in) :: N_tile_in_cell_ij - + integer, dimension(:,:,:), pointer :: tile_num_in_cell_ij ! input - + real, dimension(N_latlon), intent(in) :: lat, lon type(obs_param_type), intent(in) :: this_obs_param - + integer, dimension(N_latlon), intent(out) :: tile_num - + real, optional, intent(in) :: shift_lat real, optional, intent(in) :: shift_lon - + ! ------------------------ - + ! local variables - + real, dimension(N_latlon) :: max_dist_x ! vector [deg] real :: max_dist_y ! scalar [deg] - + character(len=*), parameter :: Iam = 'get_tile_num_for_obs' - + real, dimension(N_latlon) :: tmp_lon, tmp_lat - + ! ----------------------------------------------------------------------------- ! - ! get "max_dist" in deg lat/lon from field-of-view (FOV) + ! get "max_dist" in deg lat/lon from field-of-view (FOV) ! ! "max_dist" = Maximum distance allowed between obs lat/lon and tile com_lat/com_lon ! when searching for a tile to which the obs will be assigned. ! ! NOTE: Subroutine get_tile_num_from_latlon() computes distances in Minkowski norm. - + if ( trim(this_obs_param%FOV_units)=='deg' ) then - + max_dist_y = this_obs_param%FOV max_dist_x = this_obs_param%FOV - + elseif ( trim(this_obs_param%FOV_units)=='km' ) then ! convert from [km] (FOV) to [deg] (max_dist_*) - + call dist_km2deg( this_obs_param%FOV, N_latlon, lat, max_dist_x, max_dist_y ) - + else - + call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'unknown FOV_units') - + end if - + if (max_dist_y<0. .or. any(max_dist_x<0.)) & call ldas_abort(LDAS_GENERIC_ERROR, Iam, 'encountered negative max_dist') - + ! SMAP 36 km EASE grid readers require special accommodation: ! ! temporarily shift lat/lon of obs for computation of nearest tile to ! avoid ambiguous assignment of M09 model tile within M36 obs grid cell - ! (center of M36 grid cell is equidistant from at least two M09 model + ! (center of M36 grid cell is equidistant from at least two M09 model ! tiles) -- reichle, 23 Aug 2013 - + tmp_lat = lat tmp_lon = lon - + if (present(shift_lat)) tmp_lat = tmp_lat + shift_lat if (present(shift_lon)) tmp_lon = tmp_lon + shift_lon - + ! find tile numbers call get_tile_num_from_latlon(N_catd, tile_coord, & @@ -10388,7 +10395,7 @@ subroutine get_tile_num_for_obs(N_catd, tile_coord, & N_latlon, tmp_lat, tmp_lon, & tile_num, max_dist_x, max_dist_y ) - + end subroutine get_tile_num_for_obs ! ***************************************************************** @@ -10402,9 +10409,9 @@ end module clsm_ensupd_read_obs program test use clsm_ensupd_read_obs - + implicit none - + integer, parameter :: N_files = 3 character(200), parameter :: fpath = & @@ -10413,33 +10420,33 @@ program test 'AMSR_E_L2_Land_B01_200209050129_A.hdf', & 'AMSR_E_L2_Land_B01_200209091008_D.hdf', & 'AMSR_E_L2_Land_B01_200209092002_D.hdf' /) - + integer :: N_data - + real, dimension(:), pointer :: lon, lat, ae_l2_sm - + character(300), dimension(N_files) :: infiles - + integer :: i - + do i=1,N_files - + infiles(i) = trim(fpath) // trim(fname(i)) - + end do - + call read_ae_l2_sm_hdf(N_files, infiles, N_data, lon, lat, ae_l2_sm ) - + write (*,*) 'N_data = ', N_data - + do i=1,3 write (*,*) i, lon(i), lat(i), ae_l2_sm(i) end do - + do i=N_data-2,N_data write (*,*) i, lon(i), lat(i), ae_l2_sm(i) end do - + end program test #endif From b1fab51a5a5d4cef453ba0145a895232a79581fd Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 2 Oct 2024 08:38:17 -0400 Subject: [PATCH 2/2] if0 out more hdf code --- .../clsm_ensupd_read_obs.F90 | 25 ++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 index f7eade5..17640e6 100644 --- a/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 +++ b/GEOSlandassim_GridComp/clsm_ensupd_read_obs.F90 @@ -5125,6 +5125,11 @@ end subroutine read_obs_SMOS ! ***************************************************************** +! NOTE: This code uses HDF4. HDF4 currently does not support +! building Fortran interfaces due to bugs. For now we +! comment out the code that uses HDF4 but keep it in +! the code base for future reference. +#if 0 subroutine read_obs_MODIS_SCF( & date_time, dtstep_assim, N_catd, tile_coord, & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & @@ -5576,6 +5581,7 @@ subroutine read_obs_MODIS_SCF( & end subroutine read_obs_MODIS_SCF +#endif ! ******************************************************************************************************* @@ -5650,6 +5656,11 @@ end subroutine localtime2longitude ! ***************************************************************** +! NOTE: This code uses HDF4. HDF4 currently does not support +! building Fortran interfaces due to bugs. For now we +! comment out the code that uses HDF4 but keep it in +! the code base for future reference. +#if 0 subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & N_good_data, CMG_lon, CMG_lat, CMG_SCF ) @@ -6031,6 +6042,7 @@ subroutine read_MODIS_SCF_hdf( fname, lon_min, lon_max, lat_min, lat_max, & deallocate(Snow_Spatial_QA) end subroutine read_MODIS_SCF_hdf +#endif ! ***************************************************************** @@ -8459,6 +8471,11 @@ subroutine read_obs( & select case (trim(this_obs_param%descr)) +! NOTE: This code uses HDF4. HDF4 currently does not support +! building Fortran interfaces due to bugs. For now we +! comment out the code that uses HDF4 but keep it in +! the code base for future reference. +#if 0 case ('ae_l2_sm_a', 'ae_l2_sm_d') call read_obs_ae_l2_sm( & @@ -8478,7 +8495,7 @@ subroutine read_obs( & tmp_obs, tmp_std_obs ) end if - +#endif case ('ae_sm_LPRM_a_C', 'ae_sm_LPRM_d_C', 'ae_sm_LPRM_a_X', 'ae_sm_LPRM_d_X' ) @@ -8706,6 +8723,11 @@ subroutine read_obs( & end if +! NOTE: This code uses HDF4. HDF4 currently does not support +! building Fortran interfaces due to bugs. For now we +! comment out the code that uses HDF4 but keep it in +! the code base for future reference. +#if 0 case ('MOD10C1','MYD10C1') call read_obs_MODIS_SCF( & @@ -8713,6 +8735,7 @@ subroutine read_obs( & tile_grid_d, N_tile_in_cell_ij, tile_num_in_cell_ij, & this_obs_param, & found_obs, tmp_obs, tmp_std_obs, tmp_lon, tmp_lat ) +#endif case('SMAP_L1C_Tbh_A', 'SMAP_L1C_Tbv_A', & 'SMAP_L1C_Tbh_D', 'SMAP_L1C_Tbv_D', &