ecearth3/sources/tm5mp/base/tmm_cf.F90

!###############################################################################
!
! Converstion from TM5 names/units to CF standard names/units 
!
! USAGE
!
!   use TMM_CF
!
!   ! initialize module:
!   !   o read name of CF standard table from rcfile key 'cf-standard-name-table'
!   call TMM_CF_Init( rcf, status )
!     type(TRcFile), intent(in)           ::  rcf
!     integer, intent(out)                ::  status
!
!   ! get standard units for standard name:
!   call TMM_CF_Standard_Units( cf_standard_name, cf_units, status )
!     character(len=*), intent(in)                ::  cf_standard_name
!     character(len=*), intent(out)               ::  cf_units
!     integer, intent(out)                        ::  status
!
!   ! convert from TM5 variable name to CF standard name:
!   call TMM_CF_Convert_Name( tm5_name, cf_standard_name, status )
!     character(len=*), intent(in)                ::  tm5_name
!     character(len=*), intent(out)               ::  cf_standard_name
!     integer, intent(out)                        ::  status
!
!   ! get conversion factor from TM5 units to CF standard unit:
!   call subroutine TMM_CF_Convert_Units( tm5_units, cf_units, tm5_to_cf_factor, status )
!     character(len=*), intent(in)                ::  tm5_units
!     character(len=*), intent(in)                ::  cf_units
!     real, intent(out)                           ::  tm5_to_cf_factor
!     integer, intent(out)                        ::  status
!
!   ! done with module:
!   call TMM_CF_Done( status )
!     integer, intent(out)                        ::  status
!
!
! EXTERNAL LIBRARIES
!
!   Module uses the 'UDUnits' library (either old FORTRAN version 1.x [with_udunits1], or C version
!    2.x [with_udunits2]).
!
!###############################################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#ifdef with_udunits1
#define IF_UDUNITS_NOTOK_RETURN(action) if (status/=UDUNITS_NOERR) then; gol=trim(UDUnits_StrError(status)); call goErr; TRACEBACK; action; return; end if
#endif  
!
#include "tmm.inc"
!
!###############################################################################

module TMM_CF

  use GO     , only : gol, goErr, goPr
#ifdef with_udunits1
  use UDUnits, only : UDUNITS_NOERR, UDUnits_StrError
#endif
#ifdef with_udunits2
  use ISO_C_BINDING
#endif
  use Standard_Name_Table, only : T_Standard_Name_Table

  implicit none
  private
  
  public  ::  TMM_CF_Init, TMM_CF_Done
  public  ::  TMM_CF_Standard_Units
  public  ::  TMM_CF_Convert_Name
  public  ::  TMM_CF_Convert_Units
  
  ! --- const ------------------------------------
  
  character(len=*), parameter  ::  mname = 'TMM_CF'

  ! rcfile key with name of standard table:
  character(len=*), parameter  ::  rckey_cf_table = 'cf-standard-name-table'

  ! --- local ------------------------------------

  logical                               ::  with_cf_table
  type(T_Standard_Name_Table), save     ::  cf_table

contains

  ! ==============================================================

  subroutine TMM_CF_Init( rcf, status )
  
    use GO                 , only : TrcFile, ReadRc
    use Standard_Name_Table, only : Standard_Name_Table_Init
#ifdef with_udunits1
    use UDUnits            , only : UDUnits_Init
#endif

    ! --- in/out ---------------------------------
    
    type(TRcFile), intent(in)           ::  rcf
    integer, intent(out)                ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_CF_Init'
    
    ! --- local ----------------------------------
    
    character(len=512)      ::  cf_standard_name_table
    
    ! --- begin ----------------------------------
    
    ! read full path of cf table:
    call ReadRc( rcF, rckey_cf_table, cf_standard_name_table, &
                           status, default='None' )
    IF_ERROR_RETURN(status=1)
    ! set flag to false if key was not found:
    with_cf_table = trim(cf_standard_name_table) /= 'None'
    ! key provided ?
    if ( with_cf_table ) then
       ! standard names and units:
       call Standard_Name_Table_Init( cf_table, trim(cf_standard_name_table), status )
       if (status/=0) then
          write(gol,*) "problem with file '"//trim(cf_standard_name_table)//"'"; call goErr
          write(gol,*) "this file is defined through the '"//trim(rckey_cf_table)//"' rc key in your meteo file"; call goErr
          IF_NOTOK_RETURN(status=1)
       endif
    end if

