diff --git a/ldt/USAFSI/LDT_bratsethMod.F90 b/ldt/USAFSI/LDT_bratsethMod.F90 index 60046b249..dc2ccba0b 100644 --- a/ldt/USAFSI/LDT_bratsethMod.F90 +++ b/ldt/USAFSI/LDT_bratsethMod.F90 @@ -19,6 +19,8 @@ ! errors. Plus, added two new QC tests based on ! Brasnett. Plus, other bug fixes. ! 02 Nov 2020 Eric Kemp Removed blacklist at request of 557WW. +! 28 Jul 2023 Eric Kemp Expanded station ID to 31 characters. +! 24 Aug 2023 Eric Kemp Expanded station ID to 32 characters. ! ! DESCRIPTION: ! Source code for snow depth analysis using the Bratseth objective analysis @@ -78,7 +80,7 @@ module LDT_bratsethMod real :: back_err_v_corr_len ! Background error vert correlation length integer :: nobs ! Number of observations character*10, allocatable :: networks(:) ! Network name - character*10, allocatable :: platforms(:) ! Observation station ID + character*32, allocatable :: platforms(:) ! Observation station ID real, allocatable :: obs(:) ! Observed values real, allocatable :: lats(:) ! Latitude of observation (deg N) real, allocatable :: lons(:) ! Longitude of observation (deg E) @@ -229,8 +231,8 @@ end function LDT_bratseth_count_good_obs ! Append a single observation into a LDT_bratseth_t object. Value of ! interpolated background value is optional (useful for adding ! "superobservations") - subroutine LDT_bratseth_append_ob(this, network, platform, ob, lat, lon, & - elev, ob_err_var, back) + subroutine LDT_bratseth_append_ob(this, network, platform, ob, & + lat, lon, elev, ob_err_var, back) ! Imports use LDT_logmod, only : LDT_logunit @@ -241,7 +243,7 @@ subroutine LDT_bratseth_append_ob(this, network, platform, ob, lat, lon, & ! Arguments class(LDT_bratseth_t),intent(inout) :: this character(len=10), intent(in) :: network - character(len=10), intent(in) :: platform + character(len=32), intent(in) :: platform real, intent(in) :: ob real, intent(in) :: lat real, intent(in) :: lon @@ -838,7 +840,8 @@ subroutine LDT_bratseth_run_dup_qc(this) integer :: count_dups real :: mean, back, newlat, newlon, newelev real :: ob_err_var, ob_err_corr_len - character(len=10) :: network, platform + character(len=10) :: network + character(len=32) :: platform real :: diff integer :: i,j logical :: reject_all @@ -1239,8 +1242,8 @@ end subroutine LDT_bratseth_run_water_qc ! obs, and rejected if deviation is too high. Superobs are considered ! "close" if they are in the same LIS grid box. Based on Lespinas et al ! (2015). - subroutine LDT_bratseth_run_superstat_qc(this, n, new_name, ncols, nrows, & - silent_rejects) + subroutine LDT_bratseth_run_superstat_qc(this, n, new_name, & + ncols, nrows, silent_rejects) ! Imports use LDT_coreMod, only: LDT_domain, LDT_rc @@ -1253,7 +1256,7 @@ subroutine LDT_bratseth_run_superstat_qc(this, n, new_name, ncols, nrows, & ! Arguments class(LDT_bratseth_t), intent(inout) :: this integer, intent(in) :: n - character(len=10), intent(in) :: new_name + character(len=32), intent(in) :: new_name integer, intent(in) :: ncols integer, intent(in) :: nrows logical, optional, intent(in) :: silent_rejects @@ -1262,7 +1265,8 @@ subroutine LDT_bratseth_run_superstat_qc(this, n, new_name, ncols, nrows, & integer :: num_good_obs integer :: num_rejected_obs, num_merged_obs, num_superobs logical :: silent_rejects_local - character(len=10) :: network_new, platform_new + character(len=10) :: network_new + character(len=32) :: platform_new integer, allocatable :: actions(:) real, allocatable :: means(:,:) real, allocatable :: superobs(:,:), superbacks(:,:) @@ -1817,7 +1821,7 @@ subroutine LDT_bratseth_resort_bad_obs(this) integer :: iob,job logical :: found_replacement character*10 :: network - character*10 :: platform + character*32 :: platform real :: ob, lat, lon, elev, ob_err_var, back integer :: qc character*80 :: failed_qc_test @@ -2042,6 +2046,9 @@ end subroutine LDT_bratseth_run_missing_elev_qc ! See https://en.wikipedia.org/wiki/Comb_sort subroutine LDT_bratseth_sort_obs_by_id(this) + ! Imports + use LDT_logmod, only : LDT_logunit + ! Defaults implicit none @@ -2052,6 +2059,7 @@ subroutine LDT_bratseth_sort_obs_by_id(this) integer :: gap integer :: switch character*10 :: ctemp10 + character*32 :: ctemp32 character*80 :: ctemp80 real :: rtemp integer :: itemp @@ -2078,9 +2086,9 @@ subroutine LDT_bratseth_sort_obs_by_id(this) this%networks(i) = this%networks(j) this%networks(j) = ctemp10 - ctemp10 = this%platforms(i) + ctemp32 = this%platforms(i) this%platforms(i) = this%platforms(j) - this%platforms(j) = ctemp10 + this%platforms(j) = ctemp32 rtemp = this%obs(i) this%obs(i) = this%obs(j) @@ -2152,7 +2160,7 @@ subroutine LDT_bratseth_print_snowdepths(this,minprt) int(this%elevs(i)), this%obs(i), trim(this%reject_reasons(i)) end if end do ! i -7000 format (/,'[INFO] NETID = ',A5,' STATION ID = ',A9, & +7000 format (/,'[INFO] NETID = ',A5,' STATION ID = ',A32, & ' LAT = ',F6.2,' LON = ',F7.2, & ' ELEV(M) = ',I5,' DEPTH(M) = ', F7.5, & ' QC VERDICT = ',A) diff --git a/ldt/USAFSI/LDT_usafsiMod.F90 b/ldt/USAFSI/LDT_usafsiMod.F90 index 57451c4ea..66f67d849 100644 --- a/ldt/USAFSI/LDT_usafsiMod.F90 +++ b/ldt/USAFSI/LDT_usafsiMod.F90 @@ -27,6 +27,7 @@ module LDT_usafsiMod character*10 :: date10 character*255 :: fracdir character*255 :: modif + integer :: sfcobsfmt ! EMK 20230727 character*255 :: sfcobs character*255 :: ssmis character*255 :: gmi !kyh20201118 @@ -139,11 +140,20 @@ subroutine LDT_usafsiInit() call ESMF_ConfigGetAttribute(LDT_config, usafsi_settings%modif, rc=rc) call LDT_verify(rc, trim(cfg_entry)//" not specified") + ! New sfcobs format...EMK 20230728 + cfg_entry = "USAFSI surface obs data format:" + call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + call ESMF_ConfigGetAttribute(LDT_config, usafsi_settings%sfcobsfmt, & + rc=rc) + call LDT_verify(rc, trim(cfg_entry)//" not specified") + ! Get sfcobs cfg_entry = "USAFSI surface obs data directory:" call ESMF_ConfigFindLabel(LDT_config, trim(cfg_entry), rc=rc) call LDT_verify(rc, trim(cfg_entry)//" not specified") - call ESMF_ConfigGetAttribute(LDT_config, usafsi_settings%sfcobs, rc=rc) + call ESMF_ConfigGetAttribute(LDT_config, usafsi_settings%sfcobs, & + rc=rc) call LDT_verify(rc, trim(cfg_entry)//" not specified") !------------------------------------------------------------------------------kyh20201118 diff --git a/ldt/USAFSI/USAFSI_analysisMod.F90 b/ldt/USAFSI/USAFSI_analysisMod.F90 index 89fdcb6ca..d1d01ab14 100644 --- a/ldt/USAFSI/USAFSI_analysisMod.F90 +++ b/ldt/USAFSI/USAFSI_analysisMod.F90 @@ -17,6 +17,8 @@ ! 02 Nov 2020 Eric Kemp Removed blacklist code at request of 557WW. ! 22 Jan 2021 Yeosang Yoon Add subroutine for new 0.1 deg snow climatology ! 13 Jan 2022 Eric Kemp Added support for GRIB1 FNMOC SST file. +! 28 Jul 2023 Eric Kemp Added support for new sfcobs file format (longer +! station names. ! ! DESCRIPTION: ! Source code for Air Force snow depth analysis. @@ -676,7 +678,7 @@ subroutine getgeo (month, static, nc, nr, elevations) end subroutine getgeo subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & - stalat, stalon, staelv, stadep) + stalat, stalon, staelv, stadep, sfcobsfmt) !******************************************************************************* !******************************************************************************* @@ -741,6 +743,8 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & !** 21 Mar 19 Ported to LDT...Eric Kemp, NASA GSFC/SSAI !** 09 May 19 Renamed LDTSI...Eric Kemp, NASA GSFC/SSAI !** 13 Dec 19 Renamed USAFSI...Eric Kemp, NASA GSFC/SSAI + !** 27 Jul 23 Added new sfcobs file format...Eric Kemp, SSAI + !** 24 Aug 23 New global sfcsno file format...Eric Kemp, SSAI !** !******************************************************************************* !******************************************************************************* @@ -760,12 +764,14 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & integer, intent(in) :: month ! CURRENT MONTH (1-12) character*255, intent(in) :: sfcobs ! PATH TO DBPULL SNOW OBS DIRECTORY character*5, intent(out) :: netid (:) ! NETWORK ID OF AN OBSERVATION - character*9, intent(out) :: staid (:) ! STATION ID OF AN OBSERVATION + character*32, intent(out) :: staid (:) ! STATION ID OF AN OBSERVATION + integer, intent(out) :: stacnt ! TOTAL NUMBER OF OBSERVATIONS USED integer, intent(out) :: stalat (:) ! LATITUDE OF A STATION OBSERVATION integer, intent(out) :: stalon (:) ! LONGITUDE OF A STATION OBSERVATION integer, intent(out) :: staelv (:) ! ELEVATION OF A STATION OBSERVATION (METERS) real, intent(out) :: stadep (:) ! SNOW DEPTH REPORTED AT A STATION (METERS) + integer, intent(in) :: sfcobsfmt ! Format of sfcobs file ! Local variables character*7 :: access_type ! FILE ACCESS TYPE @@ -779,9 +785,10 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & character*90 :: message (msglns) ! ERROR MESSAGE character*255 :: obsfile ! NAME OF OBSERVATION TEXT FILE character*5 :: obsnet ! RETURNED OBS STATION NETWORK - character*9 :: obssta ! RETURNED OBS STATION ID + character*32 :: obssta ! RETURNED OBS STATION ID character*5, allocatable :: oldnet (:) ! ARRAY OF NETWORKS FOR OLDSTA - character*9, allocatable :: oldsta (:) ! ARRAY OF PROCESSED STATIONS WITH SNOW DEPTHS + character*32, allocatable :: oldsta (:) ! ARRAY OF PROCESSED STATIONS WITH SNOW DEPTHS + character*12 :: routine_name ! NAME OF THIS SUBROUTINE integer :: ctrgrd ! TEMP HOLDER FOR GROUND OBS INFO integer :: ctrtmp ! TEMP HOLDER FOR TOO WARM TEMPERATURE OBS @@ -855,13 +862,18 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & message = ' ' ! OPEN INPUT FILE. - obsfile = trim(sfcobs) // 'sfcsno_' // chemi(hemi) // & - interval // date10 // '.txt' + if (sfcobsfmt == 1) then + obsfile = trim(sfcobs) // 'sfcsno_' // chemi(hemi) // & + interval // date10 // '.txt' + else if (sfcobsfmt == 2) then ! Global file + obsfile = trim(sfcobs) // 'sfcsno_' // & + '06hr_' // date10 // '.txt' + end if inquire (file=obsfile, exist=isfile) file_check : if (isfile) then access_type = 'OPENING' - open (lunsrc(hemi), file=obsfile, iostat=istat, err=5000, & + open (lunsrc(hemi), file=obsfile, iostat=istat, err=5000, & form='formatted') isopen = .true. @@ -873,17 +885,33 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & ! LOOP THROUGH ALL OBSERVATIONS RETRIEVED FROM THE DATABASE. read_loop : do while (istat .eq. 0) - read (lunsrc(hemi), 6400, iostat=istat, end=3000, err=5000) & - date10_hourly, obsnet, obssta, oblat, oblon, obelev, & - itemp, depth, ground - + if (sfcobsfmt == 1) then + read (lunsrc(hemi), 6400, iostat=istat, end=3000, & + err=5000) & + date10_hourly, obsnet, obssta, oblat, oblon, & + obelev, & + itemp, depth, ground + else if (sfcobsfmt == 2) then + ! New format with longer station IDs + read (lunsrc(hemi), 6401, iostat=istat, end=3000, & + err=5000) & + date10_hourly, obsnet, obssta, oblat, oblon, & + obelev, & + itemp, depth, ground + end if good_read : if (istat == 0) then if (date10_hourly .ne. date10_prev) then if (totalobs > 1) then - write(ldt_logunit,6500) & - trim(routine_name), chemicap(hemi), & - date10_prev, obsrtn + if (sfcobsfmt == 1) then + write(ldt_logunit,6500) & + trim(routine_name), chemicap(hemi), & + date10_prev, obsrtn + else + write(ldt_logunit,6501) & + trim(routine_name), & + date10_prev, obsrtn + end if obsrtn = 0 end if end if @@ -924,7 +952,6 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & NETID(STCTP1) = OBSNET staid(stctp1) = obssta - stadep(stctp1) = (float (depth) / 1000.0) ! convert from mm to meters if (depth >= 1 .and. stadep(stctp1) < 0.001) then @@ -1057,38 +1084,64 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & if (totalobs > 0) then stacnt_h = stacnt - stacnt_h - write (ldt_logunit,6500) trim(routine_name), chemicap(hemi), & - date10_prev, obsrtn - - write (ldt_logunit,6800) trim(routine_name), chemicap(hemi), & - totalobs, & - stacnt_h, obwsno, ctrgrd, ctrtmp, ctrtrs - + if (sfcobsfmt == 1) then + write (ldt_logunit,6500) trim(routine_name), & + chemicap(hemi), & + date10_prev, obsrtn + write (ldt_logunit,6800) trim(routine_name), & + chemicap(hemi), & + totalobs, & + stacnt_h, obwsno, ctrgrd, ctrtmp, ctrtrs + else if (sfcobsfmt == 2) then + write (ldt_logunit,6501) trim(routine_name), & + date10_prev, obsrtn + write (ldt_logunit,6801) trim(routine_name), & + totalobs, & + stacnt_h, obwsno, ctrgrd, ctrtmp, ctrtrs + end if else - message(1) = & - '[WARN] NO SURFACE OBSERVATIONS READ FOR ' // date10 // & - ' ' // chemicap(hemi) + if (sfcobsfmt == 1) then + message(1) = & + '[WARN] NO SURFACE OBSERVATIONS READ FOR ' // & + date10 // & + ' ' // chemicap(hemi) + else if (sfcobsfmt == 2) then + message(1) = & + '[WARN] NO SURFACE OBSERVATIONS READ FOR ' // & + date10 + end if call error_message (program_name, routine_name, message) end if else file_check - message(1) = & - '[WARN] NO SURFACE OBSERVATIONS FILE FOR ' // date10 // & - ' ' // chemicap(hemi) + if (sfcobsfmt == 1) then + message(1) = & + '[WARN] NO SURFACE OBSERVATIONS FILE FOR ' // & + date10 // & + ' ' // chemicap(hemi) + else if (sfcobsfmt == 2) then + message(1) = & + '[WARN] NO SURFACE OBSERVATIONS FILE FOR ' & + // date10 + end if + message(2) = '[WARN] Looked for ' // trim(obsfile) call error_message (program_name, routine_name, message) end if file_check + ! New file format is global, so we don't need to loop again + if (sfcobsfmt == 2) exit + end do hemi_loop ! DEALLOCATE ARRAYS deallocate (oldsta) return - + ! ERROR-HANDLING SECTION. 5000 continue @@ -1103,10 +1156,14 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & 6000 format (/, '[INFO] ', A, ': READING ', A) !6200 format (I) 6400 format (A10, 1X, A5, 1X, A10, 6(I10)) +!6401 format (A10, 1X, A5, 1X, A31, 1X, 6(I10)) +6401 format (A10, 1X, A5, 1X, A32, 6(I10)) 6500 format (/, '[INFO] ', A6, ': SURFACE OBS READ FOR ', A2, ' DTG ', & A10, ' = ', I6) +6501 format (/, '[INFO] ', A6, ': SURFACE OBS READ FOR DTG ', & + A10, ' = ', I6) 6600 format (1X, '**', A6, ': DEPTH = ', I6, ' STADEP = ', I6) -6700 format (/, 1X, '[INFO] HIGH POLAR TEMP: NETW= ', A5, 1X, 'STN= ', A9, & +6700 format (/, 1X, '[INFO] HIGH POLAR TEMP: NETW= ', A5, 1X, 'STN= ', A31, & 1X, 'LAT= ', F8.2, 1X, 'LON= ', F8.2, & 1X, 'ELEV= ', I5, /, 6X, 'TEMP= ', F7.1, & 2X, 'ST OF GRND= ', I9, 2X, 'DEPTH(CM)= ', I6) @@ -1119,6 +1176,15 @@ subroutine getobs (date10, month, sfcobs, netid, staid, stacnt, & /, 5X, '[INFO] OBS NOT USED FOR SEASON AND ELEVATION = ', I4, & /, 5X, '[INFO] OBS NOT USED FOR EXCEEDED TEMP THRESH = ', I6, & /, 1X, 55('-')) +6801 format (/, 1X, 55('-'), & + /, 3X, '[INFO] SUBROUTINE: ', A6, & + /, 5X, '[INFO] TOTAL SURFACE OBS READ = ', I6, & + /, 5X, '[INFO] TOTAL NON-DUPLICATE OBS PROCESSED = ', I6, & + /, 5X, '[INFO] STATIONS WITH A FOUR-THREE GROUP = ', I4, & + /, 5X, '[INFO] OBS NOT USED FOR STATE OF GROUND = ', I4, & + /, 5X, '[INFO] OBS NOT USED FOR SEASON AND ELEVATION = ', I4, & + /, 5X, '[INFO] OBS NOT USED FOR EXCEEDED TEMP THRESH = ', I6, & + /, 1X, 55('-')) end subroutine getobs @@ -2780,7 +2846,7 @@ subroutine run_snow_analysis_noglacier(runcycle, nc, nr, landmask, & integer :: snomask_reject_count integer :: bad_back_count, glacier_zone_count real :: ob_value - character*10 :: new_name + character*32 :: new_name integer :: gindex real :: rlat diff --git a/ldt/USAFSI/USAFSI_run.F90 b/ldt/USAFSI/USAFSI_run.F90 index 770f57932..8d9df8ef8 100644 --- a/ldt/USAFSI/USAFSI_run.F90 +++ b/ldt/USAFSI/USAFSI_run.F90 @@ -71,6 +71,9 @@ subroutine USAFSI_run(n) !** 13 Jan 22 Added support for FNMOC SST GRIB1 file.........Eric Kemp/NASA GSFC/SSAI !** 27 Jun 23 Removed LDT_endrun for normal termination, to avoid error ! code 1.........................................Eric Kemp/SSAI + !** 28 Jun 23 Extended station names to 31 characters........Eric Kemp/SSAI + !** 24 Aug 23 Changed station names to 32 characters.........Eric Kemp/SSAI + !***************************************************************************************** !***************************************************************************************** @@ -109,8 +112,10 @@ subroutine USAFSI_run(n) character*5, allocatable :: netid (:) ! NETWORK ID OF AN OBSERVATION character*255 :: modif ! PATH TO MODIFIED DATA DIRECTORY character*255 :: sfcobs ! PATH TO DBPULL SNOW OBS DIRECTORY + integer :: sfcobsfmt ! Format of sfcobs file character*255 :: TB_product_path ! TB_based retrivals path !kyh20201118 - character*9, allocatable :: staid (:) ! STATION ID OF AN OBSERVATION + !character*9, allocatable :: staid (:) ! STATION ID OF AN OBSERVATION + character*32, allocatable :: staid (:) ! STATION ID OF AN OBSERVATION character*255 :: static ! STATIC FILE DIRECTORY PATH character*255 :: stmpdir ! SFC TEMP DIRECTORY PATH character*255 :: sstdir ! EMK 20220113 @@ -131,7 +136,8 @@ subroutine USAFSI_run(n) real, allocatable :: stadep (:) ! OBSERVATION SNOW DEPTH (METERS) character*12 :: routine_name type(LDT_bratseth_t) :: bratseth - character*10 :: network10, platform10 + character*10 :: network10 + character*32 :: platform32 real :: rob, rlat, rlon, relev integer :: nc,nr real, allocatable :: landmask(:,:) @@ -178,6 +184,8 @@ subroutine USAFSI_run(n) fracdir = trim(usafsi_settings%fracdir) modif = trim(usafsi_settings%modif) sfcobs = trim(usafsi_settings%sfcobs) + sfcobsfmt = usafsi_settings%sfcobsfmt ! EMK test + !---------------------------------------------------------kyh20201118 if (usafsi_settings%TB_option == 1) then !SSMIS TB_product_path = trim(usafsi_settings%ssmis) @@ -369,8 +377,16 @@ subroutine USAFSI_run(n) allocate (stalat (maxsobs)) allocate (stalon (maxsobs)) - call getobs (date10, month, sfcobs, netid, staid, stacnt, & - stalat, stalon, staelv, stadep) + if (sfcobsfmt == 1 .or. sfcobsfmt == 2) then + call getobs (date10, month, sfcobs, netid, staid, stacnt, & + stalat, stalon, staelv, stadep, sfcobsfmt) + else + write(LDT_logunit,*)'[ERR] Invalid sfcobs file format!' + write(LDT_logunit,*)'[ERR] Expected 1 (old) or 2 (new)' + write(LDT_logunit,*)'[ERR] Received ', sfcobsfmt + call LDT_endrun() + end if + write(LDT_logunit,*) & '[INFO] TOTAL OBSERVATIONS RETURNED FROM GETOBS: ', & stacnt @@ -378,12 +394,13 @@ subroutine USAFSI_run(n) ! EMK Copy observations to bratseth object do j = 1, stacnt network10 = trim(netid(j)) - platform10 = trim(staid(j)) + platform32 = trim(staid(j)) rob = stadep(j) rlat = real(stalat(j)) * 0.01 rlon = real(stalon(j)) * 0.01 relev = real(staelv(j)) - call bratseth%append_ob(network10, platform10, rob, rlat, rlon,& + call bratseth%append_ob(network10, platform32, rob, & + rlat, rlon,& relev, & usafsi_settings%ob_err_var, back=-1.) end do diff --git a/lis/core/LIS_domainMod.F90 b/lis/core/LIS_domainMod.F90 index e0efcd7b9..e9be57c8d 100644 --- a/lis/core/LIS_domainMod.F90 +++ b/lis/core/LIS_domainMod.F90 @@ -397,15 +397,22 @@ subroutine LIS_domain_setup(n) TRACE_ENTER("dom_setup") allocate(LIS_domain(n)%ntiles_pergrid(LIS_rc%gnc(n)*LIS_rc%gnr(n))) + LIS_domain(n)%ntiles_pergrid = 0 ! EMK TEST allocate(LIS_domain(n)%str_tind(LIS_rc%gnc(n)*LIS_rc%gnr(n))) + LIS_domain(n)%str_tind = 0 ! EMK TEST allocate(ntiles_pergrid(LIS_rc%lnc(n)*LIS_rc%lnr(n)),stat=ierr) + ntiles_pergrid = 0 allocate(ntiles_pergrid_red(LIS_rc%lnc_red(n)*LIS_rc%lnr_red(n)),stat=ierr) - + ntiles_pergrid_red = 0 allocate(npatch_pergrid(LIS_rc%lnc(n)*LIS_rc%lnr(n),LIS_rc%max_model_types)) + npatch_pergrid = 0 allocate(npatch_pergrid_red(LIS_rc%lnc(n)*LIS_rc%lnr(n),LIS_rc%max_model_types)) + npatch_pergrid_red = 0 do m=1,LIS_rc%max_model_types allocate(LIS_surface(n,m)%npatch_pergrid(LIS_rc%gnc(n)*LIS_rc%gnr(n))) + LIS_surface(n,m)%npatch_pergrid = 0 allocate(LIS_surface(n,m)%str_patch_ind(LIS_rc%gnc(n)*LIS_rc%gnr(n))) + LIS_surface(n,m)%str_patch_ind = 0 enddo do t=1,LIS_rc%ntiles(n) diff --git a/lis/core/LIS_logMod.F90 b/lis/core/LIS_logMod.F90 index 5c763b624..bfda37a9a 100644 --- a/lis/core/LIS_logMod.F90 +++ b/lis/core/LIS_logMod.F90 @@ -390,7 +390,7 @@ subroutine LIS_warning(ierr,msg) if ( ierr /= 0 ) then write(LIS_logunit,*) '[WARN] ****************WARNING*********************' - write(LIS_logunit,*) '[WARN] ',msg + write(LIS_logunit,*) '[WARN] ', trim(msg) write(LIS_logunit,*) '[WARN] ****************WARNING*********************' endif diff --git a/lis/core/LIS_tbotAdjustMod.F90 b/lis/core/LIS_tbotAdjustMod.F90 index b99164e6a..c6610e358 100644 --- a/lis/core/LIS_tbotAdjustMod.F90 +++ b/lis/core/LIS_tbotAdjustMod.F90 @@ -454,9 +454,9 @@ subroutine LIS_readTmnUpdateRestart(n,ftn,wformat) real, allocatable :: tmptilen(:) if ( LIS_rc%tbot_update_lag == 0 ) then - write(LIS_logunit,*) '[WARN] dynamic deep soil temperature updating '//& + write(LIS_logunit,*) '[INFO] dynamic deep soil temperature updating '//& 'was disabled in lis.config file.' - write(LIS_logunit,*) '[WARN] returning from LIS_readTmnUpdateRestart' + write(LIS_logunit,*) '[INFO] returning from LIS_readTmnUpdateRestart' return endif diff --git a/lis/metforcing/usaf/AGRMET_cdfs2_est.F90 b/lis/metforcing/usaf/AGRMET_cdfs2_est.F90 index 0833ee06c..83634ba09 100644 --- a/lis/metforcing/usaf/AGRMET_cdfs2_est.F90 +++ b/lis/metforcing/usaf/AGRMET_cdfs2_est.F90 @@ -259,7 +259,7 @@ subroutine AGRMET_cdfs2_est( n,k, cliprc, clippd,& agrmet_struc(n)%clouddir,agrmet_struc(n)%use_timestamp,hemi,& yr,mo,da,hr) - write(LIS_logunit,*)'- OPENING ', trim(ifil) + write(LIS_logunit,*)'[INFO] OPENING ', trim(ifil) open(9, file=trim(ifil), access='direct', & recl=icdfs2*jcdfs2*1, iostat=istat, status="old") @@ -268,7 +268,7 @@ subroutine AGRMET_cdfs2_est( n,k, cliprc, clippd,& ! ------------------------------------------------------------------ if (istat /= 0) then - write(LIS_logunit,*)'- ERROR OPENING ',trim(ifil),istat + write(LIS_logunit,*)'[ERR] ERROR OPENING ',trim(ifil),istat error_flag = .true. cycle READ_DATA end if @@ -277,12 +277,12 @@ subroutine AGRMET_cdfs2_est( n,k, cliprc, clippd,& ! bad read, cycle loop and try to use previous hour's data. ! ------------------------------------------------------------------ - write(LIS_logunit,*)'- READING ',trim(ifil) + write(LIS_logunit,*)'[INFO] READING ',trim(ifil) read(9, rec=1, iostat=istat) totalc(hemi,:,:) close(9) if (istat /= 0) then - write(LIS_logunit,*)'- ERROR READING ',trim(ifil) + write(LIS_logunit,*)'[ERR] ERROR READING ',trim(ifil) error_flag = .true. cycle READ_DATA end if @@ -294,7 +294,7 @@ subroutine AGRMET_cdfs2_est( n,k, cliprc, clippd,& call agrmet_cdfs_pixltime_filename(ifil,agrmet_struc(n)%agrmetdir,& agrmet_struc(n)%clouddir,agrmet_struc(n)%use_timestamp,hemi,& yr,mo,da,hr) - write(LIS_logunit,*)'- OPENING ', trim(ifil) + write(LIS_logunit,*)'[INFO] OPENING ', trim(ifil) open(9, file=trim(ifil), access='direct', & recl=icdfs2*jcdfs2*4, iostat=istat, status="old") @@ -304,12 +304,12 @@ subroutine AGRMET_cdfs2_est( n,k, cliprc, clippd,& ! ------------------------------------------------------------------ if (istat /= 0) then - write(LIS_logunit,*)'- ERROR OPENING ',trim(ifil) + write(LIS_logunit,*)'[ERR] Cannot open ',trim(ifil) error_flag = .true. cycle READ_DATA end if - write(LIS_logunit,*)'- READING ', trim(ifil) + write(LIS_logunit,*)'[INFO] READING ', trim(ifil) read(9, rec=1, iostat=istat) times(hemi, :,:) close(9) @@ -320,7 +320,7 @@ subroutine AGRMET_cdfs2_est( n,k, cliprc, clippd,& ! ------------------------------------------------------------------ if (istat /= 0) then - write(LIS_logunit,*)'- ERROR READING ',trim(ifil) + write(LIS_logunit,*)'[ERR] ERROR READING ',trim(ifil) error_flag = .true. cycle READ_DATA else @@ -339,9 +339,9 @@ subroutine AGRMET_cdfs2_est( n,k, cliprc, clippd,& if (istat /= 0) then write(LIS_logunit,*) - write(LIS_logunit,*) '* IN ROUTINE AGRMET_CDFS2_EST: ERRORS WITH CDFS2 DATA.' - write(LIS_logunit,*) '* DATA IS CORRUPT OR DOES NOT EXIST.' - write(LIS_logunit,*) '* CDFS2 PRECIP ESTIMATE WILL NOT BE PERFORMED.' + write(LIS_logunit,*) '[WARN] IN ROUTINE AGRMET_CDFS2_EST: ERRORS WITH CDFS2 DATA.' + write(LIS_logunit,*) '[WARN] DATA IS CORRUPT OR DOES NOT EXIST.' + write(LIS_logunit,*) '[WARN] CDFS2 PRECIP ESTIMATE WILL NOT BE PERFORMED.' write(LIS_logunit,*) call AGRMET_julhr_date10 ( j3hr, date10 ) diff --git a/lis/metforcing/usaf/AGRMET_fldbld_galwem.F90 b/lis/metforcing/usaf/AGRMET_fldbld_galwem.F90 index 43e717a43..41f51ac97 100644 --- a/lis/metforcing/usaf/AGRMET_fldbld_galwem.F90 +++ b/lis/metforcing/usaf/AGRMET_fldbld_galwem.F90 @@ -1127,7 +1127,7 @@ subroutine galwem_reset_interp_input(n, findex, gridDesci) agrmet_struc(n)%fg_galwem_interp = LIS_rc%met_interp(findex) - write(LIS_logunit,*) 'MSG: The GALWEM forcing resolution is coarser ' // & + write(LIS_logunit,*) '[INFO] The GALWEM forcing resolution is coarser ' // & 'than the running domain.' write(LIS_logunit,*) ' Interpolating with the ' // & trim(agrmet_struc(n)%fg_galwem_interp) // ' method.' @@ -1184,7 +1184,7 @@ subroutine galwem_reset_interp_input(n, findex, gridDesci) elseif ( howtoTransform == 'neighbor') then agrmet_struc(n)%fg_galwem_interp = 'neighbor' - write(LIS_logunit,*) 'MSG: The GALWEM forcing resolution is comparable ' // & + write(LIS_logunit,*) '[INFO] The GALWEM forcing resolution is comparable ' // & 'to the running domain.' write(LIS_logunit,*) ' Interpolating with the ' // & trim(agrmet_struc(n)%fg_galwem_interp) // ' method.' @@ -1195,7 +1195,7 @@ subroutine galwem_reset_interp_input(n, findex, gridDesci) elseif ( howtoTransform == 'upscale' ) then agrmet_struc(n)%fg_galwem_interp = LIS_rc%met_upscale(findex) - write(LIS_logunit,*) 'MSG: The GALWEM forcing resolution is finer ' // & + write(LIS_logunit,*) '[INFO] The GALWEM forcing resolution is finer ' // & 'than the running domain.' write(LIS_logunit,*) ' Upscaling with the ' // & trim(agrmet_struc(n)%fg_galwem_interp) // ' method.' diff --git a/lis/metforcing/usaf/AGRMET_fldbld_gfs.F90 b/lis/metforcing/usaf/AGRMET_fldbld_gfs.F90 index 2667559f8..eb7704c16 100644 --- a/lis/metforcing/usaf/AGRMET_fldbld_gfs.F90 +++ b/lis/metforcing/usaf/AGRMET_fldbld_gfs.F90 @@ -1587,7 +1587,7 @@ subroutine gfs_reset_interp_input(n, findex, gridDesci) agrmet_struc(n)%fg_gfs_interp = LIS_rc%met_interp(findex) - write(LIS_logunit,*) 'MSG: The GFS forcing resolution is coarser ' // & + write(LIS_logunit,*) '[INFO] The GFS forcing resolution is coarser ' // & 'than the running domain.' write(LIS_logunit,*) ' Interpolating with the ' // & trim(agrmet_struc(n)%fg_gfs_interp) // ' method.' @@ -1644,7 +1644,7 @@ subroutine gfs_reset_interp_input(n, findex, gridDesci) elseif ( howtoTransform == 'neighbor') then agrmet_struc(n)%fg_gfs_interp = 'neighbor' - write(LIS_logunit,*) 'MSG: The GFS forcing resolution is comparable ' // & + write(LIS_logunit,*) '[INFO] The GFS forcing resolution is comparable ' // & 'to the running domain.' write(LIS_logunit,*) ' Interpolating with the ' // & trim(agrmet_struc(n)%fg_gfs_interp) // ' method.' @@ -1655,7 +1655,7 @@ subroutine gfs_reset_interp_input(n, findex, gridDesci) elseif ( howtoTransform == 'upscale' ) then agrmet_struc(n)%fg_gfs_interp = LIS_rc%met_upscale(findex) - write(LIS_logunit,*) 'MSG: The GFS forcing resolution is finer ' // & + write(LIS_logunit,*) '[INFO] The GFS forcing resolution is finer ' // & 'than the running domain.' write(LIS_logunit,*) ' Upscaling with the ' // & trim(agrmet_struc(n)%fg_gfs_interp) // ' method.' diff --git a/lis/metforcing/usaf/AGRMET_fldbld_precip_galwem.F90 b/lis/metforcing/usaf/AGRMET_fldbld_precip_galwem.F90 index 6882775eb..c641846ac 100644 --- a/lis/metforcing/usaf/AGRMET_fldbld_precip_galwem.F90 +++ b/lis/metforcing/usaf/AGRMET_fldbld_precip_galwem.F90 @@ -550,7 +550,7 @@ subroutine AGRMET_fldbld_read_precip_galwem(fg_filename, ifguess, jfguess,& count_prec = 0 write(LIS_logunit,*)' ' - write(LIS_logunit,*)'[MSG] Reading first guess precip ', trim(fg_filename) + write(LIS_logunit,*)'[INFO] Reading first guess precip ', trim(fg_filename) call grib_count_in_file(ftn,nvars,ierr) call LIS_verify(ierr, 'error in grib_count_in_file in ' // & diff --git a/lis/metforcing/usaf/AGRMET_fldbld_precip_gfs.F90 b/lis/metforcing/usaf/AGRMET_fldbld_precip_gfs.F90 index 474212953..9953b8622 100644 --- a/lis/metforcing/usaf/AGRMET_fldbld_precip_gfs.F90 +++ b/lis/metforcing/usaf/AGRMET_fldbld_precip_gfs.F90 @@ -635,7 +635,7 @@ subroutine AGRMET_fldbld_read_precip_gfs( fg_filename, ifguess, jfguess,& allocate ( dum1d (ifguess*jfguess) ) count_prec = 0 - write(LIS_logunit,*)'[MSG] Reading first guess precip ', trim(fg_filename) + write(LIS_logunit,*)'[INFO] Reading first guess precip ', trim(fg_filename) call grib_count_in_file(ftn,nvars,ierr) call LIS_verify(ierr, 'error in grib_count_in_file in AGRMET_fldbld_read') diff --git a/lis/metforcing/usaf/AGRMET_forcingMod.F90 b/lis/metforcing/usaf/AGRMET_forcingMod.F90 index 220ca1563..3f242e639 100644 --- a/lis/metforcing/usaf/AGRMET_forcingMod.F90 +++ b/lis/metforcing/usaf/AGRMET_forcingMod.F90 @@ -354,6 +354,8 @@ module AGRMET_forcingMod real*8 :: agrmetpcptime1, agrmetpcptime2 real*8 :: cmortime integer :: pcpobswch + integer :: pcpobsfmt ! EMK...File format for precip obs + integer :: sfcobsfmt ! EMK...File format for sfcobs integer :: pwswch integer :: raswch integer :: cdfs2swch @@ -2075,26 +2077,41 @@ subroutine init_AGRMET(findex) agrmet_struc(n)%pcpclimoAlarmTime = 0.0 allocate(agrmet_struc(n)%cliprc1(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%cliprc1 = 0 allocate(agrmet_struc(n)%clippd1(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%clippd1 = 0 allocate(agrmet_struc(n)%clirtn1(LIS_rc%lnc(n), LIS_rc%lnr(n))) - + agrmet_struc(n)%clirtn1 = 0 allocate(agrmet_struc(n)%cliprc2(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%cliprc2 = 0 allocate(agrmet_struc(n)%clippd2(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%clippd2 = 0 allocate(agrmet_struc(n)%clirtn2(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%clirtn2 = 0 allocate(agrmet_struc(n)%cliprc(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%cliprc = 0 allocate(agrmet_struc(n)%clippd(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%clippd = 0 allocate(agrmet_struc(n)%clirtn(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%clirtn = 0 allocate(agrmet_struc(n)%sfcprs(6,LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%sfcprs = 0 allocate(agrmet_struc(n)%sfcrlh(6,LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%sfcrlh = 0 allocate(agrmet_struc(n)%sfcspd(6,LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%sfcspd = 0 allocate(agrmet_struc(n)%sfctmp(6,LIS_rc%lnc(n), LIS_rc%lnr(n))) - + agrmet_struc(n)%sfctmp = 0 allocate(agrmet_struc(n)%lasprs(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%lasprs = 0 allocate(agrmet_struc(n)%lasrlh(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%lasrlh = 0 allocate(agrmet_struc(n)%lasspd(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%lasspd = 0 allocate(agrmet_struc(n)%lastmp(LIS_rc%lnc(n), LIS_rc%lnr(n))) + agrmet_struc(n)%lastmp = 0 ! EMK BEGIN ! diff --git a/lis/metforcing/usaf/AGRMET_geoest.F90 b/lis/metforcing/usaf/AGRMET_geoest.F90 index 1ffff203a..051726fc0 100644 --- a/lis/metforcing/usaf/AGRMET_geoest.F90 +++ b/lis/metforcing/usaf/AGRMET_geoest.F90 @@ -189,9 +189,9 @@ subroutine AGRMET_geoest( n, j3hr, land ,gest, grnk, quad9r, & inquire( file = trim(ifil), exist = exists) if ( .not. exists ) then - write(LIS_logunit,*) 'AGRMET_geoest/precip: error opening file ',trim(ifil) - write(LIS_logunit,*) ' file does not exist' - write(LIS_logunit,*) ' geo precip estimate will not be performed' + write(LIS_logunit,*) '[WARN] AGRMET_geoest/precip: error opening file ',trim(ifil) + write(LIS_logunit,*) '[WARN] file does not exist' + write(LIS_logunit,*) '[WARN] geo precip estimate will not be performed' message = ' ' message(1) = 'program: LIS' message(2) = ' routine: AGRMET_geoest' @@ -205,7 +205,7 @@ subroutine AGRMET_geoest( n, j3hr, land ,gest, grnk, quad9r, & endif return endif - write(LIS_logunit,*) '- READING ', trim(ifil) + write(LIS_logunit,*) '[INFO] READING ', trim(ifil) call LIS_putget( geoprc, 'r', ifil, routine_name, & imax, jmax ) @@ -235,14 +235,14 @@ subroutine AGRMET_geoest( n, j3hr, land ,gest, grnk, quad9r, & if ( exists ) then gdgeornk = .true. - write(LIS_logunit,*) '- READING RANK FILE', trim(ifil) + write(LIS_logunit,*) '[INFO] READING RANK FILE', trim(ifil) call LIS_putget( geornk, 'r', ifil, routine_name, & imax, jmax ) else - write(LIS_logunit,*)'AGRMET_GEOEST/PRECIP: ERROR OPENING FILE ',trim(ifil) - write(LIS_logunit,*)' FILE DOES NOT EXIST' - write(LIS_logunit,*)' DEFAULT VALUES OF 4 WILL BE USED FOR GEO RANKS' + write(LIS_logunit,*)'[WARN] AGRMET_GEOEST/PRECIP: ERROR OPENING FILE ',trim(ifil) + write(LIS_logunit,*)'[WARN] FILE DOES NOT EXIST' + write(LIS_logunit,*)'[WARN] DEFAULT VALUES OF 4 WILL BE USED FOR GEO RANKS' message = ' ' message(1)='program: LIS' message(2)=' routine: AGRMET_geoest' @@ -273,8 +273,8 @@ subroutine AGRMET_geoest( n, j3hr, land ,gest, grnk, quad9r, & ! of these cases. if ( is_geo_corrupted(geoprc, imax, jmax, mo, hemi) ) then - write(LIS_logunit,*) 'AGRMET_geoest/precip: data corrupted - ',trim(ifil) - write(LIS_logunit,*) ' geo precip estimate will not be performed' + write(LIS_logunit,*) '[WARN] AGRMET_geoest/precip: data corrupted - ',trim(ifil) + write(LIS_logunit,*) '[WARN] geo precip estimate will not be performed' message = ' ' message(1) = 'program: LIS' message(2) = ' routine: AGRMET_geoest' @@ -330,10 +330,10 @@ subroutine AGRMET_geoest( n, j3hr, land ,gest, grnk, quad9r, & grnk_tmp(hemi,i,j) = 4 write(LIS_logunit,*)' ' write(LIS_logunit,*)'--------------------------------------------' - write(LIS_logunit,*)'Bad geo precip rank value' - write(LIS_logunit,*)'This value will be set to a default of 4' - write(LIS_logunit,*)'geornk(',i,',',j,') = ',geornk(i,j) - write(LIS_logunit,*)'from file :' + write(LIS_logunit,*)'[WARN] Bad geo precip rank value' + write(LIS_logunit,*)'[WARN] This value will be set to a default of 4' + write(LIS_logunit,*)'[WARN] geornk(',i,',',j,') = ',geornk(i,j) + write(LIS_logunit,*)'[WARN] from file :' write(LIS_logunit,*) trim(ifil) write(LIS_logunit,*)'--------------------------------------------' write(LIS_logunit,*)' ' diff --git a/lis/metforcing/usaf/AGRMET_getcli.F90 b/lis/metforcing/usaf/AGRMET_getcli.F90 index 2e4c8e6d5..a9eb8d0a8 100644 --- a/lis/metforcing/usaf/AGRMET_getcli.F90 +++ b/lis/metforcing/usaf/AGRMET_getcli.F90 @@ -74,7 +74,7 @@ subroutine AGRMET_getcli(n, filename,rtn,clidat) inquire( file = trim(filename), exist = exists) if ( exists ) then - write(LIS_logunit,*) '- READING ', trim(filename) + write(LIS_logunit,*) '[INFO] READING ', trim(filename) ! ------------------------------------------------------------------ ! the rtneph climo files are real valued files, while the precip diff --git a/lis/metforcing/usaf/AGRMET_getpcpobs.F90 b/lis/metforcing/usaf/AGRMET_getpcpobs.F90 index 5ff482c9f..bd60dd8d6 100644 --- a/lis/metforcing/usaf/AGRMET_getpcpobs.F90 +++ b/lis/metforcing/usaf/AGRMET_getpcpobs.F90 @@ -297,9 +297,9 @@ subroutine AGRMET_getpcpobs(n, j6hr, month, prcpwe, & message(6) = ' Observations beyond array size will be ignored' message(7) = ' Increase number of AGRMET maximum precip obs in lis.config file!' if (LIS_masterproc) then + alert_number = alert_number + 1 call LIS_alert('LIS.AGRMET_getpcpobs', & alert_number, message) - alert_number = alert_number + 1 end if nsize = agrmet_struc(n)%max_pcpobs @@ -397,9 +397,9 @@ subroutine AGRMET_getpcpobs(n, j6hr, month, prcpwe, & message(2) = ' Routine: AGRMET_getpcpobs' message(3) = ' Problem reading '// trim(filename) if (LIS_masterproc) then + alert_number = alert_number + 1 call LIS_alert('LIS.AGRMET_getpcpobs', & alert_number, message) - alert_number = alert_number + 1 end if end if else @@ -416,9 +416,9 @@ subroutine AGRMET_getpcpobs(n, j6hr, month, prcpwe, & message(2) = ' Routine: AGRMET_getpcpobs' message(3) = ' Missing rain gage file '// trim(filename) if (LIS_masterproc) then + alert_number = alert_number + 1 call LIS_alert('LIS.AGRMET_getpcpobs', & alert_number, message) - alert_number = alert_number + 1 end if endif diff --git a/lis/metforcing/usaf/AGRMET_getsfc.F90 b/lis/metforcing/usaf/AGRMET_getsfc.F90 index 641c25506..8a6eb9560 100644 --- a/lis/metforcing/usaf/AGRMET_getsfc.F90 +++ b/lis/metforcing/usaf/AGRMET_getsfc.F90 @@ -40,20 +40,23 @@ ! 25 Jun 20 Modified to check valid times of surface obs. ! ..............................Eric Kemp/NASA/SSAI ! 16 Dec 21 Replaced julhr with YYYYMMDDHH in log...Eric Kemp/NASA/SSAI +! 07 Nov 23 Added new sfcobs format.......Eric Kemp/NASA/SSAI ! ! !INTERFACE: subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & ri, rj, obstmp, obsrlh, obsspd, obscnt, & isize, minwnd, alert_number, imax, jmax, rootdir,cdmsdir,& - use_timestamp) + use_timestamp, use_wigos_sfcobs) ! !USES: + use ESMF ! EMK Patch for DTG check use AGRMET_forcingMod, only : agrmet_struc use LIS_coreMod, only : LIS_domain, LIS_masterproc use LIS_timeMgrMod, only : LIS_julhr_date - use LIS_logMod, only : LIS_logunit, LIS_alert + use LIS_logMod, only : LIS_logunit, LIS_alert, & + LIS_getNextUnitNumber, LIS_releaseUnitNumber use map_utils, only : latlon_to_ij use USAF_bratsethMod, only: USAF_ObsData, USAF_assignObsData - use ESMF ! EMK Patch for DTG check + implicit none ! !ARGUMENTS: @@ -76,6 +79,7 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & character(len=*) :: rootdir character(len=*) :: cdmsdir integer, intent(in) :: use_timestamp + logical, intent(in) :: use_wigos_sfcobs ! ! !DESCRIPTION: ! @@ -220,11 +224,14 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & character*10 :: date10 ! EMK replace cjulhr in log character*14, allocatable :: dtg ( : ) character*100 :: message ( 20 ) - character*8, allocatable :: netyp ( : ) + !character*8, allocatable :: netyp ( : ) + character*9, allocatable :: netyp ( : ) character*8 :: norsou ( 2 ) ! character*8, allocatable :: platform ( : ) - character*9, allocatable :: platform ( : ) ! EMK BUG FIX - character*8, allocatable :: rptyp ( : ) + !character*9, allocatable :: platform ( : ) ! EMK BUG FIX + character*32, allocatable :: platform ( : ) ! EMK WIGOS + !character*8, allocatable :: rptyp ( : ) + character*9, allocatable :: rptyp ( : ) ! WIGOS logical, allocatable :: skip (:) ! EMK real :: rlat real :: rlon @@ -244,9 +251,19 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & integer :: rc real, external :: AGRMET_calcrh_dpt - character(len=10) :: net10, platform10 + !character(len=10) :: net10, platform10 + character(len=8) :: netyp8 + character(len=9) :: platform9 + character(len=32) :: net32, platform32 + character(len=8) :: rptyp8 + character*32, parameter :: blank32 = " " + integer :: iunit + logical :: found_file + data norsou / 'NORTHERN', 'SOUTHERN' / + message = '' + ! ------------------------------------------------------------------ ! Executable code begins here. Intialize observation counter to 0. ! ------------------------------------------------------------------ @@ -292,20 +309,78 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & skip(:) = .false. ! EMK call getsfcobsfilename(sfcobsfile, rootdir, cdmsdir, & - use_timestamp,hemi, yr, mo, da, hr) - - write(LIS_logunit,*)'Reading OBS: ',trim(sfcobsfile) - open(22,file=trim(sfcobsfile),status='old', iostat=ierr1) - if(ierr1.eq.0) then - read(22,*, iostat=ierr2) nsize - + use_timestamp,hemi, yr, mo, da, hr, use_wigos_sfcobs) + + inquire(file=trim(sfcobsfile), exist=found_file) + if (.not. found_file) then + write(LIS_logunit,*) '[WARN] Cannot find ', trim(sfcobsfile) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: AGRMET_getsfc' + message(3) = ' Cannot find file ' // trim(sfcobsfile) + if (LIS_masterproc) then + call LIS_alert('LIS.AGRMET_getsfc', & + alert_number, message) + alert_number = alert_number + 1 + end if + message = '' + if (use_wigos_sfcobs) exit ! These files are global + cycle + end if + + write(LIS_logunit,*)'[INFO] Opening: ',trim(sfcobsfile) + iunit = LIS_getNextUnitNumber() + open(iunit,file=trim(sfcobsfile),status='old', iostat=ierr1) + + if (ierr1 .ne. 0) then + write(date10,'(i4, i2.2, i2.2, i2.2)', iostat=istat1) & + yr, mo, da, hr + write(LIS_logunit,*)' ' + write(LIS_logunit,*) & + '[WARN] ROUTINE AGRMET_GETSFC: ERROR OPENING ', & + trim(sfcobsfile) + write(LIS_logunit,*)'[WARN] ISTAT IS ', ierr1 + message(1) = 'program: LIS' + message(2) = ' routine: AGRMET_getsfc' + message(3) = ' error opening file ' // trim(sfcobsfile) + alert_number = alert_number + 1 + if(LIS_masterproc) then + call lis_alert( 'LIS.AGRMET_getsfc', alert_number, message ) + endif + message = '' + if (use_wigos_sfcobs) exit ! New WIGOS version is global + cycle + endif + + if(ierr1.eq.0) then + read(iunit,*, iostat=ierr2) nsize + + if (ierr2 .ne. 0) then + write(LIS_logunit,*) & + '[WARN] Problem reading total report count from ', & + trim(sfcobsfile) + message(1) = 'program: LIS' + message(2) = ' routine: AGRMET_getsfc' + message(3) = ' Cannot read total report count from ' // & + trim(sfcobsfile) + alert_number = alert_number + 1 + if(LIS_masterproc) then + call lis_alert( 'LIS.AGRMET_getsfc', alert_number, & + message ) + endif + message = '' + close(iunit) + call LIS_releaseUnitNumber(iunit) + if (use_wigos_sfcobs) exit + cycle + end if + ! ------------------------------------------------------------------ ! If the number of obs in the file is greater than the array size ! write an alert to the log and set back the number to read to ! prevent a segfault. ! ------------------------------------------------------------------ - if ( nsize .GT. isize ) then + if ( nsize .GT. isize ) then write(LIS_logunit,*)' ' write(LIS_logunit,*)"******************************************************" @@ -321,15 +396,32 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & end if do i=1,nsize - read(22,*, iostat=ierr3) platform(i), dtg(i), & - itmp(i), idpt(i), irelh(i), ilat(i), ilon(i), & - ispd(i), rptyp(i), netyp(i) + if (.not. use_wigos_sfcobs) then + !read(iunit, *, iostat=ierr3) platform(i), dtg(i), & + ! itmp(i), idpt(i), irelh(i), ilat(i), ilon(i), & + ! ispd(i), rptyp(i), netyp(i) + read(iunit, *, iostat=ierr3) platform9, dtg(i), & + itmp(i), idpt(i), irelh(i), ilat(i), ilon(i), & + ispd(i), rptyp8, netyp8 + if (ierr3 == 0) then + platform(i) = platform9 + rptyp(i) = rptyp8 + netyp(i) = netyp8 + end if + else + read(iunit, 1000, iostat=ierr3) platform(i), dtg(i), & + itmp(i), idpt(i), irelh(i), ilat(i), ilon(i), & + ispd(i), rptyp(i), netyp(i) +1000 format(a32, 1x, a14, 1x, i9, 1x, i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, i9, 1x, a9, 1x, a9) + end if !if (ierr3 .ne. 0) skip(i) = .true. ! EMK ! EMK Patch Skip report if problem occurred reading it if (ierr3 .ne. 0) then write(LIS_logunit,*) & - '[WARN] Problem reading report ', i, ' from file' + '[WARN] Problem reading report ', i, & + ' from sfcobs file, skipping line' skip(i) = .true. cycle end if @@ -366,7 +458,8 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & end if enddo - close(22) + close(iunit) + call LIS_releaseUnitNumber(iunit) !if(ierr2.eq.0.and.ierr3.eq.0) then ! EMK Patch...Allow observation storage even if problem occurred @@ -399,20 +492,21 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & ! ------------------------------------------------------------------ ! Only process the observations for the current hemisphere. ! ------------------------------------------------------------------ - + + ! EMK...Use all obs if we are using a WIGOS global file if( ((hemi .eq. 1) .and. (rlat .ge. 0.0)) .or. & - ((hemi .eq. 2) .and. (rlat .lt. 0.0)) ) then + ((hemi .eq. 2) .and. (rlat .lt. 0.0)) .or. & + use_wigos_sfcobs) then ! ------------------------------------------------------------------ ! Convert point's lat/lon to i/j coordinates. ! check for values of rigrid and rjgrid outside of ! the AGRMET model grid, or in the other hemisphere. ! ------------------------------------------------------------------ - call latlon_to_ij(LIS_domain(n)%lisproj, rlat, rlon, & - rigrid, rjgrid) - -! if(rigrid.ge.1.and.rigrid.le.imax.and. & + rigrid, rjgrid) + + !if(rigrid.ge.1.and.rigrid.le.imax.and. & ! rjgrid.ge.1.and.rjgrid.le.jmax) then ! EMK TEST if (.true.) then @@ -423,16 +517,16 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & ! observations had to be discared for this reason. ! ------------------------------------------------------------------ - if ( ( (itmp(irecord) .gt. -99999998) .and. & - (ispd(irecord) .gt. -99999998) ) .and. & - ( (irelh(irecord) .gt. -99999998) .or. & - (idpt(irecord) .gt. -99999998) ) ) then - - ! EMK...Make sure dew point .le. temperature - if ( (itmp(irecord) .gt. -99999998) .and. & - (idpt(irecord) .gt. -99999998) ) then - if (idpt(irecord) .gt. itmp(irecord)) cycle - end if + if ( ( (itmp(irecord) .gt. -99999998) .and. & + (ispd(irecord) .gt. -99999998) ) .and. & + ( (irelh(irecord) .gt. -99999998) .or. & + (idpt(irecord) .gt. -99999998) ) ) then + + ! EMK...Make sure dew point .le. temperature + if ( (itmp(irecord) .gt. -99999998) .and. & + (idpt(irecord) .gt. -99999998) ) then + if (idpt(irecord) .gt. itmp(irecord)) cycle + end if ! ------------------------------------------------------------------ ! Gross error check and set temperature. @@ -442,14 +536,14 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & ! so the Barnes scheme will ignore it. ! ------------------------------------------------------------------ - if (itmp(irecord) .gt. -99999998) then - rtmp = float (itmp(irecord)) / 100.0 - if ( (rtmp .lt. 200.0) .or. (rtmp .gt. 350.0) ) then + if (itmp(irecord) .gt. -99999998) then + rtmp = float (itmp(irecord)) / 100.0 + if ( (rtmp .lt. 200.0) .or. (rtmp .gt. 350.0) ) then + rtmp = -1.0 + end if + else rtmp = -1.0 end if - else - rtmp = -1.0 - end if ! ------------------------------------------------------------------ ! Gross error check and set relative humidity values. @@ -458,21 +552,21 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & ! value of -1, so the Barnes scheme will ignore it. ! ------------------------------------------------------------------ - if (irelh(irecord) .gt. -99999998) then - rrelh = float(irelh(irecord)) / 10000.0 - if ( (rrelh .gt. 1.0) .or. (rrelh .lt. 0.0) ) then - rrelh = -1.0 - end if - else if ( (idpt(irecord) .gt. -99999998) .and. & - (rtmp .ne. -1.0) ) then - rdpt = float(idpt(irecord)) / 100.0 - rrelh = AGRMET_calcrh_dpt(rtmp, rdpt) - if ( (rrelh .gt. 1.0) .or. (rrelh .lt. 0.0) ) then + if (irelh(irecord) .gt. -99999998) then + rrelh = float(irelh(irecord)) / 10000.0 + if ( (rrelh .gt. 1.0) .or. (rrelh .lt. 0.0) ) then + rrelh = -1.0 + end if + else if ( (idpt(irecord) .gt. -99999998) .and. & + (rtmp .ne. -1.0) ) then + rdpt = float(idpt(irecord)) / 100.0 + rrelh = AGRMET_calcrh_dpt(rtmp, rdpt) + if ( (rrelh .gt. 1.0) .or. (rrelh .lt. 0.0) ) then + rrelh = -1.0 + end if + else rrelh = -1.0 end if - else - rrelh = -1.0 - end if ! ------------------------------------------------------------------ ! Constrain surface wind speed to a value between @@ -481,149 +575,116 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & ! of -1, so the Barnes scheme will ignore it. ! ------------------------------------------------------------------ - if (ispd(irecord) .gt. -99999998) then - rspd = float(ispd(irecord)) / 10.0 - rspd = max( min( rspd, 75.0 ), minwnd ) - else - rspd = -1.0 - end if - + if (ispd(irecord) .gt. -99999998) then + rspd = float(ispd(irecord)) / 10.0 + rspd = max( min( rspd, 75.0 ), minwnd ) + else + rspd = -1.0 + end if + ! ------------------------------------------------------------------ ! If all the data is out of range and/or missing, don't store ! this observation. ! ------------------------------------------------------------------ - if ( (nint(rtmp) .gt. 0) .and. & - (nint(rrelh) .gt. 0) .and.& - (nint(rspd) .gt. 0) ) then - + if ( (nint(rtmp) .gt. 0) .and. & + (nint(rrelh) .gt. 0) .and.& + (nint(rspd) .gt. 0) ) then + ! ------------------------------------------------------------------ ! If we make it here, the observation probably has some good ! data, so increment the observation counter, store the ! data into the observation arrays. ! ------------------------------------------------------------------ - ! EMK...Add to data structures. Handle reformated - ! CDMS data that is missing platform and network - if (rtmp .gt. 0) then - net10 = trim(netyp(irecord)) - platform10 = trim(platform(irecord)) - if (trim(net10) .eq. 'NULL') then - net10 = 'CDMS' - end if - if (trim(platform10) .eq. '-99999999') then - platform10 = '00000000' - end if - call USAF_assignObsData(t2mObs,net10, & - platform10,rtmp,rlat,rlon, & - agrmet_struc(n)%bratseth_t2m_stn_sigma_o_sqr,& - 0.) + ! EMK...Add to data structures. Handle reformated + ! CDMS data that is missing platform and network + if (rtmp .gt. 0) then + net32 = blank32 + net32 = netyp(irecord) + platform32 = blank32 + platform32 = platform(irecord) + if (trim(net32) .eq. 'NULL') then + net32 = 'CDMS' + end if + if (trim(platform32) .eq. '-99999999') then + platform32 = '00000000' + end if + call USAF_assignObsData(t2mObs,net32, & + platform32,rtmp,rlat,rlon, & + agrmet_struc(n)%bratseth_t2m_stn_sigma_o_sqr,& + 0.) - end if - if (rrelh .gt. 0) then - net10 = trim(netyp(irecord)) - platform10 = trim(platform(irecord)) - if (trim(net10) .eq. 'NULL') then - net10 = 'CDMS' end if - if (trim(platform10) .eq. '-99999999') then - platform10 = '00000000' + if (rrelh .gt. 0) then + net32 = blank32 + net32 = netyp(irecord) + platform32 = blank32 + platform32 = platform(irecord) + if (trim(net32) .eq. 'NULL') then + net32 = 'CDMS' + end if + if (trim(platform32) .eq. '-99999999') then + platform32 = '00000000' + end if + call USAF_assignObsData(rh2mObs,net32, & + platform32,rrelh,rlat,rlon, & + agrmet_struc(n)%bratseth_t2m_stn_sigma_o_sqr, & + 0.) end if - call USAF_assignObsData(rh2mObs,net10, & - platform10,rrelh,rlat,rlon, & - agrmet_struc(n)%bratseth_t2m_stn_sigma_o_sqr, & - 0.) - end if - if (rspd .gt. 0) then - net10 = trim(netyp(irecord)) - platform10 = trim(platform(irecord)) - if (trim(net10) .eq. 'NULL') then - net10 = 'CDMS' + if (rspd .gt. 0) then + net32 = blank32 + net32 = netyp(irecord) + platform32 = blank32 + platform32 = platform(irecord) + if (trim(net32) .eq. 'NULL') then + net32 = 'CDMS' + end if + if (trim(platform32) .eq. '-99999999') then + platform32 = '00000000' + end if + call USAF_assignObsData(spd10mObs,net32, & + platform32,rspd,rlat,rlon, & + agrmet_struc(n)%bratseth_spd10m_stn_sigma_o_sqr, & + 0.) end if - if (trim(platform10) .eq. '-99999999') then - platform10 = '00000000' - end if - call USAF_assignObsData(spd10mObs,net10, & - platform10,rspd,rlat,rlon, & - agrmet_struc(n)%bratseth_spd10m_stn_sigma_o_sqr, & - 0.) - end if + obscnt = obscnt + 1 - obscnt = obscnt + 1 - - ri(obscnt) = rigrid - rj(obscnt) = rjgrid - obstmp(obscnt) = rtmp - obsrlh(obscnt) = rrelh - obsspd(obscnt) = rspd - if(rspd.lt.0) then - print*, 'probl ',rlat,rlon, rspd + ri(obscnt) = rigrid + rj(obscnt) = rjgrid + obstmp(obscnt) = rtmp + obsrlh(obscnt) = rrelh + obsspd(obscnt) = rspd endif - endif ! ------------------------------------------------------------------ ! If we have reached the number of obs that our hard-wired ! arrays can hold, then exit the loop. ! ------------------------------------------------------------------ - if ( obscnt .eq. isize ) exit + if ( obscnt .eq. isize ) then + write(LIS_logunit,*) & + '[WARN] ROUTINE GETSFC: REACHED MAXIMUM NUMBER OF SFC OBS TO SAVE IN MEMORY' + message(1) = 'program: LIS' + message(2) = ' routine: AGRMET_getsfc' + message(3) = ' reached maximum number of sfc obs to save memory in memory' + alert_number = alert_number + 1 + if(LIS_masterproc) then + call lis_alert( 'LIS.AGRMET_getsfc', & + alert_number, & + message ) + endif + message = '' + exit ! Jump out of hemi do loop + end if + endif endif endif - endif - enddo - else -! ------------------------------------------------------------------ -! There was an error retrieving obs for this Julhr and hemi. -! Send an alert message, but don't abort. -! ------------------------------------------------------------------ - - ! EMK...Replace julhr with YYYYMMDD - !write(cjulhr,'(i6)',iostat=istat1) julhr - write(date10,'(i4, i2.2, i2.2, i2.2)', iostat=istat1) yr, mo, da, hr - write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- ROUTINE GETSFC: ERROR RETRIEVING SFC OBS FOR' - write(LIS_logunit,*)'- THE '//norsou(hemi)//' HEMISPHERE.' - write(LIS_logunit,*)'- ISTAT IS ', ierr1 - message(1) = 'program: LIS' - message(2) = ' routine: AGRMET_getsfc' - message(3) = ' error retrieving sfc obs from database for' - message(4) = ' the '//norsou(hemi)//' hemisphere.' - if( istat1 .eq. 0 ) then - ! EMK...Replace julhr with YYYYMMDD - !write(LIS_logunit,*)'- JULHR IS ' // cjulhr - !message(5) = ' julhr is ' // trim(cjulhr) // '.' - write(LIS_logunit,*)' - YYYYMMDDHH is ' // date10 - message(5) = ' yyyymmddhh is ' // trim(date10) // '.' - endif - alert_number = alert_number + 1 - if(LIS_masterproc) then - call lis_alert( 'sfcalc ', alert_number, message ) - endif - endif - else - ! EMK...Replace julhr with YYYYMMDD - !write(cjulhr,'(i6)',iostat=istat1) julhr - write(date10,'(i4, i2.2, i2.2, i2.2)', iostat=istat1) yr, mo, da, hr - write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- ROUTINE AGRMET_GETSFC: ERROR RETRIEVING SFC OBS FOR' - write(LIS_logunit,*)'- THE '//norsou(hemi)//' HEMISPHERE.' - write(LIS_logunit,*)'- ISTAT IS ', ierr1 - message(1) = 'program: LIS' - message(2) = ' routine: AGRMET_getsfc' - message(3) = ' error retrieving sfc obs from database for' - message(4) = ' the '//norsou(hemi)//' hemisphere.' - if( istat1 .eq. 0 ) then - ! EMK...Replace julhr with YYYYMMDDHH - !write(LIS_logunit,*)'- JULHR IS ' // cjulhr - !message(5) = ' julhr is ' // trim(cjulhr) // '.' - write(LIS_logunit,*)' - YYYYMMDDHH is ' // date10 - message(5) = ' yyyymmddhh is ' // trim(date10) // '.' - endif - alert_number = alert_number + 1 - if(LIS_masterproc) then - call lis_alert( 'sfcalc ', alert_number, message ) + enddo endif endif + + if (use_wigos_sfcobs) exit ! New sfcobs file is global enddo deallocate ( idpt ) @@ -638,6 +699,8 @@ subroutine AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & deallocate ( rptyp ) deallocate(skip) ! EMK + + write(LIS_logunit,*)'[INFO] Stored ', obscnt, ' sfc obs' end subroutine AGRMET_getsfc !BOP @@ -647,15 +710,18 @@ end subroutine AGRMET_getsfc ! ! !INTERFACE: subroutine getsfcobsfilename(filename, rootdir, dir, & - use_timestamp, hemi, yr,mo,da,hr) + use_timestamp, hemi, yr, mo, da, hr, use_wigos_sfcobs) implicit none + ! !ARGUMENTS: - character(len=*) :: filename - character(len=*) :: rootdir - character(len=*) :: dir - integer, intent(in) :: yr,mo,da,hr, hemi + character(len=*), intent(inout) :: filename + character(len=*), intent(in) :: rootdir + character(len=*), intent(in) :: dir integer, intent(in) :: use_timestamp + integer, intent(in) :: hemi + integer, intent(in) :: yr, mo, da, hr + logical, intent(in) :: use_wigos_sfcobs ! ! !DESCRIPTION: ! This routines generates the name of the surface obs file @@ -699,14 +765,30 @@ subroutine getsfcobsfilename(filename, rootdir, dir, & write(unit=fhemi,fmt='(a2)') 'sh' endif - if(use_timestamp.eq.1) then - filename = trim(rootdir)//'/'//trim(fyr)//trim(fmo)//trim(fda)//& - '/'//trim(dir)//'/sfcobs_'& - //trim(fhemi)//'.01hr.'//trim(fyr)//trim(fmo)//trim(fda)//trim(fhr) + if(use_timestamp.eq.1) then + if (use_wigos_sfcobs) then + filename = trim(rootdir) // '/' // & + trim(fyr) // trim(fmo) // trim(fda) // & + '/' // trim(dir) // '/sfcobs_01hr_' // & + trim(fyr) // trim(fmo) // trim(fda) // trim(fhr) // ".txt" + else + filename = trim(rootdir) // '/' // & + trim(fyr) // trim(fmo) // trim(fda) // & + '/' // trim(dir) // '/sfcobs_' // & + trim(fhemi) // '.01hr.' // & + trim(fyr) // trim(fmo) // trim(fda) // trim(fhr) + end if else - filename = trim(rootdir)//& - '/'//trim(dir)//'/sfcobs_'& - //trim(fhemi)//'.01hr.'//trim(fyr)//trim(fmo)//trim(fda)//trim(fhr) + if (use_wigos_sfcobs) then + filename = trim(rootdir) // & + '/' // trim(dir) // '/sfcobs_01hr_' // & + trim(fyr) // trim(fmo) // trim(fda) // trim(fhr) // ".txt" + else + filename = trim(rootdir) // & + '/' // trim(dir) // '/sfcobs_' // & + trim(fhemi) // '.01hr.' // & + trim(fyr) // trim(fmo) // trim(fda) // trim(fhr) + end if endif end subroutine getsfcobsfilename diff --git a/lis/metforcing/usaf/AGRMET_processobs.F90 b/lis/metforcing/usaf/AGRMET_processobs.F90 index 714a65d51..eba6bf5e5 100644 --- a/lis/metforcing/usaf/AGRMET_processobs.F90 +++ b/lis/metforcing/usaf/AGRMET_processobs.F90 @@ -289,7 +289,9 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & type(rain_obs), allocatable :: obs_cur(:) type(rain_obs), allocatable :: obs_6(:) type(rain_obs), allocatable :: obs_12(:) - + + character*32 :: net32, platform32 + data chemi / '_nh.', '_sh.' / sumsqr(a,b,c,d) = ((a-b)**2) + ((c-d)**2) @@ -486,7 +488,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & if (sixyes) then - write(LIS_logunit,*)'- READING ', trim(filename_min6) + write(LIS_logunit,*)'[INFO] READING ', trim(filename_min6) iofunc = "OPEN " open(8, file=trim(filename_min6), iostat=istat, err=100) @@ -533,9 +535,9 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & write(LIS_logunit,*)' ' write(LIS_logunit,*)"******************************************************" - write(LIS_logunit,*)"* PRECIP SAVE FILE FROM 6 HOURS AGO DOES NOT EXIST." - write(LIS_logunit,*)"* FILE NAME IS ", trim(filename_min6) - write(LIS_logunit,*)"* OBSERVATION COUNT WILL BE REDUCED." + write(LIS_logunit,*)"[WARN] PRECIP SAVE FILE FROM 6 HOURS AGO DOES NOT EXIST." + write(LIS_logunit,*)"[WARN] FILE NAME IS ", trim(filename_min6) + write(LIS_logunit,*)"[WARN] OBSERVATION COUNT WILL BE REDUCED." write(LIS_logunit,*)"******************************************************" write(LIS_logunit,*)' ' @@ -554,9 +556,9 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & write(LIS_logunit,*)' ' write(LIS_logunit,*)"*******************************************************" - write(LIS_logunit,*)"* BAD ", trim(iofunc), " ON FILE ",trim(filename_min6) - write(LIS_logunit,*)"* ISTAT = ", istat - write(LIS_logunit,*)"* OBSERVATION COUNT WILL BE REDUCED." + write(LIS_logunit,*)"[WARN] BAD ", trim(iofunc), " ON FILE ",trim(filename_min6) + write(LIS_logunit,*)"[WARN] ISTAT = ", istat + write(LIS_logunit,*)"[WARN] OBSERVATION COUNT WILL BE REDUCED." write(LIS_logunit,*)"*******************************************************" write(LIS_logunit,*)' ' @@ -598,7 +600,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & iofunc = "OPEN " open(9, file=trim(filename_min12), iostat=istat, err=200) - write(LIS_logunit,*)'- READING ', trim(filename_min12) + write(LIS_logunit,*)'[INFO] READING ', trim(filename_min12) iofunc = "READ " read(9, *, iostat=istat, err=200, end=200) count12 @@ -643,9 +645,9 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & write(LIS_logunit,*)' ' write(LIS_logunit,*)"*******************************************************" - write(LIS_logunit,*)"* PRECIP SAVE FILE FROM 12 HOURS AGO DOES NOT EXIST." - write(LIS_logunit,*)"* FILE NAME IS ", trim(filename_min12) - write(LIS_logunit,*)"* OBSERVATION COUNT WILL BE REDUCED." + write(LIS_logunit,*)"[WARN] PRECIP SAVE FILE FROM 12 HOURS AGO DOES NOT EXIST." + write(LIS_logunit,*)"[WARN] FILE NAME IS ", trim(filename_min12) + write(LIS_logunit,*)"[WARN] OBSERVATION COUNT WILL BE REDUCED." write(LIS_logunit,*)"*******************************************************" write(LIS_logunit,*)' ' @@ -664,9 +666,9 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & write(LIS_logunit,*)' ' write(LIS_logunit,*)"*******************************************************" - write(LIS_logunit,*)"* BAD ", trim(iofunc), " ON FILE ",trim(filename_min12) - write(LIS_logunit,*)"* ISTAT = ", istat - write(LIS_logunit,*)"* OBSERVATION COUNT WILL BE REDUCED." + write(LIS_logunit,*)"[WARN] BAD ", trim(iofunc), " ON FILE ",trim(filename_min12) + write(LIS_logunit,*)"[WARN] ISTAT = ", istat + write(LIS_logunit,*)"[WARN] OBSERVATION COUNT WILL BE REDUCED." write(LIS_logunit,*)"*******************************************************" write(LIS_logunit,*)' ' @@ -698,7 +700,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & if (sixyes) then write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- MAKING 12 HRLY AMTS FROM PREVIOUS 6 HRLY AMTS' + write(LIS_logunit,*)'[INFO] MAKING 12 HRLY AMTS FROM PREVIOUS 6 HRLY AMTS' MAKE12_FROM_6 : do i = 1, fltrcnt @@ -757,7 +759,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & if (twelveyes) then - write(LIS_logunit,*)'- MAKING 12 HRLY AMTS FROM CURRENT 24 AND OLD 12 HRLY AM TS' + write(LIS_logunit,*)'[INFO] MAKING 12 HRLY AMTS FROM CURRENT 24 AND OLD 12 HRLY AM TS' MAKE12_FROM_24 : do i = 1, fltrcnt @@ -1046,7 +1048,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & open (11, file=trim(filename), iostat=istat, err=300) write(LIS_logunit,*)' ' - write(LIS_logunit,*)"- WRITING ",trim(filename) + write(LIS_logunit,*)"[INFO] WRITING ",trim(filename) iofunc = "WRITE" write(11,*, iostat=istat, err=300) fltrcnt @@ -1068,7 +1070,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & close (11) write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- NUMBER OF 12 AND 6 HOURLY OBS IS ',count12obs, count6obs + write(LIS_logunit,*)'[INFO] NUMBER OF 12 AND 6 HOURLY OBS IS ',count12obs, count6obs write(LIS_logunit,*)' ' 300 continue @@ -1080,8 +1082,8 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & write(LIS_logunit,*)' ' write(LIS_logunit,*)"*******************************************************" - write(LIS_logunit,*)"* BAD ", trim(iofunc), " ON FILE ",trim(filename) - write(LIS_logunit,*)"* ISTAT = ", istat + write(LIS_logunit,*)"[WARN] BAD ", trim(iofunc), " ON FILE ",trim(filename) + write(LIS_logunit,*)"[WARN] ISTAT = ", istat write(LIS_logunit,*)"*******************************************************" write(LIS_logunit,*)' ' @@ -1116,7 +1118,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & USE_6 : if ( .not. use_twelve ) then - write(LIS_logunit,*)"- PUTTING 6 HOURLY RAIN GAUGE OBSERVATIONS ON GRID." + write(LIS_logunit,*)"[INFO] PUTTING 6 HOURLY RAIN GAUGE OBSERVATIONS ON GRID." oldd = 999.0 @@ -1142,8 +1144,10 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & obs_cur(i)%lat, obs_cur(i)%lon,ri,rj) ! EMK...Add observation + net32 = obs_cur(i)%net + platform32 = obs_cur(i)%platform call USAF_assignObsData(precip6, & - obs_cur(i)%net, obs_cur(i)%platform, & + net32, platform32, & float(obs_cur(i)%amt6) * 0.1, & obs_cur(i)%lat, obs_cur(i)%lon,& agrmet_struc(n)%bratseth_precip_gauge_sigma_o_sqr, & @@ -1192,7 +1196,7 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & USE_12 : if (use_twelve) then - write(LIS_logunit,*)"- PUTTING 12 HOURLY RAIN GAUGE OBSERVATIONS ON GRID." + write(LIS_logunit,*)"[INFO] PUTTING 12 HOURLY RAIN GAUGE OBSERVATIONS ON GRID." oldd = 999.0 @@ -1209,8 +1213,10 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & obs_cur(i)%lon,ri,rj) ! EMK...Add observation + net32 = obs_cur(i)%net + platform32 = obs_cur(i)%platform call USAF_assignObsData(precip12, & - obs_cur(i)%net, obs_cur(i)%platform, & + net32, platform32, & float(obs_cur(i)%amt12) * 0.1, & obs_cur(i)%lat, obs_cur(i)%lon, & agrmet_struc(n)%bratseth_precip_gauge_sigma_o_sqr, & @@ -1275,8 +1281,10 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & obs_cur(i)%lon,ri,rj) ! EMK...Add observation + net32 = obs_cur(i)%net + platform32 = obs_cur(i)%platform call USAF_assignObsData(precip12, & - obs_cur(i)%net, obs_cur(i)%platform, & + net32, platform32, & 0.0, & obs_cur(i)%lat, obs_cur(i)%lon, & agrmet_struc(n)%bratseth_precip_gauge_sigma_o_sqr, & @@ -1326,8 +1334,10 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & ri,rj) ! EMK...Add observation + net32 = obs_cur(i)%net + platform32 = obs_cur(i)%platform call USAF_assignObsData(precip12, & - obs_cur(i)%net, obs_cur(i)%platform, & + net32, platform32, & float(obs_cur(i)%amtmsc) * 0.1, & obs_cur(i)%lat, obs_cur(i)%lon, & agrmet_struc(n)%bratseth_precip_gauge_sigma_o_sqr, & @@ -1382,8 +1392,10 @@ subroutine AGRMET_processobs(n, obs, isize, stncnt, hemi, julhr, & ri,rj) ! EMK...Add observation + net32 = obs_6(i)%net + platform32 = obs_6(i)%platform call USAF_assignObsData(precip12, & - obs_6(i)%net, obs_6(i)%platform, & + net32, platform32, & float(obs_6(i)%amtmsc) * 0.1, & obs_6(i)%lat, obs_6(i)%lon, & agrmet_struc(n)%bratseth_precip_gauge_sigma_o_sqr, & diff --git a/lis/metforcing/usaf/AGRMET_read_sfcalccntm.F90 b/lis/metforcing/usaf/AGRMET_read_sfcalccntm.F90 index 0e6b488d8..b8700de71 100644 --- a/lis/metforcing/usaf/AGRMET_read_sfcalccntm.F90 +++ b/lis/metforcing/usaf/AGRMET_read_sfcalccntm.F90 @@ -56,8 +56,8 @@ subroutine AGRMET_read_sfcalccntm(n) inquire( file = trim(agrmet_struc(n)%sfcntmfile), exist = exists) if ( exists ) then - write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- READING ', trim(agrmet_struc(n)%sfcntmfile) + !write(LIS_logunit,*)' ' + write(LIS_logunit,*)'[INFO] READING ', trim(agrmet_struc(n)%sfcntmfile) ftn= LIS_getNextUnitNumber() open(ftn, file=trim(agrmet_struc(n)%sfcntmfile), access='direct',& status='old', form="unformatted", recl=LIS_rc%gnr(n)*LIS_rc%gnc(n)*4) @@ -73,10 +73,10 @@ subroutine AGRMET_read_sfcalccntm(n) else write(LIS_logunit,*) write(LIS_logunit,*) "*****************************************************" - write(LIS_logunit,*) "* LIS: ERROR OPENING FILE:" - write(LIS_logunit,*) "* ", trim(agrmet_struc(n)%sfcntmfile) - write(LIS_logunit,*) "* FILE DOES NOT EXIST." - write(LIS_logunit,*) "* LIS WILL ABORT." + write(LIS_logunit,*) "[ERR] LIS: ERROR OPENING FILE:" + write(LIS_logunit,*) "[ERR] ", trim(agrmet_struc(n)%sfcntmfile) + write(LIS_logunit,*) "[ERR] FILE DOES NOT EXIST." + write(LIS_logunit,*) "[ERR] LIS WILL ABORT." write(LIS_logunit,*) "*****************************************************" message = ' ' message(1) = 'program: LIS' diff --git a/lis/metforcing/usaf/AGRMET_readmask.F90 b/lis/metforcing/usaf/AGRMET_readmask.F90 index 380d2e7d4..63e3b18b7 100644 --- a/lis/metforcing/usaf/AGRMET_readmask.F90 +++ b/lis/metforcing/usaf/AGRMET_readmask.F90 @@ -75,8 +75,8 @@ subroutine AGRMET_readmask(n) endif inquire( file = trim(name), exist = exists) if ( exists ) then - write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- READING ', trim(name) + !write(LIS_logunit,*)' ' + write(LIS_logunit,*)'[INFO] READING ', trim(name) if(agrmet_struc(n)%global_or_hemi .eq. 0) then call LIS_putget( agrmet_struc(n)%land(:,:,hemi), 'r', name, routine_name, & agrmet_struc(n)%imax, agrmet_struc(n)%jmax ) @@ -87,10 +87,10 @@ subroutine AGRMET_readmask(n) else write(LIS_logunit,*) write(LIS_logunit,*) "*****************************************************" - write(LIS_logunit,*) "* LIS: ERROR OPENING FILE:" - write(LIS_logunit,*) "* ", trim(name) - write(LIS_logunit,*) "* FILE DOES NOT EXIST." - write(LIS_logunit,*) "* LIS WILL ABORT." + write(LIS_logunit,*) "[ERR] LIS: ERROR OPENING FILE:" + write(LIS_logunit,*) "[ERR] ", trim(name) + write(LIS_logunit,*) "[ERR] FILE DOES NOT EXIST." + write(LIS_logunit,*) "[ERR] LIS WILL ABORT." write(LIS_logunit,*) "*****************************************************" message = ' ' message(1) = 'program: LIS' diff --git a/lis/metforcing/usaf/AGRMET_readpcpcntm.F90 b/lis/metforcing/usaf/AGRMET_readpcpcntm.F90 index e450f2e4d..0ac8572e9 100644 --- a/lis/metforcing/usaf/AGRMET_readpcpcntm.F90 +++ b/lis/metforcing/usaf/AGRMET_readpcpcntm.F90 @@ -60,8 +60,8 @@ subroutine AGRMET_readpcpcntm(n) name = trim(agrmet_struc(n)%climodir)//'/pcp_spread_radii.1gd4r' inquire( file = trim(name) , exist = exists) if ( exists ) then - write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- READING ', trim(name) + !write(LIS_logunit,*)' ' + write(LIS_logunit,*)'[INFO] READING ', trim(name) ftn= LIS_getNextUnitNumber() open(ftn, file=trim(name), access='direct',& status='old', form="unformatted", recl=LIS_rc%gnr(n)*LIS_rc%gnc(n)*4) @@ -76,10 +76,10 @@ subroutine AGRMET_readpcpcntm(n) call LIS_releaseUnitNumber(ftn) else - write(LIS_logunit,*) 'LIS: error opening file ',trim(name) - write(LIS_logunit,*) ' file does not exist' - write(LIS_logunit,*) ' this could indicate serious data discrepancies' - write(LIS_logunit,*) ' LIS will be aborted' + write(LIS_logunit,*) '[ERR] LIS: error opening file ',trim(name) + write(LIS_logunit,*) '[ERR] file does not exist' + write(LIS_logunit,*) '[ERR] this could indicate serious data discrepancies' + write(LIS_logunit,*) '[ERR] LIS will be aborted' message(1) = 'program: LIS' message(2) = ' routine: AGRMET_readpcpcntm' message(3) = ' error opening file '//trim(name) diff --git a/lis/metforcing/usaf/AGRMET_readterrain.F90 b/lis/metforcing/usaf/AGRMET_readterrain.F90 index 859676106..358acc6d7 100644 --- a/lis/metforcing/usaf/AGRMET_readterrain.F90 +++ b/lis/metforcing/usaf/AGRMET_readterrain.F90 @@ -61,17 +61,17 @@ subroutine AGRMET_readterrain(n) call get_agrmetterrain_filename(name,agrmet_struc(n)%terrainfile,hemi) inquire( file = trim(name), exist = exists) if ( exists ) then - write(LIS_logunit,*)' ' - write(LIS_logunit,*)'- READING ', trim(name) + !write(LIS_logunit,*)' ' + write(LIS_logunit,*)'[INFO] READING ', trim(name) call LIS_putget( agrmet_struc(n)%alt(:,:,hemi), 'r', name, routine_name, & agrmet_struc(n)%imax, agrmet_struc(n)%jmax ) else write(LIS_logunit,*) write(LIS_logunit,*) "*****************************************************" - write(LIS_logunit,*) "* LIS: ERROR OPENING FILE:" - write(LIS_logunit,*) "* ", trim(name) - write(LIS_logunit,*) "* FILE DOES NOT EXIST." - write(LIS_logunit,*) "* LIS WILL ABORT." + write(LIS_logunit,*) "[ERR] LIS: ERROR OPENING FILE:" + write(LIS_logunit,*) "[ERR] ", trim(name) + write(LIS_logunit,*) "[ERR] FILE DOES NOT EXIST." + write(LIS_logunit,*) "[ERR] LIS WILL ABORT." write(LIS_logunit,*) "*****************************************************" message = ' ' message(1) = 'program: LIS' diff --git a/lis/metforcing/usaf/AGRMET_sfcalc.F90 b/lis/metforcing/usaf/AGRMET_sfcalc.F90 index cd45f7a3e..b79426223 100644 --- a/lis/metforcing/usaf/AGRMET_sfcalc.F90 +++ b/lis/metforcing/usaf/AGRMET_sfcalc.F90 @@ -174,12 +174,13 @@ subroutine AGRMET_sfcalc(n) character(len=10) :: yyyymmddhh integer :: ierr integer :: r,c, L - character(len=10) :: type + character(len=32) :: type integer :: gdeltas, gid, ntiles integer :: count1 logical :: found_inq integer :: rc integer, external :: LIS_create_subdirs + logical :: use_wigos_sfcobs data lokspd / 15, 25, 30, 40, 50 / data lokrlh / 10, 15, 25, 35, 40 / @@ -367,7 +368,7 @@ subroutine AGRMET_sfcalc(n) write(LIS_logunit,*)' ' write(LIS_logunit,*)'---------------------------- ' !write(LIS_logunit,*)'- PROCESSING-SFC JULHR ', julhr - write(LIS_logunit,*)'- PROCESSING-SFC YYYYMMDDHH ', yyyymmddhh + write(LIS_logunit,*)'[INFO] PROCESSING-SFC YYYYMMDDHH ', yyyymmddhh write(LIS_logunit,*)'---------------------------- ' ! ------------------------------------------------------------------ @@ -486,12 +487,17 @@ subroutine AGRMET_sfcalc(n) call USAF_createObsData(spd10mObs, n, & maxobs=agrmet_struc(n)%max_sfcobs) + if (agrmet_struc(n)%sfcobsfmt == 1) then + use_wigos_sfcobs = .false. + else + use_wigos_sfcobs = .true. + end if call AGRMET_getsfc( n, julhr, t2mObs, rh2mObs, spd10mObs, & ri, rj, obstmp, obsrlh, obsspd, & obscnt, agrmet_struc(n)%max_sfcobs, agrmet_struc(n)%minwnd, & alert_number, LIS_rc%lnc(n), LIS_rc%lnr(n),& agrmet_struc(n)%agrmetdir,agrmet_struc(n)%cdmsdir,& - agrmet_struc(n)%use_timestamp) + agrmet_struc(n)%use_timestamp, use_wigos_sfcobs) ! call MPI_Barrier(LIS_mpi_comm, ierr) ! stop diff --git a/lis/metforcing/usaf/AGRMET_smiest.F90 b/lis/metforcing/usaf/AGRMET_smiest.F90 index 6fc4088d6..1879e3987 100644 --- a/lis/metforcing/usaf/AGRMET_smiest.F90 +++ b/lis/metforcing/usaf/AGRMET_smiest.F90 @@ -136,9 +136,9 @@ subroutine AGRMET_smiest( n,j3hr, quad9r, ra, razero,alert_number,imax,jmax) inquire(file = trim(ifil), exist = exists) if (.not. exists) then write(LIS_logunit,*) ' ' - write(LIS_logunit,*)'precip/smiedr: error opening file' - write(LIS_logunit,*)' SSMI data file ', trim(ifil),' does not exist.' - write(LIS_logunit,*)' SSMI estimates will not be used in ',& + write(LIS_logunit,*)'[WARN] precip/smiedr: error opening file' + write(LIS_logunit,*)'[WARN] SSMI data file ', trim(ifil),' does not exist.' + write(LIS_logunit,*)'[WARN] SSMI estimates will not be used in ',& 'precip analysis.' write(LIS_logunit,*) ' ' message =' ' @@ -153,7 +153,7 @@ subroutine AGRMET_smiest( n,j3hr, quad9r, ra, razero,alert_number,imax,jmax) ra_tmp(hemi,:,:) = udef else - write(LIS_logunit,*) '- READING ',trim(ifil) + write(LIS_logunit,*) '[INFO] READING ',trim(ifil) call LIS_putget( ra_tmp(hemi,:,:), 'r', ifil, routine_name, & imax, jmax) endif @@ -161,7 +161,7 @@ subroutine AGRMET_smiest( n,j3hr, quad9r, ra, razero,alert_number,imax,jmax) enddo if ( .not. use_zeros ) then - write(LIS_logunit,*)'- SSMI ZEROS NOT USED' + write(LIS_logunit,*)'[INFO] SSMI ZEROS NOT USED' where ( ra_tmp .eq. 0.0 ) ra_tmp = udef endif diff --git a/lis/metforcing/usaf/USAF_GagesMod.F90 b/lis/metforcing/usaf/USAF_GagesMod.F90 new file mode 100644 index 000000000..c926e8916 --- /dev/null +++ b/lis/metforcing/usaf/USAF_GagesMod.F90 @@ -0,0 +1,2639 @@ +! +! MODULE: USAF_GagesMod +! +! DESCRIPTION: Contains structure and related routines for storing and +! updating new database of rain gauge reports from USAF. Includes +! routines to reconcile reports by intercomparison. +! +! AUTHOR: Eric Kemp, SSAI/NASA GSFC +! +! REFERENCES: +! Durre, I, M J Menne, B E Gleason, T G Houston, and R S Vose, 2010: +! Comprehensive automated quality assurance of daily surface observations. +! J Appl Meteor Climatol, 49, 1615-1633, +! https://doi.org/10.1175/2010JAMC2375.1 +! Qi, Y, S Martinaitis, J Zhang, and S Cocks, 2016: A real-time automated +! quality control of hourly rain gauge data based on multiple sensors in +! MRMS system. J Hydrometeor, 17, 1675-1691, +! https://doi.org/10.1175/JHM-D-15-0188.1 +! WMO, 2019: Manual on Codes, International Codes, Volume I.1, Part A -- +! Alphanumeric Codes. WMO No. 306, +! https://library.wmo.int/doc_num.php?explnum_id=10235 +! WMO, 2021: Manual on Codes, International Codes, Volume I.2, Part B -- +! Binary Codes. WMO No. 306, +! https://library.wmo.int/doc_num.php?explnum_id=10722 +! WMO, 2018: Manual on Codes, Regional Codes and National Coding Practices, +! Volume II. WMO No. 306, +! https://library.wmo.int/doc_num.php?explnum_id=5730 +! WMO, 2012: Weather Reporting, Volume A -- Observing Stations. WMO No. 9, +! https://library.wmo.int/doc_num.php?explnum_id=9896 +! +! REVISION HISTORY: +! 30 Aug 2021...First developmental version committed to git. +! 03 Aug 2023...Added 1-hr and 2-hr accumulations, date/times of each +! observation, and ob-specific past weather durations. +! 23 Aug 2023...Added both WMO and FIPS country codes. +! 07 Sep 2023...Fixed WMO block check (values in preobs files are typically +! <= 100, so these are multipied by 1000 or 10000 to check for country +! ranges). +!------------------------------------------------------------------------------ + +module USAF_GagesMod + + ! Defaults + implicit none + private + + type, public :: USAF_Gages_t + private + logical :: cdms_flag + character(10) :: date10 + integer :: nobs + character(14), allocatable :: YYYYMMDDhhmmss(:) + character(32), allocatable :: networks(:) + character(32), allocatable :: platforms(:) + character(2), allocatable :: wmocode_id(:) + character(2), allocatable :: fipscode_id(:) + integer, allocatable :: wmonumbers(:) + integer, allocatable :: bsn(:) + ! NOTE: Latitudes and longitudes are in hundreths of degrees + ! (100 = 1.00 deg) + integer, allocatable :: lats(:) + integer, allocatable :: lons(:) + ! NOTE: Accumulations are in tenths of mm (10 = 1 mm) + integer, allocatable :: amts24(:) + integer, allocatable :: amts21(:) ! Non-standard + integer, allocatable :: amts18(:) + integer, allocatable :: amts15(:) + integer, allocatable :: amts12(:) + ! "Original", non-filled or reconciled 12hr amount, needed for + ! some 12Z South America sanity checks. These are not written + ! to the presav2 file. + integer, allocatable :: amts12_orig(:) + integer, allocatable :: amts09(:) + integer, allocatable :: amts06(:) + ! "Original", non-filled or reconciled 06hr amount, needed for + ! some 12Z South America sanity checks. These are not written + ! to the presav2 file. + integer, allocatable :: amts06_orig(:) + integer, allocatable :: amts03(:) + integer, allocatable :: amts02(:) + integer, allocatable :: amts01(:) + integer, allocatable :: amts00(:) ! Non-standard accumulations + integer, allocatable :: durations(:) ! Durations of non-std accums + integer, allocatable :: preswx(:) + integer, allocatable :: pastwx(:) + integer, allocatable :: pastwx_durations(:) + character(32), allocatable :: unique_networks(:) + integer, allocatable :: firsts(:) ! Starting indices for each network + integer, allocatable :: lasts(:) ! Ending indices for each network + integer :: num_unique_networks + contains + procedure :: new => USAF_gages_new + procedure :: read_data => USAF_gages_read_data + procedure :: delete => USAF_gages_delete + procedure :: check_gross_errors => & + USAF_gages_check_gross_errors + procedure :: use_misc_precip => USAF_gages_use_misc_precip + procedure :: reconcile_self => USAF_gages_reconcile_self + procedure :: reconcile_gages01 => USAF_gages_reconcile_gages01 + procedure :: reconcile_gages02 => USAF_gages_reconcile_gages02 + procedure :: reconcile_gages03 => USAF_gages_reconcile_gages03 + procedure :: reconcile_gages06 => USAF_gages_reconcile_gages06 + procedure :: reconcile_gages09 => USAF_gages_reconcile_gages09 + procedure :: reconcile_gages12 => USAF_gages_reconcile_gages12 + procedure :: reconcile_gages15 => USAF_gages_reconcile_gages15 + procedure :: reconcile_gages18 => USAF_gages_reconcile_gages18 + procedure :: reconcile_gages21 => USAF_gages_reconcile_gages21 + procedure :: correct_region3_12Z => USAF_gages_correct_region3_12Z + procedure :: use_preswx_pastwx => USAF_gages_use_preswx_pastwx + procedure :: fill_gaps => USAF_gages_fill_gaps + procedure :: write_data => USAF_gages_write_data + procedure :: set_pastwx_durations => USAF_gages_set_pastwx_durations + procedure :: copy_to_usaf_obsdata => USAF_copy_to_usaf_obsdata + end type USAF_Gages_t + + ! Local parameters + integer, parameter :: MAX_UNIQUE_NETWORKS = 25 + integer, parameter :: MISSING = -99999999 + character(4), parameter :: MISSING_NAME = "NULL" + + ! Old CDMS replacements for WMO numbers + ! The following two are guesses + integer, parameter :: CDMS_SWEDEN_LOWLIMIT = 020010 + integer, parameter :: CDMS_SWEDEN_HIGHLIMIT = 026999 + ! The following two are guesses + integer, parameter :: CDMS_DENMARK_LOWLIMIT = 060010 + integer, parameter :: CDMS_DENMARK_HIGHLIMIT = 061999 + integer, parameter :: CDMS_RUSSIA_LOWLIMIT = 200000 + integer, parameter :: CDMS_RUSSIA_HIGHLIMIT = 390009 + integer, parameter :: CDMS_INDIA_LOWLIMIT = 420010 + integer, parameter :: CDMS_INDIA_HIGHLIMIT = 433999 + integer, parameter :: CDMS_SRILANKA_LOWLIMIT = 434000 + integer, parameter :: CDMS_SRILANKA_HIGHLIMIT = 434979 + integer, parameter :: CDMS_S_AMER_LOWLIMIT = 800000 + integer, parameter :: CDMS_S_AMER_HIGHLIMIT = 889999 + + ! JMOBS WMO Numbers. + integer, parameter, public :: JMOBS_SWEDEN_LOWLIMIT = 02001 + integer, parameter, public :: JMOBS_SWEDEN_HIGHLIMIT = 02699 + integer, parameter, public :: JMOBS_DENMARK_LOWLIMIT = 06001 + integer, parameter, public :: JMOBS_DENMARK_HIGHLIMIT = 06199 + integer, parameter, public :: JMOBS_RUSSIA_LOWLIMIT = 20000 + integer, parameter, public :: JMOBS_RUSSIA_HIGHLIMIT = 39000 + integer, parameter, public :: JMOBS_INDIA_LOWLIMIT = 42000 + integer, parameter, public :: JMOBS_INDIA_HIGHLIMIT = 43399 + integer, parameter, public :: JMOBS_SRILANKA_LOWLIMIT = 43400 + integer, parameter, public :: JMOBS_SRILANKA_HIGHLIMIT = 43497 + integer, parameter, public :: JMOBS_S_AMER_LOWLIMIT = 80000 + integer, parameter, public :: JMOBS_S_AMER_HIGHLIMIT = 88999 + + ! Maximum allowed accumulations + ! Values below are in tenths of mm (10 = 1 mm) + ! Qi et al. (2016) use 50.8 mm (2 in) for 1 hour accum, based on + ! subjective evaluation in CONUS regions with poor radar coverage; + ! they judged approximately 97% of observations above this amount were + ! erroneous. + integer, parameter :: MAX_PCP_1HR = 508 ! 2 inches + integer, parameter :: MAX_PCP_2HR = 1016 ! 4 inches + integer, parameter :: MAX_PCP_3HR = 1524 ! 6 inches + integer, parameter :: MAX_PCP_6HR = 3048 ! 12 inches + integer, parameter :: MAX_PCP_9HR = 4572 ! 18 inches + integer, parameter :: MAX_PCP_12HR = 6096 ! 24 inches + integer, parameter :: MAX_PCP_15HR = 7620 ! 30 inches + integer, parameter :: MAX_PCP_18HR = 9144 ! 36 inches + integer, parameter :: MAX_PCP_21HR = 10668 ! 42 inches + integer, parameter :: MAX_PCP_24HR = 12192 ! 48 inches + + ! Durre et al (2010) use 1828.8 mm (72 in) for 24-hr limit, based on WMO + ! world record for 24-hr accumulation. + ! integer, parameter :: MAX_PCP_3HR = 2286 ! 9 inches + ! integer, parameter :: MAX_PCP_6HR = 4572 ! 18 inches + ! integer, parameter :: MAX_PCP_9HR = 6858 ! 27 inches + ! integer, parameter :: MAX_PCP_12HR = 9144 ! 36 inches + ! integer, parameter :: MAX_PCP_15HR = 11430 ! 45 inches + ! integer, parameter :: MAX_PCP_18HR = 13716 ! 54 inches + ! integer, parameter :: MAX_PCP_21HR = 16002 ! 63 inches + ! integer, parameter :: MAX_PCP_24HR = 18288 ! 72 inches + + ! Note: WMO world records. See https://wmo.asu.edu + ! 1-hr rainfall: 305 mm (12 inch), Holt MO, USA, June 1947 + ! 12-hr rainfall: 1144 mm (45 inch), Foc-Foc, La Reunion, Jan 1966 + ! 24-hr rainfall: 1825 mm (71.8 inch), Foc-Foc, La Reunion, Jan 1966 + +contains + + ! Constructor for USAF_gages_t type + subroutine USAF_gages_new(this, date10, & + nobs, YYYYMMDDhhmmss, networks, platforms, & + wmocode_id, fipscode_id, & + wmonumbers, bsn, lats, lons, & + amts24, amts21, amts18, amts15, amts12, amts09, amts06, amts03, & + amts02, amts01, & + amts00, durations, preswx, pastwx) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(out) :: this + character(10), intent(in) :: date10 + integer, intent(in) :: nobs + character(14), intent(in) :: YYYYMMDDhhmmss(nobs) + character(32), intent(in) :: networks(nobs) + character(32), intent(in) :: platforms(nobs) + character(2), intent(in) :: wmocode_id(nobs) + character(2), intent(in) :: fipscode_id(nobs) + integer, intent(in) :: wmonumbers(nobs) + integer, intent(in) :: bsn(nobs) + integer, intent(in) :: lats(nobs) + integer, intent(in) :: lons(nobs) + integer, intent(in) :: amts24(nobs) + integer, intent(in) :: amts21(nobs) + integer, intent(in) :: amts18(nobs) + integer, intent(in) :: amts15(nobs) + integer, intent(in) :: amts12(nobs) + integer, intent(in) :: amts09(nobs) + integer, intent(in) :: amts06(nobs) + integer, intent(in) :: amts03(nobs) + integer, intent(in) :: amts02(nobs) + integer, intent(in) :: amts01(nobs) + integer, intent(in) :: amts00(nobs) + integer, intent(in) :: durations(nobs) + integer, intent(in) :: preswx(nobs) + integer, intent(in) :: pastwx(nobs) + + this%date10 = date10 + + this%nobs = nobs + + allocate(this%YYYYMMDDhhmmss(nobs)) + this%YYYYMMDDhhmmss = YYYYMMDDhhmmss + + allocate(this%networks(nobs)) + this%networks = networks + + allocate(this%platforms(nobs)) + this%platforms = platforms + + allocate(this%wmocode_id(nobs)) + this%wmocode_id = wmocode_id + + allocate(this%fipscode_id(nobs)) + this%fipscode_id = fipscode_id + + allocate(this%wmonumbers(nobs)) + this%wmonumbers = wmonumbers + + allocate(this%bsn(nobs)) + this%bsn = bsn + + allocate(this%lats(nobs)) + this%lats = lats + + allocate(this%lons(nobs)) + this%lons = lons + + allocate(this%amts24(nobs)) + this%amts24 = amts24 + + allocate(this%amts21(nobs)) + this%amts21 = amts21 + + allocate(this%amts18(nobs)) + this%amts18 = amts18 + + allocate(this%amts15(nobs)) + this%amts15 = amts15 + + allocate(this%amts12(nobs)) + this%amts12 = amts12 + + allocate(this%amts12_orig(nobs)) + this%amts12_orig = amts12 + + allocate(this%amts09(nobs)) + this%amts09 = amts09 + + allocate(this%amts06(nobs)) + this%amts06 = amts06 + + allocate(this%amts06_orig(nobs)) + this%amts06_orig = amts06 + + allocate(this%amts03(nobs)) + this%amts03 = amts03 + + allocate(this%amts02(nobs)) + this%amts02 = amts02 + + allocate(this%amts01(nobs)) + this%amts01 = amts01 + + allocate(this%amts00(nobs)) + this%amts00 = amts00 + + allocate(this%durations(nobs)) + this%durations = durations + + allocate(this%preswx(nobs)) + this%preswx = preswx + + allocate(this%pastwx(nobs)) + this%pastwx = pastwx + + !this%pastwx_duration = set_pastwx_duration(date10) + allocate(this%pastwx_durations(nobs)) + call this%set_pastwx_durations() + + ! Make sure lats/lons are valid. + call check_latlons(this) + + ! Handle arrays for unique networks in private method. Also set & + ! cdms_flag. + call set_unique_networks(this) + + end subroutine USAF_gages_new + + ! Destructor + subroutine USAF_gages_delete(this) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + this%date10 = "NULL" + + this%nobs = 0 + if (allocated(this%YYYYMMDDhhmmss)) deallocate(this%YYYYMMDDhhmmss) + if (allocated(this%networks)) deallocate(this%networks) + if (allocated(this%platforms)) deallocate(this%platforms) + if (allocated(this%wmocode_id)) deallocate(this%wmocode_id) + if (allocated(this%fipscode_id)) deallocate(this%fipscode_id) + if (allocated(this%wmonumbers)) deallocate(this%wmonumbers) + if (allocated(this%bsn)) deallocate(this%bsn) + if (allocated(this%lats)) deallocate(this%lats) + if (allocated(this%lons)) deallocate(this%lons) + if (allocated(this%amts24)) deallocate(this%amts24) + if (allocated(this%amts21)) deallocate(this%amts21) + if (allocated(this%amts18)) deallocate(this%amts18) + if (allocated(this%amts15)) deallocate(this%amts15) + if (allocated(this%amts12)) deallocate(this%amts12) + if (allocated(this%amts12_orig)) deallocate(this%amts12_orig) + if (allocated(this%amts09)) deallocate(this%amts09) + if (allocated(this%amts06)) deallocate(this%amts06) + if (allocated(this%amts06_orig)) deallocate(this%amts06_orig) + if (allocated(this%amts03)) deallocate(this%amts03) + if (allocated(this%amts02)) deallocate(this%amts02) + if (allocated(this%amts01)) deallocate(this%amts01) + if (allocated(this%amts00)) deallocate(this%amts00) + if (allocated(this%durations)) deallocate(this%durations) + if (allocated(this%preswx)) deallocate(this%preswx) + if (allocated(this%pastwx)) deallocate(this%pastwx) + if (allocated(this%unique_networks)) deallocate(this%unique_networks) + if (allocated(this%firsts)) deallocate(this%firsts) + if (allocated(this%lasts)) deallocate(this%lasts) + if (allocated(this%pastwx_durations)) & + deallocate(this%pastwx_durations) + this%num_unique_networks = 0 + end subroutine USAF_gages_delete + + ! Use the miscellaneous precip amount. + subroutine USAF_gages_use_misc_precip(this) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + integer :: russia_lowlimit + integer :: russia_highlimit + integer :: india_lowlimit + integer :: india_highlimit + integer :: srilanka_lowlimit + integer :: srilanka_highlimit + integer :: wmo_block + character(10) :: date10 + integer :: nobs + integer :: i + + if (this%cdms_flag) then + india_lowlimit = CDMS_INDIA_LOWLIMIT + india_highlimit = CDMS_INDIA_HIGHLIMIT + russia_lowlimit = CDMS_RUSSIA_LOWLIMIT + russia_highlimit = CDMS_RUSSIA_HIGHLIMIT + srilanka_lowlimit = CDMS_SRILANKA_LOWLIMIT + srilanka_highlimit = CDMS_SRILANKA_HIGHLIMIT + else + india_lowlimit = JMOBS_INDIA_LOWLIMIT + india_highlimit = JMOBS_INDIA_HIGHLIMIT + russia_lowlimit = JMOBS_RUSSIA_LOWLIMIT + russia_highlimit = JMOBS_RUSSIA_HIGHLIMIT + srilanka_lowlimit = JMOBS_SRILANKA_LOWLIMIT + srilanka_highlimit = JMOBS_SRILANKA_HIGHLIMIT + end if + + date10 = this%date10 + + nobs = this%nobs + do i = 1, nobs + + ! WMO block number in preobs files usually range from 1 to 100. + ! Multiply so they are directly comparable to the + ! appropriate WMO block ID ranges. + wmo_block = this%wmonumbers(i) + if (wmo_block < 1000) then + if (this%cdms_flag) then + wmo_block = wmo_block*10000 + else + wmo_block = wmo_block*1000 + end if + end if + + ! India rule...Obs are usually relative to 03Z with an invalid + ! duration flag. So we may need to copy values from miscellaneous + ! (amts00) to a more standard duration. + if (this%durations(i) .eq. 0 .or. & + this%durations(i) .eq. MISSING) then + + if ((wmo_block .ge. india_lowlimit .and. & + wmo_block .le. india_highlimit) .or. & + (this%bsn(i) .ge. india_lowlimit .and. & + this%bsn(i) .le. india_highlimit)) then + + if (this%amts00(i) .ne. MISSING) then + + if (date10(9:10) == "03" .and. & + this%amts24(i) .eq. MISSING) then + this%amts24(i) = this%amts00(i) + else if (date10(9:10) == "06" .and. & + this%amts03(i) .eq. MISSING) then + this%amts03(i) = this%amts00(i) + else if (date10(9:10) == "09" .and. & + this%amts06(i) .eq. MISSING) then + this%amts06(i) = this%amts00(i) + else if (date10(9:10) == "12" .and. & + this%amts09(i) .eq. MISSING) then + this%amts09(i) = this%amts00(i) + else if (date10(9:10) == "15" .and. & + this%amts12(i) .eq. MISSING) then + this%amts12(i) = this%amts00(i) + else if (date10(9:10) == "18" .and. & + this%amts15(i) .eq. MISSING) then + this%amts15(i) = this%amts00(i) + else if (date10(9:10) == "21" .and. & + this%amts18(i) .eq. MISSING) then + this%amts18(i) = this%amts00(i) + else if (date10(9:10) == "00" .and. & + this%amts21(i) .eq. MISSING) then + this%amts21(i) = this%amts00(i) + end if + end if + end if + end if + + ! Sri Lanka rule. + if ((wmo_block .ge. srilanka_lowlimit .and. & + wmo_block .le. srilanka_highlimit) .or. & + (this%bsn(i) .ge. srilanka_lowlimit .and. & + this%bsn(i) .le. srilanka_highlimit) ) then + + ! At present, it doesn't appear a special Sri Lanka rule is + ! required. We can process 3, 9, and 15 hour reports. + continue + end if + + ! Russia rule (actually, post-Soviet rule). In the past, USAF + ! personnel noted many "Russian" obs were only reported at 00Z and + ! 12Z, and did not use a valid duration flag. These are decoded in + ! the miscellaneous accumulation (amts00), and the duration is set + ! to 0. In AGRMET they were treated as 12-hr accumulations. More + ! recently, a spot check in July 2021 showed similar invalid + ! durations in Kazakhstan, Kyrgyzstan, Georgia, Tajiskistan, and + ! Turkenistan at 06Z and 18Z. Since WMO Region VI (Europe) expects + ! 12-hr accumulations at 06Z and 18Z, we assume these rogue + ! reports are 12-hr accumulations. + if (this%durations(i) .eq. 0 .or. & + this%durations(i) .eq. MISSING) then + + if ((wmo_block .ge. russia_lowlimit .and. & + wmo_block .le. russia_highlimit) .or. & + (this%bsn(i) .ge. russia_lowlimit .and. & + this%bsn(i) .le. russia_highlimit) ) then + + if (date10(9:10) .eq. "00" .or. date10(9:10) .eq. "06" .or. & + date10(9:10) .eq. "12" .or. date10(9:10) .eq. "18") then + + if (this%amts00(i) .ne. MISSING) then + if (this%amts12(i) .eq. MISSING) then + this%amts12(i) = this%amts00(i) + end if + end if + end if + end if + end if ! Russia rule + + ! Copy from miscellaneous (amts00) to a more standard duration, + ! according to the reported duration. + if (this%amts00(i) .ne. MISSING) then + if (this%amts24(i) .eq. MISSING .and. & + this%durations(i) .eq. 24) then + this%amts24(i) = this%amts00(i) + else if (this%amts21(i) .eq. MISSING .and. & + this%durations(i) .eq. 21) then + this%amts21(i) = this%amts00(i) + else if (this%amts18(i) .eq. MISSING .and. & + this%durations(i) .eq. 18) then + this%amts18(i) = this%amts00(i) + else if (this%amts15(i) .eq. MISSING .and. & + this%durations(i) .eq. 15) then + this%amts15(i) = this%amts00(i) + else if (this%amts12(i) .eq. MISSING .and. & + this%durations(i) .eq. 12) then + this%amts12(i) = this%amts00(i) + else if (this%amts09(i) .eq. MISSING .and. & + this%durations(i) .eq. 9) then + this%amts09(i) = this%amts00(i) + else if (this%amts06(i) .eq. MISSING .and. & + this%durations(i) .eq. 6) then + this%amts06(i) = this%amts00(i) + else if (this%amts03(i) .eq. MISSING .and. & + this%durations(i) .eq. 3) then + this%amts03(i) = this%amts00(i) + ! EMK...Disabled saving 2-hr and 1-hr amounts from misc. + ! Old preobs files have a mix of different hours, so we + ! can't expect this to work. We will revisit with the + ! new preobs files, likely requiring a new argument to + ! this subroutine to activate. + ! else if (this%amts02(i) .eq. MISSING .and. & + ! this%durations(i) .eq. 2) then + ! this%amts02(i) = this%amts00(i) + ! else if (this%amts01(i) .eq. MISSING .and. & + ! this%durations(i) .eq. 1) then + ! this%amts01(i) = this%amts00(i) + end if + end if + + end do ! i + + ! Make sure shorter durations don't exceed larger durations. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_use_misc_precip + + ! Reconciles different accumulations from the same report. + subroutine USAF_gages_reconcile_self(this) + + use LIS_logMod, only: LIS_logunit + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + integer :: nobs + integer :: i + integer :: amts_tmp(10) ! 1, 2, 3, 6, 9, 12, 15, 18, 21, 24 + integer :: j, jj + + nobs = this%nobs + do i = 1, nobs + + ! Zero out smaller-duration accumulation if larger one is zero. + if (this%amts24(i) .eq. 0) then + this%amts21(i) = 0 + end if + if (this%amts21(i) .eq. 0) then + this%amts18(i) = 0 + end if + if (this%amts18(i) .eq. 0) then + this%amts15(i) = 0 + end if + if (this%amts15(i) .eq. 0) then + this%amts12(i) = 0 + end if + if (this%amts12(i) .eq. 0) then + this%amts09(i) = 0 + end if + if (this%amts09(i) .eq. 0) then + this%amts06(i) = 0 + end if + if (this%amts06(i) .eq. 0) then + this%amts03(i) = 0 + end if + if (this%amts03(i) .eq. 0) then + this%amts02(i) = 0 + end if + if (this%amts02(i) .eq. 0) then + this%amts01(i) = 0 + end if + + ! Sometimes smaller-duration accumulations exceed larger-durations, + ! due to precision differences in WMO reporting standard. We + ! correct that here. + amts_tmp(1) = this%amts01(i) + amts_tmp(2) = this%amts02(i) + amts_tmp(3) = this%amts03(i) + amts_tmp(4) = this%amts06(i) + amts_tmp(5) = this%amts09(i) + amts_tmp(6) = this%amts12(i) + amts_tmp(7) = this%amts15(i) + amts_tmp(8) = this%amts18(i) + amts_tmp(9) = this%amts21(i) + amts_tmp(10) = this%amts24(i) + + do j = 10, 2, -1 ! Longer duration limit, 24-hr to 02-hr + if (amts_tmp(j) .ne. MISSING) then + do jj = j-1, 1, -1 ! Shorter duration limit, 21-hr to 01-hr + if (amts_tmp(jj) .ne. MISSING) then + amts_tmp(jj) = min(amts_tmp(jj), amts_tmp(j)) + end if + end do ! jj + end if + end do ! j + + this%amts01(i) = amts_tmp(1) + this%amts02(i) = amts_tmp(2) + this%amts03(i) = amts_tmp(3) + this%amts06(i) = amts_tmp(4) + this%amts09(i) = amts_tmp(5) + this%amts12(i) = amts_tmp(6) + this%amts15(i) = amts_tmp(7) + this%amts18(i) = amts_tmp(8) + this%amts21(i) = amts_tmp(9) + this%amts24(i) = amts_tmp(10) + + end do ! i + + end subroutine USAF_gages_reconcile_self + + ! Fill gaps in precip record if bookend accumulations are identical. + subroutine USAF_gages_fill_gaps(this) + + use LIS_logmod, only: LIS_logunit + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + integer :: nobs + integer :: i + integer :: amts_tmp(10) ! 1, 2, 3, 6, 9, 12, 15, 18, 21, 24 + integer :: j, jj, jjj + logical :: is_data_gap + + nobs = this%nobs + do i = 1, nobs + + amts_tmp(1) = this%amts01(i) + amts_tmp(2) = this%amts02(i) + amts_tmp(3) = this%amts03(i) + amts_tmp(4) = this%amts06(i) + amts_tmp(5) = this%amts09(i) + amts_tmp(6) = this%amts12(i) + amts_tmp(7) = this%amts15(i) + amts_tmp(8) = this%amts18(i) + amts_tmp(9) = this%amts21(i) + amts_tmp(10) = this%amts24(i) + + ! Fill in gaps if shorter and longer accumulations are identical. + do j = 1, 8 ! Shorter duration limit, from 01-hr to 18-hr + do jj = j+2, 10 ! Longer duration limit, from 03-hr to 24-hr + if (amts_tmp(j) .ne. MISSING .and. & + amts_tmp(jj) .ne. MISSING .and. & + amts_tmp(j) == amts_tmp(jj)) then + ! We have two bookends with identical values. See + ! if all the durations in-between are missing. + is_data_gap = .true. + do jjj = j+1, jj-1 + if (amts_tmp(jjj) .ne. MISSING) then + is_data_gap = .false. + exit + end if + end do ! jjj + ! If all the durations in between are missing, just + ! copy the value (we can assume no additional precip + ! occurred during these periods, at least from this + ! particular report). + if (is_data_gap) then + do jjj = j+1, jj-1 + amts_tmp(jjj) = amts_tmp(j) + end do ! jjj + end if + end if + end do ! jj + end do ! j + + this%amts01(i) = amts_tmp(1) + this%amts02(i) = amts_tmp(2) + this%amts03(i) = amts_tmp(3) + this%amts06(i) = amts_tmp(4) + this%amts09(i) = amts_tmp(5) + this%amts12(i) = amts_tmp(6) + this%amts15(i) = amts_tmp(7) + this%amts18(i) = amts_tmp(8) + this%amts21(i) = amts_tmp(9) + this%amts24(i) = amts_tmp(10) + + end do ! i + + end subroutine USAF_gages_fill_gaps + + ! Reconcile current obs with obs from 1 hour ago + subroutine USAF_gages_reconcile_gages01(this, gages01) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages01 + + ! Locals + integer :: nobs + integer tmp + integer :: i, i1 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i1 = search_pcpobs(gages01, this%networks(i), this%platforms(i), & + this%wmocode_id(i), this%fipscode_id(i)) + + if (i1 .eq. MISSING) cycle + + ! Use current 1-hr accumulation, if available. + if (this%amts01(i) .ne. MISSING) then + + ! Update 3-hr accumulation if missing + if (this%amts03(i) .eq. MISSING .and. & + gages01%amts02(i1) .ne. MISSING ) then + tmp = this%amts01(i) + gages01%amts02(i1) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) this%amts03(i) = tmp + end if + + ! Update 2-hr accumulation if missing + if (this%amts02(i) .eq. MISSING .and. & + gages01%amts01(i1) .ne. MISSING ) then + tmp = this%amts01(i) + gages01%amts01(i1) + if (tmp >= 0 .and. tmp <= MAX_PCP_2HR) this%amts02(i) = tmp + end if + + end if + + ! Update the 01-hr accumulation + if (this%amts01(i) .eq. MISSING) then + if (this%amts02(i) .ne. MISSING .and. & + gages01%amts01(i1) .ne. MISSING) then + tmp = this%amts02(i) - gages01%amts01(i1) + if (tmp >= 0 .and. tmp <= MAX_PCP_1HR) then + this%amts01(i) = tmp + cycle + end if + end if + if (this%amts03(i) .ne. MISSING .and. & + gages01%amts02(i1) .ne. MISSING) then + tmp = this%amts03(i) - gages01%amts02(i1) + if (tmp >= 0 .and. tmp <= MAX_PCP_1HR) then + this%amts01(i) = tmp + cycle + end if + end if + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration accumulations + ! don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages01 + + ! Reconcile current obs with obs from 2 hours ago + subroutine USAF_gages_reconcile_gages02(this, gages02) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages02 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i2 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i2 = search_pcpobs(gages02, this%networks(i), this%platforms(i), & + this%wmocode_id(i), this%fipscode_id(i)) + + if (i2 .eq. MISSING) cycle + + ! Use current 2-hr accumulation, if available. + if (this%amts02(i) .ne. MISSING) then + + ! Update 3-hr accumulation if missing + if (this%amts03(i) .eq. MISSING .and. & + gages02%amts01(i2) .ne. MISSING ) then + tmp = this%amts02(i) + gages02%amts01(i2) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) this%amts03(i) = tmp + end if + end if + + ! Update the 2-hr accumulation + if (this%amts02(i) .eq. MISSING) then + if (this%amts03(i) .ne. MISSING .and. & + gages02%amts01(i2) .ne. MISSING) then + tmp = this%amts03(i) - gages02%amts01(i2) + if (tmp >= 0 .and. tmp <= MAX_PCP_2HR) this%amts02(i) = tmp + end if + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration accumulations + ! don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages02 + + ! Reconcile current obs with obs from 3 hours ago + subroutine USAF_gages_reconcile_gages03(this, gages03) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages03 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i3 + + nobs = this%nobs + do i = 1, nobs + + ! Find current station in prior report list + i3 = search_pcpobs(gages03, this%networks(i), this%platforms(i), & + this%wmocode_id(i), this%fipscode_id(i)) + + if (i3 .eq. MISSING) cycle + + ! Use current 3-hr accumulation, if available. + if (this%amts03(i) .ne. MISSING) then + + ! Update 24-hr accumulation if missing + if (this%amts24(i) .eq. MISSING .and. & + gages03%amts21(i3) .ne. MISSING ) then + tmp = this%amts03(i) + gages03%amts21(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_24HR) this%amts24(i) = tmp + end if + + ! Update 21-hr accumulation if missing + if (this%amts21(i) .eq. MISSING .and. & + gages03%amts18(i3) .ne. MISSING) then + tmp = this%amts03(i) + gages03%amts18(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_21HR) this%amts21(i) = tmp + end if + + ! Update 18-hr accumulation if missing + if (this%amts18(i) .eq. MISSING .and. & + gages03%amts15(i3) .ne. MISSING) then + tmp = this%amts03(i) + gages03%amts15(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_18HR) this%amts18(i) = tmp + end if + + ! Update 15-hr accumulation if missing + if (this%amts15(i) .eq. MISSING .and. & + gages03%amts12(i3) .ne. MISSING) then + tmp = this%amts03(i) + gages03%amts12(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_15HR) this%amts15(i) = tmp + end if + + ! Update 12-hr accumulation if missing + if (this%amts12(i) .eq. MISSING .and. & + gages03%amts09(i3) .ne. MISSING) then + tmp = this%amts03(i) + gages03%amts09(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_12HR) this%amts12(i) = tmp + end if + + ! Update 09-hr accumulation if missing + if (this%amts09(i) .eq. MISSING .and. & + gages03%amts06(i3) .ne. MISSING) then + tmp = this%amts03(i) + gages03%amts06(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_9HR) this%amts09(i) = tmp + end if + + ! Update 06-hr accumulation if missing + if (this%amts06(i) .eq. MISSING .and. & + gages03%amts03(i3) .ne. MISSING) then + tmp = this%amts03(i) + gages03%amts03(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_6HR) this%amts06(i) = tmp + end if + end if + + ! Multiple attempts to update 03-hr accumulation. + if (this%amts03(i) .eq. MISSING) then + if (this%amts06(i) .ne. MISSING .and. & + gages03%amts03(i3) .ne. MISSING) then + tmp = this%amts06(i) - gages03%amts03(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) then + this%amts03(i) = tmp + cycle + end if + end if + + if (this%amts09(i) .ne. MISSING .and. & + gages03%amts06(i3) .ne. MISSING) then + tmp = this%amts09(i) - gages03%amts06(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) then + this%amts03(i) = tmp + cycle + end if + end if + + if (this%amts12(i) .ne. MISSING .and. & + gages03%amts09(i3) .ne. MISSING) then + tmp = this%amts12(i) - gages03%amts09(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) then + this%amts03(i) = tmp + cycle + end if + end if + + if (this%amts15(i) .ne. MISSING .and. & + gages03%amts12(i3) .ne. MISSING) then + tmp = this%amts15(i) - gages03%amts12(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) then + this%amts03(i) = tmp + cycle + end if + end if + + if (this%amts18(i) .ne. MISSING .and. & + gages03%amts15(i3) .ne. MISSING) then + tmp = this%amts18(i) - gages03%amts15(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) then + this%amts03(i) = tmp + cycle + end if + end if + + if (this%amts21(i) .ne. MISSING .and. & + gages03%amts18(i3) .ne. MISSING) then + tmp = this%amts21(i) - gages03%amts18(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) then + this%amts03(i) = tmp + cycle + end if + end if + + if (this%amts24(i) .ne. MISSING .and. & + gages03%amts21(i3) .ne. MISSING ) then + tmp = this%amts24(i) - gages03%amts21(i3) + if (tmp >= 0 .and. tmp <= MAX_PCP_3HR) then + this%amts03(i) = tmp + cycle + end if + end if + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration accumulations + ! don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages03 + + ! Reconcile current obs with obs from 6 hours ago. + subroutine USAF_gages_reconcile_gages06(this, gages06) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages06 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i6 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i6 = search_pcpobs(gages06, this%networks(i), this%platforms(i), & + this%wmocode_id(i), this%fipscode_id(i)) + + if (i6 .eq. MISSING) cycle + + ! Use current 6-hr accumulation, if available. + if (this%amts06(i) .ne. MISSING) then + + ! Update 24-hr accumulation if missing + if (this%amts24(i) .eq. MISSING .and. & + gages06%amts18(i6) .ne. MISSING) then + tmp = this%amts06(i) + gages06%amts18(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_24HR) this%amts24(i) = tmp + end if + + ! Update 21-hr accumulation if missing + if (this%amts21(i) .eq. MISSING .and. & + gages06%amts15(i6) .ne. MISSING) then + tmp = this%amts06(i) + gages06%amts15(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_21HR) this%amts21(i) = tmp + end if + + ! Update 18-hr accumulation if missing + if (this%amts18(i) .eq. MISSING .and. & + gages06%amts12(i6) .ne. MISSING) then + tmp = this%amts06(i) + gages06%amts12(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_18HR) this%amts18(i) = tmp + end if + + ! Update 15-hr accumulation if missing + if (this%amts15(i) .eq. MISSING .and. & + gages06%amts09(i6) .ne. MISSING) then + tmp = this%amts06(i) + gages06%amts09(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_15HR) this%amts15(i) = tmp + end if + + ! Update 12-hr accumulation if missing + if (this%amts12(i) .eq. MISSING .and. & + gages06%amts06(i6) .ne. MISSING) then + tmp = this%amts06(i) + gages06%amts06(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_12HR) this%amts12(i) = tmp + end if + + ! Update 09-hr accumulation if missing + if (this%amts09(i) .eq. MISSING .and. & + gages06%amts03(i6) .ne. MISSING) then + tmp = this%amts06(i) + gages06%amts03(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_9HR) this%amts09(i) = tmp + end if + + end if + + ! Multiple attempts to update 06-hr accumulation. + if (this%amts06(i) .eq. MISSING) then + if ( (this%amts12(i) .ne. MISSING) .and. & + (gages06%amts06(i6) .ne. MISSING) ) then + tmp = this%amts12(i) - gages06%amts06(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_6HR) then + this%amts06(i) = tmp + cycle + end if + end if + + if ( (this%amts15(i) .ne. MISSING) .and. & + (gages06%amts09(i6) .ne. MISSING) ) then + tmp = this%amts15(i) - gages06%amts09(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_6HR) then + this%amts06(i) = tmp + cycle + end if + end if + + if ( (this%amts18(i) .ne. MISSING) .and. & + (gages06%amts12(i6) .ne. MISSING) ) then + tmp = this%amts18(i) - gages06%amts12(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_6HR) then + this%amts06(i) = tmp + cycle + end if + end if + + if ( (this%amts21(i) .ne. MISSING) .and. & + (gages06%amts15(i6) .ne. MISSING) ) then + tmp = this%amts21(i) - gages06%amts15(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_6HR) then + this%amts06(i) = tmp + cycle + end if + end if + + if ( (this%amts24(i) .ne. MISSING) .and. & + (gages06%amts18(i6) .ne. MISSING) ) then + tmp = this%amts24(i) - gages06%amts18(i6) + if (tmp >= 0 .and. tmp <= MAX_PCP_6HR) then + this%amts06(i) = tmp + cycle + end if + end if + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration + ! accumulations don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages06 + + ! Reconcile current obs with obs from 9 hours ago. + subroutine USAF_gages_reconcile_gages09(this, gages09) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages09 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i9 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i9 = search_pcpobs(gages09, this%networks(i), this%platforms(i), & + this%wmocode_id(i), this%fipscode_id(i)) + + if (i9 .eq. MISSING) cycle + + ! Use current 9-hr accumulation, if available. + if (this%amts09(i) .ne. MISSING) then + + ! Update 24-hr accumulation if missing + if (this%amts24(i) .eq. MISSING .and. & + gages09%amts15(i9) .ne. MISSING ) then + tmp = this%amts09(i) + gages09%amts15(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_24HR) this%amts24(i) = tmp + end if + + ! Update 21-hr accumulation if missing + if (this%amts21(i) .eq. MISSING .and. & + gages09%amts12(i9) .ne. MISSING ) then + tmp = this%amts09(i) + gages09%amts12(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_21HR) this%amts21(i) = tmp + end if + + ! Update 18-hr accumulation if missing + if (this%amts18(i) .eq. MISSING .and. & + gages09%amts09(i9) .ne. MISSING) then + tmp = this%amts09(i) + gages09%amts09(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_18HR) this%amts18(i) = tmp + end if + + ! Update 15-hr accumulation if missing + if (this%amts15(i) .eq. MISSING .and. & + gages09%amts06(i9) .ne. MISSING) then + tmp = this%amts09(i) + gages09%amts06(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_15HR) this%amts15(i) = tmp + end if + + ! Update 12-hr accumulation if missing + if (this%amts12(i) .eq. MISSING .and. & + gages09%amts03(i9) .ne. MISSING) then + tmp = this%amts09(i) + gages09%amts03(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_12HR) this%amts12(i) = tmp + end if + + end if + + ! Multiple attempts to update 09-hr accumulation. + if (this%amts09(i) .eq. MISSING) then + if ( (this%amts12(i) .ne. MISSING) .and. & + (gages09%amts03(i9) .ne. MISSING) ) then + tmp = this%amts12(i) - gages09%amts03(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_9HR) then + this%amts09(i) = tmp + cycle + end if + end if + + if ( (this%amts15(i) .ne. MISSING) .and. & + (gages09%amts06(i9) .ne. MISSING) ) then + this%amts09(i) = this%amts15(i) - gages09%amts06(i9) + this%amts09(i) = max(this%amts09(i), 0) + tmp = this%amts15(i) - gages09%amts06(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_9HR) then + this%amts09(i) = tmp + cycle + end if + end if + + if ( (this%amts18(i) .ne. MISSING) .and. & + (gages09%amts09(i9) .ne. MISSING) ) then + tmp = this%amts18(i) - gages09%amts09(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_9HR) then + this%amts09(i) = tmp + cycle + end if + end if + + if ( (this%amts21(i) .ne. MISSING) .and. & + (gages09%amts12(i9) .ne. MISSING) ) then + tmp = this%amts21(i) - gages09%amts12(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_9HR) then + this%amts09(i) = tmp + cycle + end if + end if + + if ( (this%amts24(i) .ne. MISSING) .and. & + (gages09%amts15(i9) .ne. MISSING) ) then + tmp = this%amts24(i) - gages09%amts15(i9) + if (tmp >= 0 .and. tmp <= MAX_PCP_9HR) then + this%amts09(i) = tmp + cycle + end if + end if + end if + end do + + ! Reconcile current obs again to ensure smaller-duration + ! accumulations don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages09 + + ! Reconcile current obs with obs from 12 hours ago. + subroutine USAF_gages_reconcile_gages12(this, gages12) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages12 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i12 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i12 = search_pcpobs(gages12, this%networks(i), & + this%platforms(i), this%wmocode_id(i), this%fipscode_id(i)) + + if (i12 .eq. MISSING) cycle + + ! Use current 12-hr accumulation, if available + if (this%amts12(i) .ne. MISSING) then + + ! Update 24-hr accumulation if missing + if (this%amts24(i) .eq. MISSING .and. & + gages12%amts12(i12) .ne. MISSING) then + tmp = this%amts12(i) + gages12%amts12(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_24HR) this%amts24(i) = tmp + end if + + ! Update 21-hr accumulation if missing + if (this%amts21(i) .eq. MISSING .and. & + gages12%amts09(i12) .ne. MISSING) then + tmp = this%amts12(i) + gages12%amts09(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_21HR) this%amts21(i) = tmp + end if + + ! Update 18-hr accumulation if missing + if (this%amts18(i) .eq. MISSING .and. & + gages12%amts06(i12) .ne. MISSING) then + tmp = this%amts12(i) + gages12%amts06(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_18HR) this%amts18(i) = tmp + end if + + ! Update 15-hr accumulation if missing + if (this%amts15(i) .eq. MISSING .and. & + gages12%amts03(i12) .ne. MISSING) then + tmp = this%amts12(i) + gages12%amts03(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_15HR) this%amts15(i) = tmp + end if + end if + + ! Multiple attempts to update 12-hr accumulation. + if (this%amts12(i) .eq. MISSING) then + if ( (this%amts15(i) .ne. MISSING) .and. & + (gages12%amts03(i12) .ne. MISSING) ) then + tmp = this%amts15(i) - gages12%amts03(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_12HR) then + this%amts12(i) = tmp + cycle + end if + end if + + if( (this%amts18(i) .ne. MISSING) .and. & + (gages12%amts06(i12) .ne. MISSING) ) then + tmp = this%amts18(i) - gages12%amts06(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_12HR) then + this%amts12(i) = tmp + cycle + end if + end if + + if( (this%amts21(i) .ne. MISSING) .and. & + (gages12%amts09(i12) .ne. MISSING) ) then + tmp = this%amts21(i) - gages12%amts09(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_12HR) then + this%amts12(i) = tmp + cycle + end if + end if + + if ( (this%amts24(i) .ne. MISSING) .and. & + (gages12%amts12(i12) .ne. MISSING) ) then + tmp = this%amts24(i) - gages12%amts12(i12) + if (tmp >= 0 .and. tmp <= MAX_PCP_12HR) then + this%amts12(i) = tmp + cycle + end if + end if + + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration + ! accumulations don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages12 + + ! Reconcile current obs with obs from 15 hours ago. + subroutine USAF_gages_reconcile_gages15(this, gages15) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages15 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i15 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i15 = search_pcpobs(gages15, this%networks(i), & + this%platforms(i), this%wmocode_id(i), this%fipscode_id(i)) + + if (i15 .eq. MISSING) cycle + + ! Use current 15-hr accumulation, if available + if (this%amts15(i) .ne. MISSING) then + + ! Update 24-hr accumulation if missing + if (this%amts24(i) .eq. MISSING .and. & + gages15%amts09(i15) .ne. MISSING) then + tmp = this%amts15(i) + gages15%amts09(i15) + if (tmp >= 0 .and. tmp <= MAX_PCP_24HR) this%amts24(i) = tmp + end if + + ! Update 21-hr accumulation if missing + if (this%amts21(i) .eq. MISSING .and. & + gages15%amts06(i15) .ne. MISSING) then + tmp = this%amts15(i) + gages15%amts06(i15) + if (tmp >= 0 .and. tmp <= MAX_PCP_21HR) this%amts21(i) = tmp + end if + + ! Update 18-hr accumulation if missing + if (this%amts18(i) .eq. MISSING .and. & + gages15%amts03(i15) .ne. MISSING) then + tmp = this%amts15(i) + gages15%amts03(i15) + if (tmp >= 0 .and. tmp <= MAX_PCP_18HR) this%amts18(i) = tmp + end if + + end if + + ! Multiple attempts to update 15-hr accumulation. + if (this%amts15(i) .eq. MISSING) then + if ( (this%amts18(i) .ne. MISSING) .and. & + (gages15%amts03(i15) .ne. MISSING) ) then + tmp = this%amts18(i) - gages15%amts03(i15) + if (tmp >= 0 .and. tmp <= MAX_PCP_15HR) then + this%amts15(i) = tmp + cycle + end if + end if + + if ( (this%amts21(i) .ne. MISSING) .and. & + (gages15%amts06(i15) .ne. MISSING) ) then + tmp = this%amts21(i) - gages15%amts06(i15) + if (tmp >= 0 .and. tmp <= MAX_PCP_15HR) then + this%amts15(i) = tmp + cycle + end if + end if + + if ( (this%amts24(i) .ne. MISSING) .and. & + (gages15%amts09(i15) .ne. MISSING) ) then + tmp = this%amts24(i) - gages15%amts09(i15) + if (tmp >= 0 .and. tmp <= MAX_PCP_15HR) then + this%amts15(i) = tmp + cycle + end if + end if + + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration & + ! accumulations don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages15 + + ! Reconcile current obs with obs from 18 hours ago. + subroutine USAF_gages_reconcile_gages18(this, gages18) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages18 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i18 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i18 = search_pcpobs(gages18, this%networks(i), & + this%platforms(i), this%wmocode_id(i), this%fipscode_id(i)) + + if (i18 .eq. MISSING) cycle + + ! Use current 18-hr accumulation, if available + if (this%amts18(i) .ne. MISSING) then + + ! Update 24-hr accumulation if missing + if (this%amts24(i) .eq. MISSING .and. & + gages18%amts06(i18) .ne. MISSING) then + tmp = this%amts18(i) + gages18%amts06(i18) + if (tmp >= 0 .and. tmp <= MAX_PCP_24HR) this%amts24(i) = tmp + end if + + ! Update 21-hr accumulation if missing + if (this%amts21(i) .eq. MISSING .and. & + gages18%amts03(i18) .ne. MISSING) then + tmp = this%amts18(i) + gages18%amts03(i18) + if (tmp >= 0 .and. tmp <= MAX_PCP_21HR) this%amts21(i) = tmp + end if + + end if + + ! Multiple attempts to update 18-hr accumulation. + if (this%amts18(i) .eq. MISSING) then + if (this%amts21(i) .ne. MISSING .and. & + gages18%amts03(i18) .ne. MISSING) then + tmp = this%amts21(i) - gages18%amts03(i18) + if (tmp >= 0 .and. tmp <= MAX_PCP_18HR) then + this%amts18(i) = tmp + cycle + end if + end if + + if (this%amts24(i) .ne. MISSING .and. & + gages18%amts06(i18) .ne. MISSING) then + tmp = this%amts24(i) - gages18%amts06(i18) + if (tmp >= 0 .and. tmp <= MAX_PCP_18HR) then + this%amts18(i) = tmp + cycle + end if + end if + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration + ! accumulations don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages18 + + ! Reconcile current obs with obs from 21 hours ago. + subroutine USAF_gages_reconcile_gages21(this, gages21) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + class(USAF_gages_t), intent(in) :: gages21 + + ! Locals + integer :: nobs + integer :: tmp + integer :: i, i21 + + nobs = this%nobs + + do i = 1, nobs + + ! Find current station in prior report list + i21 = search_pcpobs(gages21, this%networks(i), & + this%platforms(i), this%wmocode_id(i), this%fipscode_id(i)) + + if (i21 .eq. MISSING) cycle + + ! Use current 21-hr accumulation, if available + if (this%amts21(i) .ne. MISSING) then + + ! Update 24-hr accumulation if missing + if (this%amts24(i) .eq. MISSING .and. & + gages21%amts03(i21) .ne. MISSING) then + tmp = this%amts21(i) + gages21%amts03(i21) + if (tmp >= 0 .and. tmp <= MAX_PCP_24HR) this%amts24(i) = tmp + end if + + end if + + ! Update 21-hr accumulation + if (this%amts21(i) .eq. MISSING) then + if (this%amts24(i) .ne. MISSING .and. & + gages21%amts03(i21) .ne. MISSING) then + tmp = this%amts24(i) - gages21%amts03(i21) + if (tmp >= 0 .and. tmp <= MAX_PCP_21HR) this%amts21(i) = tmp + end if + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration + ! accumulations don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_reconcile_gages21 + + ! Correct for missing "zero-precip" 12-hr reports at 12Z for Region III + ! (South America). Based on logic in AGRMET_processobs. + subroutine USAF_gages_correct_region3_12Z(this) + + use LIS_logmod, only: LIS_logunit + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + integer :: s_amer_lowlimit + integer :: s_amer_highlimit + integer :: wmo_block + character(10) :: date10 + integer :: nobs + integer :: i + + ! Only work with 12Z observations. + date10 = this%date10 + if (date10(9:10) .ne. "12") return + + if (this%cdms_flag) then + s_amer_lowlimit = CDMS_S_AMER_LOWLIMIT + s_amer_highlimit = CDMS_S_AMER_HIGHLIMIT + else + s_amer_lowlimit = JMOBS_S_AMER_LOWLIMIT + s_amer_highlimit = JMOBS_S_AMER_HIGHLIMIT + end if + + nobs = this%nobs + + do i = 1, nobs + + ! WMO block number in preobs files usually range from 1 to 100. + ! Multiply so they are directly comparable to the + ! appropriate WMO block ID ranges. + wmo_block = this%wmonumbers(i) + if (wmo_block < 1000) then + if (this%cdms_flag) then + wmo_block = wmo_block*10000 + else + wmo_block = wmo_block*1000 + end if + end if + + ! In the past, USAF personnel noticed a lack of "zero-precip" + ! 12-hr accumulations at 12Z over South America, which they + ! attributed to two causes: + ! (1) Many stations fail to report at 06Z; and + ! (2) Many 12Z reports indicate zero precip without indicating + ! duration. (Speculation: SYNOP reports excluded a 6RRRt_r + ! precip group in the ob, and instead mark the "I_R" + ! indicator in Section 0 as "3" (meaning precipitation amount + ! is zero) without indicating the duration.) + ! Apparently this is decoded by default as "no precip over last + ! six hours". + ! To reconcile, South American obs at 12Z with missing 12-hr + ! accumulations but with "zero" 6-hr accumulations will also + ! be treated as "zero" 12-hr accumulations. + ! EXCEPTION: If a rare 9-hour non-zero report exists, don't + ! change the 12-hr. + if (this%amts12(i) .ne. MISSING) cycle + if (this%amts12_orig(i) .ne. MISSING) cycle + if (this%amts09(i) .gt. 0) cycle + + if ((wmo_block .ge. s_amer_lowlimit .and. & + wmo_block .le. s_amer_highlimit) .or. & + (this%bsn(i) .ge. s_amer_lowlimit .and. & + this%bsn(i) .le. s_amer_highlimit)) then + + ! We want to make sure the 6-hr report is "original", and + ! not just filled in or derived from comparing other gage + ! reports or present/past weather. That is, the station + ! really made a "zero precip" ob at 12Z, and nothing is + ! available from 06Z. + if (this%amts06(i) .eq. 0 .and. & + this%amts06_orig(i) .eq. 0) then + this%amts12(i) = 0 + this%amts09(i) = 0 + this%amts06(i) = 0 + this%amts03(i) = 0 + this%amts02(i) = 0 + this%amts01(i) = 0 + end if + end if + + end do + + ! Reconcile current obs again to ensure smaller-duration + ! accumulations don't exceed longer-duration. + call this%reconcile_self() + call this%fill_gaps() + + end subroutine USAF_gages_correct_region3_12Z + + ! Write gage data to output file. + subroutine USAF_gages_write_data(this, filename) + + ! Imports + use LIS_logMod, only: LIS_getNextUnitNumber, LIS_releaseUnitNumber, LIS_logunit + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(in) :: this + character(*), intent(in) :: filename + + ! Locals + integer :: istat + integer :: nobs, nobs_good + integer :: i + integer :: iunit + + nobs = this%nobs + + ! First, count number of good obs in the structure, after quality + ! control. + nobs_good = 0 + do i = 1, nobs + if (this%amts24(i) .eq. MISSING .and. & + this%amts21(i) .eq. MISSING .and. & + this%amts18(i) .eq. MISSING .and. & + this%amts15(i) .eq. MISSING .and. & + this%amts12(i) .eq. MISSING .and. & + this%amts09(i) .eq. MISSING .and. & + this%amts06(i) .eq. MISSING .and. & + this%amts03(i) .eq. MISSING .and. & + this%amts02(i) .eq. MISSING .and. & + this%amts01(i) .eq. MISSING) cycle + nobs_good = nobs_good + 1 + end do + + iunit = LIS_getNextUnitNumber() + open(iunit, file=trim(filename), iostat=istat, err=300) + write(iunit, *, iostat=istat, err=300) nobs_good + do i = 1, nobs + + if (this%amts24(i) .eq. MISSING .and. & + this%amts21(i) .eq. MISSING .and. & + this%amts18(i) .eq. MISSING .and. & + this%amts15(i) .eq. MISSING .and. & + this%amts12(i) .eq. MISSING .and. & + this%amts09(i) .eq. MISSING .and. & + this%amts06(i) .eq. MISSING .and. & + this%amts03(i) .eq. MISSING .and. & + this%amts02(i) .eq. MISSING .and. & + this%amts01(i) .eq. MISSING) cycle + write(iunit, 6000, iostat=istat, err=300) & + this%YYYYMMDDhhmmss(i), & + this%networks(i), this%platforms(i), & + this%wmocode_id(i), this%fipscode_id(i), & + this%wmonumbers(i), this%bsn(i), & + this%lats(i), this%lons(i), & + this%amts24(i), this%amts21(i), this%amts18(i), & + this%amts15(i), this%amts12(i), this%amts09(i), & + this%amts06(i), this%amts03(i), & + this%amts02(i), this%amts01(i), & + this%amts00(i), this%durations(i), & + this%preswx(i), this%pastwx(i) +6000 format (a14, 1x, & + a32, 1x, a32, 1x, & + a2, 1x, a2, 1x, & + i6, 1x, i6, 1x, & + i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, & + i9, 1x, i9, 1x & + i9, 1x, i9) + + end do + + 300 continue + close(iunit) + call LIS_releaseUnitNumber(iunit) + + end subroutine USAF_gages_write_data + + ! Read gage data from file. Acts as an alternative constructor. + subroutine USAF_gages_read_data(this, filename, date10, alert_number) + + ! Imports + use LIS_coreMod, only: LIS_masterproc + use LIS_logMod, only: LIS_getNextUnitNumber, LIS_releaseUnitNumber, & + LIS_alert, LIS_logunit + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(out) :: this + character(*), intent(in) :: filename + character(10), intent(in) :: date10 + integer, intent(inout) :: alert_number + + ! Locals + integer :: nobs + character(14), allocatable :: YYYYMMDDhhmmss(:) + character(32), allocatable :: networks(:) + character(32), allocatable :: platforms(:) + character(2), allocatable :: wmocode_id(:) + character(2), allocatable :: fipscode_id(:) + integer, allocatable :: wmonumbers(:) + integer, allocatable :: bsn(:) + integer, allocatable :: lats(:) + integer, allocatable :: lons(:) + integer, allocatable :: amts24(:) + integer, allocatable :: amts21(:) + integer, allocatable :: amts18(:) + integer, allocatable :: amts15(:) + integer, allocatable :: amts12(:) + integer, allocatable :: amts09(:) + integer, allocatable :: amts06(:) + integer, allocatable :: amts03(:) + integer, allocatable :: amts02(:) + integer, allocatable :: amts01(:) + integer, allocatable :: amts00(:) + integer, allocatable :: durations(:) + integer, allocatable :: preswx(:) + integer, allocatable :: pastwx(:) + integer :: istat + logical :: found + integer :: i + integer :: iunit + character(255) :: message(20) + + message = '' + call this%delete() ! Make sure structure is empty + + inquire(file=trim(filename), exist=found) + if (.not. found) then + write(LIS_logunit,*)'[WARN] Cannot find ', trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_read_data' + message(3) = ' Cannot find presav2 file ' // & + trim(filename) + message(4) = ' Observation count will be reduced' + if (LIS_masterproc) then + call LIS_alert('LIS.USAF_read_data', & + alert_number, message) + alert_number = alert_number + 1 + end if + return + end if + iunit = LIS_getNextUnitNumber() + + open(iunit, file=trim(filename), iostat=istat) + if (istat .ne. 0) then + write(LIS_logunit,*)'[WARN] Problem opening ', trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_gages_read_data' + message(3) = ' Cannot open file ' // trim(filename) + if (LIS_masterproc) then + call LIS_alert('LIS.USAF_gages_read_data', & + alert_number, message) + alert_number = alert_number + 1 + end if + return + end if + + write(LIS_logunit,*) '[INFO] Reading ', trim(filename) + read(iunit, *, iostat=istat) nobs + if (istat .ne. 0) then + write(LIS_logunit,*)'[WARN] Problem reading ', trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_gages_read_data' + message(3) = ' Problem reading file ' // trim(filename) + if (LIS_masterproc) then + call LIS_alert('LIS.USAF_gages_read_data', & + alert_number, message) + alert_number = alert_number + 1 + end if + close(iunit) + call LIS_releaseUnitNumber(iunit) + return + end if + + if (nobs .le. 0) then + write(LIS_logunit,*)'[WARN] No precip obs found in ', & + trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_gages_read_data' + message(3) = ' No precip obs found in file ' // trim(filename) + if (LIS_masterproc) then + call LIS_alert('LIS.USAF_gages_read_data', & + alert_number, message) + alert_number = alert_number + 1 + end if + close(iunit) + call LIS_releaseUnitNumber(iunit) + return + end if + + allocate(YYYYMMDDhhmmss(nobs)) + allocate(networks(nobs)) + allocate(platforms(nobs)) + allocate(wmocode_id(nobs)) + allocate(fipscode_id(nobs)) + allocate(wmonumbers(nobs)) + allocate(bsn(nobs)) + allocate(lats(nobs)) + allocate(lons(nobs)) + allocate(amts24(nobs)) + allocate(amts21(nobs)) + allocate(amts18(nobs)) + allocate(amts15(nobs)) + allocate(amts12(nobs)) + allocate(amts09(nobs)) + allocate(amts06(nobs)) + allocate(amts03(nobs)) + allocate(amts02(nobs)) + allocate(amts01(nobs)) + allocate(amts00(nobs)) + allocate(durations(nobs)) + allocate(preswx(nobs)) + allocate(pastwx(nobs)) + + do i = 1, nobs + read(iunit, 6000, iostat=istat, err=300, end=300) & + YYYYMMDDhhmmss(i), & + networks(i), platforms(i), & + wmocode_id(i), fipscode_id(i), & + wmonumbers(i), bsn(i), & + lats(i), lons(i), & + amts24(i), amts21(i), amts18(i), & + amts15(i), amts12(i), amts09(i), & + amts06(i), amts03(i), & + amts02(i), amts01(i), & + amts00(i), durations(i), & + preswx(i), pastwx(i) +6000 format (a14, 1x, & + a32, 1x, a32, 1x, & + a2, 1x, a2, 1x, & + i6, 1x, i6, 1x, & + i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, & + i9, 1x, i9, 1x, & + i9, 1x, i9) + end do + + close(iunit) + call LIS_releaseUnitNumber(iunit) +300 continue + + ! If read was successful, copy to USAF_gages_t structure. + if (istat .eq. 0) then + call this%new(date10, nobs, YYYYMMDDhhmmss, & + networks, platforms, & + wmocode_id, fipscode_id, & + wmonumbers, bsn, lats, lons, & + amts24, amts21, amts18, amts15, amts12, amts09, amts06, & + amts03, amts02, amts01, & + amts00, durations, preswx, pastwx) + else + write(LIS_logunit,*)'[WARN] Problem reading from ', & + trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_gages_read_data' + message(3) = ' Problem reading from file ' // trim(filename) + if (LIS_masterproc) then + call LIS_alert('LIS.USAF_gages_read_data', & + alert_number, message) + alert_number = alert_number + 1 + end if + end if + + ! Clean up + deallocate(YYYYMMDDhhmmss) + deallocate(networks) + deallocate(platforms) + deallocate(wmocode_id) + deallocate(fipscode_id) + deallocate(wmonumbers) + deallocate(bsn) + deallocate(lats) + deallocate(lons) + deallocate(amts24) + deallocate(amts21) + deallocate(amts18) + deallocate(amts15) + deallocate(amts12) + deallocate(amts09) + deallocate(amts06) + deallocate(amts03) + deallocate(amts02) + deallocate(amts01) + deallocate(amts00) + deallocate(durations) + deallocate(preswx) + deallocate(pastwx) + + ! Handle arrays for unique networks in private method. Also set + ! cdms_flag. + if (istat .eq. 0) then + call set_unique_networks(this) + end if + + !write(LIS_logunit,*)'[INFO] Read in ', nobs, ' gage reports' + end subroutine USAF_gages_read_data + + ! Search USAF_gages_t type for presumed erroneous accumulations. + subroutine USAF_gages_check_gross_errors(this) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + character(10) :: date10 + integer :: india_lowlimit + integer :: india_highlimit + integer :: russia_lowlimit + integer :: russia_highlimit + integer :: threshold + integer :: threshold_india + integer :: nobs + integer :: i + + if (this%cdms_flag) then + india_lowlimit = CDMS_INDIA_LOWLIMIT + india_highlimit = CDMS_INDIA_HIGHLIMIT + russia_lowlimit = CDMS_RUSSIA_LOWLIMIT + russia_highlimit = CDMS_RUSSIA_HIGHLIMIT + else + india_lowlimit = JMOBS_INDIA_LOWLIMIT + india_highlimit = JMOBS_INDIA_HIGHLIMIT + russia_lowlimit = JMOBS_RUSSIA_LOWLIMIT + russia_highlimit = JMOBS_RUSSIA_HIGHLIMIT + end if + date10 = this%date10 + + nobs = this%nobs + do i = 1, nobs + + ! Check for negatives + if (this%amts24(i) < 0) this%amts24(i) = MISSING + if (this%amts21(i) < 0) this%amts21(i) = MISSING + if (this%amts18(i) < 0) this%amts18(i) = MISSING + if (this%amts15(i) < 0) this%amts15(i) = MISSING + if (this%amts12(i) < 0) this%amts12(i) = MISSING + if (this%amts09(i) < 0) this%amts09(i) = MISSING + if (this%amts06(i) < 0) this%amts06(i) = MISSING + if (this%amts03(i) < 0) this%amts03(i) = MISSING + if (this%amts02(i) < 0) this%amts02(i) = MISSING + if (this%amts01(i) < 0) this%amts01(i) = MISSING + if (this%amts00(i) < 0) this%amts00(i) = MISSING + + ! Check for large errors + if (this%amts24(i) > MAX_PCP_24HR) this%amts24(i) = MISSING + if (this%amts21(i) > MAX_PCP_21HR) this%amts21(i) = MISSING + if (this%amts18(i) > MAX_PCP_18HR) this%amts18(i) = MISSING + if (this%amts15(i) > MAX_PCP_15HR) this%amts15(i) = MISSING + if (this%amts12(i) > MAX_PCP_12HR) this%amts12(i) = MISSING + if (this%amts09(i) > MAX_PCP_9HR) this%amts09(i) = MISSING + if (this%amts06(i) > MAX_PCP_6HR) this%amts06(i) = MISSING + if (this%amts03(i) > MAX_PCP_3HR) this%amts03(i) = MISSING + if (this%amts02(i) > MAX_PCP_2HR) this%amts02(i) = MISSING + if (this%amts01(i) > MAX_PCP_1HR) this%amts01(i) = MISSING + + ! Checking miscellaneous accumulations is more complicated. + ! First, see if we have a WMO established duration. + threshold = set_precip_duration_threshold(this%durations(i)) + if (threshold .ne. MISSING) then + if (this%amts00(i) > threshold) then + this%amts00(i) = MISSING + this%durations(i) = MISSING + end if + cycle + end if + + ! Try the "Russia" rule. Three sets of tests before invoking. + if (this%durations(i) .eq. 0 .or. & + this%durations(i) .eq. MISSING) then + + if ((this%wmonumbers(i) .ge. russia_lowlimit .and. & + this%wmonumbers(i) .le. russia_highlimit) .or. & + (this%bsn(i) .ge. russia_lowlimit .and. & + this%bsn(i) .le. russia_highlimit) ) then + + if (date10(9:10) .eq. "00" .or. & + date10(9:10) .eq. "06" .or. & + date10(9:10) .eq. "12" .or. & + date10(9:10) .eq. "18") then + + ! Russia rule (actually, post-Soviet rule). In the + ! past, USAF personnel noted many "Russian" obs were + ! only reported at 00Z and 12Z, and did not use a valid + ! duration flag. These are decoded in the miscellaneous + ! accumulation (amts00). In AGRMET they were treated as + ! 12-hr accumulations. More recently, a spot check in + ! July 2021 showed similar invalid durations in + ! Kazakhstan, Kyrgyzstan, Georgia, Tajiskistan, and + ! Turkenistan at 06Z and 18Z. Since WMO Region VI + ! (Europe) expects 12-hr accumulations at 06Z and 18Z, we + ! assume these rogue reports are also 12-hr + ! accumulations. + + if (this%amts00(i) > MAX_PCP_12HR) then + this%amts00(i) = MISSING + this%durations(i) = MISSING + end if + cycle + end if + end if + end if ! Russia rule + + ! Try the India rule. Two sets of tests. + if (this%durations(i) .eq. 0 .or. & + this%durations(i) .eq. MISSING) then + + if ((this%wmonumbers(i) .ge. india_lowlimit .and. & + this%wmonumbers(i) .le. india_highlimit) .or. & + (this%bsn(i) .ge. india_lowlimit .or. & + this%bsn(i) .le. india_highlimit)) then + + ! If these are from India, we assume duration is from 03Z. + ! But, we will only accept reports at three-hour intervals. + threshold_india = set_india_precip_threshold(date10(9:10)) + if (threshold_india .ne. MISSING) then + if (this%amts00(i) > threshold_india) then + this%amts00(i) = MISSING + this%durations(i) = MISSING + end if + cycle + end if + end if + end if ! India rule + + ! If all else fails, assume 24-hr accumulation + if (this%amts00(i) > MAX_PCP_24HR) then + this%amts00(i) = MISSING + this%durations(i) = MISSING + end if + + end do + end subroutine USAF_gages_check_gross_errors + + ! Check the observations for bad lat/lon values. If found, reset the + ! data to missing. + subroutine check_latlons(this) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + logical :: bad_latlon + integer :: nobs + integer :: i + + nobs = this%nobs + do i = 1, nobs + bad_latlon = .false. + + ! Sometimes a station report is listed with a lat/lon of 0, + ! putting it in the Atlantic Ocean. We check and flag. + if (this%lats(i) == 0 .and. this%lons(i) == 0) bad_latlon = .true. + + ! If we don't have a good lat/lon, the information from the + ! station is useless. So remove it. + if (bad_latlon) then + this%amts24(i) = MISSING + this%amts21(i) = MISSING + this%amts18(i) = MISSING + this%amts15(i) = MISSING + this%amts12(i) = MISSING + this%amts09(i) = MISSING + this%amts06(i) = MISSING + this%amts03(i) = MISSING + this%amts02(i) = MISSING + this%amts01(i) = MISSING + this%amts00(i) = MISSING + this%durations(i) = MISSING + this%preswx(i) = MISSING + this%pastwx(i) = MISSING + end if + end do + end subroutine check_latlons + + ! Search USAF_gages_t type for unique networks, and establish index + ! bounds for each network. Assumes data are sorted. Borrows logic + ! from AGRMET_processobs. + subroutine set_unique_networks(this) + + ! Imports + use LIS_logMod, only: LIS_logunit, LIS_endrun + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + character(32) :: prior_net + integer :: net_count + integer :: nobs + integer :: i + + if (.not. allocated(this%unique_networks)) then + allocate(this%unique_networks(MAX_UNIQUE_NETWORKS)) + end if + this%unique_networks = MISSING_NAME + if (.not. allocated(this%firsts)) then + allocate(this%firsts(MAX_UNIQUE_NETWORKS)) + end if + this%firsts = MISSING + if (.not. allocated(this%lasts)) then + allocate(this%lasts(MAX_UNIQUE_NETWORKS)) + end if + this%lasts = MISSING + + prior_net = MISSING_NAME + + nobs = this%nobs + do i = 1, nobs + if (prior_net .eq. MISSING_NAME) then + this%firsts(1) = i + this%unique_networks(i) = this%networks(i) + net_count = 1 + prior_net = this%networks(i) + else if (this%networks(i) .ne. prior_net) then + this%lasts(net_count) = i - 1 + net_count = net_count + 1 + if (net_count .gt. MAX_UNIQUE_NETWORKS) then + write(LIS_logunit,*) '[ERR] Too many unique networks!' + call LIS_endrun + end if + this%firsts(net_count) = i + this%unique_networks(net_count) = this%networks(i) + prior_net = this%networks(i) + end if + end do + this%lasts(net_count) = this%nobs + + this%num_unique_networks = net_count + + ! Mark if these obs are all legacy CDMS output + this%cdms_flag = .false. + if (net_count .eq. 1) then + if (this%unique_networks(net_count) .eq. "CDMS") then + this%cdms_flag = .true. + end if + end if + end subroutine set_unique_networks + + ! Search USAF_gages_t type for given network and platform, and + ! return index. Based on AGRMET_pcpobs_search + function search_pcpobs(this, network, plat_id, & + wmocode_id, fipscode_id) result(index) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(in) :: this + character(32), intent(in) :: network + character(32), intent(in) :: plat_id + character(2), intent(in) :: wmocode_id + character(2), intent(in) :: fipscode_id + + ! Result + integer :: index + + ! Locals + integer :: first + integer :: middle + integer :: last + integer :: net_index + logical :: found + + index = MISSING + first = 1 + last = this%nobs + found = .false. + + net_index = 1 + + ! Select the appropriate 'first' and 'last' search indices based on + ! requested network + do while ((net_index .le. this%num_unique_networks) .and. & + (.not. found)) + if (network .eq. this%unique_networks(net_index)) then + first = this%firsts(net_index) + last = this%lasts(net_index) + found = .true. + end if + net_index = net_index + 1 + end do + if (.not. found) return + + found = .false. ! Reuse for finding station ID below + do + if ((first .gt. last) .or. found) return + ! Use binary search + middle = (first + last) / 2 + if (plat_id .lt. this%platforms(middle)) then + last = middle - 1 + elseif (plat_id .gt. this%platforms(middle)) then + first = middle + 1 + else + ! Make sure country code matches + if (this%wmocode_id(middle) .eq. wmocode_id .and. & + this%fipscode_id(middle) .eq. fipscode_id) then + found = .true. + index = middle + end if + exit + end if + end do + + end function search_pcpobs + + ! Infer zero precip for certain durations if the accumulations are + ! missing and if present and past weather codes support no precip. + subroutine USAF_gages_use_preswx_pastwx(this) + + ! Defaults + implicit none + + ! Arguments + class(USAF_gages_t), intent(inout) :: this + + ! Locals + integer :: sweden_lowlimit + integer :: sweden_highlimit + integer :: denmark_lowlimit + integer :: denmark_highlimit + integer :: pastwx_duration + integer :: nobs + integer :: i + + if (this%cdms_flag) then + sweden_lowlimit = CDMS_SWEDEN_LOWLIMIT + sweden_highlimit = CDMS_SWEDEN_HIGHLIMIT + denmark_lowlimit = CDMS_DENMARK_LOWLIMIT + denmark_highlimit = CDMS_DENMARK_HIGHLIMIT + else + sweden_lowlimit = JMOBS_SWEDEN_LOWLIMIT + sweden_highlimit = JMOBS_SWEDEN_HIGHLIMIT + denmark_lowlimit = JMOBS_DENMARK_LOWLIMIT + denmark_highlimit = JMOBS_DENMARK_HIGHLIMIT + end if + + nobs = this%nobs + do i = 1, nobs + + ! If we don't have a valid past weather duration, we can't know + ! which precip accumulation to adjust. + pastwx_duration = this%pastwx_durations(i) + if (pastwx_duration .eq. MISSING) cycle + + if (this%preswx(i) .eq. MISSING) cycle + if (this%pastwx(i) .eq. MISSING) cycle + + ! We will only use past and present weather if the appropriate + ! accumulations are missing. + if (pastwx_duration .eq. 6) then + if (this%amts06(i) .ne. MISSING) cycle + if (this%amts03(i) .ne. MISSING) cycle + if (this%amts02(i) .ne. MISSING) cycle + if (this%amts01(i) .ne. MISSING) cycle + end if + + if (pastwx_duration .eq. 3) then + if (this%amts03(i) .ne. MISSING) cycle + if (this%amts02(i) .ne. MISSING) cycle + if (this%amts01(i) .ne. MISSING) cycle + end if + + ! Exclude Danish stations, since reports at non-standard + ! international hours have past weather of only 1 hour. + if ((this%wmonumbers(i) .ge. denmark_lowlimit .and. & + this%wmonumbers(i) .le. denmark_highlimit) .or. & + (this%bsn(i) .ge. denmark_lowlimit .and. & + this%bsn(i) .le. denmark_highlimit)) then + cycle + end if + + ! Exclude Swedish stations, since durations for past weather + ! vary from one to six hours depending on when the report was + ! made. + if ((this%wmonumbers(i) .ge. sweden_lowlimit .and. & + this%wmonumbers(i) .le. sweden_highlimit) .or. & + (this%bsn(i) .ge. sweden_lowlimit .and. & + this%bsn(i) .le. sweden_highlimit)) then + cycle + end if + + ! Exclude Antarctic stations. This is because (1) Australian + ! Antarctic stations are known to have varying past weather + ! durations due to staffing; (2) isolating the Australian + ! stations is difficult due to WMO station ID assignments; and + ! (3) Antarctic stations will probably be subfreezing and have + ! their precip accumulations rejected anyway. For simplicity, we + ! skip any report south of 60S. + if (this%lons(i) .le. -6000) cycle + + ! WMO SYNOP reports use different code tables depending on if the + ! station is manned or automatic. Unfortunately, the indicator + ! for manned vs automatic is not included in the USAF decoded + ! files, so we must be conservative and check for precipitation + ! codes from either scenario. The decision criteria are kept + ! separate in case the USAF file format is changed in the future + ! to include this information. In addition, BUFR has it's own + ! tables, and it is not clear from the USAF decoded files if BUFR + ! or SYNOP was decoded. + + ! First, check present weather for manned station (WMO Code Table + ! 4677; also WMO BUFR Code Table 0 20 003). + select case (this%preswx(i)) + case (50:99) ! Precipitation at station at time of observation + cycle + case (20:27, 29) ! Precip at station in last hour, but not at + cycle ! time of observation + end select + + ! Next, check present weather for automatic station (WMO Code + ! Table 4680; see also WMO BUFR Code Table 0 20 003) + select case (this%preswx(i)) + case (40:48, 50:58, 60:68, 70:78, 89, 90:96) + cycle ! Precipitation at station at time of observation + case (21:26) + cycle ! Precip at station in last hour, but not at time of + ! observation + case (140:148, 150:158, 160:168, 170:178, 180:188, 190:196) + cycle ! BUFR codes for precip at automated station at time of + ! observation + case (121:126) + cycle ! BUFR codes for precip at station in last hour, but + ! not at time of observation + case (250:257, 259, 260:267, 270:279, 280:291) + cycle ! More descriptive BUFR codes for precip at station. + case (510:511) + cycle ! BUFR codes for missing data + end select + + ! Next, check past weather for manned station (WMO Code Table + ! 4561) + select case (this%pastwx(i)) + case (5:9) + cycle + end select + + ! Next, check past weather for automatic station (WMO Code Table + ! 4531, and WMO BUFR Code Table 0 20 004 / 0 20 005) + select case (this%pastwx(i)) + case (4:9) ! WMO Code Table 4531 + cycle + case (14:19) ! BUFR Code Tables + cycle + case (31) ! BUFR code for missing + cycle + end select + + ! At this point, we can infer zero precipitation for the + ! appropriate duration + if (pastwx_duration .eq. 6) then + this%amts06(i) = 0 + this%amts03(i) = 0 + this%amts02(i) = 0 + this%amts01(i) = 0 + else if (pastwx_duration .eq. 3) then + this%amts03(i) = 0 + this%amts02(i) = 0 + this%amts01(i) = 0 + end if + end do + + call this%reconcile_self() + call this%fill_gaps() + end subroutine USAF_gages_use_preswx_pastwx + + ! Method for setting pastwx duration for each report. + subroutine USAF_gages_set_pastwx_durations(this) + implicit none + class(USAF_gages_t), intent(inout) :: this + integer :: nobs + character(10) :: date10 + integer :: i + nobs = this%nobs + do i = 1, nobs + if (this%YYYYMMDDhhmmss(i) .eq. "NULL") then + date10 = this%date10 + else + date10 = this%YYYYMMDDhhmmss(i)(1:10) + end if + this%pastwx_durations(i) = set_pastwx_duration(date10) + end do + end subroutine USAF_gages_set_pastwx_durations + + ! Function for setting duration of past weather report. This is based + ! on WMO SYNOP definition. (Two and one hour durations are also + ! allowed by the WMO, but these require SYNOP reports every two or one + ! hours; and it is not clear from the USAF decoded files when that + ! occurs.) + function set_pastwx_duration(date10) result (duration) + implicit none + character(10), intent(in) :: date10 + integer :: duration + select case (date10(9:10)) + case ('00', '06', '12', '18') + duration = 6 + case ('03', '09', '15', '21') + duration = 3 + case default + duration = MISSING + end select + end function set_pastwx_duration + + ! Return appropriate precipitation threshold based on duration. + ! NOTE: We currently skip 1 hour and 2-hour accumulations. + function set_precip_duration_threshold(duration) result (threshold) + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: duration + + ! Return variable + integer :: threshold + + ! Select appropriate threshold + select case (duration) + case (1) + threshold = MAX_PCP_1HR + case (2) + threshold = MAX_PCP_2HR + case (3) + threshold = MAX_PCP_3HR + case (6) + threshold = MAX_PCP_6HR + case (9) + threshold = MAX_PCP_9HR + case (12) + threshold = MAX_PCP_12HR + case (15) + threshold = MAX_PCP_15HR + case (18) + threshold = MAX_PCP_18HR + case (21) ! Not a WMO established duration, but we'll check anyway. + threshold = MAX_PCP_21HR + case (24) + threshold = MAX_PCP_24HR + case default + threshold = MISSING + end select + + end function set_precip_duration_threshold + + ! Return appropriate precipitation threshold based on report time in + ! India + function set_india_precip_threshold(utc_hour) result (threshold) + + ! Defaults + implicit none + + ! Arguments + character(2), intent(in) :: utc_hour + + ! Return variable + integer :: threshold + + ! Select appropriate threshold. India obs are assumed to accumulate + ! from 03Z + select case (utc_hour) + case ("04") + threshold = MAX_PCP_1HR + case ("05") + threshold = MAX_PCP_2HR + case ("06") + threshold = MAX_PCP_3HR + case ("09") + threshold = MAX_PCP_6HR + case ("12") + threshold = MAX_PCP_9HR + case ("15") + threshold = MAX_PCP_12HR + case ("18") + threshold = MAX_PCP_15HR + case ("21") + threshold = MAX_PCP_18HR + case ("00") + threshold = MAX_PCP_21HR + case ("03") + threshold = MAX_PCP_24HR + case default + threshold = MISSING + end select + + end function set_india_precip_threshold + + ! Method for copying appropriate data to a USAF_ObsData structure + ! for use in data assimilation. + subroutine USAF_copy_to_usaf_obsdata(this, hr, gage_sigma_o_sqr, & + precipObs) + + ! Imports + use LIS_logMod, only: LIS_logunit, LIS_endrun + use USAF_bratsethMod, only: USAF_ObsData, USAF_assignObsData + + ! Defaults + implicit none + + ! Arguments + class(USAF_Gages_t), intent(in) :: this + integer, intent(in) :: hr + real, intent(in) :: gage_sigma_o_sqr + type(USAF_ObsData), intent(inout) :: precipObs + + ! Locals + integer :: i + integer :: num_obs_copied + + ! Sanity checks + if (this%nobs == 0) return + if (hr .ne. 6 .and. hr .ne. 12) then + write(LIS_logunit,*) & + '[ERR] Invalid hour passed to USAF_copy_to_usaf_obsdata' + write(LIS_logunit,*) & + '[ERR] Should be 6 or 12, received ', hr + call LIS_endrun() + end if + + num_obs_copied = 0 + if (hr == 6) then + do i = 1, this%nobs + if (this%amts06(i) < 0) cycle + call USAF_assignObsData(precipObs, & + this%networks(i), & + this%platforms(i), & + real(this%amts06(i)) * 0.1, & + real(this%lats(i)) * 0.01, & + real(this%lons(i)) * 0.01, & + gage_sigma_o_sqr, 0.) + num_obs_copied = num_obs_copied + 1 + end do + write(LIS_logunit,*)'[INFO] Copied ', num_obs_copied, & + ' 6-hr gage reports' + else if (hr == 12) then + do i = 1, this%nobs + if (this%amts12(i) < 0) cycle + call USAF_assignObsData(precipObs, & + this%networks(i), & + this%platforms(i), & + real(this%amts12(i)) * 0.1, & + real(this%lats(i)) * 0.01, & + real(this%lons(i)) * 0.01, & + gage_sigma_o_sqr, 0.) + num_obs_copied = num_obs_copied + 1 + end do + write(LIS_logunit,*)'[INFO] Copied ', num_obs_copied, & + ' 12-hr gage reports' + end if + + + end subroutine USAF_copy_to_usaf_obsdata +end module USAF_GagesMod diff --git a/lis/metforcing/usaf/USAF_ImergHHMod.F90 b/lis/metforcing/usaf/USAF_ImergHHMod.F90 index f84c0b986..7aa7d9c52 100644 --- a/lis/metforcing/usaf/USAF_ImergHHMod.F90 +++ b/lis/metforcing/usaf/USAF_ImergHHMod.F90 @@ -111,8 +111,8 @@ subroutine copyToObsDataImergHHPrecip(this, sigmaOSqr, & type(ImergHHPrecip), intent(in) :: this real, intent(in) :: sigmaOSqr real, intent(in) :: oErrScaleLength - character*10, intent(in) :: net - character*10, intent(in) :: platform + character*32, intent(in) :: net + character*32, intent(in) :: platform type(USAF_ObsData), intent(inout) :: obsData_struc ! Local variables diff --git a/lis/metforcing/usaf/USAF_PreobsReaderMod.F90 b/lis/metforcing/usaf/USAF_PreobsReaderMod.F90 new file mode 100644 index 000000000..cf2f362df --- /dev/null +++ b/lis/metforcing/usaf/USAF_PreobsReaderMod.F90 @@ -0,0 +1,1128 @@ +! +! MODULE: USAF_PreobsReaderMod +! +! DESCRIPTION: Contains code for reading USAF preobs files, performing simple +! preprocessing, and adding to gages database. +! +! AUTHOR: Eric Kemp, SSAI/NASA GSFC +! +!------------------------------------------------------------------------------ + +module USAF_PreobsReaderMod + + ! Defaults + implicit none + private + + ! Public routines + public :: USAF_read_preobs + + integer, parameter :: MISSING = -99999999 + +contains + + ! Read preobs files, perform simple preprocessing, and store + ! in database. + subroutine USAF_read_preobs(preobsdir, presavdir, & + use_timestamp, & + year, month, day, hour, use_expanded_station_ids, & + alert_number) + + ! Imports + use ESMF + use LIS_coreMod, only: LIS_masterproc + use LIS_logMod, only: LIS_logunit, LIS_alert, LIS_getNextUnitNumber, & + LIS_releaseUnitNumber + use LIS_mpiMod, only: LIS_mpi_comm + use USAF_GagesMod, only: USAF_Gages_t + + ! Defaults + implicit none + + ! Arguments + character(*), intent(in) :: preobsdir + character(*), intent(in) :: presavdir + integer, intent(in) :: use_timestamp + integer, intent(in) :: year + integer, intent(in) :: month + integer, intent(in) :: day + integer, intent(in) :: hour + integer, intent(in) :: use_expanded_station_ids + integer, intent(inout) :: alert_number + + ! Locals + character(255) :: filename + integer, allocatable :: twfprc(:) + integer, allocatable :: duration(:) + integer, allocatable :: sixprc(:) + integer, allocatable :: mscprc(:) + integer, allocatable :: ilat(:) + integer, allocatable :: ilon(:) + integer, allocatable :: bsn(:) + character(32), allocatable :: network(:) + character(32), allocatable :: plat_id(:) + character(2), allocatable :: wmocode_id(:) + character(2), allocatable :: fipscode_id(:) + integer, allocatable :: pastwx(:) + integer, allocatable :: preswx(:) + integer, allocatable :: wmoblk(:) + integer :: twfprc_tmp + integer :: duration_tmp + integer :: sixprc_tmp + integer :: mscprc_tmp + integer :: ilat_tmp + integer :: ilon_tmp + character(14) :: YYYYMMDDhhmmss_tmp + character(32) :: network_tmp + character(32) :: plat_id_tmp + character(2) :: wmocode_id_tmp, fipscode_id_tmp + integer :: pastwx_tmp + integer :: preswx_tmp + integer :: wmoblk_tmp + integer :: bsn_tmp + integer :: nsize + integer :: nsize_total + integer :: ihemi + logical :: found_file + integer :: ierr + integer :: stncnt + integer :: i + character(10) :: date10 + character(14), allocatable :: YYYYMMDDhhmmss(:) + integer, allocatable :: amts24(:) + integer, allocatable :: amts21(:) + integer, allocatable :: amts18(:) + integer, allocatable :: amts15(:) + integer, allocatable :: amts12(:) + integer, allocatable :: amts09(:) + integer, allocatable :: amts06(:) + integer, allocatable :: amts03(:) + integer, allocatable :: amts02(:) + integer, allocatable :: amts01(:) + integer, allocatable :: amts00(:) + type(USAF_Gages_t) :: obscur, obsprev + integer :: rc + character(255) :: presav_filename + integer :: ipass, j + logical :: exchanges + type(ESMF_Time) :: curtime, prevtime, reporttime + type(ESMF_TimeInterval) :: deltatime, maxdeltatime + integer :: prevyear, prevmonth, prevday, prevhour + logical :: file_exists + integer :: deltahr + character(10) :: prevdate10 + integer :: yyyy, mm, dd, h, m, s + character(255) :: timestring + integer :: iunit + character(255) :: message(20) + + message = '' + + write(LIS_logunit,*)'---------------------------------------------' + + ! Set time limit allowed for report + call esmf_timeintervalset(maxdeltatime, m=5, rc=rc) ! 5 min after + + ! Find the total number of entries in the two hemispheric preobs + ! files for this date/time. This will be the upper-limit for how + ! much memory we allocate for temporary storage. If the newer + ! global preobs file is read, logic in the loop will break out + ! of the loop before the second pass. + nsize_total = 0 + do ihemi = 1, 2 + + call get_preobs_filename(filename, preobsdir, & + use_timestamp, & + ihemi, year, month, day, hour, use_expanded_station_ids) + + inquire(file=trim(filename), exist=found_file) + if (.not. found_file) then + write(LIS_logunit,*) '[WARN] Cannot find ', trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_read_preobs' + message(3) = ' Cannot find file ' // trim(filename) + if (LIS_masterproc) then + alert_number = alert_number + 1 + call LIS_alert('LIS.USAF_read_preobs', & + alert_number, message) + end if + if (use_expanded_station_ids == 1) exit ! These files are global + cycle + end if + + iunit = LIS_getNextUnitNumber() + open(iunit, file=trim(filename), status='old', iostat=ierr) + if (ierr .ne. 0) then + write(LIS_logunit,*) '[WARN] Problem opening ', trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_read_preobs' + message(3) = ' Cannot open file ' // trim(filename) + if (LIS_masterproc) then + alert_number = alert_number + 1 + call LIS_alert('LIS.USAF_read_preobs', & + alert_number, message) + end if + if (use_expanded_station_ids == 1) exit ! These files are global + cycle + end if + + nsize = 0 + read(iunit, *, iostat=ierr) nsize + if (ierr .ne. 0) then + write(LIS_logunit,*) '[WARN] Problem reading ', trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_read_preobs' + message(3) = ' Problem reading file ' // trim(filename) + if (LIS_masterproc) then + alert_number = alert_number + 1 + call LIS_alert('LIS.USAF_read_preobs', & + alert_number, message) + end if + close(iunit) + call LIS_releaseUnitNumber(iunit) + if (use_expanded_station_ids == 1) exit ! These files are global + cycle + end if + + if (nsize == 0) then + write(LIS_logunit,*)'[WARN] No precip obs found in ', & + trim(filename) + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_read_preobs' + message(3) = ' No precip obs found in ' // trim(filename) + if (LIS_masterproc) then + alert_number = alert_number + 1 + call LIS_alert('LIS.USAF_read_preobs', & + alert_number, message) + end if + else + write(LIS_logunit,*) '[INFO] Will process ', trim(filename) + end if + + nsize_total = nsize_total + nsize + close(iunit) + call LIS_releaseUnitNumber(iunit) + + if (use_expanded_station_ids == 1) exit ! These files are global + end do + + if (nsize_total == 0) then + write(LIS_logunit,*) '[WARN] No precip obs available!' + return + end if + + ! We now have an upper limit of how many gage reports to save. + allocate(twfprc(nsize_total)) + twfprc = MISSING + allocate(duration(nsize_total)) + duration = MISSING + allocate(sixprc(nsize_total)) + sixprc = MISSING + allocate(mscprc(nsize_total)) + mscprc = MISSING + allocate(ilat(nsize_total)) + ilat = MISSING + allocate(ilon(nsize_total)) + ilon = MISSING + allocate(bsn(nsize_total)) + bsn = MISSING + allocate(network(nsize_total)) + network = "NULL" + allocate(plat_id(nsize_total)) + plat_id = "NULL" + allocate(wmocode_id(nsize_total)) + wmocode_id = "??" + allocate(fipscode_id(nsize_total)) + fipscode_id = "??" + allocate(pastwx(nsize_total)) + pastwx = MISSING + allocate(preswx(nsize_total)) + preswx = MISSING + allocate(wmoblk(nsize_total)) + wmoblk = MISSING + allocate(YYYYMMDDhhmmss(nsize_total)) + YYYYMMDDhhmmss = "NULL" + + write(date10, "(i4.4,i2.2,i2.2,i2.2)") year, month, day, hour + + ! Need to set current time in ESMF for calculations + call esmf_timeset(curtime, yy=year, mm=month, dd=day, h=hour, & + m=0, s=0, rc=rc) + + ! Now begin reading the data in the files, performing simple checks + ! along the way. + stncnt = 0 + do ihemi = 1, 2 + call get_preobs_filename(filename, preobsdir, & + use_timestamp, & + ihemi, year, month, day, hour, use_expanded_station_ids) + + inquire(file=trim(filename), exist=found_file) + if (.not. found_file) then + if (use_expanded_station_ids == 1) exit + cycle + end if + + iunit = LIS_getNextUnitNumber() + open(iunit, file=trim(filename), status='old', iostat=ierr) + if (ierr .ne. 0) cycle + + nsize = 0 + read(iunit, *, iostat=ierr) nsize + if (ierr .ne. 0) then + close(iunit) + call LIS_releaseUnitNumber(iunit) + if (use_expanded_station_ids == 1) exit + cycle + end if + + ! Start reading the actual obs + do i = 1, nsize + if (use_expanded_station_ids == 1) then + twfprc_tmp = MISSING + sixprc_tmp = MISSING + read(iunit, 6001, iostat=ierr) YYYYMMDDhhmmss_tmp, & + network_tmp, plat_id_tmp, ilat_tmp, ilon_tmp, & + wmocode_id_tmp, fipscode_id_tmp, wmoblk_tmp, & + mscprc_tmp, duration_tmp, pastwx_tmp, preswx_tmp +6001 format(1x, a14, 1x, a10, 1x, a32, 1x, i9, 1x, i9, 1x, & + a2, 1x, a2, 1x, i9, 1x, i9, 1x, i9, 1x, i9, 1x, i9) + if (wmocode_id_tmp == " ") wmocode_id_tmp = "??" + if (fipscode_id_tmp == " ") fipscode_id_tmp = "??" + else + read(iunit, 6000, iostat=ierr) twfprc_tmp, duration_tmp, & + sixprc_tmp, & + mscprc_tmp, ilat_tmp, ilon_tmp, network_tmp, & + plat_id_tmp, & + pastwx_tmp, preswx_tmp, wmoblk_tmp +6000 format(i10, i10, i10, i10, i10, i10, 1x, a10, 2x, a10 & + ,i10, i10, i10) + wmocode_id_tmp = "??" + fipscode_id_tmp = "??" + YYYYMMDDhhmmss_tmp = date10 // "0000" + end if + + if (ierr .ne. 0) then + write(LIS_logunit,*) & + '[WARN] Problem reading report, skipping line...' + cycle + end if + + ! Skip if lat/lon is 0 (this is interpreted as missing). + if (ilat_tmp == 0 .and. ilon_tmp == 0) cycle + + ! Skip reports that are too much after the analysis time + ! (but allow earlier reports). This is a crude way of + ! allowing for Australian reports that are sometimes one or + ! two hours behind the sub-synoptic times. + if (use_expanded_station_ids == 1) then + read(YYYYMMDDhhmmss_tmp, & + '(I4.4, I2.2, I2.2, I2.2, I2.2, I2.2)') yyyy, mm, dd, & + h, m, s + call esmf_timeset(reporttime, yy=yyyy, mm=mm, dd=dd, & + h=h, m=m, s=s, rc=rc) + deltatime = reporttime - curtime + + call esmf_timeintervalget(deltatime, timestring=timestring, & + rc=rc) + + ! if (deltatime < mindeltatime) cycle + if (deltatime > maxdeltatime) cycle + end if + + ! Don't save 1-hr precip. + ! FIXME...Add special handling in new preobs files to + ! estimate 3-hr accums if enough 1-hr reports are available. + if (duration_tmp .eq. 1) then + mscprc_tmp = MISSING + duration_tmp = MISSING + end if + + ! Skip if no useable data in report. + if (twfprc_tmp == MISSING .and. & + sixprc_tmp == MISSING .and. & + mscprc_tmp == MISSING .and. & + (pastwx_tmp == MISSING .or. & + preswx_tmp == MISSING)) then + cycle + end if + + ! If we don't have FIPS or WMO codes available, try guessing + ! the FIPS code for countries with special rules based on + ! the WIGOS issuer number. + if (fipscode_id_tmp == "??" .and. & + wmocode_id_tmp == "??") then + fipscode_id_tmp = get_fips_from_wigos_issuer(plat_id_tmp) + end if + + ! Set the numeric bsn field. + if (network_tmp .eq. "WMO") then + bsn_tmp = set_bsn_wmo(use_expanded_station_ids, & + plat_id_tmp, fipscode_id_tmp, wmocode_id_tmp) + else if (network_tmp .eq. "CDMS") then + read (plat_id_tmp, '(i32)') bsn_tmp + else + bsn_tmp = 0 + end if + + ! Set the country codes if possible. + if (network_tmp .eq. "CANA") then + wmocode_id_tmp = "CA" + fipscode_id_tmp = "CA" + else if (network_tmp .eq. "FAA") then + wmocode_id_tmp = "US" + fipscode_id_tmp = "US" + end if + + ! Sometimes 6 and 24-hr data are stored in misc. Copy + ! over if necessary, then erase. + if (duration_tmp == 24) then + if (twfprc_tmp == MISSING) then + twfprc_tmp = mscprc_tmp + end if + mscprc_tmp = MISSING + duration_tmp = MISSING + else if (duration_tmp == 6) then + if (sixprc_tmp == MISSING) then + sixprc_tmp = mscprc_tmp + end if + mscprc_tmp = MISSING + duration_tmp = MISSING + end if + + ! If this is the first report, save it. + if (stncnt == 0) then + stncnt = stncnt + 1 + twfprc(stncnt) = twfprc_tmp + duration(stncnt) = duration_tmp + sixprc(stncnt) = sixprc_tmp + mscprc(stncnt) = mscprc_tmp + ilat(stncnt) = ilat_tmp + ilon(stncnt) = ilon_tmp + bsn(stncnt) = bsn_tmp + network(stncnt) = network_tmp + plat_id(stncnt) = plat_id_tmp + wmocode_id(stncnt) = wmocode_id_tmp + fipscode_id(stncnt) = fipscode_id_tmp + pastwx(stncnt) = pastwx_tmp + preswx(stncnt) = preswx_tmp + wmoblk(stncnt) = wmoblk_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + cycle + else + ! Look at the most recent saved gage report, and see if this + ! new entry is from the same station. + if (network(stncnt) .eq. network_tmp .and. & + trim(plat_id(stncnt)) .eq. trim(plat_id_tmp)) then + + ! For newer file format with date/time with each report, + ! pick the more recent date/time for duration other + ! than 1 hour. + if (YYYYMMDDhhmmss(stncnt) < YYYYMMDDhhmmss_tmp .and. & + (twfprc_tmp .ne. MISSING .or. & + sixprc_tmp .ne. MISSING .or. & + mscprc_tmp .ne. MISSING .or. & + (preswx_tmp .ne. MISSING .and. & + pastwx_tmp .ne. MISSING))) then + twfprc(stncnt) = twfprc_tmp + duration(stncnt) = duration_tmp + sixprc(stncnt) = sixprc_tmp + mscprc(stncnt) = mscprc_tmp + ilat(stncnt) = ilat_tmp + ilon(stncnt) = ilon_tmp + bsn(stncnt) = bsn_tmp + network(stncnt) = network_tmp + plat_id(stncnt) = plat_id_tmp + wmocode_id(stncnt) = wmocode_id_tmp + fipscode_id(stncnt) = fipscode_id_tmp + pastwx(stncnt) = pastwx_tmp + preswx(stncnt) = preswx_tmp + wmoblk(stncnt) = wmoblk_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + cycle + end if + + ! Logic for older file format. Code above assigns + ! the file date/time as the report date/time. + ! See if there is any additional information worth + ! saving. The AGRMET logic overwrites the older data + ! with the newer if the newer data is a higher amount. + ! An exception is that "zero" six-hour precip is ignored + ! if 24-hr or 12-hr precip is positive, apparently + ! because "zero" six-hour precip reports in this case are + ! typically from an asynoptic ob. We follow that logic + ! here. + if (YYYYMMDDhhmmss(stncnt) == YYYYMMDDhhmmss_tmp) then + if (twfprc(stncnt) < twfprc_tmp) then + twfprc(stncnt) = twfprc_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + end if + + if (mscprc(stncnt) .eq. MISSING .and. & + mscprc_tmp .ne. MISSING) then + ! We don't have a misc precip stored yet, so do it + ! here. + mscprc(stncnt) = mscprc_tmp + duration(stncnt) = duration_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + else if (duration_tmp > duration(stncnt) .and. & + duration_tmp > 1 .and. & + mscprc_tmp .ne. MISSING) then + ! Try to capture 3-hr or other longer duration + ! precip. + mscprc(stncnt) = mscprc_tmp + duration(stncnt) = duration_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + else if (mscprc_tmp > mscprc(stncnt) .and. & + duration_tmp == duration(stncnt)) then + ! Somewhat mimics AGRMET. Save larger accum for + ! current duration. + mscprc(stncnt) = mscprc_tmp + duration(stncnt) = duration_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + else if ((mscprc_tmp > mscprc(stncnt)) .and. & + (duration_tmp == MISSING .or. & + duration_tmp == 0) .and. & + (duration(stncnt) == MISSING .or. & + duration(stncnt) == 0)) then + ! Somewhat mimics AGRMET. Save larger accum for + ! undefined duration, for use in certain regions + ! (e.g., Russia) + mscprc(stncnt) = mscprc_tmp + duration(stncnt) = duration_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + end if + if ( (twfprc(stncnt) > 0) .or. & + (duration(stncnt) == 12 .and. & + mscprc(stncnt) > 0) ) then + if (sixprc(stncnt) == MISSING .and. & + sixprc_tmp == 0) then + sixprc(stncnt) = MISSING + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + else if (sixprc(stncnt) == 0 .and. & + sixprc_tmp == MISSING) then + sixprc(stncnt) = MISSING + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + else + if (sixprc(stncnt) < sixprc_tmp) then + sixprc(stncnt) = sixprc_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + end if + end if + else + if (sixprc(stncnt) < sixprc_tmp) then + sixprc(stncnt) = sixprc_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + end if + end if + end if + else + ! This isn't the same station, so save it. + !if (trim(network_tmp) == "AMIL" .and. & + ! trim(plat_id_tmp) == "KQTC") then + ! write(LIS_logunit,*)'EMK: KQTC, twfprc = ', twfprc_tmp + !end if + + stncnt = stncnt + 1 + twfprc(stncnt) = twfprc_tmp + duration(stncnt) = duration_tmp + sixprc(stncnt) = sixprc_tmp + mscprc(stncnt) = mscprc_tmp + ilat(stncnt) = ilat_tmp + ilon(stncnt) = ilon_tmp + bsn(stncnt) = bsn_tmp + network(stncnt) = network_tmp + plat_id(stncnt) = plat_id_tmp + wmocode_id(stncnt) = wmocode_id_tmp + fipscode_id(stncnt) = fipscode_id_tmp + pastwx(stncnt) = pastwx_tmp + preswx(stncnt) = preswx_tmp + wmoblk(stncnt) = wmoblk_tmp + YYYYMMDDhhmmss(stncnt) = YYYYMMDDhhmmss_tmp + cycle + end if + end if + end do + + if (use_expanded_station_ids == 1) cycle ! These files are global + end do ! ihemi + + ! Since we combined both the NH and SH files, the resulting list + ! is unsorted. + ! Bubble sort the reports by network. + if (use_expanded_station_ids == 0) then + ipass = 1 + exchanges = .true. + do while ((ipass < stncnt) .and. exchanges) + exchanges = .false. + do j = 1, stncnt - ipass + if (network(j) > network(j+1)) then + call swap_int(twfprc(j), twfprc(j+1)) + call swap_int(duration(j), duration(j+1)) + call swap_int(sixprc(j), sixprc(j+1)) + call swap_int(mscprc(j), mscprc(j+1)) + call swap_int(ilat(j), ilat(j+1)) + call swap_int(ilon(j), ilon(j+1)) + call swap_int(bsn(j), bsn(j+1)) + call swap_char32(network(j), network(j+1)) + call swap_char32(plat_id(j), plat_id(j+1)) + call swap_char2(wmocode_id(j), wmocode_id(j+1)) + call swap_char2(fipscode_id(j), fipscode_id(j+1)) + call swap_int(pastwx(j), pastwx(j+1)) + call swap_int(preswx(j), preswx(j+1)) + call swap_int(wmoblk(j), wmoblk(j+1)) + call swap_char14(YYYYMMDDhhmmss(j), YYYYMMDDhhmmss(j+1)) + exchanges = .true. + end if + end do + ipass = ipass + 1 + end do + + ! Now bubble sort by station ID within the same network + ipass = 1 + exchanges = .true. + do while ((ipass < stncnt) .and. exchanges) + exchanges = .false. + do j = 1, stncnt - ipass + if ((network(j) == network(j+1)) .and. & + (plat_id(j) > plat_id(j+1))) then + call swap_int(twfprc(j), twfprc(j+1)) + call swap_int(duration(j), duration(j+1)) + call swap_int(sixprc(j), sixprc(j+1)) + call swap_int(mscprc(j), mscprc(j+1)) + call swap_int(ilat(j), ilat(j+1)) + call swap_int(ilon(j), ilon(j+1)) + call swap_int(bsn(j), bsn(j+1)) + call swap_char32(network(j), network(j+1)) + call swap_char32(plat_id(j), plat_id(j+1)) + call swap_char2(wmocode_id(j), wmocode_id(j+1)) + call swap_char2(fipscode_id(j), fipscode_id(j+1)) + call swap_int(pastwx(j), pastwx(j+1)) + call swap_int(preswx(j), preswx(j+1)) + call swap_int(wmoblk(j), wmoblk(j+1)) + call swap_char14(YYYYMMDDhhmmss(j), YYYYMMDDhhmmss(j+1)) + exchanges = .true. + end if + end do + ipass = ipass + 1 + end do + end if + + ! Report if any WMO obs are not legacy 5-digit SYNOP + do j = 1, stncnt + if (network(j) .ne. "WMO") cycle + if (.not. len_trim(plat_id(j)) == 5) then + write(LIS_logunit,*) '[INFO] Found non SYNOP WMO station ', & + trim(plat_id(j)) + else if (verify(trim(plat_id(j)), "0123456789") .ne. 0) then + write(LIS_logunit,*) '[INFO] Found non SYNOP WMO station ', & + trim(plat_id(j)) + end if + end do + + ! Copy into Gage_type object + allocate(amts24(nsize_total)) ; amts24 = twfprc + allocate(amts21(nsize_total)) ; amts21 = MISSING + allocate(amts18(nsize_total)) ; amts18 = MISSING + allocate(amts15(nsize_total)) ; amts15 = MISSING + allocate(amts12(nsize_total)) ; amts12 = MISSING + allocate(amts09(nsize_total)) ; amts09 = MISSING + allocate(amts06(nsize_total)) ; amts06 = sixprc + allocate(amts03(nsize_total)) ; amts03 = MISSING + allocate(amts02(nsize_total)) ; amts02 = MISSING + allocate(amts01(nsize_total)) ; amts01 = MISSING + allocate(amts00(nsize_total)) ; amts00 = mscprc + + call obscur%new(date10, stncnt, YYYYMMDDhhmmss, network, plat_id, & + wmocode_id, fipscode_id, wmoblk, bsn, ilat, ilon, & + amts24, amts21, amts18, amts15, amts12, amts09, amts06, & + amts03, amts02, amts01, amts00, duration, preswx, pastwx) + + ! Clean up the temporary arrays + deallocate(twfprc) + deallocate(duration) + deallocate(sixprc) + deallocate(mscprc) + deallocate(ilat) + deallocate(ilon) + deallocate(bsn) + deallocate(network) + deallocate(plat_id) + deallocate(wmocode_id) + deallocate(fipscode_id) + deallocate(pastwx) + deallocate(preswx) + deallocate(wmoblk) + deallocate(YYYYMMDDhhmmss) + deallocate(amts24) + deallocate(amts21) + deallocate(amts18) + deallocate(amts15) + deallocate(amts12) + deallocate(amts09) + deallocate(amts06) + deallocate(amts03) + deallocate(amts02) + deallocate(amts01) + deallocate(amts00) + + ! Check for unphysical values. + call obscur%check_gross_errors() + + ! Reconcile values in same report (shorter durations should not be + ! larger than longer durations) + call obscur%reconcile_self() + + ! Fill gaps if possible + call obscur%fill_gaps() + + ! Use the miscellaneous precip, with special rules for some + ! countries. + call obscur%use_misc_precip() + + ! Leverage present and past weather information to identify + ! zero precip durations. + !EMK...Disable until more accurate inspection of time in new + !preobs files is implemented. + !call obscur%use_preswx_pastwx() + + ! Try fetching presav files from earlier hours and reconcile. + write(LIS_logunit,*)'[INFO] Will compare with earlier gage reports.' + do deltahr = 3, 21, 3 + call esmf_timeintervalset(deltatime, yy=0, mm=0, d=0, & + h=deltahr, m=0, s=0, rc=rc) + prevtime = curtime - deltatime + call esmf_timeget(prevtime, yy=prevyear, mm=prevmonth, & + dd=prevday, h=prevhour, rc=rc) + write(presav_filename,'(A,A,i4.4,i2.2,i2.2,i2.2)') & + trim(presavdir), '/presav2.03hr.', & + prevyear, prevmonth, prevday, prevhour + inquire(file=presav_filename, exist=file_exists) + if (file_exists) then + write(LIS_logunit,*)'[INFO] Comparing against data in ', & + trim(presav_filename) + write(prevdate10,'(i4.4,i2.2,i2.2,i2.2)') & + prevyear, prevmonth, prevday, prevhour + call obsprev%read_data(presav_filename, prevdate10, & + alert_number) + if (deltahr == 3) then + call obscur%reconcile_gages03(obsprev) + else if (deltahr == 6) then + call obscur%reconcile_gages06(obsprev) + else if (deltahr == 9) then + call obscur%reconcile_gages09(obsprev) + else if (deltahr == 12) then + call obscur%reconcile_gages12(obsprev) + else if (deltahr == 15) then + call obscur%reconcile_gages15(obsprev) + else if (deltahr == 18) then + call obscur%reconcile_gages18(obsprev) + else if (deltahr == 21) then + call obscur%reconcile_gages21(obsprev) + end if + call obsprev%delete() + else + write(LIS_logunit,*)'[WARN] Cannot find ', trim(presav_filename) + write(LIS_logunit,*) & + '[WARN] Will skip reconciling with obs from ', & + abs(deltahr),' hours ago' + message(1) = '[WARN] Program: LIS' + message(2) = ' Routine: USAF_read_preobs' + message(3) = ' Cannot find earlier presav2 file ' // & + trim(presav_filename) + message(4) = ' Observation count will be reduced' + if (LIS_masterproc) then + alert_number = alert_number + 1 + call LIS_alert('LIS.USAF_read_preobs', & + alert_number, message) + end if + end if + end do + + ! Correct for overnight reporting issues in South America + write(LIS_logunit,*) & + '[INFO] Correcting overnight reporting in South America' + call obscur%correct_region3_12z() + + ! Have the master process write the data to file. + if (LIS_masterproc) then + write(presav_filename,'(A,A,i4.4,i2.2,i2.2,i2.2)') & + trim(presavdir), '/presav2.03hr.', & + year, month, day, hour + write(LIS_logunit,*)'[INFO] Writing to ', trim(presav_filename) + call obscur%write_data(presav_filename) + end if +#if (defined SPMD) + call MPI_Barrier(LIS_mpi_comm, ierr) +#endif + + ! Clean up + call obscur%delete() + + end subroutine USAF_read_preobs + + ! Generate name of preobs file. Based on code in AGRMET_getpcpobs.F90 + subroutine get_preobs_filename(filename, preobsdir, & + use_timestamp, & + ihemi, year, month, day, hour, use_expanded_station_ids) + + ! Defaults + implicit none + + ! Arguments + character(*), intent(out) :: filename + character(*), intent(in) :: preobsdir + integer, intent(in) :: use_timestamp + integer, intent(in) :: ihemi + integer, intent(in) :: year + integer, intent(in) :: month + integer, intent(in) :: day + integer, intent(in) :: hour + integer, intent(in) :: use_expanded_station_ids + + ! Locals + character(8) :: ftime1 + character(10) :: ftime2 + character(2), parameter :: FHEMI(2) = (/'nh', 'sh'/) + + write(ftime2, '(i4.4, i2.2, i2.2, i2.2)') year, month, day, hour + + if (use_expanded_station_ids == 0) then + !if (use_timestamp == 1) then + ! write(ftime1, '(i4.4, i2.2, i2.2)') & + ! year, month, day + ! filename = ftime1 // '/' // trim(preobsdir) // & + ! '/preobs_' // FHEMI(ihemi) // '.03hr.' // ftime2 + !else + filename = trim(preobsdir) // & + '/preobs_' // FHEMI(ihemi) // '.03hr.' // ftime2 + !endif + else if (use_expanded_station_ids == 1) then + !if (use_timestamp == 1) then + ! write(ftime1, '(i4.4, i2.2, i2.2)') & + ! year, month, day + ! filename = ftime1 // '/' // trim(preobsdir) // & + ! '/preobs_03hr_' // ftime2 // ".txt" + !else + filename = trim(preobsdir) // & + '/preobs_03hr_' // ftime2 // ".txt" + !endif + end if + end subroutine get_preobs_filename + + ! Set block number for a WMO station + function set_bsn_wmo(use_expanded_station_id, plat_id, fipscode_id, & + wmocode_id) result(bsn) + + ! Imports + use USAF_GagesMod, only: JMOBS_SWEDEN_LOWLIMIT, & + JMOBS_DENMARK_LOWLIMIT, JMOBS_RUSSIA_LOWLIMIT, & + JMOBS_INDIA_LOWLIMIT, JMOBS_SRILANKA_LOWLIMIT, & + JMOBS_S_AMER_LOWLIMIT + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: use_expanded_station_id + character(32), intent(in) :: plat_id + character(2), intent(in) :: fipscode_id + character(2), intent(in) :: wmocode_id + + ! Return variable + integer :: bsn + + bsn = 0 + + ! Old preobs files don't have any country codes. Just copy the + ! 5-digit ID and use as is. + if (use_expanded_station_id == 0) then + read(plat_id, '(i5.5)') bsn + return + end if + + ! If this is a legacy 5-digit WMO ID, just copy it. + if (verify(trim(plat_id), "0123456789") == 0 & + .and. len_trim(plat_id) == 5) then + read(plat_id, '(i5.5)') bsn + return + end if + + ! WIGOS supports legacy 5-digit WMO IDs with a special prefix + if (len_trim(plat_id) == 15) then + if (plat_id(1:10) == "0-20000-0-" .and. & + verify(trim(plat_id(11:15)), "0123456789") == 0) then + read(plat_id(11:15), '(i5.5)') bsn + return + end if + end if + + ! Try using the FIPS country code. This code seems to be more + ! accurate than the WMO code. + select case(fipscode_id) + case('SW') + bsn = JMOBS_SWEDEN_LOWLIMIT + return + case('DA') + bsn = JMOBS_DENMARK_LOWLIMIT + return + case('AM', 'AJ', 'BO', 'EN', 'GG', 'LG', 'LH', 'KZ', & + 'KG', 'MD', 'RS', 'TI', 'TX', 'UP', 'UZ') + ! AM is Armenia + ! AJ is Azerbaijan + ! BO is Belarus + ! EN is Estonia + ! GG is Georgia + ! LG is Latvia + ! LH is Lituania + ! KZ is Kazakhstan + ! KG is Kyrgyzstan + ! MD is Moldova + ! RS is Russia + ! TI is Tajikistan + ! TX is Turkmenistan + ! UP is Ukraine + ! UZ is Uzbekistan + bsn = JMOBS_RUSSIA_LOWLIMIT + return + case ("IN") + bsn = JMOBS_INDIA_LOWLIMIT + return + case ("CE") + bsn = JMOBS_SRILANKA_LOWLIMIT + return + case ('AR', 'BL', 'BR', 'CI', 'CO', 'EC', 'FK', 'FG', 'GY', & + 'PA', 'PE', 'SX', 'NS', 'UY', 'VE') + ! AR is Argentina + ! BL is Bolivia + ! BR is Brazil + ! CI is Chile + ! CO is Colombia + ! EC is Ecuador + ! FK is Falkland Islands + ! FG is French Guiana + ! GY is Guyana + ! PA is Paraguay + ! PE is Peru + ! SX is South Georgia and the South Sandwich Islands + ! NS is Suriname + ! UY is Uruguay + ! VE is Venezuela + bsn = JMOBS_S_AMER_LOWLIMIT + return + case default + bsn = 0 + end select + + ! Try using the WMO country code + select case(wmocode_id) + case('SE') + bsn = JMOBS_SWEDEN_LOWLIMIT + return + case('DK') + bsn = JMOBS_DENMARK_LOWLIMIT + return + case('AM', 'AZ', 'BY', 'EE', 'GE', 'LV', 'LT', 'KZ', 'KG', & + 'MD', 'RU', 'TJ', 'TM', 'UA', 'UZ') + ! AM is Armenia + ! AZ is Azerbaijan + ! BY is Belarus + ! EE is Estonia + ! GE is Georgia + ! LV is Latvia + ! LT is Lithuania + ! KZ is Kazakhstan + ! KG is Kyrgyzstan + ! MD is Moldova + ! RU is Russia + ! TJ is Tajikistan + ! TM is Turkmenistan + ! UA is Ukraine + ! UZ is Uzbekistan + bsn = JMOBS_RUSSIA_LOWLIMIT + return + case ("IN") + bsn = JMOBS_INDIA_LOWLIMIT + return + case ("LK") + bsn = JMOBS_SRILANKA_LOWLIMIT + return + case ('AR', 'BO', 'BR', 'CL', 'CO', 'EC', 'FK', 'GF', 'GY', & + 'PY', 'PE', 'GS', 'SR', 'UG', 'VE') + ! AR is Argentina + ! BO is Bolivia + ! BR is Brazil + ! CL is Chile + ! CO is Colombia + ! EC is Ecuador + ! FK is Falkland Islands + ! GF is French Guiana + ! GY is Guyana + ! PY is Paraguay + ! PE is Peru + ! GS is South Georgia and the South Sandwich Islands + ! SR is Suriname + ! UG is Uruguay + ! VE is Venezuela + bsn = JMOBS_S_AMER_LOWLIMIT + return + case default + bsn = 0 + end select + + end function set_bsn_wmo + + subroutine swap_int(var1, var2) + implicit none + integer, intent(inout) :: var1 + integer, intent(inout) :: var2 + integer :: tmp + tmp = var1 + var1 = var2 + var2 = tmp + end subroutine swap_int + + subroutine swap_char2(var1, var2) + implicit none + character(2), intent(inout) :: var1 + character(2), intent(inout) :: var2 + character(2) :: tmp + tmp = var1 + var1 = var2 + var2 = tmp + end subroutine swap_char2 + + subroutine swap_char10(var1, var2) + implicit none + character(10), intent(inout) :: var1 + character(10), intent(inout) :: var2 + character(10) :: tmp + tmp = var1 + var1 = var2 + var2 = tmp + end subroutine swap_char10 + + subroutine swap_char14(var1, var2) + implicit none + character(14), intent(inout) :: var1 + character(14), intent(inout) :: var2 + character(14) :: tmp + tmp = var1 + var1 = var2 + var2 = tmp + end subroutine swap_char14 + + subroutine swap_char32(var1, var2) + implicit none + character(32), intent(inout) :: var1 + character(32), intent(inout) :: var2 + character(32) :: tmp + tmp = var1 + var1 = var2 + var2 = tmp + end subroutine swap_char32 + + ! Try to set FIPS country code based on WIGOS issuer number. Only + ! certain specific countries are checked for (for special accumulation + ! rules). + function get_fips_from_wigos_issuer(stn) result (fips) + + ! Defaults + implicit none + + ! Arguments + character(*),intent(in) :: stn + + ! Return variable + character(2) :: fips + + ! Locals + integer :: idx1, idx2 + integer :: id + + fips = '??' ! First guess + + ! See if this is a WIGOS station ID; if so, the issuer code is the + ! integer string between the first two dashes. If this is not + ! WIGOS, just return. + idx1 = scan(stn,'-') + if (idx1 == 0) return + idx2 = scan(stn(idx1+1:len(stn)),'-') + if (idx2 == 0) return + if (verify(trim(stn(idx1+1:idx2+1)), "0123456789") .ne. 0) return + read(stn(idx1+1:idx2+1), '(i5.5)') id + + ! Select FIPS code for specific countries where we have special rules. + select case(id) + case(752) + fips = 'SW' ! Sweden + case(208) + fips = 'DA' ! Denmark + case(051) ! Post-Soviet countries start here + fips = 'AM' ! Armenia + case(031) + fips = 'AJ' ! Azerbaijan + case(112) + fips = 'BO' ! Belarus + case(233) + fips = 'EN' ! Estonia + case(268) + fips = 'GG' ! Georgia + case(428) + fips = 'LG' ! Latvia + case(440) + fips = 'LH' ! Lithuania + case(398) + fips = 'KZ' ! Kazakhstan + case(417) + fips = 'KG' ! Kyrgyzstan + case(498) + fips = 'MD' ! Moldova + case(643) + fips = 'RS' ! Russia + case(762) + fips = 'TI' ! Tajikistan + case(795) + fips = 'TX' ! Turkmenistan + case(804) + fips = 'UP' ! Ukraine + case(860) + fips = 'UZ' ! Uzbekistan + case(356) ! Post-Soviet countries end here + fips = 'IN' ! India + case(144) + fips = 'CE' ! Sri Lanka + case(032) ! South American countries start here + fips = 'AR' ! Argentina + case(068) + fips = 'BL' ! Bolivia + case(076) + fips = 'BR' ! Brazil + case(152) + fips = 'CI' ! Chile + case(170) + fips = 'CO' ! Columbia + case(218) + fips = 'EC' ! Ecuador + case(238) + fips = 'FK' ! Falkland Islands + case(254) + fips = 'FG' ! French Guiana + case(328) + fips = 'GY' ! Guyana + case(600) + fips = 'PA' ! Paraguay + case(604) + fips = 'PE' ! Peru + case(239) + fips = 'SX' ! South Georgia and the South Sandwich Islands + case(740) + fips = 'NS' ! Suriname + case(858) + fips = 'UY' ! Uruguay + case(862) + fips = 'VE' ! Venezuela + end select + end function get_fips_from_wigos_issuer + +end module USAF_PreobsReaderMod diff --git a/lis/metforcing/usaf/USAF_bratsethMod.F90 b/lis/metforcing/usaf/USAF_bratsethMod.F90 index 883e79acf..46421128a 100644 --- a/lis/metforcing/usaf/USAF_bratsethMod.F90 +++ b/lis/metforcing/usaf/USAF_bratsethMod.F90 @@ -96,8 +96,8 @@ module USAF_bratsethMod type USAF_obsData private integer :: nobs - character*10, allocatable :: net(:) - character*10, allocatable :: platform(:) + character*32, allocatable :: net(:) + character*32, allocatable :: platform(:) real, allocatable :: obs(:) ! Observed variable real, allocatable :: lat(:) ! Latitude of observation (deg N) real, allocatable :: lon(:) ! Longitude of observation (deg E) @@ -145,10 +145,10 @@ module USAF_bratsethMod real, parameter :: MISSING = -9999 ! Quality control flags. Should these be public? - real, parameter :: QC_UNKNOWN = 0 - real, parameter :: QC_GOOD = 1 - real, parameter :: QC_SUSPECT = 2 - real, parameter :: QC_REJECT = 3 + integer, parameter :: QC_UNKNOWN = 0 + integer, parameter :: QC_GOOD = 1 + integer, parameter :: QC_SUSPECT = 2 + integer, parameter :: QC_REJECT = 3 contains @@ -280,8 +280,8 @@ subroutine USAF_assignObsData(this,net,platform,ob,lat,lon,sigmaOSqr, & ! Arguments type(USAF_ObsData),intent(inout) :: this - character(len=10), intent(in) :: net - character(len=10), intent(in) :: platform + character(len=32), intent(in) :: net + character(len=32), intent(in) :: platform real, intent(in) :: ob real, intent(in) :: lat real, intent(in) :: lon @@ -622,7 +622,8 @@ subroutine USAF_addSSMIObsData(this,imax,jmax,ra_tmp,nest) real, allocatable :: xpts(:), ypts(:), rlat(:), rlon(:) real :: sigmaOSqr, ob, xi1, xj1, oErrScaleLength real :: xpnmcaf, ypnmcaf, orient, xmesh, xmeshl - character(len=10) :: net, platform + character(len=32) :: net + character(len=32) :: platform integer :: icount integer :: i,j integer :: count_good_ssmi @@ -658,8 +659,8 @@ subroutine USAF_addSSMIObsData(this,imax,jmax,ra_tmp,nest) continue else write(LIS_logunit,*)'[ERR] Invalid imax dimension for SSM/I!' - write(LIS_logunit,*)'Received ', imax - write(LIS_logunit,*)'Only support 512, 1024, or 1400' + write(LIS_logunit,*)'[ERR] Received ', imax + write(LIS_logunit,*)'[ERR] Only support 512, 1024, or 1400' flush(LIS_logunit) call LIS_endrun() end if @@ -726,7 +727,7 @@ subroutine USAF_addSSMIObsData(this,imax,jmax,ra_tmp,nest) if (imax .eq. 1440) then write(LIS_logunit,*)'[ERR] Lat/lon SSM/I data not supported yet!' - write(LIS_logunit,*)'Modify USAF_addSSMIObsData and recompile!' + write(LIS_logunit,*)'[ERR] Modify USAF_addSSMIObsData and recompile!' flush(LIS_logunit) call LIS_endrun() @@ -908,7 +909,7 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex) call AGRMET_julhr_date10(j3hr, yyyymmddhh) write(LIS_logunit,*) & '[ERR] No NWP background precipitation found for ',yyyymmddhh - write(LIS_logunit,*) ' ABORTING!' + write(LIS_logunit,*) '[ERR] ABORTING!' flush(LIS_logunit) message(:) = '' message(1) = '[ERR] Program: LIS' @@ -918,6 +919,8 @@ subroutine USAF_getBackNWP(nest,back4,pcp_src, use_twelve, j6hr, findex) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_getBackNWP') #endif if(LIS_masterproc) then call LIS_alert( 'LIS.USAF_getBackNWP ', 1, & @@ -1050,9 +1053,9 @@ subroutine USAF_getSSMIObsData(nest,j6hr,use_twelve,precip3,precip6, & write(LIS_logunit,*) ' ' write(LIS_logunit,*) & '[WARN] precip/smiedr: error opening file' - write(LIS_logunit,*)' SSMI data file ', trim(ifil), & + write(LIS_logunit,*)'[WARN] SSMI data file ', trim(ifil), & ' does not exist.' - write(LIS_logunit,*)' SSMI estimates will not be used ',& + write(LIS_logunit,*)'[WARN] SSMI estimates will not be used ',& 'in precip analysis.' write(LIS_logunit,*) ' ' message =' ' @@ -1069,7 +1072,7 @@ subroutine USAF_getSSMIObsData(nest,j6hr,use_twelve,precip3,precip6, & ra(hemi,:,:) = MISSING else - write(LIS_logunit,*) '- READING ',trim(ifil) + write(LIS_logunit,*) '[INFO] READING ',trim(ifil) call LIS_putget( ra(hemi,:,:), 'r', ifil, routine_name, & imax, jmax) end if ! .not. exists @@ -1077,8 +1080,9 @@ subroutine USAF_getSSMIObsData(nest,j6hr,use_twelve,precip3,precip6, & ! Honor option to reset SSMI zero precip values to missing if (.not. use_zeros) then - write(LIS_logunit,*)'- SSMI ZEROS NOT USED' - where ( ra(:,:,:) .eq. 0.0 ) + write(LIS_logunit,*)'[INFO] SSMI ZEROS NOT USED' + where ( .not. ra(:,:,:) > 0.0 .and. & + .not. ra(:,:,:) < 0) ra(:,:,:) = MISSING end where end if @@ -1164,7 +1168,8 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& real, allocatable :: xpts(:), ypts(:), rlat(:), rlon(:) real :: sigmaOSqr, ob, xi1, xj1, oErrScaleLength real :: xpnmcaf, ypnmcaf, orient, xmesh, xmeshl - character(len=10) :: net, platform + character(len=32) :: net + character(len=32) :: platform integer :: count_good_geo_precip, icount integer :: npts integer :: i,j,jj @@ -1223,7 +1228,7 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& else write(LIS_logunit,*) & '[ERR] Invalid dimension for GEO_PRECIP data!' - write(LIS_logunit,*)'Read ', agrmet_struc(nest)%imax + write(LIS_logunit,*)'[ERR] Read ', agrmet_struc(nest)%imax call LIS_endrun() end if end if @@ -1253,8 +1258,8 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& else write(LIS_logunit,*) & '[ERR] Invalid imax dimension for GEO_PRECIP!' - write(LIS_logunit,*)'Received ',imax - write(LIS_logunit,*)'Only support 512, 1024, or 4096!' + write(LIS_logunit,*)'[ERR] Received ',imax + write(LIS_logunit,*)'[ERR] Only support 512, 1024, or 4096!' flush(LIS_logunit) call LIS_endrun() end if @@ -1281,11 +1286,11 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& inquire( file = trim(ifil), exist = exists) if ( .not. exists ) then write(LIS_logunit,*) & - 'USAF_getGeoPrecipObsData: error opening file ', & + '[WARN] USAF_getGeoPrecipObsData: error opening file ', & trim(ifil) - write(LIS_logunit,*) ' file does not exist' + write(LIS_logunit,*) '[WARN] file does not exist' write(LIS_logunit,*) & - ' geo precip estimate will not be performed' + '[WARN] geo precip estimate will not be performed' message = ' ' message(1) = 'program: LIS' message(2) = ' routine: USAF_getGeoPrecipObsData' @@ -1300,7 +1305,7 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& TRACE_EXIT("bratseth_getGeopPrcp") return endif - write(LIS_logunit,*) '- READING ', trim(ifil) + write(LIS_logunit,*) '[INFO] READING ', trim(ifil) allocate(geoprc(imax,jmax)) @@ -1319,10 +1324,10 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& ! Check for anomalous geoprecip files if (is_geo_corrupted(geoprc, imax, jmax, mo, hemi)) then write(LIS_logunit,*) & - 'USAF_getGeoPrecipObsData: data corrupted - ', & + '[WARN] USAF_getGeoPrecipObsData: data corrupted - ', & trim(ifil) write(LIS_logunit,*) & - ' geo precip estimate will not be performed' + '[WARN] geo precip estimate will not be performed' message = ' ' message(1) = 'program: LIS' message(2) = ' routine: USAF_getGeoPrecipObsData' @@ -1460,7 +1465,7 @@ subroutine USAF_getGeoPrecipObsData(nest,j6hr,use_twelve,precip3,precip6,& write(LIS_logunit,*)& '[ERR] Lat/lon GEO_PRECIP data not supported yet!' write(LIS_logunit,*)& - 'Modify USAF_addGeoPrecipObsData and recompile!' + '[ERR] Modify USAF_addGeoPrecipObsData and recompile!' flush(LIS_logunit) call LIS_endrun() @@ -1499,7 +1504,7 @@ subroutine USAF_interpBackToTypeObsData(this,nest,imax,jmax,back,type) integer,intent(in) :: imax integer,intent(in) :: jmax real, intent(in) :: back(imax,jmax) - character(len=10),intent(in) :: type + character(len=32),intent(in) :: type ! Local variables integer :: nobs @@ -1885,6 +1890,8 @@ subroutine calc_invDataDensities(this,sigmaBSqr,nest,max_dist, & #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_invDataDensities') t1 = MPI_Wtime() #endif @@ -1949,7 +1956,8 @@ subroutine calc_invDataDensities(this,sigmaBSqr,nest,max_dist, & else if (trim(this%net(iob)) .eq. trim(this%net(job))) then ! Satellite observations have correlated errors. if (.not. isUncorrObType(this%net(job))) then - if (this%oErrScaleLength(job) .eq. 0) then + if (.not. this%oErrScaleLength(job) > 0 .and. & + .not. this%oErrScaleLength(job) < 0) then write(LIS_logunit,*) & '[ERR]: job, network, oErrScaleLength: ', & job, trim(this%net(job)), & @@ -1991,8 +1999,12 @@ subroutine calc_invDataDensities(this,sigmaBSqr,nest,max_dist, & allocate(invDataDensities(nobs)) invDataDensities(:) = 0 call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_invDataDensities') call MPI_ALLREDUCE(dataDensities_pet,invDataDensities,nobs,MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_ALLREDUCE call in calc_invDataDensities') #endif ! Clean up @@ -2010,6 +2022,8 @@ subroutine calc_invDataDensities(this,sigmaBSqr,nest,max_dist, & #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_invDataDensities') t2 = MPI_Wtime() if (verbose) then write(LIS_logunit,*) & @@ -2082,7 +2096,7 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& integer :: imaxabsdiff real :: maxabsdiff, y_prev, y_new, normdev integer :: icount, num_high_dev - integer :: c,r,i,j,iob,job + integer :: c,r,i,j,iob,job,ii integer :: ierr double precision :: t0,t1, t2 logical :: verbose @@ -2115,6 +2129,12 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& call LIS_endrun() end if + ! EMK TEST + !do ii = 1, nobs + ! write(LIS_logunit,*)'EMK: ii, invDataDensity: ', ii, & + ! invDataDensities(ii) + !end do + ! Here we create a 2d hash table storing the index values of each ob ! in linked lists for each LIS grid box. This can help us screen ! out obviously unnecessary ob comparisons later. @@ -2150,6 +2170,8 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_obsAnalysis') t0 = MPI_Wtime() #endif @@ -2157,6 +2179,8 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_obsAnalysis') t1 = MPI_Wtime() #endif pnew_est_pet(:) = 0 @@ -2259,11 +2283,19 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& pnew_est(:) = 0 pnew_ana(:) = 0 call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_obsAnalysis') call MPI_ALLREDUCE(pnew_est_pet, pnew_est, nobs, MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_ALLREDUCE call in calc_obsAnalysis') call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_obsAnalysis') call MPI_ALLREDUCE(pnew_ana_pet, pnew_ana, nobs, MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_ALLREDUCE call in calc_obsAnalysis') #endif ! Finish analysis and observation estimates for this iteration @@ -2287,8 +2319,12 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& ! Share sumObsEstimates across all processors. sumObsEstimates(:) = 0 call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_obsAnalysis') call MPI_ALLREDUCE(sumObsEstimates_pet, sumObsEstimates, nobs, & MPI_REAL, MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_ALLREDUCE call in calc_obsAnalysis') #else do j = 1, nobs sumObsEstimates(j) = sum(sumObsEstimates_pet) @@ -2357,6 +2393,20 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& end if end if + ! !EMK TEST + ! do ii = 1, nobs + ! write(LIS_logunit,*) & + ! '[INFO] ii,net,platform,obs, back, ana, est, dataDensity: ', & + ! ii, & + ! trim(this%net(ii)), ' ',& + ! trim(this%platform(ii)), & + ! ' ',this%obs(ii),& + ! ' ',this%back(ii),& + ! ' ',pnew_ana(ii),& + ! ' ',pnew_est(ii),& + ! ' ',1./invDataDensities(ii) + ! end do + if (done) exit ! No more iterations! if (verbose) then @@ -2387,6 +2437,8 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_obsAnalysis') t2 = MPI_Wtime() if (verbose) then write(LIS_logunit,*) & @@ -2399,24 +2451,24 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& ! production, we will allow the program to continue instead of ! aborting. !if (npasses .eq. 100) then - if (npasses .eq. 5) then ! For Ops - - write(LIS_logunit,*) & - '[WARN] Bratseth failed to converge after ',npasses, & - ' iterations!' - write(LIS_logunit,*) & - '[WARN] Will stop iterating' - flush(LIS_logunit) - exit - end if + if (npasses .eq. 5) then ! For Ops + + write(LIS_logunit,*) & + '[WARN] Bratseth failed to converge after ',npasses, & + ' iterations!' + write(LIS_logunit,*) & + '[WARN] Will stop iterating' + flush(LIS_logunit) + exit + end if end do ! Iterate until convergence if (verbose) then if (done) then - write(LIS_logunit,*) & + write(LIS_logunit,*) & '[INFO] Bratseth analysis converged after ',npasses, & ' iterations' - endif + endif write(LIS_logunit,*) & '[INFO] Mean absolute difference against obs: ana: ',mad_ana, & ' est: ',mad_est @@ -2424,6 +2476,8 @@ subroutine calc_obsAnalysis(this,sigmaBSqr,nobs,invDataDensities,nest,& #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_obsAnalysis') t2 = MPI_Wtime() if (verbose) then write(LIS_logunit,*) & @@ -2550,12 +2604,14 @@ subroutine calc_gridAnalysis(this,nest,sigmaBSqr,nobs,invDataDensities,& if (nobs .ne. this%nobs) then write(LIS_logunit,*)'[ERR] Array size mismatch in calc_gridAnalysis!' - write(LIS_logunit,*)'nobs, this%nobs = ',nobs, this%nobs + write(LIS_logunit,*)'[ERR] nobs, this%nobs = ',nobs, this%nobs call LIS_endrun() end if #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_gridAnalysis') t1 = MPI_Wtime() #endif @@ -2641,9 +2697,13 @@ subroutine calc_gridAnalysis(this,nest,sigmaBSqr,nobs,invDataDensities,& allocate(mrgp_1d(LIS_rc%gnc(nest)*LIS_rc%gnr(nest))) mrgp_1d(:) = 0 call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_gridAnalysis') call MPI_ALLREDUCE(mrgp_1d_pet,mrgp_1d, & LIS_rc%gnc(nest)*LIS_rc%gnr(nest),MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_ALLREDUCE call in calc_gridAnalysis') deallocate(mrgp_1d_pet) #endif @@ -2670,6 +2730,8 @@ subroutine calc_gridAnalysis(this,nest,sigmaBSqr,nobs,invDataDensities,& #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in calc_gridAnalysis') t2 = MPI_Wtime() write(LIS_logunit,*)'[INFO] Elapsed time for grid analysis is ', & t2 - t1,' seconds' @@ -3018,6 +3080,7 @@ subroutine check_grib_file(gribfile,yr1,mo1,da1,hr1,found, & #if (defined USE_GRIBAPI) use grib_api #endif + use LIS_coreMod, only: LIS_masterproc use LIS_logMod, only : LIS_logunit, LIS_abort, LIS_alert, & LIS_verify, LIS_endrun use LIS_mpiMod @@ -3182,6 +3245,8 @@ subroutine check_grib_file(gribfile,yr1,mo1,da1,hr1,found, & message(3) = ' LIS was not compiled with GRIBAPI or ECCODES support!' #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in check_grib_file') #endif if(LIS_masterproc) then call LIS_alert( 'LIS.check_grib_file', 1, & @@ -3219,7 +3284,7 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) ! Arguments type(USAF_ObsData), intent(inout) :: this integer,intent(in) :: nest - character(len=10), intent(in) :: new_name + character(len=32), intent(in) :: new_name character(len=*), optional :: network logical,optional,intent(in) :: silent_rejects @@ -3229,7 +3294,8 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) real :: dlat, dlon integer :: c,r,j real :: ctrlat, ctrlon - character(len=10) :: net_new, platform_new + character(len=32) :: net_new + character(len=32) :: platform_new integer, allocatable :: actions(:), actions_pet(:) real, allocatable :: superobs_pet(:),superlat_pet(:),superlon_pet(:) real, allocatable :: means(:) @@ -3274,6 +3340,8 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_superstatQC') t1 = MPI_Wtime() #endif @@ -3435,6 +3503,8 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) #if (defined SPMD) call MPI_Allreduce(superobs_pet,means,glbcr,MPI_REAL, MPI_SUM, & LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') #else do j = 1, glbcr means(j) = sum(superobs_pet) @@ -3445,6 +3515,8 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) #if (defined SPMD) call MPI_Allreduce(superob_count_pet,superob_count,glbcr,& MPI_INTEGER, MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') #else do j = 1, glbcr superob_count(j) = sum(superob_count_pet) @@ -3469,35 +3541,49 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_superstatQC') superob_count(:) = 0 call MPI_Allreduce(superob_count_pet, superob_count, glbcr,& MPI_INTEGER, MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') superob_count_pet(:) = 0 superobs(:) = 0 call MPI_Allreduce(superobs_pet, superobs, glbcr, MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') superobs_pet(:) = 0 superlat(:) = 0 call MPI_Allreduce(superlat_pet, superlat, glbcr, MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') superlat_pet(:) = 0 superlon(:) = 0 call MPI_Allreduce(superlon_pet, superlon, glbcr, MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') superlon_pet(:) = 0 superSigmaOSqr(:) = 0 call MPI_Allreduce(superSigmaOSqr_pet, superSigmaOSqr, glbcr, & MPI_REAL, MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') superSigmaOSqr_pet(:) = 0 superOErrScaleLength(:) = 0 call MPI_Allreduce(superOErrScaleLength_pet, & superOErrScaleLength, glbcr, & MPI_REAL, MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') superOErrScaleLength_pet(:) = 0 #endif @@ -3530,9 +3616,13 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_superstatQC') actions(:) = 0 call MPI_Allreduce(actions_pet, actions, nobs,& MPI_INTEGER, MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_superstatQC') actions_pet(:) = 0 #endif end if ! ipass .eq. 3 @@ -3609,6 +3699,8 @@ subroutine USAF_superstatQC(this,nest,new_name,network,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_superstatQC') t2 = MPI_Wtime() write(LIS_logunit,*) & '[INFO] Elapsed time in superstatQC is ',t2 - t1,' seconds' @@ -3651,7 +3743,8 @@ subroutine USAF_dupQC(this) integer :: count_dups integer :: total_reject_count, total_merge_count, total_create_count real :: mean,back,newlat,newlon,sigmaOSqr,oErrScaleLength - character(len=10) :: net,platform + character(len=32) :: net + character(len=32) :: platform real :: diff integer :: r,c,i integer :: nobs @@ -3669,6 +3762,8 @@ subroutine USAF_dupQC(this) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_dupQC') t1 = MPI_Wtime() #endif @@ -3735,9 +3830,11 @@ subroutine USAF_dupQC(this) ptr => head do i = 1, count_dups diff = this%lat(ptr%ob_index) - this%lat(r) - if (diff .ne. 0) location_issue = .true. + if (.not. diff > 0 .and. & + .not. diff < 0) location_issue = .true. diff = this%lon(ptr%ob_index) - this%lon(r) - if (diff .ne. 0) location_issue = .true. + if (.not. diff > 0 .and. & + .not. diff < 0) location_issue = .true. if (location_issue) exit ! Get out of loop ptr => ptr%next end do ! i @@ -3773,7 +3870,8 @@ subroutine USAF_dupQC(this) reject_all = .false. do i = 1, count_dups diff = this%obs(ptr%ob_index) - this%obs(r) - if (diff .eq. 0) then + if (.not. diff > 0 .and. & + .not. diff < 0) then this%qc(ptr%ob_index) = QC_REJECT total_reject_count = total_reject_count + 1 write(LIS_logunit,*) & @@ -3834,7 +3932,8 @@ subroutine USAF_dupQC(this) if (count_dups .eq. 1 .and. .not. location_issue) then ptr => head diff = this%obs(ptr%ob_index) - this%obs(r) - if (diff .eq. 0) then + if (.not. diff < 0 .and. & + .not. diff > 0) then this%qc(ptr%ob_index) = QC_REJECT total_reject_count = total_reject_count + 1 write(LIS_logunit,*) & @@ -3957,6 +4056,8 @@ subroutine USAF_dupQC(this) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_dupQC') t2 = MPI_Wtime() write(LIS_logunit,*) & '[INFO] Elapsed time in dupQC is ',t2 - t1,' seconds' @@ -4013,7 +4114,7 @@ subroutine reset_negative_values(nest,mrgp) if (ifix > 0) then write(LIS_logunit,6000) ifix 6000 format (/, 1x, 55('-'), & - /, 3x, 'routine reset_negative_values:',& + /, 3x, '[INFO] routine reset_negative_values:',& /, 5x, '# of pts to which negative values were set to zero = ', & i6, /, 1x, 55('-')) end if @@ -4063,9 +4164,14 @@ subroutine USAF_backQC(this,sigmaBSqr,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_backQC') t1 = MPI_Wtime() #endif + !write(LIS_logunit,*)'EMK: nobs, sigmaBSqr = ', nobs, sigmaBSqr + !flush(LIS_Logunit) + reject_count = 0 do r = 1,nobs @@ -4101,6 +4207,8 @@ subroutine USAF_backQC(this,sigmaBSqr,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_backQC') t2 = MPI_Wtime() write(LIS_logunit,*) & '[INFO] Elapsed time in backQC is ',t2 - t1,' seconds' @@ -4112,7 +4220,7 @@ end subroutine USAF_backQC ! Checks if observation network is recognized as a gauge. logical function is_gauge(net) implicit none - character(len=10), intent(in) :: net + character(len=32), intent(in) :: net logical :: answer answer = .false. if (trim(net) .eq. "AMIL") answer = .true. @@ -4135,10 +4243,8 @@ end function is_gauge ! surface stations, all observations should have uncorrelated errors. logical function is_stn(net) implicit none - character(len=10), intent(in) :: net - character(len=10) :: net_local + character(len=32), intent(in) :: net logical :: answer - net_local = net answer = .true. is_stn = answer end function is_stn @@ -4147,7 +4253,7 @@ end function is_stn ! Checks if observation "network" is recognized as SSMI retrievals. logical function is_ssmi(net) implicit none - character(len=10), intent(in) :: net + character(len=32), intent(in) :: net logical :: answer answer = .false. if (trim(net) .eq. "SSMI") answer = .true. @@ -4158,7 +4264,7 @@ end function is_ssmi ! Checks if observation "network" is recognized as GEOPRECIP retrievals. logical function is_geoprecip(net) implicit none - character(len=10), intent(in) :: net + character(len=32), intent(in) :: net logical :: answer answer = .false. if (trim(net) .eq. "GEOPRECIP") answer = .true. @@ -4169,7 +4275,7 @@ end function is_geoprecip ! Checks if observation "network" is recognized as CMORPH estimate. logical function is_cmorph(net) implicit none - character(len=10), intent(in) :: net + character(len=32), intent(in) :: net logical :: answer answer = .false. if (trim(net) .eq. "CMORPH") answer = .true. @@ -4180,7 +4286,7 @@ end function is_cmorph ! Checks if observation "network" is recognized as IMERG retrievals. logical function is_imerg(net) implicit none - character(len=10), intent(in) :: net + character(len=32), intent(in) :: net logical :: answer answer = .false. if (trim(net) .eq. "IMERG") answer = .true. @@ -4843,6 +4949,8 @@ subroutine find_LIS_cols_rows(this,nest,cols,rows) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in find_LIS_cols_rows') #endif do j = 1, nobs @@ -4923,8 +5031,12 @@ subroutine find_LIS_cols_rows(this,nest,cols,rows) cols = 0 #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in find_LIS_cols_rows') call MPI_ALLREDUCE(cols_pet,cols,nobs,MPI_INTEGER, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in find_LIS_cols_rows') #else cols(:) = cols_pet(:) #endif @@ -4935,8 +5047,12 @@ subroutine find_LIS_cols_rows(this,nest,cols,rows) rows = 0 #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in find_LIS_cols_rows') call MPI_ALLREDUCE(rows_pet,rows,nobs,MPI_INTEGER, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in find_LIS_cols_rows') #else rows(:) = rows_pet(:) #endif @@ -4986,6 +5102,8 @@ subroutine USAF_waterQC(this,nest,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_waterQC') t1 = MPI_Wtime() #endif @@ -5023,6 +5141,8 @@ subroutine USAF_waterQC(this,nest,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_waterQC') t2 = MPI_Wtime() write(LIS_logunit,*) & '[INFO] Elapsed time in waterQC is ',t2-t1,' seconds' @@ -5083,6 +5203,8 @@ subroutine USAF_snowQC(this,nest,hourindex,threshold,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_snowQC') t1 = MPI_Wtime() #endif @@ -5126,9 +5248,13 @@ subroutine USAF_snowQC(this,nest,hourindex,threshold,silent_rejects) allocate(sfctmp_1d(LIS_rc%gnc(nest)*LIS_rc%gnr(nest))) sfctmp_1d(:) = 0 call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_snowQC') call MPI_ALLREDUCE(sfctmp_1d_pet,sfctmp_1d, & LIS_rc%gnc(nest)*LIS_rc%gnr(nest),MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_snowQC') deallocate(sfctmp_1d_pet) #endif allocate(sfctmp(LIS_rc%gnc(nest), LIS_rc%gnr(nest))) @@ -5205,6 +5331,8 @@ subroutine USAF_snowQC(this,nest,hourindex,threshold,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_snowQC') t2 = MPI_Wtime() write(LIS_logunit,*) & '[INFO] Elapsed time in snowQC is ',t2-t1,' seconds' @@ -5266,6 +5394,8 @@ subroutine USAF_snowDepthQC(this,nest,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_snowDepthQC') t1 = MPI_Wtime() #endif @@ -5317,9 +5447,13 @@ subroutine USAF_snowDepthQC(this,nest,silent_rejects) allocate(snowdepth_1d(LIS_rc%gnc(nest)*LIS_rc%gnr(nest))) snowdepth_1d(:) = 0 call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_snowDepthQC') call MPI_ALLREDUCE(snowdepth_1d_pet,snowdepth_1d, & LIS_rc%gnc(nest)*LIS_rc%gnr(nest),MPI_REAL, & MPI_SUM, LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Allreduce call in USAF_snowDepthQC') deallocate(snowdepth_1d_pet) #endif allocate(snowdepth(LIS_rc%gnc(nest), LIS_rc%gnr(nest))) @@ -5387,6 +5521,8 @@ subroutine USAF_snowDepthQC(this,nest,silent_rejects) #if (defined SPMD) call MPI_Barrier(LIS_MPI_COMM, ierr) + call handle_mpi_error(ierr, & + 'MPI_Barrier call in USAF_snowDepthQC') t2 = MPI_Wtime() write(LIS_logunit,*) & '[INFO] Elapsed time in snowDepthQC is ',t2-t1,' seconds' @@ -5477,7 +5613,7 @@ subroutine USAF_analyzeScreen(screenObs,nest,back,sigmaBSqr, & real, allocatable :: invDataDensities(:) real, allocatable :: sumObsEstimates(:) integer :: npasses - character(len=10) :: new_name, type + character(len=32) :: new_name,type real :: convergeThresh integer :: j @@ -5640,7 +5776,7 @@ subroutine USAF_getCMORPHObsData(nest,j6hr,use_twelve, & integer, parameter :: YD = 1649 integer, parameter :: NCMOR = XD*YD real :: precip(XD,YD) - character(len=10) :: net, platform + character(len=32) :: net, platform real :: sigmaOSqr, oErrScaleLength integer :: count_good_obs character(len=120) :: fname @@ -6338,9 +6474,9 @@ subroutine USAF_setBratsethPrecipStats(src,nest) write(LIS_logunit,*) & '[ERR] Unknown source of background precipitation!' write(LIS_logunit,*) & - 'Source is ',trim(src) + '[ERR] Source is ',trim(src) write(LIS_logunit, *) & - 'ABORTING....' + '[ERR] ABORTING....' flush(LIS_logunit) call LIS_endrun() end if @@ -6371,6 +6507,11 @@ subroutine USAF_setBratsethScreenStats(src,n) agrmet_struc(n)%galwem_t2m_back_err_scale_length agrmet_struc(n)%bratseth_t2m_back_sigma_b_sqr = & agrmet_struc(n)%galwem_t2m_back_sigma_b_sqr + + !write(LIS_logunit,*)'EMK: n, GALWEM T2 err = ', n, & + ! agrmet_struc(n)%bratseth_t2m_back_sigma_b_sqr + !flush(LIS_logunit) + agrmet_struc(n)%bratseth_t2m_stn_sigma_o_sqr = & agrmet_struc(n)%galwem_t2m_stn_sigma_o_sqr agrmet_struc(n)%bratseth_t2m_max_dist = & @@ -6399,6 +6540,11 @@ subroutine USAF_setBratsethScreenStats(src,n) agrmet_struc(n)%gfs_t2m_back_err_scale_length agrmet_struc(n)%bratseth_t2m_back_sigma_b_sqr = & agrmet_struc(n)%gfs_t2m_back_sigma_b_sqr + + !write(LIS_logunit,*)'EMK: n, GFS T2 err = ', n, & + ! agrmet_struc(n)%bratseth_t2m_back_sigma_b_sqr + !flush(LIS_Logunit) + agrmet_struc(n)%bratseth_t2m_stn_sigma_o_sqr = & agrmet_struc(n)%gfs_t2m_stn_sigma_o_sqr agrmet_struc(n)%bratseth_t2m_max_dist = & @@ -6426,9 +6572,9 @@ subroutine USAF_setBratsethScreenStats(src,n) write(LIS_logunit,*) & '[ERR] Unknown source of background data!' write(LIS_logunit,*) & - 'Source is ',trim(src) + '[ERR] Source is ',trim(src) write(LIS_logunit, *) & - 'ABORTING....' + '[ERR] ABORTING....' flush(LIS_logunit) call LIS_endrun() end if @@ -6502,7 +6648,6 @@ subroutine USAF_pcpBackBiasRatio_s2s(n, yyyymmddhh) 'n90_get_var failed for PPT_ratio in LIS param file') call LIS_verify(nf90_close(ncid),& 'nf90_close failed in USAF_pcpBackBiasRatio_s2s') - #endif end subroutine USAF_pcpBackBiasRatio_s2s @@ -6609,4 +6754,20 @@ subroutine USAF_pcpBackBiasRatio_nrt(n, yyyymmddhh) end subroutine USAF_pcpBackBiasRatio_nrt + subroutine handle_mpi_error(errorcode, msg) + use LIS_logMod, only: LIS_logunit, LIS_endrun + use LIS_mpiMod + implicit none + integer, intent(in) :: errorcode + character(*), intent(in) :: msg + character*(MPI_MAX_ERROR_STRING) :: buf + integer :: resultlen, ierr + if (errorcode .ne. MPI_SUCCESS) then + write(LIS_logunit,*)'[ERR] ', trim(msg) + call MPI_error_string(errorcode, buf, resultlen, ierr) + write(LIS_logunit,*)'[ERR] MPI error: ', trim(buf) + flush(LIS_logunit) + call LIS_endrun() + end if + end subroutine handle_mpi_error end module USAF_bratsethMod diff --git a/lis/metforcing/usaf/USAF_getpcpobs.F90 b/lis/metforcing/usaf/USAF_getpcpobs.F90 new file mode 100644 index 000000000..3a9679aa9 --- /dev/null +++ b/lis/metforcing/usaf/USAF_getpcpobs.F90 @@ -0,0 +1,107 @@ +! ROUTINE: USAF_getpcpcobs +! +! REVISION HISTORY: +! 17 Oct 2023 Initial version. Eric Kemp/SSAI/NASA. +! +! DESCRIPTION: +! New routine to retrieve precip observations and process with new +! USAF_Gages library. Supports legacy hemispheric and new global +! preobs file formats, and new presav2 file format. Borrows heavily +! from older AGRMET_getpcpobs subroutine. +!-------------------------------------------------------------------------- + +subroutine USAF_getpcpobs(n, j6hr, month, use_twelve, pcp_src, & + use_expanded_station_ids, alert_number, precip6, precip12) + + ! Imports + use AGRMET_forcingMod, only: agrmet_struc + use ESMF + use LIS_coreMod, only: LIS_rc + use LIS_logMod, only: LIS_logunit + use LIS_timeMgrMod, only: LIS_tick, LIS_julhr_date + use USAF_bratsethMod, only: USAF_ObsData, USAF_setbratsethprecipstats + use USAF_GagesMod, only : USAF_Gages_t + use USAF_PreobsReaderMod, only: USAF_read_preobs + + ! Defaults + implicit none + + ! Arguments + integer, intent(in) :: n + integer, intent(in) :: j6hr + integer, intent(in) :: month + logical, intent(in) :: use_twelve + character*6, intent(in) :: pcp_src(4) + integer, intent(in) :: use_expanded_station_ids + integer, intent(inout):: alert_number + type(USAF_ObsData), intent(inout) :: precip6 + type(USAF_ObsData), intent(inout) :: precip12 + + ! Locals + integer :: j3hr + integer :: k + integer :: yr, mo, da, hr + character*255 :: preobsdir + character*8 :: yyyymmdd + character*10 :: yyyymmddhh + character*255 :: presav_filename + logical :: file_exists + type(USAF_Gages_t) :: obscur + + if (use_twelve) then + k = 3 + else + k = 1 + end if + + do j3hr = j6hr+3, j6hr+6, 3 + + ! Set Bratseth error statistics based on source of background field. + call USAF_setBratsethPrecipStats(pcp_src(k), n) + + call LIS_julhr_date(j3hr, yr, mo, da, hr) + if (agrmet_struc(n)%use_timestamp == 1) then + write(unit=yyyymmdd, fmt='(i4.4, i2.2, i2.2)') & + yr, mo, da + preobsdir = trim(agrmet_struc(n)%agrmetdir) // '/' & + // yyyymmdd // '/' & + // trim(agrmet_struc(n)%cdmsdir) // '/' + else + preobsdir = trim(agrmet_struc(n)%agrmetdir) // '/' & + // trim(agrmet_struc(n)%cdmsdir) // '/' + end if + + ! Read appropriate preobs file(s), intercompare with older presav2 + ! files, and create new presav2 file for current date/time. + call USAF_read_preobs(preobsdir, & + trim(agrmet_struc(n)%analysisdir), & + agrmet_struc(n)%use_timestamp, yr, mo, da, hr, & + use_expanded_station_ids, alert_number) + + ! If this is a synoptic time, read the presav2 file back in and + ! populate the appropriate USAF_ObsData object. + if ( mod(j3hr, 6) == 0 ) then + write(presav_filename,'(A,A,i4.4,i2.2,i2.2,i2.2)') & + trim(agrmet_struc(n)%analysisdir), '/presav2.03hr.', & + yr, mo, da, hr + inquire(file=presav_filename, exist=file_exists) + if (file_exists) then + write(yyyymmddhh,'(i4.4,i2.2,i2.2,i2.2)') & + yr, mo, da, hr + call obscur%read_data(presav_filename, yyyymmddhh, & + alert_number) + if (use_twelve) then + call obscur%copy_to_usaf_obsdata(12, & + agrmet_struc(n)%bratseth_precip_gauge_sigma_o_sqr, & + precip12) + else + call obscur%copy_to_usaf_obsdata(6, & + agrmet_struc(n)%bratseth_precip_gauge_sigma_o_sqr, & + precip6) + end if + end if + end if + k = k + 1 + end do + +end subroutine USAF_getpcpobs diff --git a/lis/metforcing/usaf/readagrmetpcpforcing.F90 b/lis/metforcing/usaf/readagrmetpcpforcing.F90 index b875ff14e..c47c1f020 100644 --- a/lis/metforcing/usaf/readagrmetpcpforcing.F90 +++ b/lis/metforcing/usaf/readagrmetpcpforcing.F90 @@ -267,7 +267,7 @@ subroutine readagrmetpcpforcing(n,findex, order) ! precip12Imerg_tmp integer :: hourindex integer :: k1,k2,k3,k4 - character(len=10) :: type + character(len=32) :: type character(len=10) :: yyyymmddhh character(len=50) :: pathOBA logical :: found_inq @@ -280,9 +280,10 @@ subroutine readagrmetpcpforcing(n,findex, order) integer*2 :: imerg_plp_thresh real :: imerg_sigmaOSqr real :: imerg_oErrScaleLength - character(len=10) :: imerg_net, imerg_platform + character(len=32) :: imerg_net, imerg_platform real :: sigmaBSqr - character(len=10) :: new_name + character(len=32) :: new_name + integer :: use_expanded_station_ids data alert_number / 0 / data srcwts /100.0,50.0,4.0,4.0,1.0,1.0,60.0,1.0/ data addrad /0, -1, -2, -4, -5, -5, 0, -2/ @@ -399,14 +400,18 @@ subroutine readagrmetpcpforcing(n,findex, order) write(LIS_logunit,*) & '[INFO] Fetching 6-hr gage data' call USAF_createObsData(precip_6hr_gage_tmp,n,maxobs=15000) - call AGRMET_getpcpobs(n, julbeg, LIS_rc%mo, prcpwe, & - use_twelve, p6, p12, alert_number, precip_6hr_gage_tmp, & - precip_12hr_gage_tmp, & - pcp_src) + !call AGRMET_getpcpobs(n, julbeg, LIS_rc%mo, prcpwe, & + ! use_twelve, p6, p12, alert_number, precip_6hr_gage_tmp, & + ! precip_12hr_gage_tmp, & + ! pcp_src) + use_expanded_station_ids = agrmet_struc(n)%pcpobsfmt - 1 + call USAF_getpcpobs(n, julbeg, LIS_rc%mo, use_twelve, & + pcp_src, use_expanded_station_ids, alert_number, & + precip_6hr_gage_tmp, precip_12hr_gage_tmp) ! Reject data over water write(LIS_logunit,*) & - '[INFO] Running waterQC on 0-6 hr gauge observations' + '[INFO] Running waterQC on 6 hr gauge observations' call USAF_waterQC(precip_6hr_gage_tmp,n) nobs_good = USAF_countGoodObs(precip_6hr_gage_tmp) nobs_good_extra = nint(nobs_good*1.10) @@ -912,20 +917,28 @@ subroutine readagrmetpcpforcing(n,findex, order) call USAF_createObsData(precip_12hr_gage_tmp,n,maxobs=15000) ! EMK...Call this twice to ensure we get obs for first two time ! levels. - call AGRMET_getpcpobs(n, julbeg, LIS_rc%mo, prcpwe, & - .false., p6, p12, alert_number,precip_6hr_gage_tmp, & - precip_12hr_gage_tmp, & - pcp_src) - call AGRMET_getpcpobs(n, julbeg+6, LIS_rc%mo, prcpwe, & - .true., p6, p12, alert_number,precip_6hr_gage_tmp, & - precip_12hr_gage_tmp, & - pcp_src) + !call AGRMET_getpcpobs(n, julbeg, LIS_rc%mo, prcpwe, & + ! .false., p6, p12, alert_number,precip_6hr_gage_tmp, & + ! precip_12hr_gage_tmp, & + ! pcp_src) + !call AGRMET_getpcpobs(n, julbeg+6, LIS_rc%mo, prcpwe, & + ! .true., p6, p12, alert_number,precip_6hr_gage_tmp, & + ! precip_12hr_gage_tmp, & + ! pcp_src) + use_expanded_station_ids = agrmet_struc(n)%pcpobsfmt - 1 + call USAF_getpcpobs(n, julbeg, LIS_rc%mo, .false., pcp_src, & + use_expanded_station_ids, alert_number, & + precip_6hr_gage_tmp, precip_12hr_gage_tmp) + call USAF_getpcpobs(n, julbeg+6, LIS_rc%mo, .true., pcp_src, & + use_expanded_station_ids, alert_number, & + precip_6hr_gage_tmp, precip_12hr_gage_tmp) + ! Not needed at this point since we have the 12hr accum call USAF_destroyObsData(precip_6hr_gage_tmp) ! Reject and filter out gage reports over water write(LIS_logunit,*) & - '[INFO] Running waterQC on 6-12 hr gauge observations' + '[INFO] Running waterQC on 12 hr gauge observations' call USAF_waterQC(precip_12hr_gage_tmp,n) nobs_good = USAF_countGoodObs(precip_12hr_gage_tmp) nobs_good_extra = nint(nobs_good*1.10) diff --git a/lis/metforcing/usaf/readagrmetpcpforcinganalysis.F90 b/lis/metforcing/usaf/readagrmetpcpforcinganalysis.F90 index ac7d4352e..1b7c050c8 100644 --- a/lis/metforcing/usaf/readagrmetpcpforcinganalysis.F90 +++ b/lis/metforcing/usaf/readagrmetpcpforcinganalysis.F90 @@ -127,13 +127,13 @@ subroutine readagrmetpcpforcinganalysis(n,findex, order) //'.03hr.'//date10_03 endif - write(LIS_logunit,*)' - READING precip ',ifil + write(LIS_logunit,*)'[INFO] READING precip ',ifil inquire(file=ifil,exist=exists) if(exists) then call LIS_putget(gi(hemi,:,:), 'r', ifil, & routine_name, agrmet_struc(n)%imax, agrmet_struc(n)%jmax ) else - write(LIS_logunit,*) 'premrg file does not exist' + write(LIS_logunit,*) '[ERR] premrg file does not exist' write(LIS_logunit,*) ifil call LIS_endrun() endif diff --git a/lis/metforcing/usaf/readcrd_agrmet.F90 b/lis/metforcing/usaf/readcrd_agrmet.F90 index 47635b67f..4b307f4b2 100644 --- a/lis/metforcing/usaf/readcrd_agrmet.F90 +++ b/lis/metforcing/usaf/readcrd_agrmet.F90 @@ -67,8 +67,8 @@ subroutine readcrd_agrmet() call ESMF_ConfigFindLabel(LIS_config,"AGRMET forcing directory:",rc=rc) do n=1,LIS_rc%nnest call ESMF_ConfigGetAttribute(LIS_config,agrmet_struc(n)%agrmetdir,rc=rc) - write(LIS_logunit,*)'Using AGRMET forcing' - write(LIS_logunit,*) 'AGRMET forcing directory :',agrmet_struc(n)%agrmetdir + write(LIS_logunit,*)'[INFO] Using AGRMET forcing' + write(LIS_logunit,*) '[INFO] AGRMET forcing directory: ', trim(agrmet_struc(n)%agrmetdir) enddo call ESMF_ConfigFindLabel(LIS_config,"AGRMET first guess source:",rc=rc) @@ -195,6 +195,65 @@ subroutine readcrd_agrmet() do n=1,LIS_rc%nnest call ESMF_ConfigGetAttribute(LIS_config,agrmet_struc(n)%pcpobswch,rc=rc) enddo + + ! EMK...Precip observation file formats + call ESMF_ConfigFindLabel(LIS_config,"AGRMET precip obs file format:", & + rc=rc) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config, & + agrmet_struc(n)%pcpobsfmt, rc=rc) + if (agrmet_struc(n)%pcpobsfmt .ne. 1 .and. & + agrmet_struc(n)%pcpobsfmt .ne. 2) then + write(LIS_logunit,*) & + "[ERR] Bad 'AGRMET precip obs file format:' option" + write(LIS_logunit,*) & + '[ERR] Expected 1 or 2, found ', agrmet_struc(n)%pcpobsfmt + write(LIS_logunit,*) '[ERR] Aborting...' + + flush(LIS_logunit) + message(1) = & + '[ERR] Illegal value for AGRMET precip obs file format' +#if (defined SPMD) + call MPI_Barrier(LIS_MPI_COMM, ierr) +#endif + if (LIS_masterproc) then + call LIS_abort(message) + else + call sleep(10) + call LIS_endrun() + end if + end if + enddo + + ! EMK...Sfc observation file formats + call ESMF_ConfigFindLabel(LIS_config,"AGRMET sfc obs file format:", & + rc=rc) + do n=1,LIS_rc%nnest + call ESMF_ConfigGetAttribute(LIS_config, & + agrmet_struc(n)%sfcobsfmt, rc=rc) + if (agrmet_struc(n)%sfcobsfmt .ne. 1 .and. & + agrmet_struc(n)%sfcobsfmt .ne. 2) then + write(LIS_logunit,*) & + "[ERR] Bad 'AGRMET sfc obs file format:' option" + write(LIS_logunit,*) & + '[ERR] Expected 1 or 2, found ', agrmet_struc(n)%sfcobsfmt + write(LIS_logunit,*) '[ERR] Aborting...' + + flush(LIS_logunit) + message(1) = & + '[ERR] Illegal value for AGRMET sfc obs file format' +#if (defined SPMD) + call MPI_Barrier(LIS_MPI_COMM, ierr) +#endif + if (LIS_masterproc) then + call LIS_abort(message) + else + call sleep(10) + call LIS_endrun() + end if + end if + enddo + call ESMF_ConfigFindLabel(LIS_config,"AGRMET native imax:",rc=rc) do n=1,LIS_rc%nnest call ESMF_ConfigGetAttribute(LIS_config,agrmet_struc(n)%imaxnative,rc=rc) @@ -1170,7 +1229,7 @@ subroutine readcrd_agrmet() if (LIS_masterproc) then ios = LIS_create_subdirs(len_trim(c_string),trim(c_string)) if (ios .ne. 0) then - write(LIS_logunit,*)'ERR creating directory ', & + write(LIS_logunit,*)'[ERR] Cannot create directory ', & trim(agrmet_struc(n)%analysisdir) flush(LIS_logunit) end if diff --git a/lis/routing/RAPID_router/RAPID_routing_run.F90 b/lis/routing/RAPID_router/RAPID_routing_run.F90 index 6264b0938..b1a57cead 100644 --- a/lis/routing/RAPID_router/RAPID_routing_run.F90 +++ b/lis/routing/RAPID_router/RAPID_routing_run.F90 @@ -198,6 +198,7 @@ subroutine RAPID_routing_run(n) writeint=RAPID_routing_struc(n)%outInterval) ! run RAPID +#ifdef PETSc call RAPID_model_main (n,RAPID_routing_struc(n)%bQinit,RAPID_routing_struc(n)%bQfinal,RAPID_routing_struc(n)%bV, & RAPID_routing_struc(n)%bhum,RAPID_routing_struc(n)%bfor,RAPID_routing_struc(n)%bdam, & RAPID_routing_struc(n)%binfluence,RAPID_routing_struc(n)%buq, & @@ -209,7 +210,7 @@ subroutine RAPID_routing_run(n) RAPID_routing_struc(n)%nmlfile,qout_filename, & LIS_rc%gnc(n),LIS_rc%gnr(n),surface_runoff,baseflow,RAPID_routing_struc(n)%initCheck, & RAPID_routing_struc(n)%dt,RAPID_routing_struc(n)%routingInterval) - +#endif deallocate(surface_runoff) deallocate(baseflow) diff --git a/lis/surfacemodels/land/noah.3.9/noah39_readrst.F90 b/lis/surfacemodels/land/noah.3.9/noah39_readrst.F90 index d5c716aa1..574f3a153 100644 --- a/lis/surfacemodels/land/noah.3.9/noah39_readrst.F90 +++ b/lis/surfacemodels/land/noah.3.9/noah39_readrst.F90 @@ -102,13 +102,13 @@ subroutine noah39_readrst() if(.not.file_exists) then write(LIS_logunit,*) '[ERR] Noah-3.9 restart file ', & - noah39_struc(n)%rfile,' does not exist ' + trim(noah39_struc(n)%rfile),' does not exist ' write(LIS_logunit,*) '[ERR] Program stopping ...' call LIS_endrun() endif write(LIS_logunit,*) & - '[INFO] Noah-3.9 restart file used: ',noah39_struc(n)%rfile + '[INFO] Noah-3.9 restart file used: ',trim(noah39_struc(n)%rfile) if(wformat.eq."binary") then ftn = LIS_getNextUnitNumber() diff --git a/lis/surfacemodels/openwater/template/readtemplateOpenWatercrd.F90 b/lis/surfacemodels/openwater/template/readtemplateOpenWatercrd.F90 index f681712b4..051d60e60 100644 --- a/lis/surfacemodels/openwater/template/readtemplateOpenWatercrd.F90 +++ b/lis/surfacemodels/openwater/template/readtemplateOpenWatercrd.F90 @@ -47,6 +47,6 @@ subroutine readtemplateOpenWatercrd() call LIS_parseTimeString(time,templateOpenWater_struc(n)%ts) enddo - write(LIS_logunit,*)'Running Template Open water Option:' + write(LIS_logunit,*)'[INFO] Running Template Open water Option:' end subroutine readtemplateOpenWatercrd