MODULE cgutils IMPLICIT NONE PRIVATE ! All publicly accessible functions and variables need to be ! declared as such here. All public names are prefixed with ! "gributils_" PUBLIC :: cgutils_check_ecmwf_grib, & & cgutils_check_ncep_grib2, & & cgutils_check_fv3_grib2 ! These codes depict the origin/model of the GRIB file INTEGER, PARAMETER :: GRIB_CENTRE_NCEP = 7, & & GRIB_CENTRE_ECMWF = 98 ! Return codes INTEGER, PARAMETER :: RETURN_CODE_NORMAL = 0, & & RETURN_CODE_GRIBFILE_ERROR = 1, & & RETURN_CODE_UTILITY_ERROR = 2 ! String lengths INTEGER, PARAMETER :: GRIB_ATTRIBUTE_STR_LEN = 64 ! Global (within this module) variables ! All global variables are prefixed with "global_" ! 3D variables of interest, T/F flag for each level present ! ! Variables names correspond more to "Flexpart" names than ! actual GRIB names ! ! ECMWF uses - u, v, t, q, etadot ! ! NCEP uses - u, v, t, rh, w LOGICAL, DIMENSION(:), ALLOCATABLE :: global_levels_u, & & global_levels_v, & & global_levels_t, & & global_levels_q, & & global_levels_etadot, & & global_levels_w, & & global_levels_rh ! 2D variables of interest, T/F flag for each ! ! Variables names correspond more to "Flexpart" names than ! actual GRIB names ! ! ECMWF uses - ps, sd, msl, tcc, u10, v10, t2, td2, lsprec, ! convprec, shf, sr, ewss, nsss, oro, excessoro, lsm ! ! NCEP uses - ps, sd, msl, u10, v10, t2, td2, tcc, lsprec, convprec, ! oro, lsm, hmix, rh2, tsig1, usig1, vsig1 LOGICAL :: global_2dvar_ps, & & global_2dvar_sd, & & global_2dvar_msl, & & global_2dvar_tcc, & & global_2dvar_u10, & & global_2dvar_v10, & & global_2dvar_t2, & & global_2dvar_td2, & & global_2dvar_lsprec, & & global_2dvar_convprec, & & global_2dvar_shf, & & global_2dvar_sr, & & global_2dvar_ewss, & & global_2dvar_nsss, & & global_2dvar_oro, & & global_2dvar_excessoro, & & global_2dvar_lsm, & & global_2dvar_hmix, & & global_2dvar_rh2, & & global_2dvar_tsig1, & & global_2dvar_usig1, & & global_2dvar_vsig1 CONTAINS SUBROUTINE allocate_and_init_flags(expected_num_levels) IMPLICIT NONE ! This is used by ECMWF, NCEP and FV3 routines. There are ! several variables used by one, that isn't used by the ! other, but I think the memory load is negligible and ! prefer to take this opportunity to present a little bit ! of clean, reusable code. INTEGER, INTENT(IN) :: expected_num_levels ! 3D variables (NCEP, FV3 and ECMWF) ALLOCATE(global_levels_u(expected_num_levels)) global_levels_u = .FALSE. ALLOCATE(global_levels_v(expected_num_levels)) global_levels_v = .FALSE. ALLOCATE(global_levels_t(expected_num_levels)) global_levels_t = .FALSE. ALLOCATE(global_levels_q(expected_num_levels)) global_levels_q = .FALSE. ALLOCATE(global_levels_etadot(expected_num_levels)) global_levels_etadot = .FALSE. ALLOCATE(global_levels_w(expected_num_levels)) global_levels_w = .FALSE. ALLOCATE(global_levels_rh(expected_num_levels)) global_levels_rh = .FALSE. ! 2D variables (both NCEP, ECMWF and FV3) global_2dvar_ps = .FALSE. global_2dvar_sd = .FALSE. global_2dvar_msl = .FALSE. global_2dvar_tcc = .FALSE. global_2dvar_u10 = .FALSE. global_2dvar_v10 = .FALSE. global_2dvar_t2 = .FALSE. global_2dvar_td2 = .FALSE. global_2dvar_lsprec = .FALSE. global_2dvar_convprec = .FALSE. global_2dvar_shf = .FALSE. global_2dvar_sr = .FALSE. global_2dvar_ewss = .FALSE. global_2dvar_nsss = .FALSE. global_2dvar_oro = .FALSE. global_2dvar_excessoro = .FALSE. global_2dvar_lsm = .FALSE. global_2dvar_hmix = .FALSE. global_2dvar_rh2 = .FALSE. global_2dvar_tsig1 = .FALSE. global_2dvar_usig1 = .FALSE. global_2dvar_vsig1 = .FALSE. END SUBROUTINE allocate_and_init_flags SUBROUTINE delete_levels_flags IMPLICIT NONE DEALLOCATE(global_levels_u) DEALLOCATE(global_levels_v) DEALLOCATE(global_levels_t) DEALLOCATE(global_levels_q) DEALLOCATE(global_levels_etadot) DEALLOCATE(global_levels_w) DEALLOCATE(global_levels_rh) END SUBROUTINE delete_levels_flags LOGICAL FUNCTION ecmwf_all_2dvars_present() IMPLICIT NONE ! Checks logical flags for each 2d variable, testing for ! all true. ! ! For sake of "situational awareness," this routine will ! go through the logical flags of each 2d array of interest, ! printing error message for each one with a problem, then ! returning the error code after all that. ecmwf_all_2dvars_present = .TRUE. ! surface levels ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_ps) IF ( .NOT. global_2dvar_ps) & & WRITE(6,*) 'ERROR: sp: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_sd) IF ( .NOT. global_2dvar_sd) & & WRITE(6,*) 'ERROR: sd (or sde): not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_tcc) IF ( .NOT. global_2dvar_tcc) & & WRITE(6,*) 'ERROR: tcc: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_lsprec) IF ( .NOT. global_2dvar_lsprec) & & WRITE(6,*) 'ERROR: lsp: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_convprec) IF ( .NOT. global_2dvar_convprec) & & WRITE(6,*) 'ERROR: acpcp: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_shf) IF ( .NOT. global_2dvar_shf) & & WRITE(6,*) 'ERROR: sshf: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_sr) IF ( .NOT. global_2dvar_sr) & & WRITE(6,*) 'ERROR: ssr: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_ewss) IF ( .NOT. global_2dvar_ewss) & & WRITE(6,*) 'ERROR: ewss: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_nsss) IF ( .NOT. global_2dvar_nsss) & & WRITE(6,*) 'ERROR: nsss: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_oro) IF ( .NOT. global_2dvar_oro) & & WRITE(6,*) 'ERROR: z: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_excessoro) IF ( .NOT. global_2dvar_excessoro) & & WRITE(6,*) 'ERROR: sdor: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_lsm) IF ( .NOT. global_2dvar_lsm) & & WRITE(6,*) 'ERROR: lsm: not found' ! heightAboveGround levels ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_u10) IF ( .NOT. global_2dvar_u10) & & WRITE(6,*) 'ERROR: 10u: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_v10) IF ( .NOT. global_2dvar_v10) & & WRITE(6,*) 'ERROR: 10v: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_t2) IF ( .NOT. global_2dvar_t2) & & WRITE(6,*) 'ERROR: 2t: not found' ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_td2) IF ( .NOT. global_2dvar_td2) & & WRITE(6,*) 'ERROR: 2d: not found' ! meanSea levels ecmwf_all_2dvars_present = & & (ecmwf_all_2dvars_present .AND. global_2dvar_msl) IF ( .NOT. global_2dvar_msl) & & WRITE(6,*) 'ERROR: msl: not found' RETURN END FUNCTION ecmwf_all_2dvars_present LOGICAL FUNCTION ncep_all_2dvars_present() IMPLICIT NONE ! Checks logical flags for each 2d variable, testing for ! all true. ! ! For sake of "situational awareness," this routine will ! go through the logical flags of each 2d array of interest, ! printing error message for each one with a problem, then ! returning the error code after all that. ncep_all_2dvars_present = .TRUE. ! surface levels ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_ps) IF ( .NOT. global_2dvar_ps) & & WRITE(6,*) 'ERROR: sp: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_sd) IF ( .NOT. global_2dvar_sd) & & WRITE(6,*) 'ERROR: sdwe: not found' !!! 2019-01-20 - These two variables are currently not part of what I find in !!! NCEP GRIB files at CTBTO. I think they used to be available. So, I'm !!! commenting out their testing right now. ! ncep_all_2dvars_present = & !& (ncep_all_2dvars_present .AND. global_2dvar_lsprec) ! IF ( .NOT. global_2dvar_lsprec) & !& WRITE(6,*) 'ERROR: tp: not found' ! ! ncep_all_2dvars_present = & !& (ncep_all_2dvars_present .AND. global_2dvar_convprec) ! IF ( .NOT. global_2dvar_convprec) & !& WRITE(6,*) 'ERROR: acpcp: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_oro) IF ( .NOT. global_2dvar_oro) & & WRITE(6,*) 'ERROR: orog: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_lsm) IF ( .NOT. global_2dvar_lsm) & & WRITE(6,*) 'ERROR: lsm: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_hmix) IF ( .NOT. global_2dvar_hmix) & & WRITE(6,*) 'ERROR: hpbl: not found' ! heightAboveGround levels ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_u10) IF ( .NOT. global_2dvar_u10) & & WRITE(6,*) 'ERROR: 10u: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_v10) IF ( .NOT. global_2dvar_v10) & & WRITE(6,*) 'ERROR: 10v: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_t2) IF ( .NOT. global_2dvar_t2) & & WRITE(6,*) 'ERROR: 2t: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_td2) IF ( .NOT. global_2dvar_td2) & & WRITE(6,*) 'ERROR: 2d: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_rh2) IF ( .NOT. global_2dvar_rh2) & & WRITE(6,*) 'ERROR: 2m r: not found' ! meanSea levels ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_msl) IF ( .NOT. global_2dvar_msl) & & WRITE(6,*) 'ERROR: prmsl: not found' ! sigma levels ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_tsig1) IF ( .NOT. global_2dvar_tsig1) & & WRITE(6,*) 'ERROR: sigma t: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_usig1) IF ( .NOT. global_2dvar_usig1) & & WRITE(6,*) 'ERROR: sigma u: not found' ncep_all_2dvars_present = & & (ncep_all_2dvars_present .AND. global_2dvar_vsig1) IF ( .NOT. global_2dvar_vsig1) & & WRITE(6,*) 'ERROR: sigma v: not found' !!! 2019-01-20 - These variable is currently not part of what I find in !!! NCEP GRIB files at CTBTO. I think it used to be available. So, I'm !!! commenting out the testing right now. ! Low cloud layer levels ! ncep_all_2dvars_present = & !& (ncep_all_2dvars_present .AND. global_2dvar_tcc) ! IF ( .NOT. global_2dvar_tcc) & !& WRITE(6,*) 'ERROR: tcc: not found' RETURN END FUNCTION ncep_all_2dvars_present LOGICAL FUNCTION fv3_all_2dvars_present() IMPLICIT NONE ! Checks logical flags for each 2d variable, testing for ! all true. ! ! For sake of "situational awareness," this routine will ! go through the logical flags of each 2d array of interest, ! printing error message for each one with a problem, then ! returning the error code after all that. fv3_all_2dvars_present = .TRUE. ! surface levels fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_ps) IF ( .NOT. global_2dvar_ps) & & WRITE(6,*) 'ERROR: sp: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_sd) IF ( .NOT. global_2dvar_sd) & & WRITE(6,*) 'ERROR: sdwe: not found' !!! 2019-01-20 - These two variables are currently not part of what I find in !!! NCEP GRIB files at CTBTO. I think they used to be available. So, I'm !!! commenting out their testing right now. ! ncep_all_2dvars_present = & !& (ncep_all_2dvars_present .AND. global_2dvar_lsprec) ! IF ( .NOT. global_2dvar_lsprec) & !& WRITE(6,*) 'ERROR: tp: not found' ! ! ncep_all_2dvars_present = & !& (ncep_all_2dvars_present .AND. global_2dvar_convprec) ! IF ( .NOT. global_2dvar_convprec) & !& WRITE(6,*) 'ERROR: acpcp: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_oro) IF ( .NOT. global_2dvar_oro) & & WRITE(6,*) 'ERROR: orog: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_lsm) IF ( .NOT. global_2dvar_lsm) & & WRITE(6,*) 'ERROR: lsm: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_hmix) IF ( .NOT. global_2dvar_hmix) & & WRITE(6,*) 'ERROR: hpbl: not found' ! heightAboveGround levels fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_u10) IF ( .NOT. global_2dvar_u10) & & WRITE(6,*) 'ERROR: 10u: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_v10) IF ( .NOT. global_2dvar_v10) & & WRITE(6,*) 'ERROR: 10v: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_t2) IF ( .NOT. global_2dvar_t2) & & WRITE(6,*) 'ERROR: 2t: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_td2) IF ( .NOT. global_2dvar_td2) & & WRITE(6,*) 'ERROR: 2d: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_rh2) IF ( .NOT. global_2dvar_rh2) & & WRITE(6,*) 'ERROR: 2m r: not found' ! meanSea levels fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_msl) IF ( .NOT. global_2dvar_msl) & & WRITE(6,*) 'ERROR: prmsl: not found' ! sigma levels fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_tsig1) IF ( .NOT. global_2dvar_tsig1) & & WRITE(6,*) 'ERROR: sigma t: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_usig1) IF ( .NOT. global_2dvar_usig1) & & WRITE(6,*) 'ERROR: sigma u: not found' fv3_all_2dvars_present = & & (fv3_all_2dvars_present .AND. global_2dvar_vsig1) IF ( .NOT. global_2dvar_vsig1) & & WRITE(6,*) 'ERROR: sigma v: not found' !!! 2019-01-20 - These variable is currently not part of what I find in !!! NCEP GRIB files at CTBTO. I think it used to be available. So, I'm !!! commenting out the testing right now. ! Low cloud layer levels ! ncep_all_2dvars_present = & !& (ncep_all_2dvars_present .AND. global_2dvar_tcc) ! IF ( .NOT. global_2dvar_tcc) & !& WRITE(6,*) 'ERROR: tcc: not found' RETURN END FUNCTION fv3_all_2dvars_present LOGICAL FUNCTION ecmwf_all_levels_present() IMPLICIT NONE ! Checks logical flags for each 3d variable, testing for ! all true. ! ! For sake of "situational awareness," this routine will ! go through the logical flags of each 3d array of interest, ! printing error message for each one with a problem, then ! returning the error code after all that. ecmwf_all_levels_present = .TRUE. ecmwf_all_levels_present = & & (ecmwf_all_levels_present .AND. ALL(global_levels_u)) IF ( .NOT. ALL(global_levels_u) ) & & WRITE(6,*) 'ERROR: u: not all levels found' ecmwf_all_levels_present = & & (ecmwf_all_levels_present .AND. ALL(global_levels_v)) IF ( .NOT. ALL(global_levels_v) ) & & WRITE(6,*) 'ERROR: v: not all levels found' ecmwf_all_levels_present = & & (ecmwf_all_levels_present .AND. ALL(global_levels_t)) IF ( .NOT. ALL(global_levels_t) ) & & WRITE(6,*) 'ERROR: t: not all levels found' ecmwf_all_levels_present = & & (ecmwf_all_levels_present .AND. ALL(global_levels_q)) IF ( .NOT. ALL(global_levels_q) ) & & WRITE(6,*) 'ERROR: q: not all levels found' ecmwf_all_levels_present = & & (ecmwf_all_levels_present .AND. ALL(global_levels_etadot)) IF ( .NOT. ALL(global_levels_etadot) ) & & WRITE(6,*) 'ERROR: etadot: not all levels found' END FUNCTION ecmwf_all_levels_present LOGICAL FUNCTION ncep_all_levels_present() IMPLICIT NONE ! Checks logical flags for each 3d variable, testing for ! all true. ! ! For sake of "situational awareness," this routine will ! go through the logical flags of each 3d array of interest, ! printing error message for each one with a problem, then ! returning the error code after all that. ncep_all_levels_present = .TRUE. ncep_all_levels_present = & & (ncep_all_levels_present .AND. ALL(global_levels_u)) IF ( .NOT. ALL(global_levels_u) ) & & WRITE(6,*) 'ERROR: u: not all levels found' ncep_all_levels_present = & & (ncep_all_levels_present .AND. ALL(global_levels_v)) IF ( .NOT. ALL(global_levels_v) ) & & WRITE(6,*) 'ERROR: v: not all levels found' ncep_all_levels_present = & & (ncep_all_levels_present .AND. ALL(global_levels_t)) IF ( .NOT. ALL(global_levels_t) ) & & WRITE(6,*) 'ERROR: t: not all levels found' ncep_all_levels_present = & & (ncep_all_levels_present .AND. ALL(global_levels_rh)) IF ( .NOT. ALL(global_levels_rh) ) & & WRITE(6,*) 'ERROR: r: not all levels found' ! We can't check w like the others, because it only seems to be included for messages ! with pressure levels 100 mb and higher. So, this is checked elsewhere ! ncep_all_levels_present = & !& (ncep_all_levels_present .AND. ALL(global_levels_w)) ! IF ( .NOT. ALL(global_levels_w) ) & !& WRITE(6,*) 'ERROR: w: not all levels found' RETURN END FUNCTION ncep_all_levels_present LOGICAL FUNCTION fv3_all_levels_present() IMPLICIT NONE ! Checks logical flags for each 3d variable, testing for ! all true. ! ! For sake of "situational awareness," this routine will ! go through the logical flags of each 3d array of interest, ! printing error message for each one with a problem, then ! returning the error code after all that. fv3_all_levels_present = .TRUE. fv3_all_levels_present = & & (fv3_all_levels_present .AND. ALL(global_levels_u)) IF ( .NOT. ALL(global_levels_u) ) & & WRITE(6,*) 'ERROR: u: not all levels found' fv3_all_levels_present = & & (fv3_all_levels_present .AND. ALL(global_levels_v)) IF ( .NOT. ALL(global_levels_v) ) & & WRITE(6,*) 'ERROR: v: not all levels found' fv3_all_levels_present = & & (fv3_all_levels_present .AND. ALL(global_levels_t)) IF ( .NOT. ALL(global_levels_t) ) & & WRITE(6,*) 'ERROR: t: not all levels found' fv3_all_levels_present = & & (fv3_all_levels_present .AND. ALL(global_levels_rh)) IF ( .NOT. ALL(global_levels_rh) ) & & WRITE(6,*) 'ERROR: r: not all levels found' ! We can't check w like the others, because it only seems to be included for messages ! with pressure levels 100 mb and higher. So, this is checked elsewhere ! ncep_all_levels_present = & !& (ncep_all_levels_present .AND. ALL(global_levels_w)) ! IF ( .NOT. ALL(global_levels_w) ) & !& WRITE(6,*) 'ERROR: w: not all levels found' RETURN END FUNCTION fv3_all_levels_present LOGICAL FUNCTION ncep_all_expected_w_levels_present(expected_num_levels, millibars) IMPLICIT NONE ! It seems like the w messages are only available at 100mb and below, so I can't check ! this variable easily like I can the other 3d variables. I need to explictly go through ! each value in the millibars array and, if the value is 100 or greaters, use its index to ! determine if the corresponding value in the flag array is .TRUE. or .FALSE. INTEGER, INTENT(IN) :: expected_num_levels INTEGER, INTENT(IN), DIMENSION(expected_num_levels) :: millibars INTEGER :: i, plevel LOGICAL :: all_present all_present = .TRUE. DO i=1,expected_num_levels plevel = millibars(i) IF (plevel .GE. 100) all_present = (all_present .AND. global_levels_w(i)) END DO ncep_all_expected_w_levels_present = all_present RETURN END FUNCTION ncep_all_expected_w_levels_present LOGICAL FUNCTION fv3_all_expected_w_levels_present(expected_num_levels, millibars) IMPLICIT NONE ! It seems like the w messages are only available at 100mb and below, so I can't check ! this variable easily like I can the other 3d variables. I need to explictly go through ! each value in the millibars array and, if the value is 100 or greaters, use its index to ! determine if the corresponding value in the flag array is .TRUE. or .FALSE. INTEGER, INTENT(IN) :: expected_num_levels INTEGER, INTENT(IN), DIMENSION(expected_num_levels) :: millibars INTEGER :: i, plevel LOGICAL :: all_present all_present = .TRUE. DO i=1,expected_num_levels plevel = millibars(i) IF (plevel .GE. 100) all_present = (all_present .AND. global_levels_w(i)) END DO fv3_all_expected_w_levels_present = all_present RETURN END FUNCTION fv3_all_expected_w_levels_present LOGICAL FUNCTION message_is_good(igrib, check_anl, check_fcst, check_grib2) IMPLICIT NONE ! This function is intended to drive several checks on a ! grib message and, as requirements change, this function ! may become more complex. ! ! For now, in January 2019, it checks for message being ! grib version 2, and, depending on the settings of ! check_anl and check_fcst, checks whether the message ! is of the specified type ! ! August 2019 - the checking for GRIB2 message is now conditional ! ! It is assumed that igrib refers to an already acccessed ! and active grib message INTEGER, INTENT(IN) :: igrib LOGICAL, INTENT(IN) :: check_anl, check_fcst, check_grib2 LOGICAL :: fcst_anl_good message_is_good = .TRUE. IF (check_grib2) THEN IF (.NOT. message_is_grib2(igrib)) THEN WRITE(6,*) 'Grib message not GRIB2' message_is_good = .FALSE. RETURN ENDIF ENDIF CALL anl_fcst_check(igrib, check_anl, check_fcst, fcst_anl_good) IF (.NOT. fcst_anl_good) THEN WRITE(6,*) 'Failed fcst/anl check' message_is_good = .FALSE. RETURN ENDIF END FUNCTION message_is_good LOGICAL FUNCTION message_is_forecast(igrib) USE grib_api IMPLICIT NONE INTEGER, INTENT(IN) :: igrib ! igrib is an id for a valid, already-fetched, grib message CHARACTER(LEN=8) :: data_type ! Looking for 'fc' or 'an' message_is_forecast = .TRUE. CALL grib_get(igrib, 'dataType', data_type) IF (TRIM(data_type) .NE. 'fc') message_is_forecast = .FALSE. RETURN END FUNCTION message_is_forecast LOGICAL FUNCTION message_is_analysis(igrib) USE grib_api IMPLICIT NONE INTEGER, INTENT(IN) :: igrib ! igrib is an id for a valid, already-fetched, grib message CHARACTER(LEN=8) :: data_type ! Looking for 'fc' or 'an' message_is_analysis = .TRUE. CALL grib_get(igrib, 'dataType', data_type) IF (data_type .NE. 'an') message_is_analysis = .FALSE. RETURN END FUNCTION message_is_analysis SUBROUTINE anl_fcst_check(igrib, check_for_anl, check_for_fcst, result) USE grib_api IMPLICIT NONE ! If either (but not both) the check_for_anl or check_for_fcst ! flag is true, then check that message is good for that INTEGER, INTENT(IN) :: igrib ! igrib is an id for a valid, already-fetched, grib message LOGICAL, INTENT(IN) :: check_for_anl, check_for_fcst LOGICAL, INTENT(OUT) :: result ! Only one of the two checks can be true IF (check_for_anl .AND. check_for_fcst) THEN WRITE(0,*) 'Only one of fcst and anl check can be requested' CALL EXIT(RETURN_CODE_UTILITY_ERROR) ENDIF ! It's possible and valid for this subroutine to be called where ! both flags are False. That's OK, just set result to True. No ! request was made, so the test passes IF (.NOT.(check_for_anl .OR. check_for_fcst)) THEN result = .TRUE. RETURN ENDIF ! OK, if we're here, exactly one of the checks is true IF (check_for_anl) THEN result = .TRUE. IF (.NOT. message_is_analysis(igrib)) THEN result = .FALSE. RETURN ENDIF ENDIF IF (check_for_fcst) THEN result = .TRUE. IF (.NOT. message_is_forecast(igrib)) THEN result = .FALSE. RETURN ENDIF ENDIF END SUBROUTINE anl_fcst_check LOGICAL FUNCTION message_is_grib2(igrib) USE grib_api IMPLICIT NONE INTEGER, INTENT(IN) :: igrib ! igrib is an id for a valid, already-fetched, grib message INTEGER :: edition_number message_is_grib2 = .TRUE. CALL grib_get(igrib, 'editionNumber', edition_number) IF (edition_number .NE. 2) message_is_grib2 = .FALSE. RETURN END FUNCTION message_is_grib2 SUBROUTINE build_millibars_array(ifile, expected_num_levels, & & millibars, return_code) USE grib_api IMPLICIT NONE INTEGER, INTENT(IN) :: ifile, expected_num_levels INTEGER, INTENT(OUT), DIMENSION(expected_num_levels) :: millibars INTEGER, INTENT(OUT) :: return_code ! This routine is for the old GFS (pre-FV3) GRIB files, where it is assumed ! that all isobaric-level variables will have the same number of levels. ! ! Goes through each message in the already-open GRIBFILE ! indicated by "ifile" and creates a sorted array of all of the ! pressures found in isobaric level types. When we're done, ! we have an inventory of all of the pressure levels expected by ! 3d isobaric-level variables. ! ! Note that the returned millibars array is NOT sorted. I can't ! think of a reason right now why that would be necessary. ! ! If, in creating the inventory, we end up with fewer or more ! levels than expected by "expected_num_levels" we return with ! an error INTEGER :: level_count ! Keep track of number of new levels found LOGICAL :: end_of_file INTEGER :: igrib, iret CHARACTER(LEN=GRIB_ATTRIBUTE_STR_LEN) :: gribmsg_leveltype, & & gribmsg_shortname INTEGER :: gribmsg_level return_code = RETURN_CODE_NORMAL level_count = 0 millibars = -999 !PRINT *, 'millibars: ', millibars ! Iterate through the entire file, keeping track of what we've found end_of_file = .FALSE. DO WHILE (.NOT. end_of_file) CALL grib_new_from_file(ifile, igrib, iret) IF (iret .eq. GRIB_END_OF_FILE) THEN end_of_file = .TRUE. ELSE CALL grib_get(igrib, 'shortName', gribmsg_shortname) CALL grib_get(igrib, 'typeOfLevel', gribmsg_leveltype) CALL grib_get(igrib, 'level', gribmsg_level) ! We're only interested in the isobaric level messages IF ( TRIM(gribmsg_leveltype) == 'isobaricInhPa' ) THEN CALL grib_get(igrib, 'level', gribmsg_level) !PRINT *, 'level: ', gribmsg_level ! If this level is not already in the array, update ! the counter (make sure we haven't exceeded expected ! number of levels), then add to array. IF (myfindloc(expected_num_levels, & & millibars, VALUE=gribmsg_level) .EQ. 0) THEN level_count = level_count + 1 IF (level_count .LE. expected_num_levels) THEN millibars(level_count) = gribmsg_level ELSE WRITE(6,"('ERROR: Found more than ', I3, ' levels')") & & expected_num_levels return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ENDIF ENDIF END IF END DO ! End of main loop for reading each GRIB message ! Make sure we found the number if isobaric levels that we expected IF (level_count .LT. expected_num_levels) THEN WRITE(6,"('ERROR: Only found ', I3, ' levels')") level_count return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! OK, if we made it here, we should have a full millibars array ! containing all of the isobaric levels we found. NOTE THAT I ! AM NOT BOTHERING TO SORT THIS - I DON'T SEE A REASON FOR THAT ! SO FAR. RETURN END SUBROUTINE build_millibars_array SUBROUTINE build_fv3_millibars_array(ifile, expected_num_levels, & & millibars, return_code) USE grib_api IMPLICIT NONE INTEGER, INTENT(IN) :: ifile, expected_num_levels INTEGER, INTENT(OUT), DIMENSION(expected_num_levels) :: millibars INTEGER, INTENT(OUT) :: return_code ! Goes through each message in the already-open GRIBFILE ! indicated by "ifile" and creates a sorted array of all of the ! pressures found in isobaric level types. When we're done, ! we have an inventory of all of the pressure levels expected by ! 3d isobaric-level variables. ! ! Note that the returned millibars array is NOT sorted. I can't ! think of a reason right now why that would be necessary. ! ! If, in creating the inventory, we end up with fewer or more ! levels than expected by "expected_num_levels" we return with ! an error ! ! MODIFICATION - DJM - 2019-07-27 - New NCEP FV3 files have additional ! pressure levels at 15 mb and 40 mb, but the variables at those ! levels are not exactly those at the more standard levels. The ! logic of this routine previously added these levels to the list of ! millibars, but then other routines would fail because they wouldn't ! find expected variables at these levels. So, I have modified this ! routine to simple ignore pressure levels of 15 and 40 millibars. ! This isn't necessarily the most robust way to do things, but it is ! what it is right now. INTEGER :: level_count ! Keep track of number of new levels found LOGICAL :: end_of_file INTEGER :: igrib, iret CHARACTER(LEN=GRIB_ATTRIBUTE_STR_LEN) :: gribmsg_leveltype, & & gribmsg_shortname INTEGER :: gribmsg_level return_code = RETURN_CODE_NORMAL level_count = 0 millibars = -999 !PRINT *, 'millibars: ', millibars ! Iterate through the entire file, keeping track of what we've found end_of_file = .FALSE. DO WHILE (.NOT. end_of_file) CALL grib_new_from_file(ifile, igrib, iret) IF (iret .eq. GRIB_END_OF_FILE) THEN end_of_file = .TRUE. ELSE CALL grib_get(igrib, 'shortName', gribmsg_shortname) CALL grib_get(igrib, 'typeOfLevel', gribmsg_leveltype) CALL grib_get(igrib, 'level', gribmsg_level) ! We're only interested in the isobaric level messages IF ( TRIM(gribmsg_leveltype) == 'isobaricInhPa' ) THEN CALL grib_get(igrib, 'level', gribmsg_level) !PRINT *, 'level: ', gribmsg_level IF ( (gribmsg_level .EQ. 15) .OR. (gribmsg_level .EQ. 40) ) THEN ! Bypass this message - we are ignoring 15 and 40 mb ! pressure levels. See comments above CONTINUE ELSE ! If this level is not already in the array, update ! the counter (make sure we haven't exceeded expected ! number of levels), then add to array. IF (myfindloc(expected_num_levels, & & millibars, VALUE=gribmsg_level) .EQ. 0) THEN level_count = level_count + 1 IF (level_count .LE. expected_num_levels) THEN millibars(level_count) = gribmsg_level ELSE WRITE(6,"('ERROR: Found more than ', I3, ' levels')") & & expected_num_levels return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ENDIF ENDIF ENDIF END IF END DO ! End of main loop for reading each GRIB message ! Make sure we found the number of isobaric levels that we expected IF (level_count .LT. expected_num_levels) THEN WRITE(6,"('ERROR: Only found ', I3, ' levels')") level_count return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! OK, if we made it here, we should have a full millibars array ! containing all of the isobaric levels we found. NOTE THAT I ! AM NOT BOTHERING TO SORT THIS - I DON'T SEE A REASON FOR THAT ! SO FAR. RETURN END SUBROUTINE build_fv3_millibars_array INTEGER FUNCTION myfindloc(n, array, value) INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN), DIMENSION(n) :: array INTEGER, INTENT(IN) :: value ! My own implementation of FINDLOC, which is part of the ! Fortran 2008 standard ! Returns index of first instance of "value" in the array. Returns ! zero if not found. Assumes single-dimension INTEGER array ! Default value INTEGER :: i LOGICAL :: found found = .FALSE. myfindloc = 0 DO i=1,SIZE(array) IF (array(i) .EQ. value) THEN myfindloc = i EXIT ENDIF END DO RETURN END FUNCTION myfindloc SUBROUTINE write_err(shortname, level) IMPLICIT NONE CHARACTER(LEN=GRIB_ATTRIBUTE_STR_LEN), INTENT(IN) :: shortname INTEGER, INTENT(IN) :: level WRITE(6,*) 'ERROR: Bad gribmsg_shortname: ', TRIM(shortname),& & ', gribmsg_level: ', level END SUBROUTINE write_err INTEGER FUNCTION cgutils_check_ecmwf_grib(filepath, expected_num_levels, & & check_fcst, check_anl, check_grib2) & & RESULT(return_code) USE grib_api IMPLICIT NONE CHARACTER(LEN=256), INTENT(IN) :: filepath INTEGER, INTENT(IN) :: expected_num_levels LOGICAL, INTENT(IN) :: check_fcst, check_anl, check_grib2 ! return_code: possible values defined in module header, RETURN_CODE_* INTEGER :: ifile, igrib, iret, grib_centre CHARACTER(LEN=GRIB_ATTRIBUTE_STR_LEN) :: gribmsg_leveltype, & & gribmsg_shortname INTEGER :: gribmsg_level LOGICAL :: end_of_file !!LOGICAL :: fcst_anl_good ! Start with return_code 0. Anything wrong will change it return_code = RETURN_CODE_NORMAL CALL grib_open_file(ifile, filepath, 'r', iret) IF (iret == 0) THEN ! Use first record to detect centre, which is assumed constant ! amongst all messages CALL grib_new_from_file(ifile, igrib, iret) CALL grib_get(igrib, 'centre', grib_centre) CALL grib_close_file(ifile) ELSE WRITE(0,*) 'ERROR: problem opening GRIB file: '//TRIM(filepath) return_code = RETURN_CODE_UTILITY_ERROR RETURN END IF IF (grib_centre .NE. GRIB_CENTRE_ECMWF) THEN WRITE(6,*) 'ERROR: Unexpected grib_centre: ', grib_centre return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF CALL allocate_and_init_flags(expected_num_levels) ! Iterate through the entire file, keeping track of what we've found end_of_file = .FALSE. CALL grib_open_file(ifile, filepath, 'r', iret) DO WHILE (.NOT. end_of_file) CALL grib_new_from_file(ifile, igrib, iret) IF (iret .eq. GRIB_END_OF_FILE) THEN end_of_file = .TRUE. ELSE CALL grib_get(igrib, 'shortName', gribmsg_shortname) CALL grib_get(igrib, 'typeOfLevel', gribmsg_leveltype) CALL grib_get(igrib, 'level', gribmsg_level) ! Check general attributes of this message. Note that ! all messages will be checked, even those we might not ! care about. As of 2019-01-17, the checks are for GRIB2 ! and for fcst/anl if specified by "check_anl" or "check_fcst" IF (.NOT. message_is_good(igrib, & & check_anl, check_fcst, check_grib2)) THEN CALL write_err(gribmsg_shortname, gribmsg_level) return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! Filter for 3D variables of interest IF ( TRIM(gribmsg_leveltype) == 'hybrid' ) THEN CALL grib_get(igrib, 'level', gribmsg_level) ! Check that there aren't more levels than expected IF (gribmsg_level .GT. expected_num_levels) THEN WRITE(6,*) 'ERROR: gribmsg_shortname: ', & & TRIM(gribmsg_shortname), & & ', gribmsg_level: ', gribmsg_level, & & ', expected max levels: ', & & expected_num_levels return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF SELECT CASE(gribmsg_shortname) CASE('u') global_levels_u(gribmsg_level) = .TRUE. CASE('v') global_levels_v(gribmsg_level) = .TRUE. CASE('t') global_levels_t(gribmsg_level) = .TRUE. CASE('q') global_levels_q(gribmsg_level) = .TRUE. CASE('etadot') global_levels_etadot(gribmsg_level) = .TRUE. END SELECT !!!!!! ! End of section for 3d variables of interest ! Filter for 2d variables of interest with level type ! of "surface" ELSE IF ( TRIM(gribmsg_leveltype) == 'surface' ) THEN SELECT CASE(gribmsg_shortname) CASE('sp') global_2dvar_ps = .TRUE. ! For snow depth, I've been seeing both 'sd' and ! 'sde' used in the CTBTO ECMWF files CASE('sd') global_2dvar_sd = .TRUE. CASE('sde') global_2dvar_sd = .TRUE. CASE('tcc') global_2dvar_tcc = .TRUE. CASE('lsp') global_2dvar_lsprec = .TRUE. CASE('acpcp') global_2dvar_convprec = .TRUE. CASE('sshf') global_2dvar_shf = .TRUE. CASE('ssr') global_2dvar_sr = .TRUE. CASE('ewss') global_2dvar_ewss = .TRUE. CASE('nsss') global_2dvar_nsss = .TRUE. CASE('z') global_2dvar_oro = .TRUE. CASE('sdor') global_2dvar_excessoro = .TRUE. CASE('lsm') global_2dvar_lsm = .TRUE. ! GRIB 1 variables found in the EFYYMMDDHH files ! Note that these are also found in GRIB2 files ! under different level types, except for cp, which ! is acpcp in GRIB2. If there is ever a mixture ! (which I think is unlikely) in a GRIB file, then ! the world will likely end and this won't work ! correctly. CASE('cp') global_2dvar_convprec = .TRUE. CASE('10u') global_2dvar_u10 = .TRUE. CASE('10v') global_2dvar_v10 = .TRUE. CASE('2t') global_2dvar_t2 = .TRUE. CASE('2d') global_2dvar_td2 = .TRUE. CASE('msl') global_2dvar_msl = .TRUE. END SELECT !!!!!!!! End of section for 2d "surface" variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'heightAboveGround' ) THEN SELECT CASE(gribmsg_shortname) CASE('10u') global_2dvar_u10 = .TRUE. CASE('10v') global_2dvar_v10 = .TRUE. CASE('2t') global_2dvar_t2 = .TRUE. CASE('2d') global_2dvar_td2 = .TRUE. END SELECT !!!!!!!! End of section for 2d "heightAboveGround" !!!!!!!! variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'meanSea' ) THEN SELECT CASE(gribmsg_shortname) CASE('msl') global_2dvar_msl = .TRUE. END SELECT !!!!!!!! End of section for 2d "meanSea" !!!!!!!! variables of interest END IF ! Selection of message types (3d, surface, etc.) END IF END DO ! End of main loop for reading each GRIB message CALL grib_close_file(ifile) ! This function will have the side effect of printing any ! problems to stdout. If false, we force exit IF (ecmwf_all_levels_present()) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR CALL delete_levels_flags RETURN ENDIF IF (ecmwf_all_2dvars_present()) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF RETURN END FUNCTION cgutils_check_ecmwf_grib INTEGER FUNCTION cgutils_check_ncep_grib2(filepath, expected_num_levels, & & check_fcst, check_anl, check_grib2) & & RESULT(return_code) USE grib_api IMPLICIT NONE CHARACTER(LEN=256), INTENT(IN) :: filepath INTEGER, INTENT(IN) :: expected_num_levels LOGICAL, INTENT(IN) :: check_fcst, check_anl, check_grib2 ! return_code: possible values defined in module header, RETURN_CODE_* INTEGER :: ifile, igrib, iret, grib_centre CHARACTER(LEN=GRIB_ATTRIBUTE_STR_LEN) :: gribmsg_leveltype, & & gribmsg_shortname INTEGER :: gribmsg_level, plevel_idx LOGICAL :: end_of_file INTEGER, ALLOCATABLE, DIMENSION(:) :: millibars ! Start with return_code 0. Anything wrong will change it return_code = RETURN_CODE_NORMAL CALL grib_open_file(ifile, filepath, 'r', iret) IF (iret == 0) THEN ! Use first record to detect centre, which is assumed constant ! amongst all messages CALL grib_new_from_file(ifile, igrib, iret) CALL grib_get(igrib, 'centre', grib_centre) CALL grib_close_file(ifile) ELSE WRITE(0,*) 'ERROR: problem opening GRIB file: '//TRIM(filepath) return_code = RETURN_CODE_UTILITY_ERROR RETURN END IF IF (grib_centre .NE. GRIB_CENTRE_NCEP) THEN WRITE(6,*) 'ERROR: Unexpected grib_centre: ', grib_centre return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! Pass through the file once in order to build an inventory of ! isobaric pressure levels so that later we can make sure that ! each 3D variable has each of these pressure levels. While ! we're doing the inventory we can check that the number of ! isobaric levels found is indeed what we expect. CALL grib_open_file(ifile, filepath, 'r', iret) IF (iret == 0) THEN ALLOCATE(millibars(expected_num_levels)) CALL build_millibars_array(ifile, expected_num_levels, & & millibars, return_code) CALL grib_close_file(ifile) ELSE WRITE(0,*) 'ERROR: problem opening GRIB file: '//TRIM(filepath) return_code = RETURN_CODE_UTILITY_ERROR RETURN END IF ! Check that build_millibars_array() had no errors ! We will assume that any error messages were written from the ! subroutine IF (return_code .NE. RETURN_CODE_NORMAL) THEN return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! At this point, millibars should be an array of all the pressure ! levels (isobaricInhPa) that we found. !PRINT *, 'millibars: ', millibars CALL allocate_and_init_flags(expected_num_levels) ! Now, make our primary checking pass through the GRIB file. ! For 3D isobaric-level variables, as we read a level, we will ! find its index in millibars array, and use that to index the ! correct flag element in our logical array for each 3D variable ! Iterate through the entire file, keeping track of what we've found end_of_file = .FALSE. CALL grib_open_file(ifile, filepath, 'r', iret) DO WHILE (.NOT. end_of_file) CALL grib_new_from_file(ifile, igrib, iret) IF (iret .eq. GRIB_END_OF_FILE) THEN end_of_file = .TRUE. ELSE CALL grib_get(igrib, 'shortName', gribmsg_shortname) CALL grib_get(igrib, 'typeOfLevel', gribmsg_leveltype) CALL grib_get(igrib, 'level', gribmsg_level) ! Check general attributes of this message. Note that ! all messages will be checked, even those we might not ! care about. As of 2019-01-17, the checks are for GRIB2 ! and for fcst/anl if specified by "check_anl" or "check_fcst" IF (.NOT. message_is_good(igrib, & & check_anl, check_fcst, check_grib2)) THEN CALL write_err(gribmsg_shortname, gribmsg_level) return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF !PRINT *, 'shortname, leveltype, level, plevel_idx: ', & !& gribmsg_shortname, gribmsg_leveltype, & !& gribmsg_level, plevel_idx ! Filter for 3D variables of interest IF ( TRIM(gribmsg_leveltype) == 'isobaricInhPa' ) THEN CALL grib_get(igrib, 'level', gribmsg_level) ! Unsorted array "millibars," created above, tells us what index we ! should use to access the retrieved level flag as we maintain an ! inventory of what we found, and what we didn't ! Use my function myfindloc() to find the index for a specified pressure ! level. plevel_idx = myfindloc(expected_num_levels, millibars, gribmsg_level) IF (plevel_idx .LE. 0) THEN CALL write_err(gribmsg_shortname, gribmsg_level) return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF !PRINT *, 'shortname, leveltype, level, plevel_idx: ', & !& gribmsg_shortname, gribmsg_leveltype, & !& gribmsg_level, plevel_idx SELECT CASE(gribmsg_shortname) CASE('u') global_levels_u(plevel_idx) = .TRUE. CASE('v') global_levels_v(plevel_idx) = .TRUE. CASE('t') global_levels_t(plevel_idx) = .TRUE. CASE('r') global_levels_rh(plevel_idx) = .TRUE. CASE('w') global_levels_w(plevel_idx) = .TRUE. END SELECT !!!!!! ! End of section for 3d variables of interest ! Filter for 2d variables of interest with level type ! of "surface" ELSE IF ( TRIM(gribmsg_leveltype) == 'surface' ) THEN SELECT CASE(gribmsg_shortname) CASE('sp') global_2dvar_ps = .TRUE. CASE('sdwe') global_2dvar_sd = .TRUE. CASE('tp') global_2dvar_lsprec = .TRUE. CASE('acpcp') global_2dvar_convprec = .TRUE. CASE('orog') global_2dvar_oro = .TRUE. CASE('lsm') global_2dvar_lsm = .TRUE. CASE('hpbl') global_2dvar_hmix = .TRUE. END SELECT !!!!!!!! End of section for 2d "surface" variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'heightAboveGround' ) THEN SELECT CASE(gribmsg_shortname) CASE('10u') global_2dvar_u10 = .TRUE. CASE('10v') global_2dvar_v10 = .TRUE. CASE('2t') global_2dvar_t2 = .TRUE. CASE('2d') global_2dvar_td2 = .TRUE. ! 2m RH is different depending on version of gribapi - in 1.14 it's just r, so ! one needs to also make sure that level is 2, but in later gribapi it seems to be ! 2r CASE('r') IF (gribmsg_level .EQ. 2) global_2dvar_rh2 = .TRUE. CASE('2r') global_2dvar_rh2 = .TRUE. END SELECT !!!!!!!! End of section for 2d "heightAboveGround" !!!!!!!! variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'meanSea' ) THEN SELECT CASE(gribmsg_shortname) CASE('prmsl') global_2dvar_msl = .TRUE. END SELECT !!!!!!!! End of section for 2d "meanSea" !!!!!!!! variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'Low cloud layer' ) THEN SELECT CASE(gribmsg_shortname) CASE('tcc') global_2dvar_tcc = .TRUE. END SELECT !!!!!!!! End of section for 2d "Low cloud layer" !!!!!!!! variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'sigma' ) THEN SELECT CASE(gribmsg_shortname) CASE('t') global_2dvar_tsig1 = .TRUE. CASE('u') global_2dvar_usig1 = .TRUE. CASE('v') global_2dvar_vsig1 = .TRUE. END SELECT !!!!!!!! End of section for 2d "sigma" !!!!!!!! variables of interest END IF ! Selection of message types (3d, surface, etc.) END IF END DO ! End of main loop for reading each GRIB message CALL grib_close_file(ifile) ! Check that all the 3d levels we flagged in the above loop resulted ! in all flags being true ! This function will have the side effect of printing any ! problems to stdout. If false, we force exit IF (ncep_all_levels_present()) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR CALL delete_levels_flags RETURN ENDIF ! The w (vertical wind) variable has to be handled differently because it ! is only available in the NCEP files for 100 mb and higher pressure levels IF (ncep_all_expected_w_levels_present(expected_num_levels, & & millibars)) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR WRITE(6,*) 'ERROR: Not all expected w levels are present' CALL delete_levels_flags RETURN ENDIF IF (ncep_all_2dvars_present()) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! Be sure to deallocate our allocatable arrays DEALLOCATE(millibars) CALL delete_levels_flags RETURN END FUNCTION cgutils_check_ncep_grib2 INTEGER FUNCTION cgutils_check_fv3_grib2(filepath, expected_num_levels, & & check_fcst, check_anl, check_grib2) & & RESULT(return_code) USE grib_api IMPLICIT NONE CHARACTER(LEN=256), INTENT(IN) :: filepath INTEGER, INTENT(IN) :: expected_num_levels LOGICAL, INTENT(IN) :: check_fcst, check_anl, check_grib2 ! return_code: possible values defined in module header, RETURN_CODE_* INTEGER :: ifile, igrib, iret, grib_centre CHARACTER(LEN=GRIB_ATTRIBUTE_STR_LEN) :: gribmsg_leveltype, & & gribmsg_shortname INTEGER :: gribmsg_level, plevel_idx LOGICAL :: end_of_file INTEGER, ALLOCATABLE, DIMENSION(:) :: millibars ! Start with return_code 0. Anything wrong will change it return_code = RETURN_CODE_NORMAL CALL grib_open_file(ifile, filepath, 'r', iret) IF (iret == 0) THEN ! Use first record to detect centre, which is assumed constant ! amongst all messages CALL grib_new_from_file(ifile, igrib, iret) CALL grib_get(igrib, 'centre', grib_centre) CALL grib_close_file(ifile) ELSE WRITE(0,*) 'ERROR: problem opening GRIB file: '//TRIM(filepath) return_code = RETURN_CODE_UTILITY_ERROR RETURN END IF IF (grib_centre .NE. GRIB_CENTRE_NCEP) THEN WRITE(6,*) 'ERROR: Unexpected grib_centre: ', grib_centre return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! Pass through the file once in order to build an inventory of ! isobaric pressure levels so that later we can make sure that ! each 3D variable has each of these pressure levels. While ! we're doing the inventory we can check that the number of ! isobaric levels found is indeed what we expect. CALL grib_open_file(ifile, filepath, 'r', iret) IF (iret == 0) THEN ALLOCATE(millibars(expected_num_levels)) CALL build_fv3_millibars_array(ifile, expected_num_levels, & & millibars, return_code) CALL grib_close_file(ifile) ELSE WRITE(0,*) 'ERROR: problem opening GRIB file: '//TRIM(filepath) return_code = RETURN_CODE_UTILITY_ERROR RETURN END IF ! Check that build_millibars_array() had no errors ! We will assume that any error messages were written from the ! subroutine IF (return_code .NE. RETURN_CODE_NORMAL) THEN return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! At this point, millibars should be an array of all the pressure ! levels (isobaricInhPa) that we found. !PRINT *, 'millibars: ', millibars CALL allocate_and_init_flags(expected_num_levels) ! Now, make our primary checking pass through the GRIB file. ! For 3D isobaric-level variables, as we read a level, we will ! find its index in millibars array, and use that to index the ! correct flag element in our logical array for each 3D variable ! Iterate through the entire file, keeping track of what we've found end_of_file = .FALSE. CALL grib_open_file(ifile, filepath, 'r', iret) DO WHILE (.NOT. end_of_file) CALL grib_new_from_file(ifile, igrib, iret) IF (iret .eq. GRIB_END_OF_FILE) THEN end_of_file = .TRUE. ELSE CALL grib_get(igrib, 'shortName', gribmsg_shortname) CALL grib_get(igrib, 'typeOfLevel', gribmsg_leveltype) CALL grib_get(igrib, 'level', gribmsg_level) ! Check general attributes of this message. Note that ! all messages will be checked, even those we might not ! care about. As of 2019-01-17, the checks are for GRIB2 ! and for fcst/anl if specified by "check_anl" or "check_fcst" IF (.NOT. message_is_good(igrib, & & check_anl, check_fcst, check_grib2)) THEN CALL write_err(gribmsg_shortname, gribmsg_level) return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF !PRINT *, 'shortname, leveltype, level, plevel_idx: ', & !& gribmsg_shortname, gribmsg_leveltype, & !& gribmsg_level, plevel_idx ! Filter for 3D variables of interest IF ( TRIM(gribmsg_leveltype) == 'isobaricInhPa' ) THEN CALL grib_get(igrib, 'level', gribmsg_level) ! Unsorted array "millibars," created above, tells us what index we ! should use to access the retrieved level flag as we maintain an ! inventory of what we found, and what we didn't ! Use my function myfindloc() to find the index for a specified pressure ! level. ! This is annoyingly "special case" code. Because of new FV3 files, "partial" ! levels are present in the GRIB files. Other routines compute the ! list of "valid" pressure levels, but as we go through and process ! each message, as we're doing here, if it's on one of the invalid ! pressure levels then we simply want to ignore it. If it's a valid ! pressure level, and it's a variable we're interested in, we mark it ! as present on that pressure level plevel_idx = myfindloc(expected_num_levels, millibars, gribmsg_level) IF (plevel_idx .LE. 0) THEN !CALL write_err(gribmsg_shortname, gribmsg_level) !return_code = RETURN_CODE_GRIBFILE_ERROR !RETURN CONTINUE ELSE !PRINT *, 'shortname, leveltype, level, plevel_idx: ', & !& gribmsg_shortname, gribmsg_leveltype, & !& gribmsg_level, plevel_idx SELECT CASE(gribmsg_shortname) CASE('u') global_levels_u(plevel_idx) = .TRUE. CASE('v') global_levels_v(plevel_idx) = .TRUE. CASE('t') global_levels_t(plevel_idx) = .TRUE. CASE('r') global_levels_rh(plevel_idx) = .TRUE. CASE('w') global_levels_w(plevel_idx) = .TRUE. END SELECT ENDIF !!!!!! ! End of section for 3d variables of interest ! Filter for 2d variables of interest with level type ! of "surface" ELSE IF ( TRIM(gribmsg_leveltype) == 'surface' ) THEN SELECT CASE(gribmsg_shortname) CASE('sp') global_2dvar_ps = .TRUE. CASE('sdwe') global_2dvar_sd = .TRUE. CASE('tp') global_2dvar_lsprec = .TRUE. CASE('acpcp') global_2dvar_convprec = .TRUE. CASE('orog') global_2dvar_oro = .TRUE. CASE('lsm') global_2dvar_lsm = .TRUE. CASE('hpbl') global_2dvar_hmix = .TRUE. END SELECT !!!!!!!! End of section for 2d "surface" variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'heightAboveGround' ) THEN SELECT CASE(gribmsg_shortname) CASE('10u') global_2dvar_u10 = .TRUE. CASE('10v') global_2dvar_v10 = .TRUE. CASE('2t') global_2dvar_t2 = .TRUE. CASE('2d') global_2dvar_td2 = .TRUE. ! 2m RH is different depending on version of gribapi - in 1.14 it's just r, so ! one needs to also make sure that level is 2, but in later gribapi it seems to be ! 2r CASE('r') IF (gribmsg_level .EQ. 2) global_2dvar_rh2 = .TRUE. CASE('2r') global_2dvar_rh2 = .TRUE. END SELECT !!!!!!!! End of section for 2d "heightAboveGround" !!!!!!!! variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'meanSea' ) THEN SELECT CASE(gribmsg_shortname) CASE('prmsl') global_2dvar_msl = .TRUE. END SELECT !!!!!!!! End of section for 2d "meanSea" !!!!!!!! variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'Low cloud layer' ) THEN SELECT CASE(gribmsg_shortname) CASE('tcc') global_2dvar_tcc = .TRUE. END SELECT !!!!!!!! End of section for 2d "Low cloud layer" !!!!!!!! variables of interest ELSE IF ( TRIM(gribmsg_leveltype) == 'sigma' ) THEN SELECT CASE(gribmsg_shortname) CASE('t') global_2dvar_tsig1 = .TRUE. CASE('u') global_2dvar_usig1 = .TRUE. CASE('v') global_2dvar_vsig1 = .TRUE. END SELECT !!!!!!!! End of section for 2d "sigma" !!!!!!!! variables of interest END IF ! Selection of message types (3d, surface, etc.) END IF END DO ! End of main loop for reading each GRIB message CALL grib_close_file(ifile) ! Check that all the 3d levels we flagged in the above loop resulted ! in all flags being true ! This function will have the side effect of printing any ! problems to stdout. If false, we force exit IF (ncep_all_levels_present()) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR CALL delete_levels_flags RETURN ENDIF ! The w (vertical wind) variable has to be handled differently because it ! is only available in the NCEP files for 100 mb and higher pressure levels IF (fv3_all_expected_w_levels_present(expected_num_levels, & & millibars)) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR WRITE(6,*) 'ERROR: Not all expected w levels are present' CALL delete_levels_flags RETURN ENDIF IF (fv3_all_2dvars_present()) THEN return_code = RETURN_CODE_NORMAL ELSE return_code = RETURN_CODE_GRIBFILE_ERROR RETURN ENDIF ! Be sure to deallocate our allocatable arrays DEALLOCATE(millibars) CALL delete_levels_flags RETURN END FUNCTION cgutils_check_fv3_grib2 END MODULE cgutils