#ifdef with_udunits1
    ! initialize UDUnits module:
    call UDUnits_Init( status )
    IF_UDUNITS_NOTOK_RETURN(status=1)
#endif

    ! ok
    status = 0
    
  end subroutine TMM_CF_Init

  ! ***

  subroutine TMM_CF_Done( status )
  
    use Standard_Name_Table, only : Standard_Name_Table_Done
#ifdef with_udunits1
    use UDUnits            , only : UDUnits_Done
#endif

    ! --- in/out ---------------------------------
    
    integer, intent(out)                ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_CF_Done'
    
    ! --- local ----------------------------------
    
    ! --- begin ----------------------------------

    ! standard table present ?
    if ( with_cf_table ) then
      ! done with standard names and units:
      call Standard_Name_Table_Done( cf_table, status )
      IF_NOTOK_RETURN(status=1)
    end if

#ifdef with_udunits1
    ! done with UDUnits module:
    call UDUnits_Done( status )
    IF_UDUNITS_NOTOK_RETURN(status=1)
#endif

    ! ok
    status = 0
    
  end subroutine TMM_CF_Done

  ! ***

  subroutine TMM_CF_Standard_Units( cf_standard_name, cf_units, status )

    ! --- in/out ---------------------------------
    
    character(len=*), intent(in)                ::  cf_standard_name
    character(len=*), intent(out)               ::  cf_units
    integer, intent(out)                        ::  status

    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_CF_Standard_Units'
    
    ! --- local ----------------------------------
    
    integer                         ::  i, j

    ! --- begin ----------------------------------
    
    ! no unit yet ...
    cf_units = ''
    ! check for specials ...
    select case ( trim(cf_standard_name) )
      ! TM5 special fields
      case ( 'tm5_integrated_eastward_mass_flux_of_air'  ) ; cf_units = 'kg s-1'
      case ( 'tm5_integrated_northward_mass_flux_of_air' ) ; cf_units = 'kg s-1'
      case ( 'tm5_integrated_upward_mass_flux_of_air'    ) ; cf_units = 'kg s-1'
      ! CF standard fields
      case default
        ! *
        ! check ..
        if ( .not. with_cf_table ) then
          write (gol,'("No CF table specified with rcfile key `",a,"`")') &
                  trim(rckey_cf_table); call goErr
          TRACEBACK; status=1; return
        end if
        ! *
        ! loop over all entries:
        do i = 1, size(cf_table%entry)
          ! compare ...
          if ( trim(cf_table%entry(i)%id) == trim(cf_standard_name) ) then
            ! copy values:
            cf_units = trim(cf_table%entry(i)%canonical_units)
            ! found!
            exit
          end if
        end do   ! CF table entries
        ! *
        ! not found yet ? then check alia:
        if ( len_trim(cf_units) == 0 ) then
          ! loop over all alia:
          do j = 1, size(cf_table%alias)
            ! compare ...
            if ( trim(cf_table%alias(j)%id) == trim(cf_standard_name) ) then
              ! match; now search real entries:
              do i = 1, size(cf_table%entry)
                ! compare ...
                if ( trim(cf_table%entry(i)%id) == trim(cf_table%alias(j)%entry_id) ) then
                  ! copy values:
                  cf_units = trim(cf_table%entry(i)%canonical_units)
                  ! leave:
                  exit
                end if
              end do   ! CF table entries
              ! found!
              exit
            end if
          end do   ! CF table alia
        end if
        ! *
    end select
    ! not found ?
    if ( len_trim(cf_units) == 0 ) then
      write (gol,'("id not found in cf standard name table : ",a)') trim(cf_standard_name); call goErr
      TRACEBACK; status=1; return
    end if

    ! ok
    status = 0
    
  end subroutine TMM_CF_Standard_Units


  ! ***


  subroutine TMM_CF_Convert_Name( tm5_name, cf_standard_name, status )

    ! --- in/out ---------------------------------
    
    character(len=*), intent(in)                ::  tm5_name
    character(len=*), intent(out)               ::  cf_standard_name
    integer, intent(out)                        ::  status

    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_CF_Convert_Name'
    
    ! --- local ----------------------------------
    
    ! --- begin ----------------------------------
    
    ! convert from TM5 internal name to CF standard name:
    select case ( trim(tm5_name) )
      case ( 'lon'        ) ; cf_standard_name = 'longitude'
      case ( 'lat'        ) ; cf_standard_name = 'latitude'
      case ( 'sp'         ) ; cf_standard_name = 'surface_air_pressure'
      case ( 'mfu'        ) ; cf_standard_name = 'tm5_integrated_eastward_mass_flux_of_air'
      case ( 'mfv'        ) ; cf_standard_name = 'tm5_integrated_northward_mass_flux_of_air'
      case ( 'mfw'        ) ; cf_standard_name = 'tm5_integrated_upward_mass_flux_of_air'
      case ( 'tsp'        ) ; cf_standard_name = 'tendency_of_surface_air_pressure'
      case ( 'oro'        ) ; cf_standard_name = 'geopotential'
      case ( 'T'          ) ; cf_standard_name = 'air_temperature'
      case ( 'Q'          ) ; cf_standard_name = 'specific_humidity'
      case ( 'CLWC'       ) ; cf_standard_name = 'mass_fraction_of_cloud_liquid_water_in_air'
      case ( 'CIWC'       ) ; cf_standard_name = 'mass_fraction_of_cloud_ice_in_air'
      case ( 'CC'         ) ; cf_standard_name = 'cloud_area_fraction_in_atmosphere_layer'
      case ( 'CCO'        ) ; cf_standard_name = 'cloud_area_fraction_in_atmosphere_layer'       ! dummy for overhead cloud cover
      case ( 'CCU'        ) ; cf_standard_name = 'cloud_area_fraction_in_atmosphere_layer'       ! dummy for underfeet cloud cover
      case ( 'eu'         ) ; cf_standard_name = 'atmosphere_updraft_convective_mass_flux'       ! dummy for entrainement updraft
      case ( 'ed'         ) ; cf_standard_name = 'atmosphere_downdraft_convective_mass_flux'     ! dummy for entrainement downdraft
      case ( 'du'         ) ; cf_standard_name = 'atmosphere_updraft_convective_mass_flux'       ! dummy for detrainement updraft
      case ( 'dd'         ) ; cf_standard_name = 'atmosphere_downdraft_convective_mass_flux'     ! dummy for detrainement downdraft
      case ( 'cloud_base' ) ; cf_standard_name = 'model_level_number_at_convective_cloud_base'
      case ( 'cloud_top'  ) ; cf_standard_name = 'model_level_number_at_convective_cloud_top'
      case ( 'cloud_lfs'  ) ; cf_standard_name = 'model_level_number_at_convective_cloud_top'    ! dummy for level-of-free-sinking
      case ( 'ssr'        ) ; cf_standard_name = 'surface_net_upward_shortwave_flux'
      case ( 'ssrd'       ) ; cf_standard_name = 'surface_net_downward_shortwave_flux'
      case ( 'str'        ) ; cf_standard_name = 'surface_net_upward_longwave_flux'
      case ( 'strd'       ) ; cf_standard_name = 'surface_net_downward_longwave_flux'
      case ( 'lsm'        ) ; cf_standard_name = 'land_area_fraction'
      case ( 'albedo'     ) ; cf_standard_name = 'surface_albedo'
      case ( 'sr'         ) ; cf_standard_name = 'surface_roughness_length'
      case ( 'srols'      ) ; cf_standard_name = 'surface_roughness_length'
      case ( 'ci'         ) ; cf_standard_name = 'sea_ice_area_fraction'
      case ( 'sst'        ) ; cf_standard_name = 'sea_surface_temperature'
      case ( 'u10m'       ) ; cf_standard_name = 'eastward_wind'
      case ( 'v10m'       ) ; cf_standard_name = 'northward_wind'
      case ( 'g10m'       ) ; cf_standard_name = 'wind_speed_of_gust'
      case ( 'fg10'       ) ; cf_standard_name = 'wind_speed_of_gust'
      case ( 'd2m'        ) ; cf_standard_name = 'dew_point_temperature'
      case ( 't2m'        ) ; cf_standard_name = 'air_temperature'
      case ( 'skt'        ) ; cf_standard_name = 'canopy_temperature'
      case ( 'blh'        ) ; cf_standard_name = 'atmosphere_boundary_layer_thickness'
      case ( 'sshf'       ) ; cf_standard_name = 'surface_downward_sensible_heat_flux'
      case ( 'slhf'       ) ; cf_standard_name = 'surface_downward_latent_heat_flux'
      case ( 'ewss'       ) ; cf_standard_name = 'surface_downward_eastward_stress'
      case ( 'nsss'       ) ; cf_standard_name = 'surface_downward_northward_stress'
      case ( 'cp'         ) ; cf_standard_name = 'lwe_convective_precipitation_rate'
      case ( 'lsp'        ) ; cf_standard_name = 'lwe_large_scale_precipitation_rate'
      case ( 'sf'         ) ; cf_standard_name = 'lwe_thickness_of_snowfall_amount'
      case ( 'sd'         ) ; cf_standard_name = 'lwe_thickness_of_surface_snow_amount'
      case ( 'src'        ) ; cf_standard_name = 'lwe_thickness_of_canopy_water_amount'
      case ( 'swvl1'      ) ; cf_standard_name = 'volume_fraction_of_condensed_water_in_soil'
      case ( 'stl1'       ) ; cf_standard_name = 'soil_temperature'
      case ( 'tv01', 'tv02', 'tv03', 'tv04', 'tv05', &
             'tv06', 'tv07', 'tv08', 'tv09', 'tv10', &
             'tv11', 'tv12', 'tv13', 'tv14', 'tv15', &
             'tv16', 'tv17', 'tv18', 'tv19', 'tv20' )
                              cf_standard_name = 'vegetation_area_fraction'
      case ( 'cvl', 'cvh' ) ; cf_standard_name = 'vegetation_area_fraction'
      case ( 'K'          ) ; cf_standard_name = 'atmosphere_heat_diffusivity' ! 162.109 grib code = Turbulent diffusion coefficient for heat
      case default
        !write (gol,'("do not know cf standard name for tm5 name : ",a)') trim(tm5_name); call goErr
        !TRACEBACK; status=1; return
        ! assume name follows CF already ..
        cf_standard_name = trim(tm5_name)
    end select

    ! ok
    status = 0
    
  end subroutine TMM_CF_Convert_Name


  ! ***
  

  !
  ! Return factor for conversion from "tm5_units" to "cf_units"
  !

  SUBROUTINE TMM_CF_Convert_Units( tm5_units, cf_units, tm5_to_cf_factor, status )
  
