diff --git a/CMakeLists.txt b/CMakeLists.txt index 43a83f75..2af54741 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,7 +9,7 @@ file(STRINGS "VERSION" pVersion) project( w3emc VERSION ${pVersion} - LANGUAGES C Fortran) + LANGUAGES Fortran) # Handle user options. option(ENABLE_DOCS "Enable generation of doxygen-based documentation." OFF) @@ -24,12 +24,6 @@ if(NOT CMAKE_BUILD_TYPE MATCHES "^(Debug|Release|RelWithDebInfo|MinSizeRel)$") set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo") endif() -# Find required packages. -find_package(nemsio 2.5.0 CONFIG REQUIRED) -find_package(sigio 2.3.0 CONFIG REQUIRED) - -find_package(NetCDF 4.3.3 REQUIRED Fortran) - # Compile the source code in the src directory. add_subdirectory(src) diff --git a/cmake/FindNetCDF.cmake b/cmake/FindNetCDF.cmake deleted file mode 100644 index 9e32378f..00000000 --- a/cmake/FindNetCDF.cmake +++ /dev/null @@ -1,347 +0,0 @@ -# (C) Copyright 2011- ECMWF. -# -# This software is licensed under the terms of the Apache Licence Version 2.0 -# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. -# In applying this licence, ECMWF does not waive the privileges and immunities -# granted to it by virtue of its status as an intergovernmental organisation nor -# does it submit to any jurisdiction. - -# Try to find NetCDF includes and library. -# Supports static and shared libaries and allows each component to be found in sepearte prefixes. -# -# This module defines -# -# - NetCDF_FOUND - System has NetCDF -# - NetCDF_INCLUDE_DIRS - the NetCDF include directories -# - NetCDF_VERSION - the version of NetCDF -# - NetCDF_CONFIG_EXECUTABLE - the netcdf-config executable if found -# - NetCDF_PARALLEL - Boolean True if NetCDF4 has parallel IO support via hdf5 and/or pnetcdf -# - NetCDF_HAS_PNETCDF - Boolean True if NetCDF4 has pnetcdf support -# -# Deprecated Defines -# - NetCDF_LIBRARIES - [Deprecated] Use NetCDF::NetCDF_ targets instead. -# -# -# Following components are available: -# -# - C - C interface to NetCDF (netcdf) -# - CXX - CXX4 interface to NetCDF (netcdf_c++4) -# - Fortran - Fortran interface to NetCDF (netcdff) -# -# For each component the following are defined: -# -# - NetCDF__FOUND - whether the component is found -# - NetCDF__LIBRARIES - the libraries for the component -# - NetCDF__LIBRARY_SHARED - Boolean is true if libraries for component are shared -# - NetCDF__INCLUDE_DIRS - the include directories for specified component -# - NetCDF::NetCDF_ - target of component to be used with target_link_libraries() -# -# The following paths will be searched in order if set in CMake (first priority) or environment (second priority) -# -# - NetCDF_ROOT - root of NetCDF installation -# - NetCDF_PATH - root of NetCDF installation -# -# The search process begins with locating NetCDF Include headers. If these are in a non-standard location, -# set one of the following CMake or environment variables to point to the location: -# -# - NetCDF_INCLUDE_DIR or NetCDF_${comp}_INCLUDE_DIR -# - NetCDF_INCLUDE_DIRS or NetCDF_${comp}_INCLUDE_DIR -# -# Notes: -# -# - Use "NetCDF::NetCDF_" targets only. NetCDF_LIBRARIES exists for backwards compatibility and should not be used. -# - These targets have all the knowledge of include directories and library search directories, and a single -# call to target_link_libraries will provide all these transitive properties to your target. Normally all that is -# needed to build and link against NetCDF is, e.g.: -# target_link_libraries(my_c_tgt PUBLIC NetCDF::NetCDF_C) -# - "NetCDF" is always the preferred naming for this package, its targets, variables, and environment variables -# - For compatibility, some variables are also set/checked using alternate names NetCDF4, NETCDF, or NETCDF4 -# - Environments relying on these older environment variable names should move to using a "NetCDF_ROOT" environment variable -# - Preferred component capitalization follows the CMake LANGUAGES variables: i.e., C, Fortran, CXX -# - For compatibility, alternate capitalizations are supported but should not be used. -# - If no components are defined, all components will be searched -# - -list( APPEND _possible_components C CXX Fortran ) - -## Include names for each component -set( NetCDF_C_INCLUDE_NAME netcdf.h ) -set( NetCDF_CXX_INCLUDE_NAME netcdf ) -set( NetCDF_Fortran_INCLUDE_NAME netcdf.mod ) - -## Library names for each component -set( NetCDF_C_LIBRARY_NAME netcdf ) -set( NetCDF_CXX_LIBRARY_NAME netcdf_c++4 ) -set( NetCDF_Fortran_LIBRARY_NAME netcdff ) - -## Enumerate search components -foreach( _comp ${_possible_components} ) - string( TOUPPER "${_comp}" _COMP ) - set( _arg_${_COMP} ${_comp} ) - set( _name_${_COMP} ${_comp} ) -endforeach() - -set( _search_components C) -foreach( _comp ${${CMAKE_FIND_PACKAGE_NAME}_FIND_COMPONENTS} ) - string( TOUPPER "${_comp}" _COMP ) - set( _arg_${_COMP} ${_comp} ) - list( APPEND _search_components ${_name_${_COMP}} ) - if( NOT _name_${_COMP} ) - message(SEND_ERROR "Find${CMAKE_FIND_PACKAGE_NAME}: COMPONENT ${_comp} is not a valid component. Valid components: ${_possible_components}" ) - endif() -endforeach() -list( REMOVE_DUPLICATES _search_components ) - -## Search hints for finding include directories and libraries -foreach( _comp IN ITEMS "_" "_C_" "_Fortran_" "_CXX_" ) - foreach( _name IN ITEMS NetCDF4 NetCDF NETCDF4 NETCDF ) - foreach( _var IN ITEMS ROOT PATH ) - list(APPEND _search_hints ${${_name}${_comp}${_var}} $ENV{${_name}${_comp}${_var}} ) - list(APPEND _include_search_hints - ${${_name}${_comp}INCLUDE_DIR} $ENV{${_name}${_comp}INCLUDE_DIR} - ${${_name}${_comp}INCLUDE_DIRS} $ENV{${_name}${_comp}INCLUDE_DIRS} ) - endforeach() - endforeach() -endforeach() -#Old-school HPC module env variable names -foreach( _name IN ITEMS NetCDF4 NetCDF NETCDF4 NETCDF ) - foreach( _comp IN ITEMS "_C" "_Fortran" "_CXX" ) - list(APPEND _search_hints ${${_name}} $ENV{${_name}}) - list(APPEND _search_hints ${${_name}${_comp}} $ENV{${_name}${_comp}}) - endforeach() -endforeach() - -## Find headers for each component -set(NetCDF_INCLUDE_DIRS) -set(_new_search_components) -foreach( _comp IN LISTS _search_components ) - if(NOT ${PROJECT_NAME}_NetCDF_${_comp}_FOUND) - list(APPEND _new_search_components ${_comp}) - endif() - find_file(NetCDF_${_comp}_INCLUDE_FILE - NAMES ${NetCDF_${_comp}_INCLUDE_NAME} - DOC "NetCDF ${_comp} include directory" - HINTS ${_include_search_hints} ${_search_hints} - PATH_SUFFIXES include include/netcdf - ) - mark_as_advanced(NetCDF_${_comp}_INCLUDE_FILE) - message(DEBUG "NetCDF_${_comp}_INCLUDE_FILE: ${NetCDF_${_comp}_INCLUDE_FILE}") - if( NetCDF_${_comp}_INCLUDE_FILE ) - get_filename_component(NetCDF_${_comp}_INCLUDE_FILE ${NetCDF_${_comp}_INCLUDE_FILE} ABSOLUTE) - get_filename_component(NetCDF_${_comp}_INCLUDE_DIR ${NetCDF_${_comp}_INCLUDE_FILE} DIRECTORY) - list(APPEND NetCDF_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIR}) - endif() -endforeach() -if(NetCDF_INCLUDE_DIRS) - list(REMOVE_DUPLICATES NetCDF_INCLUDE_DIRS) -endif() -set(NetCDF_INCLUDE_DIRS "${NetCDF_INCLUDE_DIRS}" CACHE STRING "NetCDF Include directory paths" FORCE) - -## Find n*-config executables for search components -foreach( _comp IN LISTS _search_components ) - if( _comp MATCHES "^(C)$" ) - set(_conf "c") - elseif( _comp MATCHES "^(Fortran)$" ) - set(_conf "f") - elseif( _comp MATCHES "^(CXX)$" ) - set(_conf "cxx4") - endif() - find_program( NetCDF_${_comp}_CONFIG_EXECUTABLE - NAMES n${_conf}-config - HINTS ${NetCDF_INCLUDE_DIRS} ${_include_search_hints} ${_search_hints} - PATH_SUFFIXES bin Bin ../bin ../../bin - DOC "NetCDF n${_conf}-config helper" ) - message(DEBUG "NetCDF_${_comp}_CONFIG_EXECUTABLE: ${NetCDF_${_comp}_CONFIG_EXECUTABLE}") -endforeach() - -set(_C_libs_flag --libs) -set(_Fortran_libs_flag --flibs) -set(_CXX_libs_flag --libs) -set(_C_includes_flag --includedir) -set(_Fortran_includes_flag --includedir) -set(_CXX_includes_flag --includedir) -function(netcdf_config exec flag output_var) - set(${output_var} False PARENT_SCOPE) - if( exec ) - execute_process( COMMAND ${exec} ${flag} RESULT_VARIABLE _ret OUTPUT_VARIABLE _val) - if( _ret EQUAL 0 ) - string( STRIP ${_val} _val ) - set( ${output_var} ${_val} PARENT_SCOPE ) - endif() - endif() -endfunction() - -## Detect additional package properties -netcdf_config(${NetCDF_C_CONFIG_EXECUTABLE} --has-parallel4 _val) -if( NOT _val MATCHES "^(yes|no)$" ) - netcdf_config(${NetCDF_C_CONFIG_EXECUTABLE} --has-parallel _val) -endif() -if( _val MATCHES "^(yes)$" ) - set(NetCDF_PARALLEL TRUE CACHE STRING "NetCDF has parallel IO capability via pnetcdf or hdf5." FORCE) -else() - set(NetCDF_PARALLEL FALSE CACHE STRING "NetCDF has no parallel IO capability." FORCE) -endif() - -if(NetCDF_PARALLEL) - find_package(MPI REQUIRED) -endif() - -## Find libraries for each component -set( NetCDF_LIBRARIES ) -foreach( _comp IN LISTS _search_components ) - string( TOUPPER "${_comp}" _COMP ) - - find_library( NetCDF_${_comp}_LIBRARY - NAMES ${NetCDF_${_comp}_LIBRARY_NAME} - DOC "NetCDF ${_comp} library" - HINTS ${NetCDF_${_comp}_INCLUDE_DIRS} ${_search_hints} - PATH_SUFFIXES lib64 lib ../lib64 ../lib ../../lib64 ../../lib ) - mark_as_advanced( NetCDF_${_comp}_LIBRARY ) - get_filename_component(NetCDF_${_comp}_LIBRARY ${NetCDF_${_comp}_LIBRARY} ABSOLUTE) - set(NetCDF_${_comp}_LIBRARY ${NetCDF_${_comp}_LIBRARY} CACHE STRING "NetCDF ${_comp} library" FORCE) - message(DEBUG "NetCDF_${_comp}_LIBRARY: ${NetCDF_${_comp}_LIBRARY}") - - if( NetCDF_${_comp}_LIBRARY ) - if( NetCDF_${_comp}_LIBRARY MATCHES ".a$" ) - set( NetCDF_${_comp}_LIBRARY_SHARED FALSE ) - set( _library_type STATIC) - else() - list( APPEND NetCDF_LIBRARIES ${NetCDF_${_comp}_LIBRARY} ) - set( NetCDF_${_comp}_LIBRARY_SHARED TRUE ) - set( _library_type SHARED) - endif() - endif() - - #Use nc-config to set per-component LIBRARIES variable if possible - netcdf_config( ${NetCDF_${_comp}_CONFIG_EXECUTABLE} ${_${_comp}_libs_flag} _val ) - if( _val ) - set( NetCDF_${_comp}_LIBRARIES ${_val} ) - if(NOT NetCDF_${_comp}_LIBRARY_SHARED AND NOT NetCDF_${_comp}_FOUND) #Static targets should use nc_config to get a proper link line with all necessary static targets. - list( APPEND NetCDF_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) - endif() - else() - set( NetCDF_${_comp}_LIBRARIES ${NetCDF_${_comp}_LIBRARY} ) - if(NOT NetCDF_${_comp}_LIBRARY_SHARED) - message(SEND_ERROR "Unable to properly find NetCDF. Found static libraries at: ${NetCDF_${_comp}_LIBRARY} but could not run nc-config: ${NetCDF_CONFIG_EXECUTABLE}") - endif() - endif() - - #Use nc-config to set per-component INCLUDE_DIRS variable if possible - netcdf_config( ${NetCDF_${_comp}_CONFIG_EXECUTABLE} ${_${_comp}_includes_flag} _val ) - if( _val ) - string( REPLACE " " ";" _val ${_val} ) - set( NetCDF_${_comp}_INCLUDE_DIRS ${_val} ) - else() - set( NetCDF_${_comp}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIR} ) - endif() - - if( NetCDF_${_comp}_LIBRARIES AND NetCDF_${_comp}_INCLUDE_DIRS ) - set( ${CMAKE_FIND_PACKAGE_NAME}_${_arg_${_COMP}}_FOUND TRUE ) - if (NOT TARGET NetCDF::NetCDF_${_comp}) - add_library(NetCDF::NetCDF_${_comp} ${_library_type} IMPORTED) - set_target_properties(NetCDF::NetCDF_${_comp} PROPERTIES - IMPORTED_LOCATION ${NetCDF_${_comp}_LIBRARY} - INTERFACE_INCLUDE_DIRECTORIES "${NetCDF_${_comp}_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) - if( NOT _comp MATCHES "^(C)$" ) - target_link_libraries(NetCDF::NetCDF_${_comp} INTERFACE NetCDF::NetCDF_C) - endif() - if(MPI_${_comp}_FOUND) - target_link_libraries(NetCDF::NetCDF_${_comp} INTERFACE MPI::MPI_${_comp}) - endif() - endif() - endif() -endforeach() -if(NetCDF_LIBRARIES AND NetCDF_${_comp}_LIBRARY_SHARED) - list(REMOVE_DUPLICATES NetCDF_LIBRARIES) -endif() -set(NetCDF_LIBRARIES "${NetCDF_LIBRARIES}" CACHE STRING "NetCDF library targets" FORCE) - -## Find version via netcdf-config if possible -if (NetCDF_INCLUDE_DIRS) - if( NetCDF_C_CONFIG_EXECUTABLE ) - netcdf_config( ${NetCDF_C_CONFIG_EXECUTABLE} --version _vers ) - if( _vers ) - string(REGEX REPLACE ".* ((([0-9]+)\\.)+([0-9]+)).*" "\\1" NetCDF_VERSION "${_vers}" ) - endif() - else() - foreach( _dir IN LISTS NetCDF_INCLUDE_DIRS) - if( EXISTS "${_dir}/netcdf_meta.h" ) - file(STRINGS "${_dir}/netcdf_meta.h" _netcdf_version_lines - REGEX "#define[ \t]+NC_VERSION_(MAJOR|MINOR|PATCH|NOTE)") - string(REGEX REPLACE ".*NC_VERSION_MAJOR *\([0-9]*\).*" "\\1" _netcdf_version_major "${_netcdf_version_lines}") - string(REGEX REPLACE ".*NC_VERSION_MINOR *\([0-9]*\).*" "\\1" _netcdf_version_minor "${_netcdf_version_lines}") - string(REGEX REPLACE ".*NC_VERSION_PATCH *\([0-9]*\).*" "\\1" _netcdf_version_patch "${_netcdf_version_lines}") - string(REGEX REPLACE ".*NC_VERSION_NOTE *\"\([^\"]*\)\".*" "\\1" _netcdf_version_note "${_netcdf_version_lines}") - set(NetCDF_VERSION "${_netcdf_version_major}.${_netcdf_version_minor}.${_netcdf_version_patch}${_netcdf_version_note}") - unset(_netcdf_version_major) - unset(_netcdf_version_minor) - unset(_netcdf_version_patch) - unset(_netcdf_version_note) - unset(_netcdf_version_lines) - endif() - endforeach() - endif() -endif () - -## Finalize find_package -include(FindPackageHandleStandardArgs) - -if(NOT NetCDF_FOUND OR _new_search_components) - find_package_handle_standard_args( ${CMAKE_FIND_PACKAGE_NAME} - REQUIRED_VARS NetCDF_INCLUDE_DIRS NetCDF_LIBRARIES - VERSION_VAR NetCDF_VERSION - HANDLE_COMPONENTS ) -endif() - -foreach( _comp IN LISTS _search_components ) - if( NetCDF_${_comp}_FOUND ) - #Record found components to avoid duplication in NetCDF_LIBRARIES for static libraries - set(NetCDF_${_comp}_FOUND ${NetCDF_${_comp}_FOUND} CACHE BOOL "NetCDF ${_comp} Found" FORCE) - #Set a per-package, per-component found variable to communicate between multiple calls to find_package() - set(${PROJECT_NAME}_NetCDF_${_comp}_FOUND True) - endif() -endforeach() - -if( ${CMAKE_FIND_PACKAGE_NAME}_FOUND AND NOT ${CMAKE_FIND_PACKAGE_NAME}_FIND_QUIETLY AND _new_search_components) - message( STATUS "Find${CMAKE_FIND_PACKAGE_NAME} defines targets:" ) - message( STATUS " - NetCDF_VERSION [${NetCDF_VERSION}]") - message( STATUS " - NetCDF_PARALLEL [${NetCDF_PARALLEL}]") - foreach( _comp IN LISTS _new_search_components ) - string( TOUPPER "${_comp}" _COMP ) - message( STATUS " - NetCDF_${_comp}_CONFIG_EXECUTABLE [${NetCDF_${_comp}_CONFIG_EXECUTABLE}]") - if( ${CMAKE_FIND_PACKAGE_NAME}_${_arg_${_COMP}}_FOUND ) - get_filename_component(_root ${NetCDF_${_comp}_INCLUDE_DIR}/.. ABSOLUTE) - if( NetCDF_${_comp}_LIBRARY_SHARED ) - message( STATUS " - NetCDF::NetCDF_${_comp} [SHARED] [Root: ${_root}] Lib: ${NetCDF_${_comp}_LIBRARY} ") - else() - message( STATUS " - NetCDF::NetCDF_${_comp} [STATIC] [Root: ${_root}] Lib: ${NetCDF_${_comp}_LIBRARY} ") - endif() - endif() - endforeach() -endif() - -foreach( _prefix NetCDF NetCDF4 NETCDF NETCDF4 ${CMAKE_FIND_PACKAGE_NAME} ) - set( ${_prefix}_INCLUDE_DIRS ${NetCDF_INCLUDE_DIRS} ) - set( ${_prefix}_LIBRARIES ${NetCDF_LIBRARIES}) - set( ${_prefix}_VERSION ${NetCDF_VERSION} ) - set( ${_prefix}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_FOUND} ) - set( ${_prefix}_CONFIG_EXECUTABLE ${NetCDF_CONFIG_EXECUTABLE} ) - set( ${_prefix}_PARALLEL ${NetCDF_PARALLEL} ) - - foreach( _comp ${_search_components} ) - string( TOUPPER "${_comp}" _COMP ) - set( _arg_comp ${_arg_${_COMP}} ) - set( ${_prefix}_${_comp}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_${_arg_comp}_FOUND} ) - set( ${_prefix}_${_COMP}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_${_arg_comp}_FOUND} ) - set( ${_prefix}_${_arg_comp}_FOUND ${${CMAKE_FIND_PACKAGE_NAME}_${_arg_comp}_FOUND} ) - - set( ${_prefix}_${_comp}_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) - set( ${_prefix}_${_COMP}_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) - set( ${_prefix}_${_arg_comp}_LIBRARIES ${NetCDF_${_comp}_LIBRARIES} ) - - set( ${_prefix}_${_comp}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIRS} ) - set( ${_prefix}_${_COMP}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIRS} ) - set( ${_prefix}_${_arg_comp}_INCLUDE_DIRS ${NetCDF_${_comp}_INCLUDE_DIRS} ) - endforeach() -endforeach() diff --git a/cmake/PackageConfig.cmake.in b/cmake/PackageConfig.cmake.in index ac576f6d..e4fda558 100644 --- a/cmake/PackageConfig.cmake.in +++ b/cmake/PackageConfig.cmake.in @@ -11,10 +11,6 @@ include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@-targets.cmake") include(CMakeFindDependencyMacro) -find_dependency(nemsio CONFIG) -find_dependency(sigio CONFIG) -find_dependency(NetCDF COMPONENTS Fortran) - # Get the build type from real32 library target get_target_property(@PROJECT_NAME@_BUILD_TYPES @PROJECT_NAME@::@PROJECT_NAME@_4 IMPORTED_CONFIGURATIONS) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4e152310..a44d9fe3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,19 +1,5 @@ include("list_of_files.cmake") -if(APPLE) - add_compile_definitions(APPLE) -elseif(UNIX) - add_compile_definitions(LINUX) -endif() - -if(CMAKE_C_COMPILER_ID MATCHES "^(Intel)$") - set(CMAKE_C_FLAGS "-g ${CMAKE_C_FLAGS}") - set(CMAKE_C_FLAGS_RELEASE "-O3") -elseif(CMAKE_C_COMPILER_ID MATCHES "^(GNU|Clang|AppleClang)$") - set(CMAKE_C_FLAGS "-ggdb ${CMAKE_C_FLAGS}") - set(CMAKE_C_FLAGS_RELEASE "-O3") -endif() - if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") set(CMAKE_Fortran_FLAGS "-g -traceback -fixed ${CMAKE_Fortran_FLAGS}") set(CMAKE_Fortran_FLAGS_RELEASE "-O2") @@ -31,7 +17,6 @@ endif() if(${CMAKE_Fortran_COMPILER_ID} MATCHES "^(GNU)$" AND ${CMAKE_Fortran_COMPILER_VERSION} VERSION_GREATER_EQUAL 10) set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -w -fallow-argument-mismatch -fallow-invalid-boz") endif() -set_source_files_properties(${c_src} PROPERTIES COMPILE_OPTIONS "${c_flags}") set(kinds "4" "8" "d") foreach(kind ${kinds}) @@ -45,16 +30,10 @@ foreach(kind ${kinds}) "${fortran_${kind}_flags}") set_target_properties(${lib_name} PROPERTIES Fortran_MODULE_DIRECTORY ${module_dir}) -add_library(${lib_name}_c OBJECT ${c_src}) - -# add_library(${lib_name} STATIC $ -# $) target_include_directories(${lib_name} INTERFACE $ $) - target_link_libraries(${lib_name} PRIVATE nemsio::nemsio sigio::sigio - NetCDF::NetCDF_Fortran) list(APPEND LIB_TARGETS ${lib_name}) install(DIRECTORY ${module_dir} DESTINATION ${CMAKE_INSTALL_PREFIX}) diff --git a/src/gblevents.f b/src/gblevents.f deleted file mode 100644 index ad2e541e..00000000 --- a/src/gblevents.f +++ /dev/null @@ -1,3662 +0,0 @@ -C> @file -C -C> SUBPROGRAM: GBLEVENTS PRE/POST PROCESSING OF PREPBUFR EVENTS -C> PRGMMR: DENNIS KEYSER ORG: EMC DATE: 2017-02-22 -C> -C> RUNS IN TWO MODES: "PREVENTS" AND "POSTEVENTS". IN THE -C> PREVENTS MODE, PREPARES OBSERVATIONAL PREPBUFR REPORTS FOR -C> SUBSEQUENT QUALITY CONTROL AND ANALYSIS PROGRAMS. THIS IS DONE -C> THROUGH THE FOLLOWING: INTERPOLATION OF GLOBAL FIRST GUESS {FROM -C> EITHER SIGIO (SIGMA OR HYBRID) OR NEMSIO INPUT} TO PREPBUFR -C> OBSERVATION LOCATIONS WITH ENCODING OF FIRST GUESS VALUES INTO -C> PREPBUFR REPORTS; ENCODING OF "PREVENT" AND/OR "VIRTMP" EVENTS INTO -C> PREPBUFR REPORTS; IN CERTAIN CASES, ENCODING OF A DERIVED PMSL INTO -C> SURFACE PREPBUFR REPORTS; AND ENCODING OF OBSERVATION ERRORS FROM -C> THE ERROR SPECIFICATION FILE INTO PREPBUFR REPORTS. IN THE -C> POSTEVENTS MODE, AFTER ALL QUALITY CONTROL AND ANALYSIS PROGRAMS -C> HAVE RUN, INTERPOLATES THE GLOBAL ANALYSIS {FROM EITHER SIGIO -C> (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT} TO PREPBUFR OBSERVATION -C> LOCATIONS AND ENCODES THESE ANALYZED VALUES INTO PREPBUFR REPORTS. -C> THE REMAINDER OF THIS ABSTRACT APPLIES ONLY TO THE PREVENTS MODE. -C> THE "PREVENT" EVENT CAN CHANGE A QUALITY MARKER TO FLAG AN -C> OBSERVATION DATUM FOR NON-USE BY SUBSEQUENT QC AND ANALYSIS -C> PROGRAMS (FILTERING). EXAMPLES WHERE THIS SUBROUTINE WILL WRITE -C> AN EVENT TO FLAG A DATUM INCLUDE: THE OBSERVATION ERROR FOR THAT -C> DATUM IS READ IN AS MISSING IN THE INPUT ERROR FILE, THE DATUM -C> ITSELF VIOLATES A GROSS OR "SANITY" CHECK, OR THE OBSERVED -C> PRESSURE DATUM IS MORE THAN 100 MB BELOW THE GUESS SURFACE -C> PRESSURE. THE "VIRTMP" EVENT CAN CHANGE THE SPECIFIC HUMIDITY -C> OBSERVATION (RE-CALCULATED) AS WELL AS THE TEMPERATURE -C> OBSERVATION (FROM SENSIBLE TO VIRTUAL TEMPERATURE, BASED ON -C> JUST-CALCULATED SPECIFIC HUMIDITY). CURRENTLY THIS APPLIES ONLY -C> TO SURFACE (LAND, MARINE AND MESONET) DATA TYPES, POSSIBLY TO -C> RAOB, DROP AND MULTI-LEVEL RECCO DATA TYPES IF THE SWITCH -C> "ADPUPA_VIRT" IS TRUE (NORMALLY, HOWEVER IT IS FALSE) [OTHER DATA -C> TYPES WITH REPORTED SENSIBLE TEMPERATURE EITHER HAVE MISSING -C> MOISTURE (E.G., ALL AIRCRAFT TYPES EXCEPT FOR SOME ACARS, SATELLITE -C> WIND TYPES), FLAGGED MOISTURE (E.G., SOME ACARS) OR CALCULATE -C> SPECIFIC HUMIDITY/VIRTUAL TEMPERATURE IN SUBSEQUENT PROGRAMS (E.G., -C> RAOBS, DROPS AND MULTI-LEVEL RECCOS WHICH CALCULATE THESE IN -C> PROGRAM "CQCBUFR", IN WHICH CASE THE SWITCH "ADPUPA_VIRT" HERE MUST -C> BE FALSE!)]. FOR CASES WHERE THE SWITCH "DOBERR" IS FALSE, THE -C> OBSERVATION ERROR FOR ALL DATA REMAINS MISSING IN THE PREPBUFR -C> FILE. IN THIS CASE, THE INPUT ERROR FILE IS USUALLY A NULL FILE -C> AND THE "PREVENT" EVENT TO FLAG THE DATUM IS NOT INVOKED. FOR -C> CASES WHERE THE SWITCH "DOFCST" IS FALSE, IF THE SWITCH "SOME_FCST" -C> IS ALSO FALSE, THEN FORECAST VALUES ARE NOT ENCODED FOR ANY MESSAGE -C> TYPE; IF "SOME_FCST" IS TRUE THEN FORECAST VALUES ARE ENCODED, BUT -C> ONLY FOR REPORTS IN THOSE MESSAGE TYPES FOR WHICH A GUESS VALUE IS -C> NEEDED BY SUBSEQUENT QC PROGRAMS. IT SHOULD BE NOTED THAT THE -C> FILTERING OF DATA ASSOCIATED WITH THE "PREVENT" EVENT PROCESSING IS -C> NOT INVOKED IF ALL THREE ARE TRUE: DOBERR= FALSE, THE FORECAST -C> VALUES ARE MISSING (DOFCST=FALSE & SOME_FCST=TRUE & MESSAGE TYPE IS -C> NOT "ADPUPA", "AIRCFT", "AIRCAR", "PROFLR", OR "VADWND" -- OR -- -C> DOFCST=FALSE & SOME_FCST=FALSE), AND "VIRTMP" EVENT PROCESSING IS -C> NOT INVOKED (EITHER MESSAGE TYPE IS NOT "ADPSFC", "SFCSHP" OR -C> "MSONET" WHEN "ADPUPA_VIRT" IS FALSE, OR MESSAGE TYPE IS NOT -C> "ADPSFC", "SFCSHP", "MSONET" OR "ADPUPA" WHEN "ADPUPA_VIRT" IS -C> TRUE). ALSO, IF VIRTUAL TEMPERATURE PROCESSING IS PERFORMED, ALL -C> SURFACE REPORTS WITH MISSING PMSL WILL ENCODE A DERIVED PMSL INTO -C> PREPBUFR IF THE SWITCH DOPMSL IS TRUE AND A VIRTUAL TEMPERATURE WAS -C> SUCCESSFULLY CALCULATED. -C> -C> PROGRAM HISTORY LOG: -C> -1999-07-01 D. A. KEYSER -- ORIGINAL AUTHOR (ADAPTED FROM PREVENTS -C> SUBROUTINE IN PREPDATA PROGRAM, BUT NOW GENERALIZED FOR -C> POSTEVENTS MODE) -C> -1999-07-12 D. A. KEYSER -- MODIFIED TO INTERPOLATE MODEL SPECIFIC -C> HUMIDITY TO OBSERVATION LOCATION WHEN OBS. SPECIFIC HUMIDITY IS -C> MISSING AS LONG AS OBS. TEMPERATURE IS NON-MISSING -C> -1999-09-09 D. A. KEYSER -- ADDED "VADWND" TO THE LIST OF MESSAGE -C> TYPES FOR WHICH FORECAST VALUES MUST BE ENCODED, EVEN WHEN -C> DOFCST=FALSE (NECESSARY BECAUSE THE NEW PROGRAM CQCVAD NEEDS THE -C> BACKGROUND DATA) -C> -1999-09-09 D. A. KEYSER -- CHANGES TO MAKE CODE MORE PORTABLE; -C> 'TFC' NOW GENERATED FOR VADWND MESSAGE TYPES EVEN THOUGH TOB IS -C> MISSING (NEEDED BY CQCVAD PROGRAM) -C> -1999-12-01 D. A. KEYSER -- SPEC. HUMIDITY AND VIRT. TEMPERATURE ARE -C> NOW CALCULATED WHEN SPEC. HUMIDITY QUAL. MARKER IS BAD (SUBJECT -C> TO A SANITY CHECK), HOWEVER THE VIRT. TEMPERATURE GETS A BAD -C> QUAL. MARKER (8) -C> -2000-09-21 D. A. KEYSER -- THE PRESSURE LEVEL ABOVE WHICH ALL SPEC. -C> HUMIDITY QUAL. MARKERS ARE "REJECTED" (Q.M. SET TO 9) IS NOW READ -C> IN AS A N-LIST SWITCH (QTOP_REJ), BEFORE IT WAS HARDWIRED TO 300 -C> MB -C> -2000-12-13 D. A. KEYSER -- WILL NO LONGER PERFORM VIRTUAL TEMPERATURE -C> PROCESSING FOR ACARS DATA SINCE MOISTURE IS FLAGGED RIGHT NOW -C> (ACARS MOISTURE ONLY WRITTEN INTO PREPBUFR FILE FOR STATISTICAL -C> REASONS) -C> -2001-02-02 D. A. KEYSER -- RESTORED LEGACY LOGIC TO FLAG CERTAIN -C> SATELLITE TEMPERATURE SOUNDINGS EITHER BELOW 100 MB (TEMP. OBS) -C> OR ON ALL LEVELS (SPEC. HUM. OBS), CONTROLLED BY NEW NAMELIST -C> SWITCH "SATMQC" -C> -2001-09-27 D. A. KEYSER -- 'TFC' AND 'QFC' NOW GENERATED FOR REPORT -C> TYPE 111 (SYNDAT REPORTS AT STORM CENTER) EVEN THOUGH "TOB" AND -C> "QOB" ARE MISSING (NEEDED BY SYNDATA PROGRAM); IN PREPARATION FOR -C> CHANGE FROM T170L42 TO T254L64 SGES, NOW MAKES COEFFICIENT ARRAYS -C> ALLOCATABLE TO ALLOW THEM TO OBTAIN MEMORY FROM "HEAP" RATHER -C> THAN FROM "STACK", ALSO HAVE INCREASED THE MAX NUMBER OF LEVELS -C> IN ARRAYS FROM 42 TO 64, FINALLY ALSO NO LONGER STOPS WITH C. -C> CODE 70 IF EVEN NUMBER OF LONGITUDES IN SIGMA GUESS (IMAX, -C> HARDWIRED TO 384) IS .LT. SPECTRAL RESOLUTION (JCAP) * 2 -C> -2001-10-10 D. A. KEYSER -- AT PREPBUFR CENTER DATES WITH AN HOUR THAT -C> IS NOT A MULTIPLE OF 3 (WHEN A GLOBAL SIGMA GUESS/ANAL FILE IS -C> NOT AVAILABLE; E.G., IN RUC2A RUNS) NOW PERFORMS A LINEAR -C> INTERPOLATION BETWEEN SPECTRAL COEFFICIENTS IN 2 SPANNING SIGMA -C> GUESS/ANAL FILES 3-HRS APART TO CENERATE A GUESS/ANAL FILE VALID -C> AT THE PREPBUFR CENTER TIME -C> -2002-05-10 D. A. KEYSER -- ADDED "AIRCAR" TO THE LIST OF TABLE A -C> MESSAGE TYPES THAT WILL STILL HAVE THE BACKGROUND ENCODED WHEN -C> DOFCST IS FALSE (BECAUSE ACARS ARE NOW Q.C.'d IN PREPOBS_ACARSQC -C> PROGRAM) -C> -2003-09-02 D. A. KEYSER -- ADDED "MSONET" TO THE LIST OF TABLE A -C> MESSAGE TYPES THAT WILL HAVE THE VIRTUAL TEMPERATURE CALCULATED; -C> DOES NOT CALL UFBINT FOR OUTPUTTING DATA IF "NLEV" (4'TH -C> ARGUMENT) IS ZERO (NOW CAN ONLY HAPPEN FOR GOESND FORECAST DATA -C> WHEN ONLY RADIANCES ARE PRESENT) -C> -2004-08-30 D. A. KEYSER -- NOW INCLUDES THE 4 LAYER PWATERS, THESE -C> GET AN OBS. ERROR (EACH THE SAME AS TOTAL PWATER) AND AN EVENT -C> IS GENERATED WITH A REJECTED Q.M. FOR THE 4 LAYER PWATERS IF THE -C> PWATER OBS. ERROR READ IN IS MISSING (THIS CHANGE ALLOWS THE ETA/ -C> GSI TO PROCESS OBS. ERRORS IN THE PREPBUFR FILE THE SAME AS THE -C> ETA/3DVAR DID WHEN READING THE OBS. ERRORS FROM AN EXTERNAL -C> FILE); FOR "RASSDA" TYPES, ENCODES A SIMPLE COPY OF THE REPORTED -C> (VIRTUAL) TEMPERATURE AS A "VIRTMP" EVENT IF DOVTMP IS TRUE, GETS -C> NEW REASON CODE 3 -C> -2004-09-10 D. T. KLEIST -- ADDED CAPABILITY TO READ GUESS FIELDS FROM -C> EITHER HYBRID OR, AS BEFORE, SIGMA GLOBAL FORECAST FILES -C> -2005-01-03 D. A. KEYSER -- FIXED ERROR READING CDAS SGES FILE WHICH -C> STILL HAS A 207-WORD HEADER (T62) {2004-09-10 CHANGE ASSUMED ALL -C> SGES FILES HAD A 226-WORD HEADER (T254), BUT THIS IS VALID ONLY -C> FOR GFS SGES) -C> -2006-05-05 R. E. TREADON -- CHANGE VERTICAL INTERPOLATION TO DIRECTLY -C> USE PRESSURE PROFILE, NOT PRESSURE PROFILE CONVERTED TO SIGMA. -C> THIS CHANGE IS IN SUBROUTINE GBLEVN03. AS A RESULT OF THIS -C> CHANGE, SUBROUTINE GBLEVN07 WAS REMOVED. -C> -2006-07-14 D. A. KEYSER -- ADDED NEW NAMELIST SWITCH "SOME_FCST" -C> WHICH APPLIES ONLY WHEN EXISTING SWITCH "DOFCST" IS FALSE: IF -C> DOFCST=F AND SOME_FCST=T THEN, JUST AS BEFORE WHEN DOFCST=F, A -C> FORECAST WILL STILL BE ENCODED FOR REPORTS IN CERTAIN MESSAGE -C> TYPES USED IN SUBSEQUENT Q.C. PROGRAMS (I.E, "ADPUPA", "AIRCFT", -C> "AIRCAR", "PROFLR" OR "VADWND") (THE DEFAULT FOR SOME_FCST IS -C> TRUE); HOWEVER IF DOFCST=F AND SOME_FCST=F THEN A FORECAST WILL -C> NOT BE ENCODED INTO REPORTS IN ANY MESSAGE TYPE (THIS ALLOWS -C> THIS PROGRAM TO ENCODE OBS ERRORS AND/OR VIRTUAL TEMPERATURE -C> EVENTS INTO A PREPBUFR FILE WITHOUT ENCODING A FORECAST); ADDED -C> NEW NAMELIST SWITCH "ADPUPA_VIRT" WHICH, WHEN TRUE, INCLUDES -C> REPORTS IN MESSAGE TYPE ADPUPA (I.E., RAOBS, DROPS, MULTI-LEVEL -C> RECCOS) IN THE "VIRTMP" PROCESSING (PROCESSING THEM WITH SAME -C> LOGIC AS IN SUBROUTINE VTPEVN OF PROGRAM PREPOBS_CQCBUFR) -C> {NORMALLY "ADPUPA_VIRT" IS FALSE (DEFAULT) BECAUSE SUBSEQUENT -C> PROGRAM PREPOBS_CQCBUFR PERFORMS THIS FUNCTION} -C> -2007-09-14 S. MOORTHI -- ADDED CAPABILITY TO READ GENERALIZED SIGMA/ -C> HYBRID FILES FROM THE GFS USING "SIGIO" UTILITY; ALSO, CLEANED UP -C> SOME CODE; NEW ERROR CONDITION CODES 70 AND 71 ADDED -C> -2007-09-14 D. A. KEYSER -- FUNCTION OEFG01, WHICH RETURNS THE OBS -C> ERROR FOR A REQUESTED VARIABLE INTERP. TO A DEFINED PRESSURE -C> LEVEL FOR A DEFINED REPORT TYPE, MODIFIED TO USE EXACT LOGIC AS -C> IN GSI (MINIMUM LIMITING VALUE FOR OBS ERROR BASED ON VARIABLE -C> TYPE, LEVEL PRESSURE LIMITED TO MAX OF 2000 MB AND MIN OF ZERO -C> MB, A FEW OTHER MINOR CHANGES) - THIS WILL ALLOW GSI TO READ OBS -C> ERROR DIRECTLY OUT OF PREPBUFR FILE RATHER THAN OUT OF AN -C> EXTERNAL FILE; FOR PW TYPES, NOW PASSES REPORTED SURFACE PRESSURE -C> (PRSS * 0.01) INTO FUNCTION OEFG01 RATHER THAN VERTICAL -C> COORDINATE PRESSURE (POB), SINCE LATTER IS ALWAYS MISSING FOR -C> THESE TYPES (DOESN'T CHANGE VALUE COMING OUT OF OEFG01 SINCE IT -C> IS CONSTANT ON ALL LEVELS ANYWAY FOR PW); IN SUBR. GBLEVN02, Q.M. -C> 9 IS NOW ASSIGNED TO A VARIABLE ONLY IF ITS OBS ERROR IS MISSING, -C> OR IN THE CASE OF MOISTURE IF THE LEVEL IS ABOVE PRESSURE LEVEL -C> "QTOP_REJ" OR IF ITS TEMPERATURE OBS ERROR IS MISSING, ALL OTHER -C> EVENT (E.G., GROSS CHECK ERRORS) ASSIGN Q.M. 8 (EVEN IF OBS ERROR -C> IS MISSING), PRIOR TO THIS ONLY REJECTION OF PRESSURE ON LEVEL -C> RESULTED IN Q.M. 8, ALL OTHER REJECTIONS GOT Q.M. 9 - THIS MEANS -C> TRULY "BAD" OBS WILL NOW ALWAYS GET Q.M. 8 AND ONLY OBS FLAGGED -C> FOR NON-USE BY ASSIMILATION (BUT STILL "GOOD") WILL NOW GET Q.M. -C> 9 (GSI MONITORS, BUT DOES NOT USE, OBS WITH Q.M. 9, BUT IT DOES -C> NOT EVEN CONSIDER OBS WITH Q.M. 8); CORRECTED ERROR WHICH -C> MISTAKENLY ASSIGNED REASON CODE OF 9 INSTEAD OF 3 TO MOISTURE -C> WITH MISSING OBS ERROR; IN SUBR. GBLEVN02, Q.M. 9 WILL NOT BE -C> ASSIGNED TO A VARIABLE IF THAT VARIABLE ALREADY HAS A "BAD" Q.M. -C> (I.E., > 3 BUT < 15), IN FACT THE "PREVENT" EVENT WHICH WOULD -C> ASSIGN Q.M. 9 IS SKIPPED ENTIRELY (DO NOT WANT THE GSI TO MONITOR -C> THE OBS WHICH REALLY ARE ARE "BAD"); IN SUBR. GBLEVN08, FOR NON- -C> "ADPUPA" TYPES, Q.M. 9 IS NOW ASSIGNED TO CALCULATED VIRT. TEMPS -C> IF THE MOISTURE Q.M. IS 9 OR 15 AND ORIG. TEMP NOT "BAD", THESE -C> "VIRTMP" EVENTS RECEIVE NEW REASON CODE 4, HAD RECEIVED Q.M. 8 -C> WITH REASON CODE 2 LIKE VIRT. TEMPS CALCULATED FROM "BAD" -C> MOISTURE - THIS MEANS ONLY TRULY "BAD" VIRT. TEMPS WILL NOW GET -C> Q.M. 8 AND VIRT. TEMPS FLAGGED FOR NON-USE BY ASSIMILATION (BUT -C> STILL "GOOD") WILL NOW GET Q.M. 9 (GSI MONITORS, BUT DOES NOT -C> USE, OBS WITH Q.M. 9, BUT IT DOES NOT EVEN CONSIDER OBS WITH Q.M. -C> 8); IN SUBR. GBLEVN08, FOR "ADPUPA" TYPES, Q.M. 3 IS NOW ASSIGNED -C> TO CALCULATED VIRT. TEMPS ONLY IF THE MOISTURE Q.M. IS TRULY BAD -C> (I.E. > 3 BUT NOT 9 OR 15) (AND, AS BEFORE, ORIG. TQM IS 1 OR 2 -C> AND POB IS BELOW 700 MB) - BEFORE, TQM SET TO 3 WHEN QQM WAS 9 OR -C> 15 AND ALL OTHER CONDITIONS MET; FOR "SATEMP" TYPES, ENCODES A -C> SIMPLE COPY OF THE REPORTED (VIRTUAL) TEMPERATURE AS A "VIRTMP" -C> EVENT IF DOVTMP IS TRUE, GETS REASON CODE 3 (SIMILAR TO WHAT IS -C> ALREADY DONE FOR "RASSDA" TYPES) -C> -2010-01-29 D. A. KEYSER -- ADDED NEW NAMELIST SWITCH "RECALC_Q" -C> WHICH APPLIES ONLY WHEN EXISTING SWITCH "DOVTMP" IS FALSE: IF -C> DOVTMP=F AND RECALC_Q=T THEN, JUST AS BEFORE WHEN DOVTMP=F, SPEC. -C> HUMIDITY IS STILL RE-CALCULATED AND THE EVENT IS ENCODED INTO THE -C> PREPBUFR FILE (BUT VIRTUAL TEMP. IS NOT ENCODED) (THE DEFAULT FOR -C> RECALC_Q IS TRUE), HOWEVER IF DOVTMP=F AND RECALC_Q=F THEN SPEC. -C> HUMIDITY IS NOT RE-CALCULATED (AND NEITHER IS VIRTUAL -C> TEMPERATURE) (THIS ALLOWS THIS PROGRAM TO BYPASS ALL "VIRTMP" -C> EVENT PROCESSING); ADDED NEW NAMELIST SWITCH "DOPREV" WHICH, -C> WHEN TRUE, WRITES "PREVENT" EVENTS INTO THE PREPBUFR FILE (IT -C> ALWAYS DID THIS BEFORE) (DEFAULT), BUT NOW ALLOWS THE PROGRAM -C> TO BYPASS "PREVENT" EVENT PROCESSING WHEN DOPREV=F; INITIALIZED -C> ARRAY IDATE AS ZERO IN SUBR. GBLEVN10, CORRECTED BUG WHICH -C> EXPOSED PREVIOUSLY HIDDEN MEMORY CLOBBERING WHEN CALLING PROGRAMS -C> WERE LINKED TO NEW BUFRLIB; RULES IN SUBROUTINE GBLEVN02 REFINED -C> TO INCLUDE FULL SFC PRESSURE SANITY CHECK FOR ALL SFC REPORTS -C> (MASS, 18x, & WIND, 28x), BEFORE ONLY DONE FOR SFC MASS REPORTS -C> (18x) AND STILL NOT DONE FOR NON-SFC WIND REPORTS SINCE LOWEST -C> LEVEL PRESSURE NOT NECESSARILY AT THE SFC), AS A RESULT 28x WINDS -C> WILL NOW GET QM=8 IF PRESSURE FAILS SANITY CHECK (OFTEN HAPPENS -C> IN MESONET REPORTS) (GSI WAS ALREADY NOT USING THESE WINDS SINCE -C> PRESSURE QM SET TO 8 ALL ALONG) -C> -2012-11-20 J. WOOLLEN INITIAL PORT TO WCOSS. ADDED CALL TO BUFRLIB -C> ROUTINE GETBMISS TO ADAPT BMISS TO LINUX ENVIRONMENT IF NEED BE -C> {I.E., OBTAINS BUFRLIB MISSING (BMISS) VIA CALL TO GETBMISS -C> RATHER THAN HARDWIRING IT TO 10E10 (10E10 CAN CAUSE INTEGER -C> OVERFLOW ON WCOSS - SEE CALLING PROGRAM FOR MORE INFO)} -C> -2013-02-13 D. A. KEYSER -- FINAL CHANGES TO RUN ON WCOSS: USE -C> FORMATTED PRINT STATEMENTS WHERE PREVIOUSLY UNFORMATTED PRINT WAS -C> > 80 CHARACTERS; RENAME ALL REAL(8) VARIABLES AS *_8 -C> -2013-04-12 D. A. KEYSER -- IN SUBROUTINE GBLEVN08, DON'T ALLOW -C> CALCULATED Q TO BE < 0 WHICH CAN OCCUR ON WCOSS FOR CASES OF -C> HORRIBLY BAD MESONET DATA -C> -2014-03-25 S. MELCHIOR -- ADDED NEW NAMELIST SWITCH "DOPMSL" WHICH, -C> WHEN TRUE, DERIVES PMSL (MNEMONIC "PMO") FROM REPORTED STATION -C> PRESSURE ("POB"), STATION HEIGHT/ELEVATION ("ZOB") AND THE JUST- -C> COMPUTED VIRTUAL TEMPERATURE FOR SURFACE REPORTS IN CASES WHEN -C> REPORTED PMSL IS MISSING (DONE IN SUBROUTINE GBLEVN08). DOVTMP -C> MUST BE TRUE AND DOANLS MUST BE FALSE ("PREVENTS" MODE). THE -C> DERIVED PMSL EITHER GETS A QUALITY MARK ("PMQ") OF 3 OR INHERITS -C> THE STATION PRESSURE QUALITY MARK ("PQM"), WHICHEVER IS GREATER. -C> THE VALUE OF THE PMSL INDICATOR (NEW MNEMONIC "PMIN") IS SET TO 1 -C> TO DENOTE PMSL WAS DERIVED RATHER THAN OBSERVED. THE DEFAULT FOR -C> "DOPMSL" IS FALSE (NORMALLY ONLY TRUE IN RTMA AND URMA RUNS). IT -C> IS FORCED TO BE FALSE IN "POSTEVENTS" MODE (DOANLS=TRUE). IN -C> SUBROUTINE GBLEVN02, SFCSHP REPORTS WITH CALM WINDS AND NON- -C> MISSING BACKGROUND U- OR V-COMPONENT WIND .GE. 5 M/SEC ARE -C> FLAGGED WITH Q.M. 8 (EVENT PGM "PREVENT", REASON CODE 8). -C> -2014-05-08 JWhiting -- altered print statement (2 format) in GBLEVN10 -C> subroutine; increased field width for spectral resolution to -C> accommodate models w/ up to 5-digit resolution (I3 to I5). -C> -2016-06-13 FANGLIN YANG AND RUSS TREADON -- HANG LEI ADDED NEMSIO -C> TO SUBROUTINE GBLEVN10 AND REMOVED ALL SIGIO CAPABILITY. -C> THIS UPDATE RESTORES GBLEVN10 FOR PROCESSING SIGIO INPUT, AND -C> ADDS A NEW GBLEVN12 FOR PROCESSING NEMSIO INPUT. THE INPUT GFS -C> FILE TYPE SIGIO VS NEMSIO IS NOW DETERMINED IN THE MAIN PROGRAM. -C> THE CODE IS ALSO UPDATED TO REMOVE BUGS. SUBROUTINE SIGIO_MODPR -C> IS USED TO COMPUTE LAYER AND INTERFACE PRESSURE FOR NEMSIO INPUT. -C> -2017-02-17 D Keyser & J Whiting -- In subroutine GBLEVN12, removed -C> references to multiple input files, since only 1 nemsio formatted -C> input is needed (no interpolation is attempted; c.f.; what is -C> done with sigio formatted inputs); GBLEVN12 routine still retains -c some structures (esp array references) left over from multi file -C> input. Updated comments and docblock to account for new NEMSIO -C> input. -C> -2017-02-22 D. Keyser -- Further changes to subr. GBLEVN12 to remove -C> array and logic references to multiple input files. -C> -2019-10-31 Hang Lei -- Add GBLEVN13 to process netcdf input. -C> -C> -C> USAGE: CALL GBLEVENTS(IDATEP,IUNITF,IUNITE,IUNITP,IUNITS,SUBSET, -C> $ NEWTYP) -C> -C> INPUT ARGUMENT LIST: -C> @param IDATEP - CENTER DATE FOR PREPBUFR FILE IN THE FORM YYYYMMDDHH -C> @param IUNITF - 2-WORD ARRAY: -C> For SIGIO input: -C> - WORD 1 - UNIT NUMBER OF FIRST INPUT SIGIO-BASED GLOBAL -C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR -C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN -C> THIS FILE IS VALID AT THE DATE IN IDATEP, IF HH IN -C> IDATEP IS NOT A MULTIPLE OF 3 THEN THIS FILE IS VALID -C> AT THE CLOSEST TIME PRIOR TO THE DATE IN IDATEP THAT -C> IS A MULTIPLE OF 3 -C> - WORD 2 - UNIT NUMBER OF SECOND INPUT SIGIO-BASED GLOBAL -C> (SIGMA OR HYBRID) FILE (EITHER FIRST GUESS OR -C> ANALYSIS); IF HH IN IDATEP IS A MULTIPLE OF 3 THEN -C> THIS FILE IS EMPTY, IF HH IN IDATEP IS NOT A MULTIPLE -C> OF 3 THEN THIS FILE IS VALID AT THE CLOSEST TIME AFTER -C> THE DATE IN IDATEP THAT IS A MULTIPLE OF 3 -C> For NEMSIO input: -C> - WORD 1 - UNIT NUMBER OF INPUT NEMSIO-BASED GLOBAL FILE -C> (EITHER FIRST GUESS OR ANALYSIS); ALWAYS VALID AT AT -C> THE DATE IN IDATEP -C> - WORD 2 - NOT USED, SHOULD BE A NULL FILE -C> @param IUNITE - UNIT NUMBER OF INPUT OBSERVATION ERROR FILE -C> - (USED ONLY IN PREVENTS MODE) -C> @param IUNITP - UNIT NUMBER OF OUTPUT PREPBUFR DATA SET -C> @param IUNITS - UNIT NUMBER OF "PREVENT" EVENTS DATA FILTERING -C> - SUMMARY PRINT FILE -C> - (USED ONLY IN PREVENTS MODE) -C> @param SUBSET - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR -C> - REPORT BEING PROCESSED -C> @param NEWTYP - INDICATOR IF THE BUFR MESSAGE TABLE A ENTRY HAS -C> - CHANGED FROM THAT OF THE PREVIOUS REPORT (=0 - NO, -C> - =1 - YES) -C> -C> -C> INPUT FILES: -C> - UNIT 05 - STANDARD INPUT (DATA CARDS - SEE NAMELIST -C> DOCUMENTATION BELOW) -C> (NOTE: IF STANDARD INPUT FILE IS NULL, THEN THIS -C> SUBROUTINE RUNS IN POSTEVENTS MODE) -C> - UNIT AA - PREPBUFR DATA SET -C> (WHERE AA IS UNIT NUMBER DEFINED AS IUNITP IN -C> INPUT ARGUMENT LIST) -C> - UNIT BB - GUESS (PREVENTS MODE) OR ANALYSIS (POSTEVENTS MODE) -C> FILE -C> (WHERE BB IS UNIT NUMBER DEFINED AS IUNITF(1) IN -C> INPUT ARGUMENT LIST) -C> - UNIT CC - GUESS (PREVENTS MODE) OR ANALYSIS (POSTEVENTS MODE) -C> FILE -C> (WHERE CC IS UNIT NUMBER DEFINED AS IUNITF(2) IN -C> INPUT ARGUMENT LIST) -C> NOTE: only valid for SIGIO input -C> - UNIT DD - OBSERVATION ERROR FILE (WHERE DD IS UNIT NUMBER -C> DEFINED AS IUNITE IN INPUT ARGUMENT LIST) -C> (USED ONLY IN PREVENTS MODE) -C> -C> OUTPUT FILES: -C> - UNIT 06 - STANDARD OUTPUT PRINT -C> - UNIT AA - PREPBUFR DATA SET -C> (WHERE AA IS UNIT NUMBER DEFINED AS IUNITP IN -C> INPUT ARGUMENT LIST) -C> - UNIT DD - "PREVENT" EVENTS DATA FILTERING SUMMARY PRINT FILE -C> (WHERE DD IS UNIT NUMBER DEFINED AS IUNITS IN -C> INPUT ARGUMENT LIST) -C> (USED ONLY IN PREVENTS MODE) -C> -C> SUBPROGRAMS CALLED: -C> UNIQUE: GBLEVN02 GBLEVN03 GBLEVN04 GBLEVN06 OEFG01 -C> GBLEVN08 GBLEVN10 GBLEVN11 GBLEVN11D GBLEVN12 -C> GBLEVN13 -C> GETLATS -C> MODULES: GBLEVN_MODULE SIGIO_MODULE SIGIO_R_MODULE -C> NEMSIO_MODULE NEMSIO_OPENCLOSE NEMSIO_READ -C> NEMSIO_WRITE -C> LIBRARY: -C> SIGIO - SIGIO_RROPEN SIGIO_RRHEAD SIGIO_SCLOSE SIGIO_ALDATS -C> SIGIO_ALDATM SIGIO_RRDATS SIGIO_RRDATM SIGIO_AXDATS -C> SIGIO_AXDATM SIGIO_MODPR SIGIO_CNVTDV -C> SPLIB - SPTEZM SPTEZMV SPLAT -C> W3NCO - W3MOVDAT ERREXIT -C> BUFRLIB - UFBINT UFBQCD GETBMISS IBFMS -C> NEMSIO - NEMSIO_OPEN NEMSIO_CLOSE NEMSIO_INIT -C> NEMSIO_GETFILEHEAD NEMSIO_READRECV NEMSIO_FINALIZE -C> NEMSIO_GETHEADVAR NEMSIO_READRECVw34 -C> -C> EXIT STATES: -C> - COND = 0 - SUCCESSFUL RUN -C> - COND = 60 - OBSERVATION ERROR TABLE EMPTY OR DOES NOT EXIST -C> - COND = 61 - VARIABLE NLTD .NE. VARIABLE NLEV -C> - COND = 62 - VARIABLE NLTQ .NE. VARIABLE NLEV -C> - COND = 63 - VARIABLE NLQQ .NE. VARIABLE NLEV -C> - COND = 68 - DATE OF FIRST GUESS/ANALYSIS FILE(S) DOES NOT MATCH, -C> OR AT LEAST SPAN, THE CENTER DATE FOR THE PREPBUFR -C> FILE -C> - COND = 69 - FOR SIGIO INPUT GLOBAL FILES, VARIABLE KMAX TOO BIG -C> - UNABLE TO TRANSFORM FIRST GUESS OR ANALYSIS FILE(S) -C> - COND = 70 - FOR SIGIO INPUT GLOBAL FILES, CALL TO SIGIO_RROPEN -C> RETURNED WITH NON-ZERO R.C. -C> - COND = 71 - FOR SIGIO INPUT GLOBAL FILES, CALL TO SIGIO_RRHEAD -C> RETURNED WITH NON-ZERO R.C. -C> -C> -C> REMARKS: THIS SUBROUTINE MAY NOT WORK CORRECTLY IN THE EIGHT BYTE -C> INTEGER W3NCO (_8) LIBRARY. PLEASE COMPILE APPLICATION CODE -C> USING A FOUR BYTE REAL W3NCO LIBRARY (_4 OR _d). -C> -C> THIS ROUTINE PROCESSES ONE REPORT AT A TIME. IT EXPECTS THAT THE -C> CALLING PROGRAM HAS ALREADY ENCODED THE REPORT INTO THE PREPBUFR -C> FILE VIA THE UFBINT OR UFBCPY ROUTINES. THE CALLING PROGRAM -C> SHOULD THEN CALL THIS ROUTINE AND, UPON ITS RETURN, THE CALLING -C> PROGRAM SHOULD CALL WRITSB TO ACTUALLY WRITE THE UPDATED SUBSET -C> (REPORT) INTO THE BUFR MESSAGE. -C>C -C> ***** VARIABLES IN NAMELIST PREVDATA READ IN BY THIS SUBROUTINE ***** -C> (NOTE: IF STANDARD INPUT FILE IS NULL, THEN THIS -C> SUBROUTINE RUNS IN POSTEVENTS MODE - DOANLS=TRUE -C> AND ALL OTHER VARIABLES ARE SET TO FALSE) -C>C -C>C -C> - DOPREV - WRITE "PREVENT" EVENT INTO THE PREPBUFR FILE? -C> DOPREV = .TRUE. ---> YES (DEFAULT) -C> DOPREV = .FALSE. ---> NO -C> - DOVTMP, ADPUPA_VIRT & RECALC_Q: -C> DOVTMP - WRITE VIRTUAL TEMPERATURE EVENT ("VIRTMP") INTO THE -C> PREPBUFR FILE (I.E., RE-CALCULATE SPECIFIC HUMIDITY -C> THEN CALCULATE VIRTUAL TEMPERATURE) FOR THE FOLLOWING -C> TYPES OF REPORTS: -C> ADPUPA_VIRT = .FALSE. ---> SURFACE LAND, MARINE, -C> MESONET AND RASS REPORTS? -C> ADPUPA_VIRT = .TRUE. ---> SURFACE LAND, MARINE, -C> MESONET RASS, RAOB, DROP -C> AND MULTI-LEVEL RECCO -C> REPORTS? -C> FOR ALL TYPES EXCEPT RASS, THIS WILL ATTEMPT TO -C> CALCULATE VIRTUAL TEMPERATURE FROM SENSIBLE TEMPERATURE -C> AND THE JUST RE-CALCULATED SPECIFIC HUMIDITY AND ENCODE -C> IT AS A STACKED EVENT IN THE PREPBUFR FILE. FOR RASS -C> REPORTS THIS WILL JUST ENCODE THE REPORTED TEMPERATURE -C> AS A STACKED EVENT IN THE PREPBUFR FILE SINCE THE -C> REPORTED TEMPERATURE IS ALREADY VIRTUAL (NO MOISTURE IS -C> PRESENT SO Q IS NOT RE-CALCULATED FOR RASS REPORTS). -C> DOVTMP = .TRUE. ---> YES (DEFAULT) -C> DOVTMP = .FALSE. -C> RECALC_Q = .TRUE. ---> RE-CALCULATE SPECIFIC -C> HUMIDITY BUT DO NOT THEN -C> CALCULATE VIRTUAL -C> TEMPERATURE (DEFAULT) -C> RECALC_Q = .FALSE. ---> NO, DO NOT RE-CALCULATE -C> SPECIFIC HUMIDITY AND DO -C> NOT CALCULATE VIRTUAL -C> TEMPERATURE -C> {NOTE1: FOR SURFACE LAND, MARINE AND MESONET REPORTS, (AND -C> RAOB, DROP AND MULTI-LEVEL RECCO REPORTS IF -C> "ADPUPA_VIRT"=TRUE) DOVTMP=FALSE WILL STILL RE-CALCULATE -C> SPECIFIC HUMIDITY AND ENCODE IT AS A STACKED EVENT IN -C> THE PREPBUFR FILE UNLESS EITHER DOANLS IS TRUE OR -C> RECALC_Q IS FALSE.) -C> (NOTE2: DOES NOT APPLY TO ANY REPORT TYPES OTHER THAN THOSE -C> MENTIONED ABOVE) -C> (NOTE3: IF DOANLS=TRUE, THEN DOVTMP IS NOT ONLY FORCED TO BE -C> FALSE, BUT ALSO SPECIFIC HUMIDITY IS NOT RE-CALCULATED.) -C> (NOTE4: ADPUPA_VIRT DEFAULTS TO FALSE.) -C> (NOTE5: IF DOVTMP=TRUE, THEN RECALC_Q IS MEANINGLESS.) -C> (NOTE6: RECALC_Q DEFAULTS TO TRUE.) -C> -C> - DOFCST & SOME_FCST: -C> DOFCST - ENCODE FORECAST (FIRST GUESS) VALUES, INTERPOLATED FROM -C> EITHER A SIGIO (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT -C> GLOBAL FILE, INTO THE PREPBUFR FILE FOR ALL MESSAGE -C> TYPES OR AT LEAST SOME MESSAGE TYPES? -C> DOFCST = .TRUE. ---> YES, ENCODE FORECST FOR ALL -C> MESSAGE TYPES (DEFAULT) -C> DOFCST = .FALSE. -C> SOME_FCST = .FALSE. ---> NO, DO NOT ENCODE FORECAST -C> FOR ANY MESSAGE TYPE -C> (VALUES REMAIN MISSING) -C> SOME_FCST = .TRUE. ---> YES, BUT ONLY FOR MESSAGE -C> TYPES "ADPUPA", "AIRCFT", -C> "AIRCAR", "PROFLR" OR -C> "VADWND" (VALUES REMAIN -C> MISSING FOR ALL OTHER -C> MESSAGE TYPES) -C> (NOTE1: THE CASE DOFCST=FALSE & SOME_FCST=TRUE WRITES THE -C> FORECAST VALUES FOR THE TYPES MENTIONED ABOVE BECAUSE -C> THEY ARE NEEDED BY SUBSEQUENT QUALITY CONTROL PROGRAMS.) -C> (NOTE2: THIS WAS ADDED AS A TIME SAVING FEATURE IN THE -C> NON-GLOBAL VERSIONS SINCE ONLY THE GLOBAL REQUIRES A -C> FIRST GUESS TO BE PRESENT FOR ALL CONVENTIONAL MESSAGE -C> TYPES IN THE PREPBUFR FILE.) -C> (NOTE3: IF DOANLS=TRUE, THEN DOFCST & SOME_FCST ARE FORCED TO BE -C> FALSE, MEANING A GUESS WILL NOT BE ENCODED FOR ANY -C> MESSAGE TYPE.) -C> (NOTE4: IF DOFCST=TRUE, THEN SOME_FCST IS MEANINGLESS.) -C> (NOTE5: SOME_FCST DEFAULTS TO TRUE.) -C> -C> - DOANLS - ENCODE ANALYZED VALUES, INTERPOLATED FROM EITHER A SIGIO -C> (SIGMA OR HYBRID) INPUT OR NEMSIO INPUT GLOBAL FILE, -C> INTO THE PREPBUFR FILE - -C> POSTEVENTS MODE - ? -C> DOANLS = .TRUE. ---> YES, FOR ALL MESSAGE TYPES -C> DOANLS = .FALSE. ---> NO, FOR ALL MESSAGE TYPES -C> - PREVENTS MODE - (DEFAULT) -C> (NOTE: DOANLS=TRUE WILL OVERRIDE AND FORCE TO FALSE ALL OTHER -C> SWITCHES. IN ADDITION, THE FORECAST VALUES WILL NOT -C> BE ENCODED FOR ANY MESSAGE TYPE AND SPECIFIC HUMIDITY -C> WILL NOT BE RE-CALCULATED.) -C> -C> - DOBERR - ENCODE OBSERVATIONAL ERROR VALUES, AS READ FROM OBS. -C> ERROR FILE, INTO THE PREPBUFR FILE? -C> DOBERR = .TRUE. ---> YES (DEFAULT) -C> DOBERR = .FALSE. ---> NO (VALUES REMAIN MISSING) -C> (NOTE1: THIS WAS ADDED AS A TIME SAVING FEATURE IN THE -C> RAP -AND PREVIOUS RUC- VERSION SINCE IT DOES NOT REQUIRE -C> OBSERVATIONAL ERRORS TO BE PRESENT IN THE PREPBUFR FILE.) -C> (NOTE2: IF DOANLS=TRUE, THEN DOBERR IS FORCED TO BE FALSE.) -C> -C> - QTOP_REJ - THE PRESSURE LEVEL (IN MB) ABOVE WHICH ALL SPECIFIC -C> HUMIDITY QUALITY MARKERS ARE "REJECTED" (THE QUALITY -C> MARKER IS SET TO 9 ON ALL PRESSURE LEVELS LESS THAN -C> THIS LEVEL) (DEFAULT=300.) -C> -C> - SATMQC - PERFORM SPECIAL QUALITY CONTROL ON SATELLITE TEMPERATURE -C> SOUNDINGS IN REPORT TYPES 160-179? -C> SATMQC = .TRUE. ---> YES -C> SATMQC = .FALSE. ---> NO (DEFAULT) -C> (NOTE: THIS APPLIES ONLY TO THE CDAS OR HISTORICAL RE-RUNS -C> WITH TEMPERATURE SOUNDINGS IN THESE REPORT TYPES) -C> -C> - DOPMSL - ENCODE DERIVED PMSL ("PMO") FOR ALL SURFACE REPORTS WHEN -C> REPORTED PMSL IS MISSING - ? -C> DOPMSL = .TRUE. ---> YES -C> DOPMSL = .FALSE. ---> NO ("PMO" REMAINS MISSING)(DEFAULT) -C> {NOTE: THIS APPLIES ONLY WHEN DOVTMP=TRUE AND DOANLS=FALSE -C> ("PREVENTS" MODE), VIRTUAL TEMPERATURE CAN BE CALCULATED, -C> AND STATION PRESSURE AND SURFACE HEIGHT/ELEVATION ARE -C> BOTH PRESENT. THE DERIVED PMSL EITHER GETS A QUALITY -C> MARK ("PMQ") OF 3 OR INHERITS THE STATION PRESSURE -C> QUALITY MARK ("PQM") PQM, WHICHEVER IS GREATER. THE VALUE -C> OF THE PMSL INDICATOR ("PMIN") IS SET TO 1 TO DENOTE PMSL -C> WAS DERIVED RATHER THAN OBSERVED.} -C> - MODULE GBLEVN_MODULE - - IMPLICIT NONE - - SAVE - - INTEGER IMAX,JMAX,KMAX,KMAXS - INTEGER*4 IDVC,IDSL,NVCOORD,SFCPRESS_ID,THERMODYN_ID - REAL (KIND=8), ALLOCATABLE :: IAR14T(:,:,:), IAR15U(:,:,:), - $ IAR16V(:,:,:), IAR17Q(:,:,:), - $ IAR12Z(:,:), IAR13P(:,:), - $ IARPSI(:,:,:), IARPSL(:,:,:), - $ IARPSD(:,:,:) - REAL (KIND=4), ALLOCATABLE :: VCOORD(:,:) - REAL DLAT,DLON - - END MODULE GBLEVN_MODULE - - SUBROUTINE GBLEVENTS(IDATEP,IUNITF,IUNITE,IUNITP,IUNITS,SUBSET, - $ NEWTYP) - - USE SIGIO_MODULE - USE SIGIO_R_MODULE - USE NEMSIO_MODULE - USE NEMSIO_OPENCLOSE - USE NEMSIO_READ - USE GBLEVN_MODULE - use netcdf - - INTEGER, PARAMETER :: IM=384, JM=IM/2+1 - integer :: idrt = 0 - - CHARACTER*80 HEADR,OBSTR,QMSTR,FCSTR,OESTR,ANSTR - CHARACTER*8 SUBSET - REAL(8) OBS_8,QMS_8,BAK_8,SID_8,HDR_8(10) - REAL(8) BMISS,GETBMISS - LOGICAL DOVTMP,DOFCST,SOME_FCST,DOBERR,FCST,VIRT,DOANLS, - $ SATMQC,ADPUPA_VIRT,RECALC_Q,DOPREV,dopmsl - INTEGER*4 IRET,IRET1 - INTEGER*4 IUNITF(2),ncid - - CHARACTER*20 CFILE1 - TYPE(SIGIO_HEAD) :: HEAD1 - TYPE(NEMSIO_GFILE) :: GFILE1 - - COMMON /GBEVAA/ SID_8,OBS_8(13,255),QMS_8(12,255),BAK_8(12,255), - $ XOB,YOB,DHR,TYP,NLEV - COMMON /GBEVBB/ PVCD,VTCD - COMMON /GBEVCC/ DOVTMP,DOFCST,SOME_FCST,DOBERR,FCST,VIRT, - $ QTOP_REJ,SATMQC,ADPUPA_VIRT,RECALC_Q,DOPREV,dopmsl - COMMON /GBEVDD/ ERRS(300,33,6) - COMMON /GBEVFF/ BMISS - - SAVE - - DATA IFIRST/0/ - - DATA HEADR / - $ 'SID XOB YOB DHR TYP '/ - DATA OBSTR / - $ 'POB QOB TOB ZOB UOB VOB PWO PW1O PW2O PW3O PW4O CAT PRSS '/ - DATA QMSTR / - $ 'PQM QQM TQM ZQM WQM PWQ PW1Q PW2Q PW3Q PW4Q NUL NUL '/ - DATA FCSTR / - $ 'PFC QFC TFC ZFC UFC VFC PWF PW1F PW2F PW3F PW4F NUL '/ - DATA ANSTR / - $ 'PAN QAN TAN ZAN UAN VAN PWA PW1A PW2A PW3A PW4A NUL '/ - DATA OESTR / - $ 'POE QOE TOE ZOE WOE PWE PW1E PW2E PW3E PW4E NUL NUL '/ - - NAMELIST /PREVDATA/DOVTMP,DOFCST,SOME_FCST,DOBERR,DOANLS, - $ QTOP_REJ,SATMQC,ADPUPA_VIRT,RECALC_Q,DOPREV,dopmsl - -C---------------------------------------------------------------------- -C---------------------------------------------------------------------- - - IF(IFIRST.EQ.0) THEN - -C ------------------------------- -C FIRST TIME IN DO A FEW THINGS... -C ------------------------------- - - IFIRST = 1 - PRINT 700 - 700 FORMAT(/1X,100('#')/' =====> SUBROUTINE GBLEVENTS INVOKED FOR ', - $ 'THE FIRST TIME - VERSION LAST UPDATED 2017-02-22'/) - - BMISS = GETBMISS() - print * - print *, 'BUFRLIB value for missing passed into GBLEVENTS ', - $ 'is: ',bmiss - print * - -C INITIALIZE NAMELIST SWITCHES TO DEFAULT VALUES -C ---------------------------------------------- - - DOVTMP = .TRUE. - DOPREV = .TRUE. - RECALC_Q = .TRUE. - DOFCST = .TRUE. - SOME_FCST = .TRUE. - DOBERR = .TRUE. - DOANLS = .FALSE. - QTOP_REJ = 300. - SATMQC = .FALSE. - ADPUPA_VIRT = .FALSE. - dopmsl = .false. - READ(5,PREVDATA,ERR=101,END=102) - GO TO 103 -C----------------------------------------------------------------------- - 101 CONTINUE - -C ERROR READING STANDARD INPUT - THIS DEFAULTS TO POSTEVENTS MODE -C --------------------------------------------------------------- - - PRINT 7013 - 7013 FORMAT(/' ##> GBLEVENTS: ERROR READING STANDARD INPUT DATA CARDS', - $ ' -- DEFAULTS TO "POSTEVENTS" MODE'/) - DOANLS = .TRUE. - GO TO 103 - -C----------------------------------------------------------------------- - 102 CONTINUE - -C STANDARD INPUT IS EMPTY - THIS DEFAULTS TO POSTEVENTS MODE -C ---------------------------------------------------------- - - PRINT 7014 - 7014 FORMAT(/' ##> GBLEVENTS: STANDARD INPUT DATA CARDS DO NOT ', - $ 'EXIST -- DEFAULTS TO "POSTEVENTS" MODE'/) - DOANLS = .TRUE. - -C----------------------------------------------------------------------- - 103 CONTINUE - IF(DOANLS) THEN - DOVTMP = .FALSE. - DOPREV = .FALSE. - DOFCST = .FALSE. - SOME_FCST = .FALSE. - DOBERR = .FALSE. - ADPUPA_VIRT = .FALSE. - dopmsl = .false. - ENDIF - IF(DOVTMP) RECALC_Q=.TRUE. ! RECALC_Q must be T if DOVTMP is T - WRITE (6,PREVDATA) - - FCST = DOFCST - VIRT = .FALSE. - -C CHECK VALID-TIME DATE OF GUESS/ANALYSIS FILE(S) AGAINST THE CENTER -C DATE FOR THE PREPBUFR FILE AND OBTAIN THE FIRST GUESS/ANALYSIS -C UNLESS ALL OF DOFCST, SOME_FCST, DOANLS ARE FALSE -C ------------------------------------------------------------------ - - IF(.NOT.DOANLS) THEN - IF(.NOT.DOFCST.AND..NOT.SOME_FCST) THEN - PRINT 901 - 901 FORMAT(/' --> GBLEVENTS: PREVENTS MODE - FIRST GUESS NOT READ ', - $ 'IN'/) - ELSE - PRINT 701 - 701 FORMAT(/' --> GBLEVENTS: PREVENTS MODE - DATE CHECK AND ', - $ 'TRANSFORM THE FIRST GUESS'/) - ENDIF - ELSE - PRINT 7701 - 7701 FORMAT(/' --> GBLEVENTS: POSTEVENTS MODE - DATE CHECK AND ', - $ 'TRANSFORM THE ANALYSIS'/) - ENDIF - - IF(DOFCST .OR. SOME_FCST .OR. DOANLS) THEN - WRITE(CFILE1,'("fort.",I2.2)') IUNITF(1) - - iret = nf90_open(trim(cfile1),nf90_nowrite,ncid) - if (iret == 0) then - print *,' ===> GFS FCST/ANAL INPUT IS NETCDF' - CALL GBLEVN13(IUNITF,IDATEP,IM,JM,IDRT) - else - CALL SIGIO_RROPEN(IUNITF(1),CFILE1,IRET) - CALL SIGIO_SRHEAD(IUNITF(1),HEAD1,IRET1) - IF(IRET == 0 .AND. IRET1 == 0) THEN - print *,' ===> GLOBAL FCST/ANAL INPUT IS SIGIO' - CALL SIGIO_SCLOSE(IUNITF(1),IRET) - CALL GBLEVN10(IUNITF,IDATEP,IM,JM,IDRT) - ELSE - CALL NEMSIO_OPEN(GFILE1,trim(CFILE1),'read',IRET=IRET) - IF(IRET == 0) THEN - CALL NEMSIO_CLOSE(GFILE1,IRET=IRET) - print *,' ===> GFS FCST/ANAL INPUT IS NEMSIO' - CALL GBLEVN12(IUNITF,IDATEP,IM,JM,IDRT) - ENDIF - ENDIF - ENDIF - ENDIF - - print*,'after returning from GBLEVN10, GBLEVN12,', - $ ' or GBLEVN13. idrt=',idrt - - IF(DOBERR) THEN - -C IF REQUESTED, READ ERROR FILES (ONLY POSSIBLE IN PREVENTS MODE) -C --------------------------------------------------------------- - - PRINT 702 - 702 FORMAT(/' --> GBLEVENTS: READ ERROR FILES'/) - - CALL GBLEVN01(IUNITE) - - ELSE - - ERRS = 0 - IF(.NOT.DOANLS) PRINT 3702 - 3702 FORMAT(/' --> GBLEVENTS: OBS. ERROR NOT ENCODED IN PREPBUFR ', - $ '(BY CHOICE)'/) - - ENDIF - -C OBTAIN NECESSARY PROGRAM CODES (ONLY USED IN PREVENTS MODE) -C ----------------------------------------------------------- - - CALL UFBQCD(IUNITP,'PREVENT',PVCD) - CALL UFBQCD(IUNITP,'VIRTMP ',VTCD) - - PRINT 703 - 703 FORMAT(/1X,100('#')/) - -C SET-UP OUTPUT "PREVENT" EVENTS DATA FILTERING SUMMARY PRINT FILE -C (ONLY USED IN PREVENTS MODE) -C ---------------------------------------------------------------- - - IF(.NOT.DOANLS) WRITE(IUNITS,1701) IDATEP - 1701 FORMAT(//130('#')//38X,'*** "PREVENT" EVENTS DATA FILTERING ', - $ 'SUMMARY ***'/35X,'--> CENTER DATE FOR PREPBUFR FILE IS: ',I10, - $ ' <--'//) - -C---------------------------------------------------------------------- -C---------------------------------------------------------------------- - - ENDIF - - IF(.NOT.DOANLS) THEN - - IF(NEWTYP.EQ.1) WRITE(IUNITS,1702) SUBSET - 1702 FORMAT(130('-')/39X,'--> SUMMARY FOR TABLE A ENTRY "',A8,'" <--'/) - - IF(.NOT.DOFCST .AND. SOME_FCST) FCST = (SUBSET.EQ.'ADPUPA ' - $ .OR.SUBSET.EQ.'PROFLR '.OR.SUBSET .EQ.'AIRCFT '.OR.SUBSET - $ .EQ.'AIRCAR '.OR.SUBSET .EQ.'VADWND ') - -C Will not subject ACARS reports to virtual temp. processing until -C spec. humidity is used in production - -ccccc VIRT = (SUBSET.EQ.'ADPSFC '.OR.SUBSET.EQ.'SFCSHP '.OR. -ccccc$ SUBSET.EQ.'MSONET '.OR.SUBSET.EQ.'AIRCAR '.OR. -ccccc$ SUBSET.EQ.'RASSDA '.OR.SUBSET.EQ.'SATEMP '.OR. -ccccc$ (SUBSET.EQ.'ADPUPA '.AND.ADPUPA_VIRT)) - VIRT = (RECALC_Q.AND.(SUBSET.EQ.'ADPSFC '.OR. - $ SUBSET.EQ.'SFCSHP '.OR. - $ SUBSET.EQ.'MSONET '.OR. - $ SUBSET.EQ.'RASSDA '.OR. - $ SUBSET.EQ.'SATEMP '.OR. - $ (SUBSET.EQ.'ADPUPA '.AND.ADPUPA_VIRT))) - - - IF(.NOT.(FCST.OR.DOBERR.OR.VIRT.OR.DOPREV)) THEN - IF(NEWTYP.EQ.1) WRITE(IUNITS,1703) - 1703 FORMAT(/' ==> DATA FILTERING NOT PERFORMED FOR THIS TABLE A ', - $ 'ENTRY -- FORECAST, OBS ERROR, "VIRTMP", "PREVENT" PROCESSING ', - $ 'NOT DONE'/) - RETURN - ENDIF - - ENDIF - -C READY TO RETRIEVE NECESSARY INFORMATION OUT OF THE NEXT REPORT WHICH -C HAS BEEN "UFB" ENCODED INTO THE PREPBUFR FILE BY THE CALLING PROGRAM -C (USE NEGATIVE UNIT NUMBER HERE SINCE FILE OPEN FOR OUTPUT) -C (NOTE: THE CALLING PROGRAM HAS NOT YET WRIITEN THE REPORT INTO -C THE PREPBUFR FILE VIA WRITSB!) -C ---------------------------------------------------------------- - - CALL UFBINT(-IUNITP,OBS_8,13,255,NLEV,OBSTR) - CALL UFBINT(-IUNITP,QMS_8,12,255,NLEV,QMSTR) - CALL UFBINT(-IUNITP,HDR_8,10, 1,IRET,HEADR) - SID_8 = HDR_8(1) - XOB = HDR_8(2) - YOB = HDR_8(3) - DHR = HDR_8(4) - TYP = HDR_8(5) - - IF(FCST.OR.DOANLS) THEN - -C PREVENTS MODE: ENCODE FIRST GUESS VALUES INTO PREPBUFR FILE -C ------------------------------------------------------------ - -C POSTEVENTS MODE: ENCODE ANALYSIS VALUES INTO REPORT AND RETURN TO -C CALLING PROGRAM TO WRITE GBL-EVENTED REPORT -C (SUBSET) INTO PREPBUFR FILE -C ----------------------------------------------------------------- - - CALL GBLEVN03(SUBSET) - IF(NLEV.GT.0) THEN - IF(FCST) THEN - CALL UFBINT(IUNITP,BAK_8,12,NLEV,IRET,FCSTR) - ELSE - CALL UFBINT(IUNITP,BAK_8,12,NLEV,IRET,ANSTR) - RETURN - ENDIF - ENDIF - ENDIF - -C -------------------------------------------------------------------- -C LOGIC FROM HERE ON PERTAINS ONLY TO PREVENTS MODE OF THIS SUBROUTINE -C -------------------------------------------------------------------- - -C ENCODE OBSERVATION ERRORS INTO REPORT -C ------------------------------------- - - IF(DOBERR) THEN - IF(NEWTYP.EQ.1) WRITE(IUNITS,1710) - 1710 FORMAT(/' ==> OBS ERROR VALUES ARE ENCODED FOR THIS TABLE A ', - $ 'ENTRY'//' ==> FILTERING VIA MISSING OBS ERROR TEST IS ', - $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE OBS ERROR VALUES ARE ', - $ 'PROCESSED/STORED'/) - CALL GBLEVN04 - IF(NLEV.GT.0) CALL UFBINT(IUNITP,BAK_8,12,NLEV,IRET,OESTR) - ELSE - IF(NEWTYP.EQ.1) WRITE(IUNITS,1705) - 1705 FORMAT(/' ==> OBS ERROR VALUES NOT ENCODED FOR THIS TABLE A ', - $ 'ENTRY'//' ==> FILTERING VIA MISSING OBS ERROR TEST NOT ', - $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE OBS ERROR VALUES NOT ', - $ 'PROCESSED/STORED'/) - ENDIF - -C MAKE THE GBLEVENTS EVENTS AND ENCODE INTO REPORT -C ------------------------------------------------ - - IF(.NOT.FCST) THEN - IF(NEWTYP.EQ.1) WRITE(IUNITS,1704) - 1704 FORMAT(/' ==> FORECAST VALUES NOT ENCODED FOR THIS TABLE A ', - $ 'ENTRY'//' ==> FILTERING VIA POB VS. GESS PSFC TEST NOT ', - $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE FORECAST VALUES NOT ', - $ 'PROCESSED/STORED'/) - ELSE - IF(NEWTYP.EQ.1) WRITE(IUNITS,1708) - 1708 FORMAT(/' ==> FORECAST VALUES ARE ENCODED FOR THIS TABLE A ', - $ 'ENTRY'//' ==> FILTERING VIA POB VS. GESS PSFC TEST IS ', - $ 'PERFORMED FOR THIS TABLE A ENTRY SINCE FORECAST VALUES ARE ', - $ 'PROCESSED/STORED'/) - ENDIF - - IF(DOPREV) THEN - IF(NEWTYP.EQ.1) WRITE(IUNITS,1807) - 1807 FORMAT(/' ==> "PREVENT" EVENT PROCESSING IS PERFORMED FOR THIS', - $ ' TABLE A ENTRY'/) - CALL GBLEVN02(IUNITP,IUNITS,NEWTYP,subset) - ELSE - IF(NEWTYP.EQ.1) WRITE(IUNITS,1806) - 1806 FORMAT(/' ==> "PREVENT" EVENT PROCESSING NOT PERFORMED FOR THIS', - $ ' TABLE A ENTRY'/) - ENDIF - -C MAKE THE VIRTUAL TEMPERATURE EVENTS AND ENCODE INTO REPORT -C ---------------------------------------------------------- - - IF(.NOT.VIRT) THEN - IF(NEWTYP.EQ.1) WRITE(IUNITS,1706) - 1706 FORMAT(/' ==> "VIRTMP" EVENT PROCESSING NOT PERFORMED FOR THIS ', - $ 'TABLE A ENTRY'/) - ELSE - IF(NEWTYP.EQ.1) WRITE(IUNITS,1707) - 1707 FORMAT(/' ==> "VIRTMP" EVENT PROCESSING IS PERFORMED FOR THIS ', - $ 'TABLE A ENTRY'/) - CALL GBLEVN08(IUNITP,SUBSET) - ENDIF - -C RETURN TO CALLING PROGRAM TO WRITE GBL-EVENTED REPORT (SUBSET) INTO -C PREPBUFR FILE -C ------------------------------------------------------------------- - - RETURN - - END -C*********************************************************************** -C*********************************************************************** - SUBROUTINE GBLEVN01(IUNITE) ! FORMERLY SUBROUTINE ETABLE - - COMMON /GBEVDD/ ERRS(300,33,6) - -C READ THE OBSERVATION ERROR TABLES -C --------------------------------- - - REWIND IUNITE - - IREC = 0 - - 10 CONTINUE - READ(IUNITE,'(1X,I3)',END=100) KX - IREC = IREC + 1 - DO K=1,33 - READ(IUNITE,'(1X,6E12.5)') (ERRS(KX,K,M),M=1,6) - ENDDO - GO TO 10 - - 100 CONTINUE - IF(IREC.LE.0) THEN - PRINT'(" ##GBLEVENTS/GBLEVN01 - OBS. ERROR TABLE EMPTY OR ", - $ "DOES NOT EXIST - STOP 60")' - CALL ERREXIT(60) - ENDIF - - RETURN - - END -C*********************************************************************** -C*********************************************************************** - SUBROUTINE GBLEVN02(IUNITP,IUNITS,NEWTYP,subset) - ! FORMERLY SUBROUTINE FILTAN - - DIMENSION NFLGRT(100:299,12),OEMIN(2:6) - CHARACTER*8 STNID,subset - CHARACTER*40 PEVN,QEVN,TEVN,WEVN,PWVN,PW1VN,PW2VN,PW3VN,PW4VN - REAL(8) PEV_8(4,255),QEV_8(4,255),TEV_8(4,255),WEV_8(5,255), - $ PWV_8(4,255),PW1V_8(4,255),PW2V_8(4,255), - $ PW3V_8(4,255),PW4V_8(4,255),OBS_8,QMS_8,BAK_8,SID_8, - $ ufc_8,vfc_8 - LOGICAL FCST,REJP_PS,REJPS,REJT,REJQ,REJW,REJPW,REJPW1, - $ REJPW2,REJPW3,REJPW4,SATMQC,SATEMP,SOLN60,SOLS60, - $ MOERR_P,MOERR_T,ADPUPA_VIRT,DOBERR,DOFCST,SOME_FCST, - $ DOVTMP,VIRT,RECALC_Q,DOPREV - REAL(8) BMISS - - COMMON /GBEVAA/ SID_8,OBS_8(13,255),QMS_8(12,255),BAK_8(12,255), - $ XOB,YOB,DHR,TYP,NLEV - COMMON /GBEVBB/ PVCD,VTCD - COMMON /GBEVCC/ DOVTMP,DOFCST,SOME_FCST,DOBERR,FCST,VIRT, - $ QTOP_REJ,SATMQC,ADPUPA_VIRT,RECALC_Q,DOPREV,dopmsl - COMMON /GBEVEE/PSG01,ZSG01,TG01(500),UG01(500),VG01(500), - x QG01(500),zint(500),pint(500),pintlog(500),plev(500), - x plevlog(500) - COMMON /GBEVFF/ BMISS - - EQUIVALENCE (SID_8,STNID) - - DATA PEVN /'POB PQM PPC PRC '/ - DATA QEVN /'QOB QQM QPC QRC '/ - DATA TEVN /'TOB TQM TPC TRC '/ - DATA WEVN /'UOB VOB WQM WPC WRC '/ - DATA PWVN /'PWO PWQ PWP PWR '/ - DATA PW1VN /'PW1O PW1Q PW1P PW1R '/ - DATA PW2VN /'PW2O PW2Q PW2P PW2R '/ - DATA PW3VN /'PW3O PW3Q PW3P PW3R '/ - DATA PW4VN /'PW4O PW4Q PW4P PW4R '/ - - DATA NFLGRT/2400*0/ - - DATA OEMIN /0.5,0.1,1.0,0.5,1.0/ - - NI = MOD((NINT(TYP)/10),10) - - IF(NEWTYP.EQ.1) NFLGRT = 0 - -C LOGICAL SWITCHES FOR OBSERVATION LOCATION FILTERING -C --------------------------------------------------- - - SATEMP = ((TYP.GE.160.AND.TYP.LE.179).AND.SATMQC) - SOLN60 = ((TYP.GE.160.AND.TYP.LE.163).AND.YOB.GE.-60.AND.SATMQC) - SOLS60 = ((TYP.EQ.160.OR.TYP.EQ.162.OR.TYP.EQ.163).AND.YOB.LT.-60 - $ .AND.SATMQC) - -C CLEAR THE EVENT ARRAYS -C ---------------------- - - PEV_8 = BMISS - QEV_8 = BMISS - TEV_8 = BMISS - WEV_8 = BMISS - PWV_8 = BMISS - PW1V_8 = BMISS - PW2V_8 = BMISS - PW3V_8 = BMISS - PW4V_8 = BMISS - - MAXPEV = 0 - MAXQEV = 0 - MAXTEV = 0 - MAXWEV = 0 - MAXPWV = 0 - MAXPW1V = 0 - MAXPW2V = 0 - MAXPW3V = 0 - MAXPW4V = 0 - -C LOOP OVER LEVELS APPLYING UNDERGROUND FILTERING AND SPECIAL RULES -C ----------------------------------------------------------------- - - IF(NLEV.GT.0) THEN - DO L=1,NLEV - - POB = OBS_8( 1,L) - QOB = OBS_8( 2,L) - TOB = OBS_8( 3,L) - UOB = OBS_8( 5,L) - VOB = OBS_8( 6,L) - PWO = OBS_8( 7,L) - PW1O = OBS_8( 8,L) - PW2O = OBS_8( 9,L) - PW3O = OBS_8(10,L) - PW4O = OBS_8(11,L) - CAT = OBS_8(12,L) - PRSS = OBS_8(13,L) - - PQM = QMS_8( 1,L) - QQM = QMS_8( 2,L) - TQM = QMS_8( 3,L) - ZQM = QMS_8( 4,L) - WQM = QMS_8( 5,L) - PWQ = QMS_8( 6,L) - PW1Q = QMS_8( 7,L) - PW2Q = QMS_8( 8,L) - PW3Q = QMS_8( 9,L) - PW4Q = QMS_8(10,L) - - REJP_PS = .FALSE. - MOERR_P = .FALSE. - MOERR_T = .FALSE. - RCD = 99999 - -C ------------------------------------------------------------------- -C RULES FOR PRESSURE (ON ANY LEVEL) -- ALL DATA (MASS AND WIND) ON -C LEVEL REJECTED IF: -C - PRESSURE MORE THAN 100 MB BELOW MODEL (GUESS) SURFACE PRESSURE -C (AND SWITCH FCST=TRUE) -- "PREVENT" PGM REASON CODE 1 -C - PRESSURE IS ZERO OR IS NEGATIVE -- "PREVENT" PGM REASON CODE 2 -C - SURFACE (MASS OR WIND) REPORT PRESSURE IS REPORTED ABOVE 450 MB -C OR BELOW 1100 MB -- "PREVENT" PGM REASON CODE 2 -C REJECTION MEANS Q.M. SET TO 8 -C ------------------------------------------------------------------- - - IF(POB.LT.BMISS) THEN - IF(.NOT.FCST) PSG01 = POB - IF(POB-PSG01.GE.100. .OR. POB.LE.0. .OR. - $ ((POB.LE.450..OR.POB.GE.1100.) .AND. NI.EQ.8)) THEN - IF(POB.LE.0..OR.POB.LE.450..OR.POB.GE.1100.) THEN - IF(NI.EQ.8) THEN - WRITE(IUNITS,302) STNID,NINT(TYP),YOB,XOB,POB - 302 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, REJECT ALL DATA ON SFC LVL - POB=',F6.1,' MB, FAILS SANITY ', - $ 'CHECK') - ELSE - WRITE(IUNITS,101) STNID,NINT(TYP),YOB,XOB,POB - 101 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $'E, REJECT ALL DATA ON LVL - POB=',F6.1,' MB, FAILS SANITY CHECK') - ENDIF - RCD = 2 - ELSE - IF(NI.EQ.8) THEN - WRITE(IUNITS,303) STNID,NINT(TYP),YOB,XOB,POB,PSG01 - 303 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, REJECT ALL DATA ON SFC LVL - POB=',F6.1,' MB, > 100 MB ', - $ 'BELOW GES PSFC(=',F6.1,'MB)') - ELSE - WRITE(IUNITS,102) STNID,NINT(TYP),YOB,XOB,POB,PSG01 - 102 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, REJECT ALL DATA ON LVL - POB=',F6.1,' MB, > 100 MB BELOW ', - $ 'GES PSFC(=',F6.1,' MB)') - ENDIF - RCD = 1 - ENDIF - REJ = 8 - REJP_PS = .TRUE. - PEV_8(1,L) = POB - PEV_8(2,L) = REJ - PEV_8(3,L) = PVCD - PEV_8(4,L) = RCD - MAXPEV = L - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR SURFACE PRESSURE -- ALL MASS DATA ON SURFACE LEVEL -C REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C - PRESSURE IS MORE THAN 100 MB ABOVE OR BELOW MODEL (GUESS) -C SURFACE PRESSURE (AND SWITCH FCST=TRUE) -- -C "PREVENT" PGM REASON CODE 4 -C - PRESSURE IS REPORTED ABOVE 450 MB OR BELOW 1100 MB -- "PREVENT" -C PGM REASON CODE 2 (NOTE: DOES NOT APPLY TO SURFACE REPORTS, -C THESE WERE TESTED FOR THIS CRITERION IN ABOVE PRESSURE TEST) -C - PRESSURE VIOLATES RULES FOR PRESSURE ON ANY LEVEL (SEE ABOVE) -C REJECTION FOR FIRST RULE MEANS Q.M. SET TO 9 UNLESS: -C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8 -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(POB.LT.BMISS .AND. CAT.EQ.0) THEN - IF(.NOT.FCST) PSG01 = POB - REJPS = OEFG01(POB,TYP,5,OEMIN(5)).GE.BMISS .OR. - $ ABS(POB-PSG01).GE.100. .OR. - $ POB.LE.450. .OR. - $ POB.GE.1100. - IF(REJPS.OR.REJP_PS) THEN - REJ = 8 - IF(.NOT.REJP_PS) THEN - IF(ABS(POB-PSG01).GE.100.) THEN - WRITE(IUNITS,104) STNID,NINT(TYP),YOB,XOB,POB,PSG01 - 104 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, REJECT ALL DATA ON SFC LVL - POB=',F6.1,' MB, > 100 MB ', - $ 'ABOVE GES PSFC(=',F6.1,'MB)') - RCD = 4 - ELSE IF(POB.LE.450..OR.POB.GE.1100.) THEN - WRITE(IUNITS,105) STNID,NINT(TYP),YOB,XOB,POB - 105 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, REJECT ALL DATA ON SFC LVL - POB=',F6.1,' MB, FAILS SANITY ', - $ 'CHECK - this should never be printed since test now made in ', - $ 'section above ') - RCD = 2 - ELSE - IF(NFLGRT(NINT(TYP),1).EQ.0) THEN - WRITE(IUNITS,201) NINT(TYP) - 201 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'ALL DATA ON SURFACE LEVEL DUE TO MISSING SFC-P OBS ERROR'/) - NFLGRT(NINT(TYP),1) = 1 - ENDIF -CDAK CDAK CDAK CDAK WRITE(IUNITS,103) STNID,NINT(TYP),YOB,XOB,POB -CD103 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -CDAK $ 'E, REJECT ALL DATA ON SFC LVL - POB=',F6.1,'MB, MISSING OBS.', -CDAK $ ' ERROR') - RCD = 3 - REJ = 9 - ENDIF - ENDIF - REJP_PS = .TRUE. - IF(RCD.EQ.3) MOERR_P = .TRUE. - IF(REJ.EQ.9.AND.(PQM.GT.3.AND.PQM.LT.15)) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1401) STNID,NINT(TYP),YOB,XOB,PQM - 1401 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT PQM =',F4.0,' -- DO NOT APPLY PSFC QM=9 EVENT') - ELSE - PEV_8(1,L) = POB - PEV_8(2,L) = REJ - PEV_8(3,L) = PVCD - PEV_8(4,L) = RCD - MAXPEV = L - ENDIF - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR TEMPERATURE -- TOB AND QOB ON LEVEL REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS -C MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3 -C - REPORT IS TYPE 160-163 (LAND TOVS/RTOVS/ATOVS TEMPERATURE -C SOUNDINGS, ALL PATHS), AND IS AT OR NORTH OF 60 DEGREES SOUTH -C LATITUDE, AND PRESSURE ON LEVEL IS AT OR BELOW 100 MB (AND -C SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 6 -C - REPORT IS TYPE 160,162,163 (LAND TOVS/RTOVS/ATOVS TEMPERATURE -C SOUNDINGS, ALL PATHS BUT CLEAR), AND IS SOUTH OF 60 DEGREES -C SOUTH LATITUDE, AND PRESSURE ON LEVEL IS BELOW 100 MB (AND -C SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 6 -C - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE -C (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE -C 2 ABOVE) (SEE ABOVE) -C - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE) -C REJECTION FOR FIRST TWO RULES MEANS Q.M. SET TO 9 UNLESS: -C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8 -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(TOB.LT.BMISS) THEN - REJT = OEFG01(POB,TYP,2,OEMIN(2)).GE.BMISS .OR. - $ (SOLN60.AND.NINT(POB*10.).GE.1000) .OR. - $ (SOLS60.AND.NINT(POB*10.).GT.1000) - IF(REJT.OR.REJP_PS) THEN - REJ = 8 - IF(.NOT.REJP_PS) THEN - IF(SOLN60.AND.NINT(POB*10.).GE.1000) THEN - IF(NFLGRT(NINT(TYP),6).EQ.0) THEN - WRITE(IUNITS,7304) NINT(TYP) - 7304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $'TOB/QOB AT AND BELOW 100 MB IF REPORT IS NORTH OF 60S LATITUDE'/) - NFLGRT(NINT(TYP),6) = 1 - ENDIF - RCD = 6 - ELSE IF(SOLS60.AND.NINT(POB*10.).GT.1000) THEN - IF(NFLGRT(NINT(TYP),7).EQ.0) THEN - WRITE(IUNITS,7305) NINT(TYP) - 7305 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'TOB/QOB BELOW 100 MB IF REPORT IS SOUTH OF 60S LATITUDE'/) - NFLGRT(NINT(TYP),7) = 1 - ENDIF - RCD = 6 - ELSE - IF(NFLGRT(NINT(TYP),2).EQ.0) THEN - IF(NI.EQ.8) THEN - WRITE(IUNITS,304) NINT(TYP) - 304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'TOB/QOB ON SFC LVL DUE TO MISSING OBS ERROR'/) - ELSE - WRITE(IUNITS,202) NINT(TYP) - 202 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'TOB/QOB ON AT LEAST ONE LVL (IF AVAILABLE ON THAT LVL) DUE TO ', - $ 'MISSING OBS ERROR'/) - ENDIF - NFLGRT(NINT(TYP),2) = 1 -cdak cdak cdak cdak cdakWRITE(IUNITS,106) STNID,NINT(TYP),YOB,XOB,TOB -cd106 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $ 'E, REJECT TOB/QOB ON LVL - TOB=',F5.1,'C, MISSING OBS. ERROR') - ENDIF - RCD = 3 - REJ = 9 - ENDIF - ELSE - IF(MOERR_P) THEN - RCD = 3 - REJ = 9 - ENDIF - ENDIF - IF(RCD.EQ.3) MOERR_T = .TRUE. - IF(REJ.EQ.9.AND.(TQM.GT.3.AND.TQM.LT.15)) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1402) STNID,NINT(TYP),YOB,XOB,TQM - 1402 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT TQM =',F4.0,' -- DO NOT APPLY TEMP QM=9 EVENT') - ELSE - TEV_8(1,L) = TOB - TEV_8(2,L) = REJ - TEV_8(3,L) = PVCD - TEV_8(4,L) = RCD - MAXTEV = L - ENDIF - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR SPECIFIC HUMIDITY -- QOB ON LEVEL REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C - PRESSURE ON LEVEL IS ABOVE "QTOP_REJ" MB {WHERE QTOP_REJ IS -C READ IN FROM NAMELIST "PREVDATA" (SEE DOCBLOCK)} -- "PREVENT" -C PGM REASON CODE 5 -C - OBSERVATION ERROR FOR TEMPERATURE ON LEVEL IS MISSING (AND -C SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3 -C - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS -C MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3 -C - TEMPERATURE ON LEVEL IS MISSING OR IS LESS THAN -150 DEG. C -- -C "PREVENT" PGM REASON CODE 2 -C - SPECIFIC HUMIDITY IS ZERO OR IS NEGATIVE -- "PREVENT" PGM REASON -C CODE 2 -C - REPORT IS TYPE 160-179 (SATELLITE TEMPERATURE SOUNDINGS, ALL -C TYPES, ALL PATHS, LAND AND SEA), ALL PRESSURE LEVELS (AND -C SWITCH SATMQC=TRUE) -- "PREVENT" PGM REASON CODE 7 -C - TEMPERATURE ON LEVEL VIOLATES RULES FOR TEMPERATURE (EXCEPT FOR -C MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE 2 ABOVE) -C (SEE ABOVE) -C - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE -C (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE -C 4 ABOVE) (SEE ABOVE) -C - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE) -C REJECTION FOR FIRST FOUR RULES MEANS Q.M. SET TO 9 UNLESS: -C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8 -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(QOB.LT.BMISS) THEN - REJQ = OEFG01(POB,TYP,3,OEMIN(3)).GE.BMISS .OR. - $ TOB.GE.BMISS .OR. - $ TOB.LE.-150. .OR. - $ NINT(POB * 10.).LT.NINT(QTOP_REJ * 10.) .OR. - $ QOB.LE.0. .OR. - $ SATEMP .OR. - $ REJT - IF(REJQ.OR.REJP_PS) THEN - REJ = 8 - IF(.NOT.REJP_PS.AND..NOT.REJT) THEN - IF(SATEMP) THEN - IF(NFLGRT(NINT(TYP),8).EQ.0) THEN - WRITE(IUNITS,7306) NINT(TYP) - 7306 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'QOB ON ALL LEVELS'/) - NFLGRT(NINT(TYP),8) = 1 - ENDIF - RCD = 7 - ELSE IF(QOB.LE.0..OR.TOB.GE.BMISS.OR.TOB.LE.-150.)THEN - WRITE(IUNITS,111) STNID,NINT(TYP),YOB,XOB, - $ QOB/1000.,TOB - 111 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, REJECT QOB ON LVL - QOB=',F6.3,' G/KG, FAILS SANITY CHECK ', - $ '(TOB=',F5.1,' C)') - RCD = 2 - ELSE IF(OEFG01(POB,TYP,3,OEMIN(3)).GE.BMISS) THEN - IF(NFLGRT(NINT(TYP),3).EQ.0) THEN - IF(NI.EQ.8) THEN - WRITE(IUNITS,305) NINT(TYP) - 305 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'QOB ON SFC LVL DUE TO MISSING OBS ERROR'/) - ELSE - WRITE(IUNITS,203) NINT(TYP) - 203 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'QOB ON AT LEAST ONE LEVEL (IF AVAILABLE ON THAT LEVEL) DUE TO ', - $ 'MISSING OBS ERROR'/) - ENDIF - NFLGRT(NINT(TYP),3) = 1 - ENDIF - RCD = 3 - REJ = 9 -cdak cdak cdak cdak WRITE(IUNITS,108) STNID,NINT(TYP),YOB,XOB,QOB/1000. -cd108 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $'E, REJECT QOB ON LVL - QOB=',F6.3,'G/KG, MISSING OBS. ERROR') - ELSE - WRITE(IUNITS,109) STNID,NINT(TYP),YOB,XOB, - $ QOB/1000.,QTOP_REJ,POB - 109 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, REJECT QOB ON LVL - QOB=',F6.3,' G/KG, ABOVE ',F6.1, - $ 'MB (POB=',F6.1,' MB)') - RCD = 5 - REJ = 9 - ENDIF - ELSE - IF(MOERR_P.OR.MOERR_T) THEN - RCD = 3 - REJ = 9 - ENDIF - ENDIF - IF(REJ.EQ.9.AND.(QQM.GT.3.AND.QQM.LT.15)) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1403) STNID,NINT(TYP),YOB,XOB,QQM - 1403 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT QQM =',F4.0,' -- DO NOT APPLY MSTR QM=9 EVENT') - ELSE - QEV_8(1,L) = QOB - QEV_8(2,L) = REJ - QEV_8(3,L) = PVCD - QEV_8(4,L) = RCD - MAXQEV = L - ENDIF - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR WIND -- UOB AND VOB ON LEVEL REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C - THIS IS SFC LEVEL AND OBSERVATION ERROR FOR SFC PRESSURE IS -C MISSING (AND SWITCH DOBERR=TRUE) -- "PREVENT" PGM REASON CODE 3 -C (Note: This can currently never happen because earlier check -C for missing obs error for sfc pressure is only done if -C "surface" level is category 0 and this is not possible -C for wind reports.) -C - THIS IS SFC LEVEL AND PRESSURE VIOLATES RULES FOR SFC PRESSURE -C (EXCEPT FOR MISSING OBSERVATION ERROR, ALREADY COVERED IN RULE -C 2 ABOVE) (SEE ABOVE) -C - This is a SFCSHP report with calm winds and non-missing -C background u- or v-component wind is .GE. 5 m/sec -- -C "PREVENT" PGM REASON CODE 8, Q.M. set to 8 -C - PRESSURE ON LEVEL VIOLATES RULES FOR PRESSURE (SEE ABOVE) -C REJECTION FOR FIRST TWO RULES MEANS Q.M. SET TO 9 UNLESS: -C - ANY OTHER RULE CAUSES REJECTION, THEN Q.M. SET TO 8 -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(MIN(UOB,VOB).LT.BMISS) THEN - REJW = OEFG01(POB,TYP,4,OEMIN(4)).GE.BMISS - IF(REJW.OR.REJP_PS) THEN - REJ = 8 - IF(.NOT.REJP_PS) THEN - IF(NFLGRT(NINT(TYP),4).EQ.0) THEN - IF(NI.EQ.8) THEN - WRITE(IUNITS,1304) NINT(TYP) - 1304 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'UOB/VOB ON SFC LVL DUE TO MISSING OBS ERROR'/) - ELSE - WRITE(IUNITS,204) NINT(TYP) - 204 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'UOB/VOB ON AT LEAST ONE LVL (IF AVAILABLE ON THAT LVL) DUE TO ', - $ 'MISSING OBS ERROR'/) - ENDIF - NFLGRT(NINT(TYP),4) = 1 - ENDIF -cdak cdak cdak WRITE(IUNITS,112) STNID,NINT(TYP),YOB,XOB -cd112 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $'E, REJECT UOB/VOB ON LVL - MISSING OBS. ERROR') - RCD = 3 - REJ = 9 - ELSE - IF(MOERR_P) THEN ! This currently can never be TRUE - ! since CAT is never 0 for "sfc" - ! level of wind reports - RCD = 3 - REJ = 9 - ENDIF - ENDIF - IF(REJ.EQ.9.AND.(WQM.GT.3.AND.WQM.LT.15)) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1404) STNID,NINT(TYP),YOB,XOB,WQM - 1404 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT WQM =',F4.0,' -- DO NOT APPLY WIND QM=9 EVENT') - ELSE - WEV_8(1,L) = UOB - WEV_8(2,L) = VOB - WEV_8(3,L) = REJ - WEV_8(4,L) = PVCD - WEV_8(5,L) = RCD - MAXWEV = L - ENDIF - ENDIF -c Reject calm sfc marine winds if background u- or v- wind .ge. 5 m/s - if(subset.eq."SFCSHP".and.uob.eq.0..and.vob.eq.0.) then - call ufbint(-iunitp,ufc_8,1,1,iret,'UFC') - call ufbint(-iunitp,vfc_8,1,1,iret,'VFC') -cccccccc if(ibfms(ufc_8).eq.0.or.ibfms(vfc_8).eq.0) then - ! DAK: changed to only allow this test if - ! both UFC and VFC are non-missing - if(ibfms(ufc_8).eq.0.and.ibfms(vfc_8).eq.0) then - if(abs(ufc_8).ge.5..or.abs(vfc_8).ge.5.) then - if(wqm.le.3) then - write(iunits,1504) stnid,nint(typ),ufc_8,vfc_8 - 1504 FORMAT(/' --> ID ',A8,', (RTP ',I3,' UFC=',F5.2,' VFC=',F5.2,') ', - $ 'NEW WIND EVENT, UOB/VOB CALM (0 M/S) WHILE |UFC| AND/OR |VFC| ', - $ '>= 5 M/S, QM SET TO 8'/) - wev_8(1,1) = uob - wev_8(2,1) = vob - wev_8(3,1) = 8 - wev_8(4,1) = pvcd - wev_8(5,1) = 8 - maxwev = 1 - end if - end if - end if - end if - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR TOTAL COLUMN PRECIPITABLE WATER -- PWO REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C REJECTION MEANS Q.M. SET TO 9 UNLESS: -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(PWO.LT.BMISS) THEN - REJPW = OEFG01(PRSS*0.01,TYP,6,OEMIN(6)).GE.BMISS - IF(REJPW) THEN - IF(NFLGRT(NINT(TYP),5).EQ.0) THEN - WRITE(IUNITS,205) NINT(TYP) - 205 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'PWO DUE TO MISSING OBS ERROR'/) - NFLGRT(NINT(TYP),5) = 1 - ENDIF -cdakcdakcdak WRITE(IUNITS,113) STNID,NINT(TYP),YOB,XOB,PWO -cd113 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $'E, REJECT PWO ON LVL - PWO=',F5.1,'MM, MISSING OBS. ERROR') - IF(PWQ.GT.3.AND.PWQ.LT.15) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1405) STNID,NINT(TYP),YOB,XOB,PWQ - 1405 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT PWQ =',F4.0,' -- DO NOT APPLY PWtO QM=9 EVENT') - ELSE - PWV_8(1,L) = PWO - PWV_8(2,L) = 9 - PWV_8(3,L) = PVCD - PWV_8(4,L) = 3 - MAXPWV = L - ENDIF - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR LAYER 1 PRECIPITABLE WATER -- PW1O REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C REJECTION MEANS Q.M. SET TO 9 UNLESS: -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(PW1O.LT.BMISS) THEN - REJPW1 = OEFG01(PRSS*0.01,TYP,6,OEMIN(6)).GE.BMISS - IF(REJPW1) THEN - IF(NFLGRT(NINT(TYP),9).EQ.0) THEN - WRITE(IUNITS,206) NINT(TYP) - 206 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'PW1O DUE TO MISSING OBS ERROR'/) - NFLGRT(NINT(TYP),9) = 1 - ENDIF -cdakcdakcdak WRITE(IUNITS,114) STNID,NINT(TYP),YOB,XOB,PW1O -cd114 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $'E, REJECT PW1O ON LVL - PW1O=',F5.1,'MM, MISSING OBS. ERROR') - IF(PW1Q.GT.3.AND.PW1Q.LT.15) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1406) STNID,NINT(TYP),YOB,XOB,PW1Q - 1406 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT PW1Q =',F4.0,' -- DO NOT APPLY PW1O QM=9 EVENT') - ELSE - PW1V_8(1,L) = PW1O - PW1V_8(2,L) = 9 - PW1V_8(3,L) = PVCD - PW1V_8(4,L) = 3 - MAXPW1V = L - ENDIF - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR LAYER 2 PRECIPITABLE WATER -- PW2O REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C REJECTION MEANS Q.M. SET TO 9 UNLESS: -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(PW2O.LT.BMISS) THEN - REJPW2 = OEFG01(PRSS*0.01,TYP,6,OEMIN(6)).GE.BMISS - IF(REJPW2) THEN - IF(NFLGRT(NINT(TYP),10).EQ.0) THEN - WRITE(IUNITS,207) NINT(TYP) - 207 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'PW2O DUE TO MISSING OBS ERROR'/) - NFLGRT(NINT(TYP),10) = 1 - ENDIF -cdakcdakcdak WRITE(IUNITS,115) STNID,NINT(TYP),YOB,XOB,PW2O -cd115 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $'E, REJECT PW2O ON LVL - PW2O=',F5.1,'MM, MISSING OBS. ERROR') - IF(PW2Q.GT.3.AND.PW2Q.LT.15) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1407) STNID,NINT(TYP),YOB,XOB,PW2Q - 1407 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT PW2Q =',F4.0,' -- DO NOT APPLY PW2O QM=9 EVENT') - ELSE - PW2V_8(1,L) = PW2O - PW2V_8(2,L) = 9 - PW2V_8(3,L) = PVCD - PW2V_8(4,L) = 3 - MAXPW2V = L - ENDIF - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR LAYER 3 PRECIPITABLE WATER -- PW3O REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C REJECTION MEANS Q.M. SET TO 9 UNLESS: -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(PW3O.LT.BMISS) THEN - REJPW3 = OEFG01(PRSS*0.01,TYP,6,OEMIN(6)).GE.BMISS - IF(REJPW3) THEN - IF(NFLGRT(NINT(TYP),11).EQ.0) THEN - WRITE(IUNITS,208) NINT(TYP) - 208 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'PW3O DUE TO MISSING OBS ERROR'/) - NFLGRT(NINT(TYP),11) = 1 - ENDIF -cdakcdakcdak WRITE(IUNITS,116) STNID,NINT(TYP),YOB,XOB,PW3O -cd116 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $'E, REJECT PW3O ON LVL - PW3O=',F5.1,'MM, MISSING OBS. ERROR') - IF(PW3Q.GT.3.AND.PW3Q.LT.15) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1408) STNID,NINT(TYP),YOB,XOB,PW3Q - 1408 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT PW3Q =',F4.0,' -- DO NOT APPLY PW3O QM=9 EVENT') - ELSE - PW3V_8(1,L) = PW3O - PW3V_8(2,L) = 9 - PW3V_8(3,L) = PVCD - PW3V_8(4,L) = 3 - MAXPW3V = L - ENDIF - ENDIF - ENDIF - -C ------------------------------------------------------------------- -C RULES FOR LAYER 4 PRECIPITABLE WATER -- PW4O REJECTED IF: -C - OBSERVATION ERROR IS MISSING (AND SWITCH DOBERR=TRUE) -- -C "PREVENT" PGM REASON CODE 3 -C REJECTION MEANS Q.M. SET TO 9 UNLESS: -C - Q.M. IS ALREADY > 3 BUT NOT 15, THEN SKIP EVENT AND LEAVE Q.M. -C AS IS -C ------------------------------------------------------------------- - - IF(PW4O.LT.BMISS) THEN - REJPW4 = OEFG01(PRSS*0.01,TYP,6,OEMIN(6)).GE.BMISS - IF(REJPW4) THEN - IF(NFLGRT(NINT(TYP),12).EQ.0) THEN - WRITE(IUNITS,209) NINT(TYP) - 209 FORMAT(/' --> FOR ALL REPORTS WITH REPORT TYPE ',I3,', REJECT ', - $ 'PW4O DUE TO MISSING OBS ERROR'/) - NFLGRT(NINT(TYP),12) = 1 - ENDIF -cdakcdakcdak WRITE(IUNITS,117) STNID,NINT(TYP),YOB,XOB,PW4O -cd117 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, -cdak $'E, REJECT PW4O ON LVL - PW4O=',F5.1,'MM, MISSING OBS. ERROR') - IF(PW4Q.GT.3.AND.PW4Q.LT.15) THEN -CDAKCDAKCDAKCDAK WRITE(IUNITS,1409) STNID,NINT(TYP),YOB,XOB,PW4Q - 1409 FORMAT(' ~~> ID ',A8,' (RTP ',I3,'), LAT=',F6.2,'N, LON=',F6.2, - $ 'E, INPUT PW4Q =',F4.0,' -- DO NOT APPLY PW4O QM=9 EVENT') - ELSE - PW4V_8(1,L) = PW4O - PW4V_8(2,L) = 9 - PW4V_8(3,L) = PVCD - PW4V_8(4,L) = 3 - MAXPW4V = L - ENDIF - ENDIF - ENDIF - - ENDDO - ENDIF - -C APPLY THE PROPER EVENTS -C ----------------------- - - IF(MAXPEV .GT.0) CALL UFBINT(IUNITP,PEV_8, 4,MAXPEV, IRET,PEVN) - IF(MAXQEV .GT.0) CALL UFBINT(IUNITP,QEV_8, 4,MAXQEV, IRET,QEVN) - IF(MAXTEV .GT.0) CALL UFBINT(IUNITP,TEV_8, 4,MAXTEV, IRET,TEVN) - IF(MAXWEV .GT.0) CALL UFBINT(IUNITP,WEV_8, 5,MAXWEV, IRET,WEVN) - IF(MAXPWV .GT.0) CALL UFBINT(IUNITP,PWV_8, 4,MAXPWV, IRET,PWVN) - IF(MAXPW1V.GT.0) CALL UFBINT(IUNITP,PW1V_8,4,MAXPW1V,IRET,PW1VN) - IF(MAXPW2V.GT.0) CALL UFBINT(IUNITP,PW2V_8,4,MAXPW2V,IRET,PW2VN) - IF(MAXPW3V.GT.0) CALL UFBINT(IUNITP,PW3V_8,4,MAXPW3V,IRET,PW3VN) - IF(MAXPW4V.GT.0) CALL UFBINT(IUNITP,PW4V_8,4,MAXPW4V,IRET,PW4VN) - - RETURN - END -C*********************************************************************** -C*********************************************************************** -C GBLEVN03 - INTERPOLATE MODEL DATA (FIRST GUESS OR ANALYSIS) TO OB -C LOCATIONS -C----------------------------------------------------------------------- - SUBROUTINE GBLEVN03(SUBSET) ! FORMERLY SUBROUTINE GETFC - - USE GBLEVN_MODULE - - REAL(8) OBS_8,QMS_8,BAK_8,SID_8 - REAL(8) BMISS - CHARACTER*8 SUBSET - - COMMON /GBEVAA/ SID_8,OBS_8(13,255),QMS_8(12,255),BAK_8(12,255), - $ XOB,YOB,DHR,TYP,NLEV - COMMON /GBEVEE/PSG01,ZSG01,TG01(500),UG01(500),VG01(500), - x QG01(500),zint(500),pint(500),pintlog(500),plev(500), - x plevlog(500) - COMMON /GBEVFF/ BMISS - - - DATA TZERO / 273.15 / - DATA BETAP / .0552 / - DATA BETA / .00650 / - DATA ROG / 29.261 / - -C CLEAR THE BACKGROUND EVENT ARRAY -C -------------------------------- - - BAK_8 = BMISS - -C GET GUESS PROFILE AT OB LOCATION -C -------------------------------- - CALL GBLEVN06(XOB,YOB) - - -C INTERPOLATE GUESS PROFILES TO OB PRESSURES -C ------------------------------------------ - - IF(NLEV.GT.0) THEN - DO 10 L=1,NLEV - - POB = OBS_8( 1,L) - QOB = OBS_8( 2,L) - TOB = OBS_8( 3,L) - ZOB = OBS_8( 4,L) - UOB = OBS_8( 5,L) - VOB = OBS_8( 6,L) - PWO = OBS_8( 7,L) - PW1O = OBS_8( 8,L) - PW2O = OBS_8( 9,L) - PW3O = OBS_8(10,L) - PW4O = OBS_8(11,L) - CAT = OBS_8(12,L) - - IF(POB.LE.0. .OR. POB.GE.BMISS) GOTO 10 - - poblog = log(pob) - - la = -999 - lb = -999 - do k=1,kmax-1 - if (poblog<=plevlog(k) .and. poblog>plevlog(k+1)) then - la = k - lb = k+1 - exit - endif - enddo - if (la > 0) then - wt = (poblog-plevlog(lb)) / (plevlog(la)-plevlog(lb)) - else - la = 1 - lb = la+1 - wt = 0.0 - endif - - li=0 - do k=1,kmax-1 - if (poblog<=pintlog(k) .and. poblog>pintlog(k+1)) then - li = k - exit - endif - enddo - -C SURFACE PRESSURE -C ---------------- - - IF(CAT.EQ.0 .AND. ZOB.LT.BMISS) THEN - TS = TG01(1) + (PSG01-PLEV(1))*BETAP - DZ = ZOB-ZSG01 - TM = TS - DZ*BETA*.5 - PFC = PSG01*EXP(-DZ/(TM*ROG)) - ELSE - PFC = BMISS - ENDIF - -C SPECIFIC HUMIDITY -C ----------------- - - IF(QOB.LT.BMISS.OR.TOB.LT.BMISS.OR.TYP.EQ.111) THEN - -C (QFC NEEDED BY SYNDATA PROGRAM BUT ONLY FOR REPORT TYPE 111) - - QOB = QG01(LB) + (QG01(LA)-QG01(LB))*WT - ENDIF - -C TEMPERATURE -C ----------- - - IF(TOB.LT.BMISS.OR.SUBSET.EQ.'VADWND '.OR.TYP.EQ.111) THEN - -C (TFC NEEDED BY CQCVAD AND SYNDATA PROGRAMS, LATTER ONLY FOR REPORT -C TYPE 111) - - IF(POB.GT.PLEV(1)) THEN - TOB = TG01(1) + (POB-PLEV(1))*BETAP - ELSE - TOB = TG01(LB) + (TG01(LA)-TG01(LB))*WT - ENDIF - TOB = TOB - TZERO - ENDIF - -C HEIGHT -C ------ - - IF(ZOB.LT.BMISS) THEN - IF(POB.GT.PLEV(1)) THEN - TM = TG01(1) + (.5*(PINT(1)+POB)-PLEV(1))*BETAP - ZOB = ZINT(1) - ROG*TM*LOG(POB/PINT(1)) - ELSE - TM = TG01(LB) + (TG01(LA)-TG01(LB))*WT - ZOB = ZINT(LI) - ROG*TM*LOG(POB/PINT(LI)) - ENDIF - ENDIF - -C U AND V COMPONENTS -C ------------------ - - IF(UOB.LT.BMISS .OR. VOB.LT.BMISS) THEN - UOB = UG01(LB) + (UG01(LA)-UG01(LB))*WT - VOB = VG01(LB) + (VG01(LA)-VG01(LB))*WT - ENDIF - - -C PRECIPITABLE WATER -C ------------------ - - PWO = BMISS - PW1O = BMISS - PW2O = BMISS - PW3O = BMISS - PW4O = BMISS - -C RELATIVE HUMIDITY -C ----------------- - - RHO = BMISS - -C SCATTER THE PROPER FIRST GUESS/ANALYSIS VALUES -C ---------------------------------------------- - - BAK_8(1,L) = PFC - BAK_8(2,L) = QOB - BAK_8(3,L) = TOB - BAK_8(4,L) = ZOB - BAK_8(5,L) = UOB - BAK_8(6,L) = VOB - BAK_8(7,L) = PWO - BAK_8(8,L) = PW1O - BAK_8(9,L) = PW2O - BAK_8(10,L) = PW3O - BAK_8(11,L) = PW4O - BAK_8(12,L) = RHO - - 10 ENDDO - ENDIF - - RETURN - END -C*********************************************************************** -C*********************************************************************** - SUBROUTINE GBLEVN04 ! FORMERLY SUBROUTINE GETOE - - DIMENSION OEMIN(2:6) - REAL(8) OBS_8,QMS_8,BAK_8,SID_8 - REAL(8) BMISS - - COMMON /GBEVAA/ SID_8,OBS_8(13,255),QMS_8(12,255),BAK_8(12,255), - $ XOB,YOB,DHR,TYP,NLEV - COMMON /GBEVFF/ BMISS - - DATA OEMIN /0.5,0.1,1.0,0.5,1.0/ - -C CLEAR THE EVENT ARRAY -C --------------------- - - BAK_8 = BMISS - -C LOOP OVER LEVELS LOOKING UP THE OBSERVATION ERROR -C ------------------------------------------------- - - IF(NLEV.GT.0) THEN - DO L=1,NLEV - - POB = OBS_8( 1,L) - QOB = OBS_8( 2,L) - TOB = OBS_8( 3,L) - WOB = MAX(OBS_8(5,L),OBS_8(6,L)) - PWO = OBS_8( 7,L) - PW1O = OBS_8( 8,L) - PW2O = OBS_8( 9,L) - PW3O = OBS_8(10,L) - PW4O = OBS_8(11,L) - CAT = OBS_8(12,L) - - IF(CAT .EQ.0 ) BAK_8( 1,L) = OEFG01(POB,TYP,5,OEMIN(5)) - IF(QOB .LT.BMISS) BAK_8( 2,L) = OEFG01(POB,TYP,3,OEMIN(3)) - IF(TOB .LT.BMISS) BAK_8( 3,L) = OEFG01(POB,TYP,2,OEMIN(2)) - IF(WOB .LT.BMISS) BAK_8( 5,L) = OEFG01(POB,TYP,4,OEMIN(4)) - IF(PWO .LT.BMISS) BAK_8( 6,L) = OEFG01(POB,TYP,6,OEMIN(6)) - IF(PW1O.LT.BMISS) BAK_8( 7,L) = OEFG01(POB,TYP,6,OEMIN(6)) - IF(PW2O.LT.BMISS) BAK_8( 8,L) = OEFG01(POB,TYP,6,OEMIN(6)) - IF(PW3O.LT.BMISS) BAK_8( 9,L) = OEFG01(POB,TYP,6,OEMIN(6)) - IF(PW4O.LT.BMISS) BAK_8(10,L) = OEFG01(POB,TYP,6,OEMIN(6)) - - ENDDO - ENDIF - - RETURN - END -C*********************************************************************** -C*********************************************************************** -C SUBROUTINE GBLEVN06 - 2D LINEAR HORIZONTAL INTERPOLATION -C----------------------------------------------------------------------- - SUBROUTINE GBLEVN06(XOB,YOB) ! FORMERLY SUBROUTINE HTERP - - USE GBLEVN_MODULE - - COMMON /GBEVEE/ PSI,ZSI,TI(500),UI(500),VI(500),QI(500), - x zint(500),pint(500),pintlog(500),plev(500),plevlog(500) - - DATA ROG / 29.261 / - - -C CALCULATE HORIZONTAL WEIGHTS AND INTERPOLATE -C -------------------------------------------- - - WX = XOB/DLON + 1.0 - I0 = WX - I1 = MOD(I0,IMAX) + 1 - WX = WX-I0 - - WY = (YOB+90.)/DLAT + 1.0 - J0 = WY - J1 = MIN(J0+1,JMAX) - WY = WY-J0 - -C HTERP FOR SURFACE HEIGHT -C ------------------------ - - P1 = iar12z(I0,J0) - P2 = iar12z(I0,J1) - P3 = iar12z(I1,J0) - P4 = iar12z(I1,J1) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - ZSI = P5+(P6-P5)*WX - -C HTERP FOR SURFACE PRESSURE -C -------------------------- - - P1 = iar13p(I0,J0) - P2 = iar13p(I0,J1) - P3 = iar13p(I1,J0) - P4 = iar13p(I1,J1) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - PSI = P5+(P6-P5)*WX - -C HTERP FOR UPA T,U,V,Q -C --------------------- - - DO K=1,KMAX - - P1 = iar14t(I0,J0,K) - P2 = iar14t(I0,J1,K) - P3 = iar14t(I1,J0,K) - P4 = iar14t(I1,J1,K) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - TI(K) = P5+(P6-P5)*WX - - P1 = iar15u(I0,J0,K) - P2 = iar15u(I0,J1,K) - P3 = iar15u(I1,J0,K) - P4 = iar15u(I1,J1,K) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - UI(K) = P5+(P6-P5)*WX - - P1 = iar16v(I0,J0,K) - P2 = iar16v(I0,J1,K) - P3 = iar16v(I1,J0,K) - P4 = iar16v(I1,J1,K) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - VI(K) = P5+(P6-P5)*WX - - P1 = iar17q(I0,J0,K) - P2 = iar17q(I0,J1,K) - P3 = iar17q(I1,J0,K) - P4 = iar17q(I1,J1,K) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - QI(K) = P5+(P6-P5)*WX - -C Layer Pressure -C -------------- - - P1 = iarpsl(I0,J0,K) - P2 = iarpsl(I0,J1,K) - P3 = iarpsl(I1,J0,K) - P4 = iarpsl(I1,J1,K) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - PLEV(K) = P5+(P6-P5)*WX - -C Interface Pressure -C ------------------ - - P1 = iarpsi(I0,J0,K+1) - P2 = iarpsi(I0,J1,K+1) - P3 = iarpsi(I1,J0,K+1) - P4 = iarpsi(I1,J1,K+1) - P5 = P1+(P2-P1)*WY - P6 = P3+(P4-P3)*WY - PINT(K+1) = P5+(P6-P5)*WX - - ENDDO - -C Compute interface heights -C ------------------------- - - zint(1) = zsi - pint(1) = psi - pintlog(1) = log(pint(1)) - do k=2,kmax - k0 = k-1 - zint(k) = zint(k0) - rog*ti(k0)*log(pint(k)/pint(k0)) - pintlog(k) = log(pint(k)) - enddo - pint(kmax+1) = 0.0 - -C Compute log(pressure) at layer midpoints -C ---------------------------------------- - - do k=1,kmax - plevlog(k) = log(plev(k)) - enddo - -ccccc print *,' pint=',pint(1:kmax) -ccccc print *,' zint=',zint(1:kmax) - - RETURN - END - -C$$$ SUBPROGRAM DOCUMENTATION BLOCK -C -C SUBPROGRAM: OEFG01 -C PRGMMR: D. A. KEYSER ORG: NP22 DATE: 2007-09-14 -C -C ABSTRACT: FUNCTION WHICH RETURNS THE OBSERVATION ERROR FOR A -C REQUESTED VARIABLE INTERPOLATED TO A DEFINED PRESSURE LEVEL FOR A -C DEFINED REPORT TYPE. IT IS OBTAINED FROM AN INPUT ARRAY CONTAINING -C OBSERVATION ERRORS ON FIXED PRESSURE LEVELS BY VARIABLE AND REPORT -C TYPE (READ EARLIER FROM THE EXTERNAL OBSERVATION ERROR TABLE) -C -C PROGRAM HISTORY LOG: -C 1995-05-17 J. WOOLLEN (NP20) - ORIGINAL AUTHOR (FUNCTION OEF) -C 2007-09-14 D. A. KEYSER -- MODIFIED TO USE EXACT LOGIC AS IN GSI -C (MINIMUM LIMITING VALUE FOR OBS ERROR BASED ON VARIABLE TYPE, -C LEVEL PRESSURE LIMITED TO MAX OF 2000 MB AND MIN OF ZERO MB, A -C FEW OTHER MINOR CHANGES) -C -C USAGE: XX = OEFG01(P,TYP,IE,OEMIN) -C INPUT ARGUMENT LIST: -C P - REAL PRESSURE LEVEL (MB) TO INTERPOLATE OBS ERROR TO -C TYP - REAL PREPBUFR REPORT TYPE -C IE - VARIABLE TYPE BEING INTERPOLATED (=2 - TEMPERATURE, -C - =3 - MOISTURE, =4 - WIND, =5 - SURFACE PRESSURE, =6 - -C - PRECIPITABLE WATER) -C - (USED ONLY IN PREVENTS MODE) -C OEMIN - REAL MINIMUM VALUE FOR OBS ERROR (FOR VARIABLE BEING -C - INTERPOLATED) -C -C REMARKS: 'OEFG01' RETURNED IS OBSERVATION ERROR FOR VARIABLE "IE" IN -C REPORT TYPE "TYP", INTERPOLATED TO PRESSURE LEVEL "P". -CC -C -C ATTRIBUTES: -C LANGUAGE: FORTRAN 90 -C MACHINE: NCEP WCOSS -C -C$$$ - FUNCTION OEFG01(P,TYP,IE,OEMIN) - - REAL(8) BMISS - - COMMON /GBEVDD/ERRS(300,33,6) - COMMON /GBEVFF/ BMISS - - OEFG01 = BMISS - -C LOOK UP ERRORS FOR PARTICULAR OB TYPES -C -------------------------------------- - - KX = TYP - - P = MAX(0.,MIN(P,2000.)) - - IF(P.GE.ERRS(KX,1,1)) K1 = 1 - - DO KL = 1,32 - IF(P.GE.ERRS(KX,KL+1,1).AND.P.LE.ERRS(KX,KL,1)) K1 = KL - ENDDO - - IF(P.LE.ERRS(KX,33,1)) K1 = 5 - - K2 = K1 + 1 - - EDIFF = ERRS(KX,K2,1) - ERRS(KX,K1,1) - IF(ABS(EDIFF).GT.0.0) THEN - DEL = (P - ERRS(KX,K1,1))/EDIFF - ELSE - DEL = BMISS - ENDIF - - DEL = MAX(0.,MIN(DEL,1.0)) - - OEFG01 = ((1.0 - DEL) * ERRS(KX,K1,IE)) + (DEL * ERRS(KX,K2,IE)) - OEFG01 = MAX(OEFG01,OEMIN) - -C SET MISSING ERROR VALUE TO "BMISS" -C ---------------------------------- - - IF(OEFG01.GE.5E5) OEFG01 = BMISS - - RETURN - - END - -C$$$ SUBPROGRAM DOCUMENTATION BLOCK -C -C SUBPROGRAM: GBLEVN08 -C PRGMMR: S. MELCHIOR ORG: NP22 DATE: 2014-03-25 -C -C ABSTRACT: CREATE VIRTUAL TEMPERATURE EVENTS WITHIN GBLEVENTS -C SUBROUTINE. FOR ALL TYPES EXCEPT RASS, THIS CONSISTS OF FIRST RE- -C CALCULATING THE SPECIFIC HUMIDITY FROM THE REPORTED DEWPOINT -C TEMPERATURE AND PRESSURE, FOLLOWED BY THE CALCULATION OF VIRTUAL -C TEMPERATURE FROM THE JUST-CALCULATED SPECIFIC HUMIDITY AND THE -C REPORTED (SENSIBLE) TEMPERATURE. THE RE-CALCULATED SPECIFIC -C HUMIDITY IS THEN ENCODED AS A STACKED EVENT TO BE LATER WRITTEN -C INTO THE PREPBUFR FILE (UNDER PROGRAM "VIRTMP", REASON CODE 0). -C IF THE NAMELIST SWITCH DOVTMP IS TRUE, THEN THE JUST-CALCULATED -C VIRTUAL TEMPERATURE IS THEN ALSO ENCODED AS A STACKED EVENT TO BE -C LATER WRITTEN INTO THE PREPBUFR FILE (UNDER PROGRAM "VIRTMP", -C REASON CODE 0, 2 OR 6). FOR RASS DATA, SPECIFIC HUMIDITY IS -C MISSING HOWEVER IF THE NAMELIST SWITCH DOVTMP IS TRUE, A SIMPLE -C COPY OF THE REPORTED (VIRTUAL) TEMPERATURE IS ENCODED AS A STACKED -C EVENT TO BE LATER WRITTEN INTO THE PREPBUFR FILE (UNDER PROGRAM -C "VIRTMP", REASON CODE 3). FOR SURFACE DATA WITH A MISSING PMSL, IF -C DOVTMP=T AND DOPMSL=T AND A VIRTUAL TEMPERATURE HAS BEEN COMPUTED, -C CALCULATE AN ESTIMATED PMSL AND ENCODE IT INTO PREPBUFR FILE ALONG -C WITH AN INDICATOR THAT IS WAS DERIVED HERE. THIS SUBROUTINE IS -C CURRENTLY ONLY CALLED FOR SURFACE LAND ("ADPSFC"), MARINE -C ("SFCSHP"), MESONET ("MSONET"), RASS ("RASSDA") OR SATELLITE -C TEMPERATURE RETRIEVAL ("SATEMP") DATA TYPES WHEN SWITCH -C "ADPUPA_VIRT" IS FALSE AND ONLY FOR SURFACE LAND ("ADPSFC"), MARINE -C ("SFCSHP"), MESONET ("MSONET"), RASS ("RASSDA"), SATELLITE -C TEMPERATURE RETRIEVAL ("SATEMP") OR RAOB/DROP/MULTI-LVL RECCO -C ("ADPUPA") DATA TYPES WHEN SWITCH "ADPUPA_VIRT" IS TRUE. IT IS -C ALSO ONLY CALLED IN THE PREVENTS MODE. THIS ROUTINE IS CALLED ONCE -C FOR EACH VALID REPORT IN THE PREPBUFR FILE. -C -C PROGRAM HISTORY LOG: -C 1995-05-17 J. WOOLLEN (NP20) - ORIGINAL AUTHOR -C 1997-06-01 D.A. KEYSER - STREAMLINED, ADDED SWITCH DOVTMP -C 1999-12-01 D. A. KEYSER -- SPEC. HUMIDITY AND VIRT. TEMPERATURE ARE -C NOW CALCULATED WHEN SPEC. HUMIDITY QUAL. MARKER IS BAD (SUBJECT -C TO A SANITY CHECK), HOWEVER THE VIRT. TEMPERATURE GETS A BAD -C QUAL. MARKER (8) -C 2004-08-30 D. A. KEYSER -- FOR "RASSDA" TYPES, ENCODES A SIMPLE COPY -C OF THE REPORTED (VIRTUAL) TEMPERATURE AS A "VIRTMP" EVENT IF -C DOVTMP IS TRUE, GETS NEW REASON CODE 3 -C 2006-07-14 D. A. KEYSER -- PROCESSES REPORTS IN MESSAGE TYPE ADPUPA -C (I.E., RAOBS, DROPS, MULTI-LEVEL RECCOS) WITH SAME LOGIC AS IN -C SUBROUTINE VTPEVN OF PROGRAM PREPOBS_CQCBUFR WHEN NEW NAMELIST -C SWITCH "ADPUPA_VIRT" IS TRUE {NORMALLY "ADPUPA_VIRT" IS FALSE -C (DEFAULT) BECAUSE SUBSEQUENT PROGRAM PREPOBS_CQCBUFR PERFORMS -C THIS FUNCTION} -C 2007-09-14 D. A. KEYSER -- FOR NON-"ADPUPA" TYPES, Q.M. 9 IS NOW -C ASSIGNED TO CALCULATED VIRT. TEMPS IF THE MOISTURE Q.M. IS 9 OR -C 15 AND ORIG. TEMP NOT "BAD", THESE "VIRTMP" EVENTS RECEIVE NEW -C REASON CODE 4, HAD RECEIVED Q.M. 8 WITH REASON CODE 2 LIKE VIRT. -C TEMPS CALCULATED FROM "BAD" MOISTURE - THIS MEANS ONLY TRULY -C "BAD" VIRT. TEMPS WILL NOW GET Q.M. 8 AND VIRT. TEMPS FLAGGED FOR -C NON-USE BY ASSIMILATION (BUT STILL "GOOD") WILL NOW GET Q.M. 9 -C (GSI MONITORS, BUT DOES NOT USE, OBS WITH Q.M. 9, BUT IT DOES NOT -C EVEN CONSIDER OBS WITH Q.M. 8); FOR "ADPUPA" TYPES, Q.M. 3 IS NOW -C ASSIGNED TO CALCULATED VIRT. TEMPS ONLY IF THE MOISTURE Q.M. IS -C TRULY BAD (I.E. > 3 BUT NOT 9 OR 15) (AND, AS BEFORE, ORIG. TQM -C IS 1 OR 2 AND POB IS BELOW 700 MB) - BEFORE, TQM SET TO 3 WHEN -C QQM WAS 9 OR 15 AND ALL OTHER CONDITIONS MET; FOR "SATEMP" TYPES, -C ENCODES A SIMPLE COPY OF THE REPORTED (VIRTUAL) TEMPERATURE AS A -C "VIRTMP" EVENT IF DOVTMP IS TRUE, GETS REASON CODE 3 (SIMILAR TO -C WHAT IS ALREADY DONE FOR "RASSDA" TYPES) -C 2013-04-12 D. A. KEYSER -- DON'T ALLOW CALCULATED Q TO BE < 0 WHICH -C CAN OCCUR ON WCOSS FOR CASES OF HORRIBLY BAD MESONET DATA -C 2014-03-25 S. MELCHIOR -- ADDED NEW NAMELIST SWITCH "DOPMSL" WHICH, -C WHEN TRUE, DERIVES PMSL (MNEMONIC "PMO") FROM REPORTED STATION -C PRESSURE ("POB"), STATION HEIGHT/ELEVATION ("ZOB") AND THE JUST- -C COMPUTED VIRTUAL TEMPERATURE FOR SURFACE REPORTS IN CASES WHEN -C REPORTED PMSL IS MISSING. DOVTMP MUST BE TRUE AND DOANLS MUST BE -C FALSE ("PREVENTS" MODE). THE DERIVED PMSL EITHER GETS A QUALITY -C MARK ("PMQ") OF 3 OR INHERITS THE STATION PRESSURE QUALITY MARK -C ("PQM"), WHICHEVER IS GREATER. THE VALUE OF THE PMSL INDICATOR -C (NEW MNEMONIC "PMIN") IS SET TO 1 TO DENOTE PMSL WAS DERIVED -C RATHER THAN OBSERVED. THE DEFAULT FOR "DOPMSL" IS FALSE -C (NORMALLY ONLY TRUE IN RTMA AND URMA RUNS). IT IS FORCED TO BE -C FALSE IN "POSTEVENTS" MODE (DOANLS=TRUE). -C -C USAGE: CALL GBLEVN08(IUNITP) -C INPUT ARGUMENT LIST: -C IUNITP - BUFR OUTPUT FILE UNIT -C SUBSET - THE BUFR MESSAGE TABLE A ENTRY FOR THE PARTICULAR -C - REPORT BEING PROCESSED -C -C REMARKS: WILL IMMEDIATELY RETURN TO CALLING PROGRAM IF ANY OF THE -C FOLLOWING CONDITIONS EXIST: THERE ARE NO LEVELS OF VALID DEWPOINT, -C OBS, TEMPERATURE Q.M. OR SPEC. HUMIDITY Q.M. IN THE INPUT PREPBUFR -C FILE FOR THE REPORT. WILL NOT ATTEMPT EITHER SPEC. HUMIDITY NOR -C VIRT. TEMP CALC. ON A GIVEN LEVEL IF ANY OF THE FOLLOWING -C CONDITIONS EXIST: REPORTED PRESSURE OBS IS MISSING, REPORTED -C (SENSIBLE) TEMPERATURE OBS IS MISSING, OR REPORTED DEWPOINT OBS IS -C MISSING. -C -C ATTRIBUTES: -C LANGUAGE: FORTRAN 90 -C MACHINE: NCEP WCOSS -C -C$$$ - SUBROUTINE GBLEVN08(IUNITP,SUBSET) ! FORMERLY SUBROUTINE VTPEVN - - parameter (Rd=287.04, g=9.81) - CHARACTER*80 EVNSTQ,EVNSTV,evnstp - CHARACTER*8 SUBSET,stnid - REAL(8) TDP_8(255),TQM_8(255),QQM_8(255),BAKQ_8(4,255), - $ BAKV_8(4,255),bakp_8(3),OBS_8,QMS_8,BAK_8, - $ SID_8,pqm_8 - real(8) pmo_8,zob_8,pmsl_8 - REAL(8) BMISS - - LOGICAL EVNQ,EVNV,DOVTMP,TROP,ADPUPA_VIRT,DOBERR,DOFCST, - $ SOME_FCST,FCST,VIRT,SATMQC,RECALC_Q,DOPREV, - $ evnp,dopmsl,surf - - COMMON /GBEVAA/ SID_8,OBS_8(13,255),QMS_8(12,255),BAK_8(12,255), - $ XOB,YOB,DHR,TYP,NLEV - COMMON /GBEVBB/ PVCD,VTCD - COMMON /GBEVCC/ DOVTMP,DOFCST,SOME_FCST,DOBERR,FCST,VIRT, - $ QTOP_REJ,SATMQC,ADPUPA_VIRT,RECALC_Q,DOPREV,dopmsl - COMMON /GBEVFF/ BMISS - - DATA EVNSTQ /'QOB QQM QPC QRC'/ - DATA EVNSTV /'TOB TQM TPC TRC'/ - data evnstp /'PMO PMQ PMIN'/ - - equivalence (sid_8,stnid) - -C----------------------------------------------------------------------- -C FCNS BELOW CONVERT TEMP/TD (K) & PRESS (MB) INTO SAT./ SPEC. HUM.(G/G) -C----------------------------------------------------------------------- - ES(T) = 6.1078*EXP((17.269*(T - 273.16))/((T - 273.16)+237.3)) - QS(T,P) = (0.622*ES(T))/(P-(0.378*ES(T))) -C----------------------------------------------------------------------- -C FCN BELOW derives mean sea-level pressure (mb)from P (station -C pressure, mb), Tv (virtual temperature, K), and Z (station height/ -C elevation, m) for surface reports -C----------------------------------------------------------------------- - PMSL_fcn(P,Tv,Z) = P*exp((g*Z)/(Rd*Tv)) -C----------------------------------------------------------------------- - -C CLEAR TEMPERATURE AND SPECIFIC HUMIDITY EVENTS -C ---------------------------------------------- - - EVNQ = .FALSE. - EVNV = .FALSE. - evnp = .false. - BAKQ_8 = BMISS - BAKV_8 = BMISS - bakp_8 = bmiss - TROP = .FALSE. - - surf = (subset.eq.'ADPSFC'.or.subset.eq.'SFCSHP'.or. - $ subset.eq.'MSONET') - -C GET DEWPOINT TEMPERATURE AND CURRENT T,Q QUALITY MARKERS -C Also get mean sea level pressure, station pressure quality mark, -C and station height (elevation) for surface reports if dopmsl=T -C ---------------------------------------------------------------- - - CALL UFBINT(-IUNITP,TDP_8,1,255,NLTD,'TDO') - CALL UFBINT(-IUNITP,QQM_8,1,255,NLQQ,'QQM') - CALL UFBINT(-IUNITP,TQM_8,1,255,NLTQ,'TQM') - if(surf.and.dopmsl) then - call ufbint(-iunitp,pmo_8,1,1,nlpm,'PMO') - call ufbint(-iunitp,zob_8,1,1,nlzo,'ZOB') - call ufbint(-iunitp,pqm_8,1,1,nlpq,'PQM') - endif - IF(SUBSET.NE.'RASSDA '.AND.SUBSET.NE.'SATEMP ') THEN - IF(NLTD.EQ.0) RETURN - IF(NLQQ.EQ.0) RETURN - ENDIF - IF(NLTQ.EQ.0) RETURN - IF(SUBSET.NE.'RASSDA '.AND.SUBSET.NE.'SATEMP ') THEN - IF(NLTD.NE.NLEV) THEN - PRINT'(" ##GBLEVENTS/GBLEVN08 - NLTD .NE. NLEV - STOP 61")' - CALL ERREXIT(61) - ENDIF - IF(NLQQ.NE.NLEV) THEN - PRINT'(" ##GBLEVENTS/GBLEVN08 - NLQQ .NE. NLEV - STOP 63")' - CALL ERREXIT(63) - ENDIF - ENDIF - IF(NLTQ.NE.NLEV) THEN - PRINT'(" ##GBLEVENTS/GBLEVN08 - NLTQ .NE. NLEV - STOP 62")' - CALL ERREXIT(62) - ENDIF - -C COMPUTE SPECIFIC HUMIDITY AND VIRTUAL TEMPERATURE USING REPORTED DEWP -C For surface reports, also calculate PMSL if it is missing, dopmsl=T, -C and a virtual temperature has been computed -C --------------------------------------------------------------------- - - IF(NLEV.GT.0) THEN - DO L=1,NLEV - POB = OBS_8(1,L) - TDO = TDP_8(L) - TOB = OBS_8(3,L) - CAT = OBS_8(12,L) - IF(DOVTMP) THEN - IF(SUBSET.EQ.'RASSDA '.OR.SUBSET.EQ.'SATEMP ') THEN - IF(TOB.LT.BMISS) THEN - BAKV_8(1,L) = TOB - BAKV_8(2,L) = TQM_8(L) - BAKV_8(3,L) = VTCD - BAKV_8(4,L) = 3 - EVNV = .TRUE. - CYCLE - ENDIF - ENDIF - ENDIF - IF(POB.LT.BMISS .AND. TOB.LT.BMISS - $ .AND. TDO.LT.BMISS) THEN - IF(QQM_8(L).GT.3) THEN -C Don't update q or calculate Tv if bad moisture obs fails sanity check -C (also don't calc. PMSL if it is missing and dopmsl=T for sfc rpts) -cdak IF(TDO.LT.-103.15 .OR. TDO.GT.46.83 .OR. POB.LT.0.1 .OR. -cdak $ POB.GT.1100.) -cdak $ print *, '&&& bad QM fails sanity check' - IF(TDO.LT.-103.15 .OR. TDO.GT.46.83 .OR. POB.LT.0.1 .OR. - $ POB.GT.1100.) CYCLE - ENDIF - QOB = QS(TDO+273.16,POB) -ccccc BAKQ_8(1,L) = QOB*1E6 ! dak fix 2/27/13: can't be > bmiss - ! else flting pt overflow in BUFRLIB -ccc IF(QOB*1E6.LT.BMISS) BAKQ_8(1,L) = QOB*1E6 - ! dak add'l fix 4/12/13: don't allow - ! calc. q to be < 0 which can occur - ! on WCOSS for cases of horribly bad - ! mesonet data - IF(QOB*1E6.LT.BMISS .AND. QOB.GT.0.) BAKQ_8(1,L) = QOB*1E6 - BAKQ_8(2,L) = QQM_8(L) ! Moist qm same as before for - ! re-calc. q - BAKQ_8(3,L) = VTCD - BAKQ_8(4,L) = 0 ! Re-calc. q gets unique reason code 0 - EVNQ = .TRUE. -C If message type ADPUPA, test this level to see if at or above trop -C (trop must be above 500 mb to pass test; if no trop level found -C assume it's at 80 mb) -C Don't calculate Tv on this level if at or above trop (doesn't affect -C q calculation) - TROP = (SUBSET.EQ.'ADPUPA ' .AND. - $ ((CAT.EQ.5 .AND. POB.LT.500.) .OR. POB.LT. 80. .OR. TROP)) - IF(DOVTMP .AND. .NOT.TROP) THEN - BAKV_8(1,L) = (TOB+273.16)*(1.+.61*QOB)-273.16 - BAKV_8(3,L) = VTCD - IF(SUBSET.EQ.'ADPUPA ') THEN -C Message type ADPUPA comes here - IF((QQM_8(L).LT.4.OR.QQM_8(L).EQ.9.OR.QQM_8(L).EQ.15) - $ .OR. TQM_8(L).EQ.0 .OR. TQM_8(L).GT.3 - $ .OR. POB.LE.700.) THEN - BAKV_8(2,L) = TQM_8(L) ! Tv qm same as for T when - ! q ok or q flagged by - ! PREPRO (but not bad) - BAKV_8(4,L) = 0 ! Tv gets unique reason code 0 - ELSE - BAKV_8(2,L) = 3 !Tv qm susp for bad moist below - ! 700 mb - BAKV_8(4,L) = 6 !Tv gets unique reason code 6 - ENDIF - ELSE -C All other message types come here - IF(QQM_8(L).LT.4) THEN - BAKV_8(2,L) = TQM_8(L) ! Tv qm same as for T when - ! q ok - BAKV_8(4,L) = 0 ! Tv gets unique reason code 0 - ELSE IF((QQM_8(L).EQ.9.OR.QQM_8(L).EQ.15).AND. - $ (TQM_8(L).LE.3.OR.TQM_8(L).GE.15.OR. - $ TQM_8(L).EQ.9)) THEN -cdak print'(" %%% process tvirt on lvl ",I0," for missing moist obs ", -cdak $ "error/high-up moist case when orig. T not ""bad"" (set TQM_8=", -cdak $ "9)")', l - BAKV_8(2,L) = 9 ! Tv qm 9 for moist w/ missing obs - ! error or moist flagged by - ! PREPRO (but not bad) and T qm - ! orig not "bad" - BAKV_8(4,L) = 4 ! Tv gets unique reason code 4 - ELSE -cdak print'(" %%% process tvirt on lvl ",I0," for ""bad"" QQM_8 case ", -cdak $ "or missing moist obs error/high-up moist w/ ""bad"" TQM_8 case", -cdak $ " (set TQM_8=8)")', l - BAKV_8(2,L) = 8 ! Tv qm 8 (bad) for "bad" moist or - ! moist w/ missing obs error or - ! moist flagged by PREPRO (but - ! not bad) and T qm orig "bad" - BAKV_8(4,L) = 2 ! Tv gets unique reason code 2 - ENDIF - if(surf) then -c For sfc rpts & dopmsl=T, derive Pmsl in cases when it is not reported -c DAK: Note right now this can only happen for sfc rpts where a Tv is -c set to be calculated and where it is successfully calculated. -c Eventually this may be able to be relaxed such that PMSL can be -c derived even if, e.g., no moisture were reported (in this case -c Ts would have to be used, still a decent estimate of PMSL could -c liekly be made). Might also be able to be derived if no T rpted. - if(dopmsl) then - if(ibfms(pmo_8).ne.0) then - tv = bakv_8(1,1) + 273.16 - zob=zob_8 - pmsl_8=PMSL_fcn(pob,tv,zob) - bakp_8(1) = pmsl_8 - bakp_8(2) = max(3,int(pqm_8))! qm>=3 b/c derived - bakp_8(3) = 1. ! pmin=1 for derived value -cccccccccc write(*,4000) stnid,bakp_8(3),bakp_8(2) - 4000 format('--> ID ',A8,' Pmsl missing - derived from Pstn; ', - $ 'PMIN = ',F4.1,' PMQ = ',F4.1,'') - evnp = .true. -c Diagnostics for Pmsl values that are suspiciously high - if(pmsl_8.gt.1060) then - write(*,4001) stnid,pob,bakp_8(1) - 4001 format('--> ID ',A8,' Derived PMSL unrealistic; FLAG; ', - $ 'POB = ',F7.2,' PMO = ',F7.2,'') - end if - end if - end if - end if - ENDIF - EVNV = .TRUE. - ENDIF - ENDIF - ENDDO - ENDIF - -C ENCODE EVENTS INTO REPORT -C ------------------------- - - IF(NLEV.GT.0) THEN - IF(EVNQ) CALL UFBINT(IUNITP,BAKQ_8,4,NLEV,IRET,EVNSTQ) - IF(EVNV) CALL UFBINT(IUNITP,BAKV_8,4,NLEV,IRET,EVNSTV) - if(nlev.eq.1.and.evnp) - $ call ufbint(iunitp,bakp_8,3,nlev,iret,evnstp) - ENDIF - - RETURN - END -C####################### -C####################### -C####################### -C####################### -C####################### -C####################### -C*********************************************************************** -C*********************************************************************** - SUBROUTINE GBLEVN10(IUNITF,IDATEP,IM,JM,IDRT) ! FORMERLY - ! SUBROUTINE GESRES - USE GBLEVN_MODULE - USE SIGIO_MODULE - USE SIGIO_R_MODULE - - IMPLICIT NONE - INTEGER IUNITF(2), IDATEP, IM, JM, IDRT - REAL, PARAMETER :: PI180=.0174532 - INTEGER*4, PARAMETER :: ONE=1, TEN=10 - - TYPE(SIGIO_HEAD) :: HEAD(2) - TYPE(SIGIO_DATS) :: DATS - TYPE(SIGIO_DATM) :: DATM - - INTEGER*4 IRET,IRET1,IRETS,IMJM4,KM4,IDVM,NTRAC,IUNIT4(2) - - INTEGER KFILES,IFILE,JFILE,IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2, - $ MDIMA,MDIMB,MDIMC,IROMB,MAXWV,IDIR,NS,I,J,K,L,II,JJ,IB,IE - - INTEGER IDATE(8,2),JDATE(8,2),KDATE(8,2),KINDX(2) - - CHARACTER*6 COORD(3) - CHARACTER*20 CFILE - - REAL FHOUR,RINC(5) - - DATA COORD /'SIGMA ','HYBRID','GENHYB'/ - - REAL, ALLOCATABLE :: cofs(:,:), cofv(:,:,:) - REAL, ALLOCATABLE :: cofs_f(:,:,:), cofv_f(:,:,:,:) - REAL (KIND=4),ALLOCATABLE :: grds(:,:,:), grdv(:,:,:,:), - $ wrk1(:,:), wrk2(:,:) - - IMAX = IM - JMAX = JM - IMJM4 = IM*JM - IUNIT4(:) = IUNITF(:) - - IF(MOD(MOD(IDATEP,100),3).EQ.0) THEN - KFILES = 1 - KINDX = 0 - PRINT 331, MOD(IDATEP,100) - 331 FORMAT(/' --> GBLEVENTS: THE PREPBUFR CENTER HOUR (',I2.2, - $ ') IS A MULTIPLE OF 3 - ONLY ONE SIGIO (SIGMA OR HYBRID) ', - $ 'INPUT GLOBAL FILE IS READ,'/16X,'NO TIME INTERPOLATION OF ', - $ 'SPECTRAL COEFFICIENTS IS PERFORMED'/) - ELSE - KFILES = 2 - KINDX(1) = MOD(MOD(IDATEP,100),3) - KINDX(2) = KINDX(1) - 3 - PRINT 332, MOD(IDATEP,100) - 332 FORMAT(/' --> GBLEVENTS: THE PREPBUFR CENTER HOUR (',I2.2, - $ ') IS NOT A MULTIPLE OF 3 - TWO SPANNING SIGIO (SIGMA OR ', - $ 'HYBRID) INPUT GLOBAL FILES'/16X,'ARE READ AND THE SPECTRAL ', - $ 'COEFFICIENTS ARE INTERPOLATED TO THE PREPBUFR CENTER TIME'/) - ENDIF - -C GET VALID-TIME DATE OF SIGMA OR HYBRID FILE(S), ALSO READ HEADERS -C ----------------------------------------------------------------- - - JFILE = 0 - RINC = 0 - DO IFILE=1,KFILES - JFILE = IFILE - - WRITE(CFILE,'("fort.",I2.2)') IUNITF(IFILE) - - CALL SIGIO_RROPEN(IUNIT4(IFILE),CFILE,IRET) - CALL SIGIO_RRHEAD(IUNIT4(IFILE),HEAD(IFILE),IRET1) - print *,' cfile_sig=',cfile,'sigio_rropen: iret=',iret, - $ 'sigio_rrhead: iret1=',iret1 - - IF(IRET.NE.0) GO TO 903 - IF(IRET1.NE.0) GO TO 904 - - IDATE(:,IFILE) = 0 - - IDATE(1,IFILE) = HEAD(IFILE)%IDATE(4) - IDATE(2:3,IFILE) = HEAD(IFILE)%IDATE(2:3) - IDATE(5,IFILE) = HEAD(IFILE)%IDATE(1) - - FHOUR = HEAD(IFILE)%FHOUR - print'(" idate=",I5,7I3.2," fhour=",F7.3)', idate(:,ifile), - $ head(ifile)%fhour - - IF(IDATE(1,IFILE).LT.100) THEN - -C IF 2-DIGIT YEAR FOUND IN GLOBAL SIMGA FILE INITIAL DATE -C (IDATE(1,IFILE)), MUST USE "WINDOWING" TECHNIQUE TO CREATE A 4-DIGIT -C YEAR (NOTE: THE T170 IMPLEMENTATION IN JUNE 1998 WAS TO INCLUDE THE -C WRITING OF A 4-DIGIT YEAR HERE. PRIOR TO THIS, THE YEAR HERE WAS -C 2-DIGIT.) - - PRINT'(" ##GBLEVENTS/GBLEVN10 - 2-DIGIT YEAR FOUND IN ", - $ "SIGIO (SIGMA OR HYBRID) INPUT GLOBAL FILE ",I0, - $ "; INITIAL DATE (YEAR IS: ",I0,")")', IFILE,idate(1,IFILE) - PRINT'(" - USE WINDOWING TECHNIQUE TO OBTAIN 4-DIGIT", - $ " YEAR")' - IF(IDATE(1,IFILE).GT.20) THEN - IDATE(1,IFILE) = 1900 + IDATE(1,IFILE) - ELSE - IDATE(1,IFILE) = 2000 + IDATE(1,IFILE) - ENDIF - PRINT'(" ##GBLEVENTS/GBLEVN10 - CORRECTED 4-DIGIT YEAR IS", - $ " NOW: ",I0)', IDATE(1,IFILE) - ENDIF - - RINC(2) = FHOUR - CALL W3MOVDAT(RINC,IDATE(:,IFILE),JDATE(:,IFILE)) - - PRINT 1, IFILE,HEAD(IFILE)%FHOUR, - $ (IDATE(II,IFILE),II=1,3),IDATE(5,IFILE),(JDATE(II,IFILE), - $ II=1,3),JDATE(5,IFILE) - 1 FORMAT(' --> GBLEVENTS: SIGIO (SIGMA OR HYBRID) INPUT GLOBAL ', - $ 'FILE',I2,' HERE IS A ',F5.1,' HOUR FORECAST FROM ',I5.4, - $ 3I3.2,' VALID AT ',I5.4,3I3.2) - - KDATE(:,IFILE) = JDATE(:,IFILE) - - IF(KFILES.EQ.2) THEN - RINC(2) = REAL(KINDX(IFILE)) - CALL W3MOVDAT(RINC,JDATE(:,IFILE),KDATE(:,IFILE)) - ENDIF - - IDATGS_COR = (KDATE(1,IFILE) * 1000000) + (KDATE(2,IFILE) * - $ 10000) + (KDATE(3,IFILE) * 100) + KDATE(5,IFILE) - -C VALID DATES MUST MATCH -C ---------------------- - - IF(IDATEP.NE.IDATGS_COR) GO TO 901 - - ENDDO - - -C EXTRACT HEADER INFO -C ------------------- - - JCAP = HEAD(1)%JCAP - KMAX = HEAD(1)%LEVS - KM4 = KMAX - IDVC = HEAD(1)%IDVC - IDSL = head(1)%IDSL - IDVM = HEAD(1)%IDVM - NTRAC = HEAD(1)%NTRAC - NVCOORD = HEAD(1)%NVCOORD - ALLOCATE (VCOORD(KMAX+1,NVCOORD)) - VCOORD = 0.0 - VCOORD = HEAD(1)%VCOORD - - SFCPRESS_ID = MOD(HEAD(1)%IDVM,TEN) - THERMODYN_ID = MOD(HEAD(1)%IDVM/TEN,TEN) - IF(IDVC == 3 .AND. THERMODYN_ID == 3) THEN - KMAXS = (NTRAC+1)*KMAX + 2 - ELSE - KMAXS = 2*KMAX + 2 - NTRAC = 1 - ENDIF - - ALLOCATE (iar12z(im,jm), iar13p(im,jm)) - ALLOCATE (iar14t(im,jm,kmax), iar15u(im,jm,kmax), - $ iar16v(im,jm,kmax), iar17q(im,jm,kmax), - $ iarpsl(im,jm,kmax), iarpsi(im,jm,kmax+1)) - - - if(idvc.eq.0) idvc = 1 ! Reset IDVC=0 to 1 (sigma coord.) - IF(IDVC < 0 .or. IDVC > 3) THEN - PRINT *, '##GBLEVENTS/GBLEVN10: INVALID VERT COORD ID (=',IDVC - ENDIF - - -C DEFINE THE OTHER RESOLUTION PARAMETERS -C -------------------------------------- - - JCAP1 = JCAP+1 - JCAP2 = JCAP+2 - JCAP1X2 = JCAP1*2 - MDIMA = JCAP1*JCAP2 - MDIMB = MDIMA/2+JCAP1 - MDIMC = MDIMB*2 - IMAX = 384 - JMAX = IMAX/2+1 - - DLAT = 180./(JMAX-1) - DLON = 360./IMAX - - PRINT 2, JCAP,IMAX,JMAX,KMAX,kmaxs,DLAT,DLON,COORD(IDVC) - 2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',I5,' ', - $ I5,' x ',I5,' GRID WITH ',I3,' LEVELS ',I3, - $' SCALARS -------> ',F5.2,' X ',F5.2,' VERT. ', - $ 'COORD: ',A) - - GO TO 902 - - 901 CONTINUE - PRINT 9901, JFILE,(JDATE(II,JFILE),II=1,3),JDATE(5,JFILE),IDATEP - 9901 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',I2,' DATE', - $ ' (',I4.4,3(I2.2),'), DOES NOT MATCH -OR SPAN- PREPBUFR FILE ', - $ 'CENTER DATE (',I10,') -STOP 68'/) - CALL ERREXIT(68) - 903 CONTINUE - PRINT 9903, JFILE,IRET - 9903 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',I2, - $ ' RETURNED FROM SIGIO_RROPEN WITH R.C.',I3,' -STOP 70'/) - CALL ERREXIT(70) - 904 CONTINUE - PRINT 9904, JFILE,IRET1 - 9904 FORMAT(/' ##GBLEVENTS/GBLEVN10 - SIGMA OR HYBRID FILE',I2, - $ ' RETURNED FROM SIGIO_RRHEAD WITH R.C.',I3,' -STOP 71'/) - CALL ERREXIT(71) - 902 CONTINUE - IF(KMAX.GT.500) then - PRINT'(" ##GBLEVENTS/GBLEVN10 - KMAX TOO BIG = ",I0, - $ " - UNABLE TO TRANSFORM GLOBAL SIGMA FILE(S) - STOP 69")', - $ KMAX - CALL ERREXIT(69) - ENDIF - -C*********************************************************************** -C*********************************************************************** - -C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -C USAGE: CALL SPTEZM(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,WAVE,GRID,IDIR) -C INPUT ARGUMENTS: -C IROMB - INTEGER SPECTRAL DOMAIN SHAPE -C (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C MAXWV - INTEGER SPECTRAL TRUNCATION -C IDRT - INTEGER GRID IDENTIFIER -C (IDRT=4 FOR GAUSSIAN GRID, -C IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES, -C IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES) -C IMAX - INTEGER EVEN NUMBER OF LONGITUDES -C JMAX - INTEGER NUMBER OF LATITUDES -C KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM -C WAVE - REAL (2*MX,KMAX) WAVE FIELD IF IDIR>0 -C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 -C GRID - REAL (IMAX,JMAX,KMAX) GRID FIELD (E->W,N->S) IF IDIR<0 -C IDIR - INTEGER TRANSFORM FLAG -C (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE) -C OUTPUT ARGUMENTS: -C WAVE - REAL (2*MX,KMAX) WAVE FIELD IF IDIR<0 -C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 -C GRID - REAL (IMAX,JMAX,KMAX) GRID FIELD (E->W,N->S) IF IDIR>0 - - -C USAGE: CALL SPTEZMV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, -C & WAVED,WAVEZ,GRIDU,GRIDV,IDIR) -C INPUT ARGUMENTS: -C WAVED - REAL (2*MX,KMAX) WAVE DIVERGENCE FIELD IF IDIR>0 -C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 -C WAVEZ - REAL (2*MX,KMAX) WAVE VORTICITY FIELD IF IDIR>0 -C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 -C GRIDU - REAL (IMAX,JMAX,KMAX) GRID U-WIND (E->W,N->S) IF IDIR<0 -C GRIDV - REAL (IMAX,JMAX,KMAX) GRID V-WIND (E->W,N->S) IF IDIR<0 -C OUTPUT ARGUMENTS: -C WAVED - REAL (2*MX,KMAX) WAVE DIVERGENCE FIELD IF IDIR<0 -C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 -C WAVEZ - REAL (2*MX,KMAX) WAVE VORTICITY FIELD IF IDIR>0 -C WHERE MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 -C GRIDU - REAL (IMAX,JMAX,KMAX) GRID U-WIND (E->W,N->S) IF IDIR>0 -C GRIDV - REAL (IMAX,JMAX,KMAX) GRID V-WIND (E->W,N->S) IF IDIR>0 - -C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - allocate (cofs_f(mdima,kmaxs,2), cofv_f(mdima,kmax,2,2)) - - IROMB = 0 - MAXWV = JCAP - if (idrt < 0 .or. idrt > 256) IDRT = 0 - IDIR = 1 - - IF(KINDX(1).EQ.0) THEN - KFILES = 1 - ELSE - KFILES = 2 - ENDIF - -C Allocate for sigio read -C ----------------------- - - SFCPRESS_ID = MOD(HEAD(1)%IDVM,TEN) - THERMODYN_ID = MOD(HEAD(1)%IDVM/TEN,TEN) - - print *,'in sig sfcpress_id=',sfcpress_id,' thermodyn_id=', - $ thermodyn_id,' ntrac=',ntrac - print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm, - $ ' nvcoord=',nvcoord - - DO IFILE=1,KFILES - CALL SIGIO_ALDATS(HEAD(IFILE),DATS,IRETS) - CALL SIGIO_ALDATM(HEAD(IFILE),ONE,KM4,DATM,IRETS) - ! Read surface fields - CALL SIGIO_RRDATS(IUNIT4(IFILE),HEAD(IFILE),DATS,IRETS) - - IF(IRETS.NE.0) THEN - print *,' irets from sigio_rrdats = ', irets - RETURN - ENDIF - - DO I=1,MDIMA - COFS_F(I,1,IFILE) = DATS%HS(I) - COFS_F(I,2,IFILE) = DATS%PS(I) - ENDDO - - ! Read fields on levels 1 through kmax - CALL SIGIO_RRDATM(IUNIT4(IFILE),HEAD(IFILE),DATM,IRETS) - IF(IRETS.NE.0) THEN - print *,' irets from sigio_rrdatm = ', irets - RETURN - ENDIF -ccccc print *,' aft sigio_rrdatm irets=',irets - - IE = KMAX + 2 - COFS_F(:,3:IE,IFILE) = DATM%T - DO I=1,NTRAC - IB = IE + 1 - IE = IB + KMAX - 1 - COFS_F(:,IB:IE,IFILE) = DATM%Q(:,1:KMAX,I) - ENDDO - COFV_F(:,:,1,IFILE) = DATM%D - COFV_F(:,:,2,IFILE) = DATM%Z - - CALL SIGIO_AXDATS(DATS,IRETS) - CALL SIGIO_AXDATM(DATM,IRETS) - ENDDO - -ccccc print *,' after sigio_axdatm' - - ALLOCATE (COFS(MDIMA,KMAXS), COFV(MDIMA,KMAX,2)) - ALLOCATE (GRDS(imax,jmax,KMAXS), GRDV(imax,jmax,KMAX,2)) - ALLOCATE (WRK1(imax*jmax,KMAX), WRK2(imax*jmax,KMAX+1)) - - IF(KFILES.EQ.1) THEN - DO I = 1,MDIMA - COFS(I,1:KMAXS) = COFS_f(I,1:KMAXS,1) - COFV(I,1:KMAX,1:2) = COFV_f(I,1:KMAX,1:2,1) - ENDDO - - ELSE - COFS= - $ ((ABS(KINDX(2))*COFS_f(:,:,1)) +(KINDX(1)*COFS_f(:,:,2)))/3. - COFV= - $ ((ABS(KINDX(2))*COFV_f(:,:,:,1))+(KINDX(1)*COFV_f(:,:,:,2)))/3. - ENDIF - - DEALLOCATE (COFS_F, COFV_F) - - CALL SPTEZM(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAXS,COFS,GRDS,IDIR) - CALL SPTEZMV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, - & COFV(1,1,1),COFV(1,1,2),GRDV(1,1,1,1),GRDV(1,1,1,2),IDIR) - - IF( SFCPRESS_ID == 2 ) THEN ! for enthalpy version - GRDS(:,:,2) = 1000.0*GRDS(:,:,2) ! Now in Pa - ELSE - GRDS(:,:,2) = 1000.0*EXP(GRDS(:,:,2)) ! Now in Pa - ENDIF - - DO NS=1, KMAXS - CALL GBLEVN11(IMAX,JMAX,GRDS(1,1,NS)) - ENDDO - DO J=1,JMAX - DO I=1,IMAX - IAR12Z(I,J) = GRDS(I,J,1) ! Orography - IAR13P(I,J) = GRDS(I,J,2) * 0.01 ! Surface Pressure in hPa - ENDDO - ENDDO - - IF(THERMODYN_ID == 3 .AND. IDVC == 3) THEN - GRDS(:,:,3:2+KMAX) = GRDS(:,:,3:2+KMAX) / HEAD(1)%CPI(1) - print *,' cpi(0)=',head(1)%cpi(1) - ENDIF - CALL SIGIO_MODPR(IMJM4,IMJM4,KM4,NVCOORD,IDVC,IDSL,VCOORD,IRET, - $ GRDS(1,1,2),GRDS(1,1,3),PM=WRK1,PD=WRK2(1,2)) - - DO J=1,JMAX - JJ = (J-1)*IMAX - DO I=1,IMAX - WRK2(I+JJ,1) = GRDS(I,J,2) ! in Pa - ENDDO - ENDDO - DO L=1,KMAX - WRK2(:,L+1) = WRK2(:,L) - WRK2(:,L+1) ! in Pa - ENDDO - -ccccc print *,' wrk1=',wrk1(1001,:) -ccccc print *,' wrk2=',wrk2(1001,:) - -ccccc CALL GBLEVN11(IMAX,JMAX,WRK2(1,KMAX+1)) -ccccc DO L=1,KMAX -ccccc CALL GBLEVN11(IMAX,JMAX,WRK1(1,L)) -ccccc CALL GBLEVN11(IMAX,JMAX,WRK2(1,L)) -ccccc ENDDO - - IF(THERMODYN_ID == 3 .AND. IDVC == 3) THEN - ! Convert from enthalpy to temperature - GRDS(:,:,3:2+KMAX) = GRDS(:,:,3:2+KMAX) * HEAD(1)%CPI(1) - CALL SIGIO_CNVTDV(IMJM4,IMJM4,KM4,IDVC,IDVM,NTRAC,IRET, - $ GRDS(1,1,3),GRDS(1,1,3+KMAX),HEAD(1)%CPI,1_4) - ! Convert back to virtual temperature - GRDS(:,:,3:KMAX+2) = GRDS(:,:,3:KMAX+2) * - $ (1.+(461.50/287.05-1)*GRDS(:,:,3+KMAX:2+KMAX*2)) - ENDIF - - DO L=1,KMAX - DO K=1,2 - CALL GBLEVN11(IMAX,JMAX,GRDV(1,1,L,K)) - ENDDO - DO J=1,JMAX - JJ = (J-1)*IMAX - DO I=1,IMAX - IAR14T(I,J,L) = GRDS(I,J,2+L) ! Temp (virtual) - IAR15U(I,J,L) = GRDV(I,J,L,1) ! U component - IAR16V(I,J,L) = GRDV(I,J,L,2) ! V component - ! specific humidity - IAR17Q(I,J,L) = MAX(0.0,GRDS(I,J,2+KMAX+L)*1.0E6) - IARPSL(I,J,L) = WRK1(I+JJ,L)*0.01 ! 3D layer pres(hPa) - ENDDO - ENDDO - ENDDO - DO L=1,KMAX+1 - DO J=1,JMAX - JJ = (J-1)*IMAX - DO I=1,IMAX - IARPSI(I,J,L) = WRK2(I+JJ,L)*0.01 ! 3D interface pressure - ! (hPa) - ENDDO - ENDDO - ENDDO - - -ccccc print *,' iar14t=',iar14t(1,80,:) -ccccc print *,' iar15u=',iar15u(1,80,:) -ccccc print *,' iar16v=',iar16v(1,80,:) -ccccc print *,' iarpsi=',iarpsi(1,80,:) -ccccc print *,' iarpsl=',iarpsl(1,80,:) - - DEALLOCATE (COFS, COFV) - DEALLOCATE (GRDS, GRDV, WRK1, WRK2) - - print *,' RETURNING from GBLENV10' - - RETURN - - END -C*********************************************************************** -C*********************************************************************** - subroutine gblevn11(imax,jmax,grid) ! formerly subroutine n_s_swap - implicit none - integer imax, jmax - real grid(imax,jmax) - real temp (imax) - integer i, j, jj - - do j=1,jmax/2 - jj = jmax-j+1 - do i=1,imax - temp(i) = grid(i,j) - grid(i,j) = grid(i,jj) - grid(i,jj) = temp(i) - enddo - enddo - return - end -C*********************************************************************** -C*********************************************************************** - subroutine gblevn11d(imax,jmax,grid) - implicit none - integer imax, jmax - real(kind=8) grid(imax,jmax) - real(kind=8) temp (imax) - integer i, j, jj - - do j=1,jmax/2 - jj = jmax-j+1 - do i=1,imax - temp(i) = grid(i,j) - grid(i,j) = grid(i,jj) - grid(i,jj) = temp(i) - enddo - enddo - return - end -!--------------------------------------------------------- - SUBROUTINE GBLEVN12(IUNITF,IDATEP,IM,JM,IDRT) -!--------------------------------------------------------- -!--INITIALLY DEVELOPED BY HANG LEI BASED ON GBLEVN10 FOR NEMSIO -!--AND SUBSEQUENTLY MODIFIED BY FANGLIN YANG AND RUSS TREADON - - USE GBLEVN_MODULE - USE NEMSIO_MODULE - USE nemsio_openclose - USE nemsio_read - USE nemsio_write - use sigio_module - - IMPLICIT NONE - INTEGER IUNITF(2), IDATEP, IM, JM, IDRT - REAL(KIND=8),PARAMETER:: CON_RV=4.6150E+2,CON_RD=2.8705E+2, - $ FV=CON_RV/CON_RD-1.,ONER=1.0,QMIN=1.0E-10 - INTEGER*4, PARAMETER :: TEN=10 - - INTEGER*4 IDRTNEMS - INTEGER*4 IRET,IMJM4,KM4,IDVM,NTRAC - - INTEGER IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,MDIMA,MDIMB,MDIMC, - $ IROMB,MAXWV,IDIR,I,J,K,L,II,JJ - - INTEGER(4) IDATE7(7),JCAP4,IDVC4,DIMZ4,K4,IM4,JM4,KINDREAL,KINDINT - INTEGER(4) NFHOUR,NFMINUTE,NFSECONDN,NFSECONDD,idsl4,idvm4 - REAL(8) tfac - REAL(4),ALLOCATABLE :: VCOORD4(:,:,:),CPI(:) - REAL,ALLOCATABLE :: tmp(:) - TYPE(NEMSIO_GFILE) :: GFILE - - INTEGER IDATE(8),JDATE(8) - CHARACTER*20 CFILE2 - REAL FHOUR,RINC(5) - REAL (KIND=4),ALLOCATABLE :: psfc(:,:), tv(:,:,:), - $ wrk1(:,:), wrk2(:,:) - - IMAX = IM - JMAX = JM - IMJM4 = IM*JM - -! no time interpolation needed for nemsio input since files are -! available every hour - original logic to perform time interpolation -! is flawed and has been removed in this revision -! -------------------------------------------------------------------- - - print 331, mod(idatep,100) - 331 format(/' --> GBLEVENTS: ONLY ONE NEMSIO INPUT GLOBAL FILE IS ', - $ 'READ SINCE FILES ARE AVAILABLE EACH HOUR (NO NEED TO ', - $ 'INTERPOLATE'/16X,'SPANNING FILES WHEN THE PREPBUFR CENTER HOUR', - $ ' (',I2.2,') IS NOT A MULTIPLE OF 3)'/) - -C GET VALID-TIME DATE OF NEMSIO INPUT FILE, ALSO READ HEADERS -C ----------------------------------------------------------- - CALL NEMSIO_INIT() - - RINC = 0 - IDATE = 0 - WRITE(CFILE2,'("fort.",I2.2)') IUNITF(1) - - CALL NEMSIO_OPEN(GFILE,trim(CFILE2),'read',iret=iret) - print *,' cfile_nems2=',cfile2,'nemsio open,iret=',iret - - IDRTNEMS=IDRT - CALL NEMSIO_GETFILEHEAD(GFILE,IDATE=IDATE7, - & JCAP=JCAP4,DIMX=IM4,DIMY=JM4, - & DIMZ=DIMZ4,IDVC=IDVC4,NTRAC=NTRAC,IDRT=IDRTNEMS, - & NFHOUR=NFHOUR,NFMINUTE=NFMINUTE,NFSECONDN=NFSECONDN, - & NFSECONDD=NFSECONDD,idsl=idsl4,idvm=idvm4,IRET=IRET) - JCAP=JCAP4 - IDVC=IDVC4 - idsl=idsl4 - idvm=idvm4 - IMAX=IM4 - JMAX=JM4 - KMAX=DIMZ4 - if(IDRT==-9999) IDRT=4 ! default for gaussian grid - IDATE(1:3)=IDATE7(1:3) - IDATE(5)=IDATE7(4) - - ALLOCATE(VCOORD(KMAX+1,3)) - VCOORD=0.0 - - IF(NFSECONDD/=0.AND. NFSECONDD/=-9999) THEN - FHOUR=REAL(NFHOUR,8)+REAL(NFMINUTE/60.,8)+ - $ REAL(NFSECONDN*1./(NFSECONDD*360.),8) - ELSE - FHOUR=REAL(NFHOUR,8)+REAL(NFMINUTE/60.,8) - ENDIF - print'(" idate=",I5,7I3.2," fhour=",F7.3)', idate,fhour - - ALLOCATE(VCOORD4(KMAX+1,3,2)) - CALL NEMSIO_GETFILEHEAD(GFILE,VCOORD=VCOORD4,IRET=IRET) - NVCOORD=3 - IF(MAXVAL(VCOORD4(:,3,1))==0..AND.MINVAL(VCOORD4(:,3,1))==0.) THEN - NVCOORD=2 - IF(MAXVAL(VCOORD4(:,2,1))==0..AND.MINVAL(VCOORD4(:,2,1))==0.) - $ NVCOORD=1 - ENDIF - - VCOORD(1:KMAX+1,1:NVCOORD)=VCOORD4(1:KMAX+1,1:NVCOORD,1) - DEALLOCATE(VCOORD4) - - ALLOCATE(CPI(NTRAC+1)) - CALL NEMSIO_GETHEADVAR(GFILE,'CPI',CPI,IRET=IRET) - - RINC(2) = FHOUR - CALL W3MOVDAT(RINC,IDATE,JDATE) - - PRINT 1, FHOUR,(IDATE(II),II=1,3),IDATE(5),(JDATE(II),II=1,3), - $ JDATE(5) - 1 FORMAT(' --> GBLEVENTS: GLOBAL NEMSIO FILE: HERE IS A ',F5.1, - $ ' HOUR FORECAST FROM ',I5.4,3I3.2,' VALID AT ',I5.4,3I3.2) - - IDATGS_COR = (JDATE(1) * 1000000) + (JDATE(2) * 10000) + - $ (JDATE(3) * 100) + JDATE(5) - -C VALID DATES MUST MATCH -C ---------------------- - - IF(IDATEP.NE.IDATGS_COR) GO TO 901 - -C------------------------------------------ - - SFCPRESS_ID = MOD(IDVM,TEN) - THERMODYN_ID = MOD(IDVM/TEN,TEN) - - IF(IDVC == 3 .AND. THERMODYN_ID == 3) THEN - KMAXS = (NTRAC+1)*KMAX + 2 - ELSE - KMAXS = 2*KMAX + 2 - NTRAC = 1 - ENDIF - - ALLOCATE (iar12z(imax,jmax), iar13p(imax,jmax)) - ALLOCATE (iar14t(imax,jmax,kmax), iar15u(imax,jmax,kmax), - $ iar16v(imax,jmax,kmax), iar17q(imax,jmax,kmax), - $ iarpsl(imax,jmax,kmax), iarpsi(imax,jmax,kmax+1), - $ iarpsd(imax,jmax,kmax)) - -!!! DAK: is this right??????? - if(idvc.eq.0) idvc = 1 ! Reset IDVC=0 to 1 (sigma coord.) - IF(IDVC < 0 .or. IDVC > 3) THEN - PRINT *, '##GBLEVENTS/GBLEVN12: INVALID VERT COORD ID (=',IDVC - ENDIF - - -C DEFINE THE OTHER RESOLUTION PARAMETERS -C -------------------------------------- - - JCAP1 = JCAP+1 - JCAP2 = JCAP+2 - JCAP1X2 = JCAP1*2 - MDIMA = JCAP1*JCAP2 - MDIMB = MDIMA/2+JCAP1 - MDIMC = MDIMB*2 - - DLAT = 180./(JMAX-1) - DLON = 360./IMAX - - PRINT 2, JCAP,IMAX,JMAX,KMAX,kmaxs,DLAT,DLON,idvc - -! DAK: is below correct?????????? - 2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',I5,' ', - $ I5,' x ',I5,' GRID WITH ',I3,' LEVELS ',I3, - $' SCALARS -------> ',F5.2,' X ',F5.2,' VERT. ', - $ 'COORD ID IS: ',I0,' (not sure what this means!') - -C*********************************************************************** - - IROMB = 0 - MAXWV = JCAP - if (idrt < 0 .or. idrt > 256) IDRT = 0 - IDIR = 1 - - print *,'in nem sfcpress_id=',sfcpress_id,' thermodyn_id=', - $ thermodyn_id,' ntrac=',ntrac - print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,' nvcoord=', - $ nvcoord - - call w3kind(kindreal,kindint) - -! print *,'in gblevent, get fields,imax,jmax - allocate(tmp(imax*jmax)) -! -! hgt - K4=1 - if(kindreal==4) then - CALL NEMSIO_READRECVw34(GFILE,'hgt','sfc',K4,tmp,iret=iret) - else - CALL NEMSIO_READRECV(GFILE,'hgt','sfc',K4,tmp,iret=iret) - endif -! print *,'in getnemsio,hgtsfc,',maxval(tmp),minval(tmp),iret - if(iret==0) THEN - IAR12Z(:,:)=reshape(tmp,(/imax,jmax/)) - call gblevn11d(imax,jmax,IAR12Z) - ENDIF - -!sfc pres, input in Pa, out in hPa - if(kindreal==4) then - CALL NEMSIO_READRECVw34(GFILE,'pres','sfc',K4,tmp,iret=iret) - else - CALL NEMSIO_READRECV(GFILE,'pres','sfc',K4,tmp,iret=iret) - endif -! print *,'in getnemsio,pressfc,',maxval(tmp),minval(tmp),iret - if(iret==0) THEN - IAR13P(:,:)=reshape(tmp*0.01,(/imax,jmax/)) - call gblevn11d(imax,jmax,IAR13P) - ENDIF -!dpres -! DO K4=1,KMAX -! if(kindreal==4) then -! CALL NEMSIO_READRECVw34(GFILE,'dpres','mid layer',K4,tmp, -! $ iret=iret) -! else -! CALL NEMSIO_READRECV(GFILE,'dpres','mid layer',K4,tmp, -! $ iret=iret) -! endif -! print *,'in getnemsio,dpres,k4=',k4,maxval(tmp),minval(tmp), -! $ iret -! if(iret==0) THEN -! IARPSI(:,:,K4+1)=reshape(tmp*0.01,(/imax,jmax/)) -! ENDIF -! ENDDO - -!pres -! DO K4=1,KMAX -! if(kindreal==4) then -! CALL NEMSIO_READRECVw34(GFILE,'pres','mid layer',K4,tmp, -! $ iret=iret) -! else -! CALL NEMSIO_READRECV(GFILE,'pres','mid layer',K4,tmp, -! $ iret=iret) -! endif -! print *,'in getnemsio,pres,k4=',k4,maxval(tmp),minval(tmp),iret -! if(iret==0) THEN -! IARPSL(:,:,k4)=reshape(tmp*0.01,(/imax,jmax/)) -! ENDIF -! ENDDO -! tmp - DO K4=1,kmax - if(kindreal==4) then - CALL NEMSIO_READRECVw34(GFILE,'tmp','mid layer',K4,tmp, - $ iret=iret) - else - CALL NEMSIO_READRECV(GFILE,'tmp','mid layer',K4,tmp, - $ iret=iret) - endif -! print *,'in getnemsio,tmp,k4=',k4,maxval(tmp),minval(tmp),iret - if(iret==0) THEN - IAR14T(:,:,K4)=reshape(tmp,(/imax,jmax/)) - call gblevn11d(imax,jmax,IAR14T(1,1,K4)) - ENDIF - ENDDO -!u - DO K4=1,kmax - if(kindreal==4) then - CALL NEMSIO_READRECVw34(GFILE,'ugrd','mid layer',K4,tmp, - $ iret=iret) - else - CALL NEMSIO_READRECV(GFILE,'ugrd','mid layer',K4,tmp, - $ iret=iret) - endif -! print *,'in getnemsio,utmp,k4=',k4,maxval(tmp),minval(tmp),iret - if(iret==0) THEN - IAR15U(:,:,k4)=reshape(tmp,(/imax,jmax/)) - call gblevn11d(imax,jmax,IAR15U(1,1,K4)) - ENDIF - ENDDO -! do k=1,kmax -! print *,'in getnemsio,i15u,k=',k,maxval(iar15u(:,:,k)), -! $ minval(iar15u(:,:,k)),iret -! end do -! print *,'in getnemsio,IAR15U.1,',maxval(IAR15U),minval(IAR15U) - -!v - DO K4=1,kmax - if(kindreal==4) then - CALL NEMSIO_READRECVw34(GFILE,'vgrd','mid layer',K4,tmp, - $ iret=iret) - else - CALL NEMSIO_READRECV(GFILE,'vgrd','mid layer',K4,tmp, - $ iret=iret) - endif - if(iret==0) THEN - IAR16V(:,:,K4)=reshape(tmp,(/imax,jmax/)) - call gblevn11d(imax,jmax,IAR16V(1,1,K4)) - ENDIF - ENDDO -! print *,'in getnemsio,IAR15U.2,',maxval(IAR15U),minval(IAR15U) -!q - DO K4=1,kmax - if(kindreal==4) then - CALL NEMSIO_READRECVw34(GFILE,'spfh','mid layer',K4,tmp, - $ iret=iret) - else - CALL NEMSIO_READRECV(GFILE,'spfh','mid layer',K4,tmp, - $ iret=iret) - endif - if(iret==0) THEN - IAR17Q(:,:,K4)=max(0.0,reshape(tmp*1.0E6,(/imax,jmax/)) ) - call gblevn11d(imax,jmax,IAR17Q(1,1,K4)) - endif - ENDDO -! print *,'in getnemsio,IAR15U.3,',maxval(IAR15U),minval(IAR15U) - -! - deallocate(tmp) - CALL NEMSIO_CLOSE(GFILE,iret=iret) -! print *,'end of nemsio' -! - -! print *,'in getnemsio,IAR15U.4,',maxval(IAR15U),minval(IAR15U) - -!--compute Tv from temp - DO K=1,KMAX - DO J=1,JMAX - DO i=1,IMAX - TFAC=ONER+FV*MAX(IAR17Q(I,J,K)*1.0E-6,QMIN) - IAR14T(I,J,K)=IAR14T(I,J,K)*TFAC - ENDDO - ENDDO - ENDDO - -! print *,'in getnemsio,IAR15U.5,',maxval(IAR15U),minval(IAR15U) - -!--COMPUTE PM AND PD - ALLOCATE (psfc(imax,jmax), tv(imax,jmax,kmax)) - ALLOCATE (WRK1(imax*jmax,KMAX), WRK2(imax*jmax,KMAX+1)) - - imjm4=imax*jmax - km4=kmax - psfc(:,:) = iar13p(:,:)*100. - tv(:,:,:) = iar14t(:,:,:) - - IF(THERMODYN_ID == 3 .AND. IDVC == 3) THEN - tv(:,:,:) = tv(:,:,:) / CPI(1) - print *,' cpi(1)=',cpi(1) - ENDIF - - CALL SIGIO_MODPR(IMJM4,IMJM4,KM4,NVCOORD,IDVC,IDSL,VCOORD,IRET, - $ psfc,tv,PM=WRK1,PD=WRK2(1,2)) - DO J=1,JMAX - JJ = (J-1)*IMAX - DO I=1,IMAX - WRK2(I+JJ,1) = psfc(i,j) ! in Pa - ENDDO - ENDDO - DO L=1,KMAX - WRK2(:,L+1) = WRK2(:,L) - WRK2(:,L+1) ! in Pa - ENDDO - - DO L=1,KMAX - DO J=1,JMAX - JJ = (J-1)*IMAX - DO I=1,IMAX - IARPSL(I,J,L) = WRK1(I+JJ,L)*0.01 ! 3D layer pres(hPa) - ENDDO - ENDDO - ENDDO - DO L=1,KMAX+1 - DO J=1,JMAX - JJ = (J-1)*IMAX - DO I=1,IMAX - IARPSI(I,J,L) = WRK2(I+JJ,L)*0.01 ! 3D interface pressure - ! (hPa) - ENDDO - ENDDO - ENDDO - -! print*,'GBLEVN12 IARPSI,',maxval(IARPSI),minval(IARPSI) -! print*,'GBLEVN12 IARPSL,',maxval(IARPSL),minval(IARPSL) - - CALL NEMSIO_FINALIZE() -! - CALL GETLATS(IDRT) - - RETURN - - 901 CONTINUE - PRINT 9901, (JDATE(II),II=1,3),JDATE(5),IDATEP - 9901 FORMAT(/' ##GBLEVENTS/GBLEVN12 - NEMSIO INPUT GLOBAL FILE DATE', - $ ' (',I4.4,3(I2.2),'), DOES NOT MATCH PREPBUFR FILE CENTER ', - $ 'DATE (',I10,') -STOP 68'/) - CALL ERREXIT(68) - - END -! -------------------------------------------------------------------- - SUBROUTINE GBLEVN13(IUNITF,IDATEP,IMX,JMX,IDRT) -!--------------------------------------------------------- -!--INITIALLY DEVELOPED BY HANG LEI BASED ON GBLEVN10 FOR NETCDF INPUT - - USE GBLEVN_MODULE - use sigio_module - use netcdf - - IMPLICIT NONE - INTEGER IUNITF(2), IDATEP, IDRT,IMX,JMX - INTEGER yyyy,mm,dd,hh - INTEGER*4 ncid - integer*4 error, id_var, dimid, len - integer*4 im,jm,km, lm,n, k,nargs - REAL(KIND=8),PARAMETER:: CON_RV=4.6150E+2,CON_RD=2.8705E+2, - $ FV=CON_RV/CON_RD-1.,ONER=1.0,QMIN=1.0E-10 - INTEGER*4, PARAMETER :: TEN=10 - - INTEGER*4 IRET,IMJM4,KM4,IDVM,NTRAC - - INTEGER IDATGS_COR,JCAP,JCAP1,JCAP2,JCAP1X2,MDIMA,MDIMB,MDIMC, - $ IROMB,MAXWV,IDIR,I,J,L,II,JJ - - INTEGER(4) IDATE7(7),JCAP4,IDVC4,DIMZ4,K4,IM4,JM4,KINDREAL,KINDINT - INTEGER(4) NFHOUR,NFMINUTE,NFSECONDN,NFSECONDD,idsl4,idvm4,kr - REAL(8) tfac, time - REAL(4),ALLOCATABLE :: VCOORD4(:,:,:),CPI(:),ak(:),bk(:) - REAL,ALLOCATABLE :: temp(:,:), temp3d(:,:,:) -! CHARACTER,ALLOCATABLE :: attone(:) - character (len = 80) :: attone - INTEGER nt1, nt2 - character(len=10) :: dim_nam, grid - CHARACTER*20 gfname - - INTEGER(4) IDATE(8),JDATE(8) - REAL(4) FHOUR,RINC(5) - - REAL (KIND=4),ALLOCATABLE :: psfc(:,:), tv(:,:,:), - $ wrk1(:,:), wrk2(:,:) - - IMAX = IM - JMAX = JM - -! no time interpolation needed for netcdf input since files are -! available every hour - original logic to perform time interpolation -! is flawed and has been removed in this revision -! -------------------------------------------------------------------- - - print 331, mod(idatep,100) - 331 format(/' --> GBLEVENTS: ONLY ONE NETCDF INPUT GLOBAL FILE IS ', - $ 'READ SINCE FILES ARE AVAILABLE EACH HOUR (NO NEED TO ', - $ 'INTERPOLATE'/16X,'SPANNING FILES WHEN THE PREPBUFR CENTER HOUR', - $ ' (',I2.2,') IS NOT A MULTIPLE OF 3)'/) - -C GET VALID-TIME DATE OF NETCDF INPUT FILE, ALSO READ HEADERS -C ----------------------------------------------------------- - - WRITE(gfname,'("fort.",I2.2)') IUNITF(1) - - error=nf90_open(trim(gfname),nf90_nowrite,ncid) - error=nf90_inq_dimid(ncid,"grid_xt",dimid) - error=nf90_inquire_dimension(ncid,dimid,dim_nam,im) - error=nf90_inq_dimid(ncid,"grid_yt",dimid) - error=nf90_inquire_dimension(ncid,dimid,dim_nam,jm) - error=nf90_inq_dimid(ncid,"pfull",dimid) - error=nf90_inquire_dimension(ncid,dimid,dim_nam,km) - error=nf90_inq_dimid(ncid,"phalf",dimid) - error=nf90_inquire_dimension(ncid,dimid,dim_nam,lm) - print*, "im,jm,kmi,lm:",im,jm,km,lm - print*, "IMX, JMX:", IMX, JMX - IF (IMX .NE. IM .OR. JMX .NE. JM) print*, "Different Resolutions" - imax = im - jmax = jm - kmax = km - -! error=nf90_inq_varid(ncid, 'pfull', id_var) -! error=nf90_get_var(ncid, id_var, ) - - -! CALL W3MOVDAT(RINC,IDATE,JDATE) - - - -C VALID DATES MUST MATCH -C ---------------------- - error=nf90_inq_varid(ncid, "time", id_var) - error=nf90_inquire_attribute(ncid, id_var, "units", len=len) -! ALLOCATE (attone(len)) - error=nf90_get_att(ncid, id_var, "units", attone) - read(attone(len-18:len-15),'(i4)') yyyy - read(attone(len-13:len-12),'(i2)') mm - read(attone(len-10:len-9),'(i2)') dd - read(attone(len-7:len-6),'(i2)') hh - IF(attone(1:5) .NE. "hours") THEN - print*, "checkout the time unit, not hourly",attone - ELSE - print*, "base time", yyyy,mm,dd,hh - ENDIF -! DEALLOCATE (attone) - error=nf90_get_var(ncid, id_var, time) - fhour=time - - idate = 0 - idate(1)=yyyy - idate(2)=mm - idate(3)=dd - idate(5)=hh - - rinc=0.0 - rinc(2)=fhour - - CALL W3MOVDAT(RINC,IDATE,JDATE) - - PRINT 1, FHOUR,(IDATE(II),II=1,3),IDATE(5),(JDATE(II),II=1,3), - $ JDATE(5) - 1 FORMAT(' --> GBLEVENTS: GLOBAL NEMSIO FILE: HERE IS A ',F5.1, - $ ' HOUR FORECAST FROM ',I5.4,3I3.2,' VALID AT ',I5.4,3I3.2) - - IDATGS_COR = (JDATE(1) * 1000000) + (JDATE(2) * 10000) + - $ (JDATE(3) * 100) + JDATE(5) - -C VALID DATES MUST MATCH -C ---------------------- - - IF(IDATEP.NE.IDATGS_COR) GO TO 901 - -C------------------------------------------ - - ALLOCATE (iar12z(im,jm), iar13p(im,jm)) - ALLOCATE (iar14t(im,jm,km), iar15u(im,jm,km), - $ iar16v(im,jm,km), iar17q(im,jm,km), - $ iarpsl(im,jm,km), iarpsi(im,jm,km+1), - $ iarpsd(im,jm,km)) - - error=nf90_get_att(ncid, NF90_GLOBAL, "grid", grid) - IF (grid == "gaussian")THEN - IDRT=4 - ENDIF - - NVCOORD=2 !for idvc=2, nvcoord=2: hybrid interface a and b - IDVC=2 !(1 for sigma and 2 for hybrid) - IDSL=1 !(1 for phillips or 2 for mean) - idvm=1 - jcap = -9999 - - SFCPRESS_ID = MOD(IDVM,TEN) - THERMODYN_ID = MOD(IDVM/TEN,TEN) - - IF(IDVC == 3 .AND. THERMODYN_ID == 3) THEN - KMAXS = (NTRAC+1)*KMAX + 2 - ELSE - KMAXS = 2*KMAX + 2 - NTRAC = 1 - ENDIF - -C DEFINE THE OTHER RESOLUTION PARAMETERS -C -------------------------------------- - - DLAT = 180./(JMAX-1) - DLON = 360./IMAX - - PRINT 2, JCAP,IMAX,JMAX,KMAX,kmaxs,DLAT,DLON,idvc - - 2 FORMAT(/' --> GBLEVENTS: GLOBAL MODEL SPECS: T',I5,' ', - $ I5,' x ',I5,' GRID WITH ',I3,' LEVELS ',I3, - $' SCALARS -------> ',F5.2,' X ',F5.2,' VERT. ', - $ 'COORD ID IS: ',I0) - - print *,'in netcdf sfcpress_id=',sfcpress_id,' thermodyn_id=', - $ thermodyn_id,' ntrac=',ntrac - print *,' idvc=',idvc,' idsl=',idsl,' idvm=',idvm,' nvcoord=', - $ nvcoord - - ALLOCATE(VCOORD(km+1,3)) - VCOORD=0.0 - allocate(ak(km+1)) - allocate(bk(km+1)) - - error=nf90_get_att(ncid, NF90_GLOBAL, "ak", ak) - error=nf90_get_att(ncid, NF90_GLOBAL, "bk", bk) - - do k=1,km+1 - kr=km+2-k - vcoord(k,1) = ak(kr) - vcoord(k,2) = bk(kr) - end do - - deallocate(ak,bk) - - allocate(temp(im,jm)) - error=nf90_inq_varid(ncid, 'hgtsfc', id_var) - error=nf90_get_var(ncid, id_var, temp) - - IAR12Z(:,:)=temp(:,:) - call gblevn11d(imax,jmax,IAR12Z) - - - error=nf90_inq_varid(ncid, 'pressfc', id_var) - error=nf90_get_var(ncid, id_var, temp) - - IAR13P(:,:)=temp(:,:)*0.01 - call gblevn11d(imax,jmax,IAR13P) - - deallocate(temp) - -!! write(6,111)'hgtsfc', 1,minval(iar12z),maxval(iar12z) -!! write(6,111)'pressfc',1,minval(iar13p),maxval(iar13p) -!! 111 format(a8,1x,i3,1x,2(g13.6,1x)) - - - allocate(temp3d(im,jm,km)) - error=nf90_inq_varid(ncid, 'tmp', id_var) - error=nf90_get_var(ncid, id_var, temp3d) - DO K4=1,km - kr=km+1-k4 - IAR14T(:,:,K4)=temp3d(:,:,kr) - call gblevn11d(imax,jmax,IAR14T(1,1,K4)) -!! write(6,111)'tmp',k4,minval(iar14t(:,:,k4)), -!! & maxval(iar14t(:,:,k4)) - ENDDO - - error=nf90_inq_varid(ncid, 'ugrd', id_var) - error=nf90_get_var(ncid, id_var, temp3d) - DO K4=1,km - kr=km+1-k4 - IAR15U(:,:,k4)=temp3d(:,:,kr) - call gblevn11d(imax,jmax,IAR15U(1,1,K4)) -!! write(6,111)'ugrd',k4,minval(iar15u(:,:,k4)), -!! & maxval(iar15u(:,:,k4)) - ENDDO - - error=nf90_inq_varid(ncid, 'vgrd', id_var) - error=nf90_get_var(ncid, id_var, temp3d) - DO K4=1,km - kr=km+1-k4 - IAR16V(:,:,k4)=temp3d(:,:,kr) - call gblevn11d(imax,jmax,IAR16V(1,1,K4)) -!! write(6,111)'vgrd',k4,minval(iar16v(:,:,k4)), -!! & maxval(iar16v(:,:,k4)) - ENDDO - - error=nf90_inq_varid(ncid, 'spfh', id_var) - error=nf90_get_var(ncid, id_var, temp3d) - DO K4=1,km - kr=km+1-k4 - IAR17Q(:,:,k4)=max(0.0,temp3d(:,:,kr)*1.e6) - call gblevn11d(imax,jmax,IAR17Q(1,1,K4)) -!! write(6,111)'spfh',k4,minval(iar17q(:,:,k4)), -!! & maxval(iar17q(:,:,k4)) - ENDDO - - deallocate(temp3d) -! - - -!--compute Tv from temp - DO K=1,km - DO J=1,jm - DO i=1,im - TFAC=ONER+FV*MAX(IAR17Q(I,J,K)*1.0E-6,QMIN) - IAR14T(I,J,K)=IAR14T(I,J,K)*TFAC - ENDDO - ENDDO - ENDDO - -! print *,'in getnemsio,IAR15U.5,',maxval(IAR15U),minval(IAR15U) - - - -!--COMPUTE PM AND PD - ALLOCATE (psfc(im,jm), tv(im,jm,km)) - ALLOCATE (WRK1(im*jm,km), WRK2(im*jm,km+1)) - - imjm4=im*jm - km4=km - psfc(:,:) = iar13p(:,:)*100. - tv(:,:,:) = iar14t(:,:,:) - - IF(THERMODYN_ID == 3 .AND. IDVC == 3) THEN - tv(:,:,:) = tv(:,:,:) / CPI(1) - print *,' cpi(1)=',cpi(1) - ENDIF - - CALL SIGIO_MODPR(imjm4,imjm4,km4,NVCOORD,IDVC,IDSL,VCOORD,IRET, - $ psfc,tv,PM=WRK1,PD=WRK2(1,2)) - DO J=1,jm - JJ = (J-1)*im - DO I=1,im - WRK2(I+JJ,1) = psfc(i,j) ! in Pa - ENDDO - ENDDO - DO L=1,km - WRK2(:,L+1) = WRK2(:,L) - WRK2(:,L+1) ! in Pa - ENDDO - - DO L=1,km - DO J=1,jm - JJ = (J-1)*im - DO I=1,im - IARPSL(I,J,L) = WRK1(I+JJ,L)*0.01 ! 3D layer pres(hPa) - ENDDO - ENDDO - ENDDO - DO L=1,km+1 - DO J=1,jm - JJ = (J-1)*im - DO I=1,im - IARPSI(I,J,L) = WRK2(I+JJ,L)*0.01 ! 3D interface pressure - ! (hPa) - ENDDO - ENDDO - ENDDO - - error=nf90_close(ncid) -! print*,'GBLEVN13 IARPSI,',maxval(IARPSI),minval(IARPSI) -! print*,'GBLEVN13 IARPSL,',maxval(IARPSL),minval(IARPSL) - -! - CALL GETLATS(IDRT) - - RETURN - - 901 CONTINUE - PRINT 9901, IDATEP,IDATGS_COR - 9901 FORMAT(/' ##GBLEVENTS/GBLEVN13 - NETCDF INPUT GLOBAL FILE DATE', - $ ' (',I4.4,3(I2.2),'), DOES NOT MATCH PREPBUFR FILE CENTER ', - $ 'DATE (',I10,') -STOP 68'/) - CALL ERREXIT(68) - - END - -!----------------------------------------------------------------- -!----------------------------------------------------------------- - subroutine getlats(idrt) - - USE GBLEVN_MODULE -! integer :: jmax -! real,allocatable :: slat(:),wlat(:) -! real :: rad2deg - real,allocatable :: slats(:) - real(4) slat(jmax),wlat(jmax),rad2deg -!get gaussian or regular latitude array based on idrt - -! print *,'in getlats,idrt=',idrt - call splat(idrt,jmax,slat,wlat) - rad2deg=180./acos(-1.) - allocate(slats(jmax)); - rad2deg=180./acos(-1.) - slats(:)=-asin(slat(:))*rad2deg - dlat=180./float(jmax-1) - - return - end diff --git a/src/list_of_files.cmake b/src/list_of_files.cmake index 3571b8e5..29ebd1fa 100644 --- a/src/list_of_files.cmake +++ b/src/list_of_files.cmake @@ -1,6 +1,5 @@ set(fortran_src ${CMAKE_CURRENT_SOURCE_DIR}/args_mod.f - ${CMAKE_CURRENT_SOURCE_DIR}/gblevents.f ${CMAKE_CURRENT_SOURCE_DIR}/getgbens.f ${CMAKE_CURRENT_SOURCE_DIR}/isrchne.f ${CMAKE_CURRENT_SOURCE_DIR}/iw3mat.f @@ -213,6 +212,3 @@ set(fortran_src ${CMAKE_CURRENT_SOURCE_DIR}/xstore.f ) -set(c_src - ${CMAKE_CURRENT_SOURCE_DIR}/mova2i.c - ) diff --git a/src/mova2i.c b/src/mova2i.c deleted file mode 100644 index 239c494c..00000000 --- a/src/mova2i.c +++ /dev/null @@ -1,68 +0,0 @@ -/*$$$ SUBPROGRAM DOCUMENTATION BLOCK -C . . . . -C SUBPROGRAM: mova2i Moves a bit string from a char*1 to int -C PRGMMR: Gilbert ORG: W/NP11 DATE: 02-08-15 -C -C ABSTRACT: This Function copies a bit string from a Character*1 variable -C to an integer variable. It is intended to replace the Fortran Intrinsic -C Function ICHAR, which only supports 0 <= ICHAR(a) <= 127 on the -C IBM SP. If "a" is greater than 127 in the collating sequence, -C ICHAR(a) does not return the expected bit value. -C This function can be used for all values 0 <= ICHAR(a) <= 255. -C -C PROGRAM HISTORY LOG: -C 98-12-15 Gilbert -C -C USAGE: I = mova2i(a) -C -C INPUT ARGUMENT : -C -C a - Character*1 variable that holds the bitstring to extract -C -C RETURN ARGUMENT : -C -C mova2i - Integer value of the bitstring in character a -C -C REMARKS: -C -C None -C -C ATTRIBUTES: -C LANGUAGE: C -C MACHINE: IBM SP - -C -C$$$i*/ - -#ifdef CRAY90 - #include - int MOVA2I(unsigned char *a) -#endif -#ifdef HP - int mova2i(unsigned char *a) -#endif -#ifdef SGI - int mova2i_(unsigned char *a) -#endif -#ifdef LINUX - int mova2i_(unsigned char *a) -#endif -#ifdef APPLE - int mova2i_(unsigned char *a) -#endif -#ifdef LINUXF90 - int MOVA2I(unsigned char *a) -#endif -#ifdef VPP5000 - int mova2i_(unsigned char *a) -#endif -#ifdef IBM4 - int mova2i(unsigned char *a) -#endif -#ifdef IBM8 - long long int mova2i(unsigned char *a) -#endif - -{ - return (int)(*a); -} diff --git a/src/mova2i.f b/src/mova2i.f index 1d9408af..873620a8 100644 --- a/src/mova2i.f +++ b/src/mova2i.f @@ -37,7 +37,7 @@ C> LANGUAGE: IBM XL FORTRAN C> MACHINE: IBM SP C> - Integer Function mova2i(a) + Integer Function mova2i(a) C integer mold character(len=1) a