#ifdef with_udunits1
    use UDUnits, only : UDUnits_ConversionFactor
#endif
#ifdef with_udunits2
    use f_udunits_2, only : UT_UNIT_PTR, CV_CONVERTER_PTR,   UT_SYSTEM_PTR, f_ut_read_xml
    use f_udunits_2, only : UT_ASCII,    f_ut_get_converter, f_ut_parse,    f_cv_convert
    use f_udunits_2, only : f_ut_free,   f_ut_free_system,   f_cv_free
#endif  
    
    ! --- in/out ---------------------------------    

    character(len=*), intent(in)  ::  tm5_units
    character(len=*), intent(in)  ::  cf_units
    real,             intent(out) ::  tm5_to_cf_factor
    integer,          intent(out) ::  status
    
    ! --- local ----------------------------------

#ifdef with_udunits2
    type(UT_SYSTEM_PTR)    :: sys
    type(CV_CONVERTER_PTR) :: conv1
    type(UT_UNIT_PTR)      :: unit1, unit2
    real                   :: a
    integer                :: charset
#endif    
    
    ! --- const ----------------------------------

    character(len=*), parameter  ::  rname = mname//'/TMM_CF_Convert'
    
    ! --- local ----------------------------------
    
    character(len=64)       ::  conversion
    
    ! --- begin ----------------------------------
    
    ! same?
    if ( trim(tm5_units) == trim(cf_units) ) then
    
      ! no conversion:
      tm5_to_cf_factor = 1.0
    
    else

      ! conversion:
      write (conversion,'(a," -> ",a)') trim(tm5_units), trim(cf_units)
      ! known conversions ...
      select case ( trim(conversion) )
      
        case ( 'kg s-1 -> kg/s', &
               'm2 s-2 -> m m/s2', &
               '1 -> kg/kg', &
               'kg m-2 s-1 -> kg/m2/s', &
               'W m-2 -> W/m2', &
               'Pa -> N/m2', &
               'm s-1 -> m/s' )

          ! no conversion needed:
          tm5_to_cf_factor = 1.0
          
        case ( '1 -> %' )
          tm5_to_cf_factor = 100.0
      
        case default

#ifdef with_udunits1

          ! unit conversion factor:
          call UDUnits_ConversionFactor( trim(tm5_units), trim(cf_units), tm5_to_cf_factor, status )
          if ( status /= UDUNITS_NOERR ) then
            write (gol,'("from conversion of TM5 units to CF units:")'); call goErr
            write (gol,'("  TM5 units    : ",a)') trim(tm5_units); call goErr
            write (gol,'("  CF units     : ",a)') trim(cf_units); call goErr
            TRACEBACK; status=1; return
          end if

#elif defined(with_udunits2)
    
          charset = UT_ASCII
          sys = f_ut_read_xml("")

          unit1 = f_ut_parse(sys,tm5_units,charset)
          if ( .not. c_associated(unit1%ptr) ) then
             write (gol,'("from parsing input units `",a,"`")') trim(tm5_units); call goErr
             TRACEBACK; status=1; return
          endif

          unit2 = f_ut_parse(sys,cf_units,charset)
          if ( .not. c_associated(unit2%ptr) ) then
             write (gol,'("from parsing target units `",a,"`")') trim(cf_units); call goErr
             TRACEBACK; status=1; return
          endif

          conv1 = f_ut_get_converter(unit1,unit2)
          if ( .not. c_associated(conv1%ptr)) then
             write (gol,'("from converting input units `",a,"` to target units `",a,"`")') &
                      trim(tm5_units), trim(cf_units); call goErr
             TRACEBACK; status=1; return
          endif

          a = 1.0
          tm5_to_cf_factor = f_cv_convert(conv1,a)

          call f_ut_free(unit1)
          call f_ut_free(unit2)
          call f_cv_free(conv1)
          call f_ut_free_system(sys)

#else
          ! dummy assignment to avoid compiler warnings:
          tm5_to_cf_factor = 1.0
          ! error ...
          write (gol,'("Unsupported unit conversion `",a,"` required.")') trim(conversion); call goPr
          write (gol,'("Either:")'); call goErr
          write (gol,'(" 1. hardcode the conversion in routine `",a,"` in `",a,"`;")') &
                              trim(rname), __FILE__; call goErr
          write (gol,'(" 2. enable conversion using UDunits library by defining one of the macros:")'); call goErr
          write (gol,'("       with_udunits1")'); call goErr
          write (gol,'("       with_udunits2")'); call goErr
          write (gol,'("    and link with appropriate library.")'); call goErr
          TRACEBACK; status=1; return
#endif

      end select  ! known conversions
      
    end if  ! same umits
        
    ! ok
    status = 0
    
  END SUBROUTINE TMM_CF_Convert_Units

END MODULE TMM_CF