Skip to content
Snippets Groups Projects
sfx_run.f90 22.9 KiB
Newer Older
    !< @brief sfx driver 

    ! directives list
    ! --------------------------------------------------------------------------------
    implicit none
    private
    ! --------------------------------------------------------------------------------

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: set_run
    public :: run_dataset
    ! --------------------------------------------------------------------------------

contains 

    !> @brief run sfx on dataset
    ! ----------------------------------------------------------------------------
    subroutine run_dataset(filename_out, dataset, model, ierr)

        ! modules used
        ! --------------------------------------------------------------------------------
        use sfx_phys_const
        use sfx_common
        use sfx_config
Evgeny Mortikov's avatar
Evgeny Mortikov committed
        use sfx_surface
        use sfx_io
        use sfx_data
    
        use sfx_esm, only: &
                get_surface_fluxes_vec_esm => get_surface_fluxes_vec, &
                numericsType_esm => numericsType
        use sfx_log, only: &
                get_surface_fluxes_vec_log => get_surface_fluxes_vec, &
                numericsType_log => numericsType
        use sfx_most, only: &
                get_surface_fluxes_vec_most => get_surface_fluxes_vec, &
                numericsType_most => numericsType
        use sfx_sheba, only: &
                get_surface_fluxes_vec_sheba => get_surface_fluxes_vec, &
                numericsType_sheba => numericsType
        use sfx_sheba_noniterative, only: &
                get_surface_fluxes_vec_sheba_noit => get_surface_fluxes_vec, &
                numericsType_sheba_noit => numericsType
        use sfx_most_snow, only: &
                get_surface_fluxes_vec_most_snow => get_surface_fluxes_vec, &
                numericsType_most_snow => numericsType
        use sfx_sheba_coare, only: &
                get_surface_fluxes_vec_sheba_coare => get_surface_fluxes_vec, &
                numericsType_sheba_coare => numericsType
        ! --------------------------------------------------------------------------------
    
        character(len=*), intent(in) :: filename_out
        type(sfxDatasetType), intent(in) :: dataset
        integer, intent(in) :: model

        integer, intent(out) :: ierr
    
        ! input/output model data
        ! --------------------------------------------------------------------------------
        type(meteoDataVecType) :: meteo         !< meteorological data (input)
        type(meteoDataType) :: meteo_cell
    
        type(sfxDataVecType) :: sfx             !< surface fluxes (output)
    
        type(numericsType_esm) :: numerics_esm      !< surface flux module (ESM) numerics parameters
        type(numericsType_log) :: numerics_log      !< surface flux module (LOG) numerics parameters
        type(numericsType_most) :: numerics_most    !< surface flux module (MOST) numerics parameters
        type(numericsType_sheba) :: numerics_sheba  !< surface flux module (SHEBA) numerics parameters
        type(numericsType_sheba_noit) :: numerics_sheba_noit  !< surface flux module (SHEBA) numerics parameters
        type(numericsType_most_snow) :: numerics_most_snow    !< surface flux module (MOST_SNOW) numerics parameters
        type(numericsType_sheba_coare) :: numerics_sheba_coare   !< surface flux module (MOST_SNOW) numerics parameters
        integer :: num                          !< number of 'cells' in input
        ! --------------------------------------------------------------------------------
    
        ! local variables
        ! --------------------------------------------------------------------------------
        integer :: i
        integer :: io, status
        ! --------------------------------------------------------------------------------
        write(*, *) ' Running SFX:'
        write(*, '(a,a)') '    model            = ', trim(get_model_tag(model))
        write(*, '(a,a)') '    dataset          = ', trim(get_dataset_tag(dataset%id))
        write(*, '(a,a)') '    filename[IN]     = ', trim(dataset%filename)
        write(*, '(a,a)') '    filename[OUT]    = ', trim(filename_out)
Evgeny Mortikov's avatar
Evgeny Mortikov committed
        write(*, '(a,g0)') '    surface type     = ', trim(get_surface_tag(dataset%surface))
        write(*, '(a,g0)') '    h                = ', dataset%h
        write(*, '(a,g0)') '    z0(m)            = ', dataset%z0_m
        write(*, '(a,g0)') '    z0(h)            = ', dataset%z0_h
        write(*, '(a,g0)') '    lai              = ', dataset%lai
        write(*, '(a,g0)') '    depth            = ', dataset%depth
        write(*, '(a,g0)') '    surface_type     = ', dataset%surface_type
        ! --- define number of elements
        open(newunit = io, file = dataset%filename, iostat = status, status ='old')
        if (status /= 0) then
            write(*, *) ' FAILURE! > unable to open file: ', trim(dataset%filename)
            return
        end if
    
        num = 0
        status = 0
        do while (status.eq.0)
            read (io, *, iostat = status) meteo_cell%U, meteo_cell%dT, meteo_cell%Tsemi, meteo_cell%dQ
            num = num + 1
        enddo
        num = num - 1
    
        close(io)
    
        ! --- print number of elements in dataset
        write(*, '(a,g0)') '    size             = ', num
        if (dataset%nmax > 0) then
            write(*, '(a,g0)') '    nmax             = ', dataset%nmax
            num = min(num, dataset%nmax)
        end if
    
    
        ! --- allocate input & output data
        call allocate_meteo_vec(meteo, num)
        call allocate_sfx_vec(sfx, num)
    
    
        ! --- setting height & roughness
        meteo_cell%h = dataset%h
        meteo_cell%z0_m = dataset%z0_m
        meteo_cell%depth = dataset%depth
        meteo_cell%lai = dataset%lai
        meteo_cell%surface_type = dataset%surface_type
        ! --- read input data
        open(newunit = io, file = dataset%filename, iostat = status, status = 'old')
        if (status /= 0) then
            write(*, *) ' FAILURE! > unable to open file: ', trim(dataset%filename)
            return
        end if
        do i = 1, num
            read(io, *) meteo_cell%U,  meteo_cell%dT,  meteo_cell%Tsemi,  meteo_cell%dQ
    
            meteo%h(i) = meteo_cell%h
            meteo%depth(i)=meteo_cell%depth 
            meteo%lai(i)=meteo_cell%lai
            meteo%surface_type(i)=meteo_cell%surface_type 

            meteo%U(i) = meteo_cell%U
            meteo%dT(i) = meteo_cell%dT
            meteo%Tsemi(i) = meteo_cell%Tsemi
            meteo%dQ(i) = meteo_cell%dQ
            meteo%z0_m(i) = meteo_cell%z0_m
        ! --- calling flux module
        if (model == model_esm) then
            call get_surface_fluxes_vec_esm(sfx, meteo, numerics_esm, num)
        else if (model == model_log) then
            call get_surface_fluxes_vec_log(sfx, meteo, numerics_log, num)
        else if (model == model_most) then
            call get_surface_fluxes_vec_most(sfx, meteo, numerics_most, num)
        else if (model == model_sheba) then
            call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, num)
        else if (model == model_sheba_noit) then
            call get_surface_fluxes_vec_sheba_noit(sfx, meteo, numerics_sheba_noit, num)
        else if (model == model_most_snow) then
            call get_surface_fluxes_vec_most_snow(sfx, meteo, numerics_most_snow, num)
        else if (model == model_sheba_coare) then
            call get_surface_fluxes_vec_sheba_coare(sfx, meteo, numerics_sheba_coare, num)
        ! --- write output data
        call write_ascii_vec11(filename_out, &
            sfx%zeta, sfx%Rib, &
            sfx%Re, sfx%B, sfx%z0_m, sfx%z0_t, &
            sfx%Rib_conv_lim, &
            sfx%Cm,sfx%Ct, sfx%Km, sfx%Pr_t_inv, num, '(11(f10.4,3x))', status)
         if (status /= 0) then
            write(*, *) ' FAILURE! > unable to write to file: ', trim(filename_out)
            ierr = status
            return
        end if
        ! --- deallocate input & output data
        call deallocate_meteo_vec(meteo)
        call deallocate_sfx_vec(sfx)
    
    end subroutine
    
    !> @brief set sfx run on dataset
    ! ----------------------------------------------------------------------------
    subroutine set_run(filename_out, dataset, model, ierr)
    
        ! modules used
        ! --------------------------------------------------------------------------------
#ifdef USE_CONFIG_PARSER
        use iso_c_binding, only: C_NULL_CHAR
        use config_parser
#endif
    
        use sfx_common
        use sfx_config
Evgeny Mortikov's avatar
Evgeny Mortikov committed
        use sfx_surface
        ! --------------------------------------------------------------------------------
    
        character(len=:), allocatable, intent(out) :: filename_out
        type(sfxDatasetType), intent(out) :: dataset
        integer, intent(out) :: model

        integer, intent(out) :: ierr
    
        ! command line arguments
        ! --------------------------------------------------------------------------------
        integer :: num_args
        character(len = 128) :: arg
    
        character(len = 128), parameter :: arg_key_model = '--model'
        character(len = 128), parameter :: arg_key_dataset = '--dataset'
        character(len = 128), parameter :: arg_key_output = '--output'
        character(len = 128), parameter :: arg_key_nmax = '--nmax'
        character(len = 128), parameter :: arg_key_help = '--help'
        character(len = 128), parameter :: arg_key_config = "--config"
    
        integer :: is_output_set
        ! --------------------------------------------------------------------------------
    
        ! local variables
        ! --------------------------------------------------------------------------------
        integer :: i, status
        integer :: id, nmax
    
#ifdef USE_CONFIG_PARSER
        character, allocatable :: config_field(:)
#endif
        ! --------------------------------------------------------------------------------
        ! --- default model & dataset
        model = model_esm                           ! default = ESM
        call set_dataset(dataset, dataset_mosaic)   ! default = MOSAiC
        is_output_set = 0                           ! default = auto define output filename
        ! --- command line arguments processing
        num_args = command_argument_count()
        do i = 1, num_args
            call get_command_argument(i, arg)
            if (trim(arg) == trim(arg_key_help)) then
                write(*, *) ' sfx model, usage:'
                write(*, *) ' --help'
                write(*, *) '    print usage options'
                write(*, *) ' --model [key]'
                write(*, *) '    key = esm (default) || log || most || sheba || sheba_noit || most_snow || sheba_coare'
                write(*, *) ' --dataset [key]'
                write(*, *) '    key = mosaic (default) || irgason || sheba'
                write(*, *) '        = lake || papa || toga || user [filename] [param]'
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                write(*, *) '    param = [surface(type)] [h] [z0(m)] [z0(h)]'
                write(*, *) '    surface(type) = ocean | land | lake | snow' 
                write(*, *) ' --output [filename]'
                write(*, *) '    set output filename'
                write(*, *) ' --config [filename]'
                write(*, *) '    use configuration file'
                write(*, *) ' --nmax [value]'
                write(*, *) '    max number of data points > 0'
            end if
            if (trim(arg) == trim(arg_key_model)) then
                if (i == num_args) then
                    write(*, *) ' FAILURE! > missing model [key] argument'
                    ierr = 1        ! signal ERROR
                    return
                end if
                
                call get_command_argument(i + 1, arg)
                model = get_model_id(arg)
                if (model == -1) then
                    write(*, *) ' FAILURE! > unknown model [key]: ', trim(arg)
                    ierr = 1        ! signal ERROR
                    return
                end if
            else if (trim(arg) == trim(arg_key_dataset)) then
                if (i == num_args) then
                    write(*, *) ' FAILURE! > missing dataset [key] argument'
                    ierr = 1        ! signal ERROR
                    return
                end if
                
                call get_command_argument(i + 1, arg)
                id = get_dataset_id(arg)
                if (id == -1) then
                    write(*, *) ' FAILURE! > unknown dataset [key]: ', trim(arg)
                    ierr = 1        ! signal ERROR
                    return
                end if
                !< save nmax if previously set
                nmax = dataset%nmax
                call set_dataset(dataset, id)
                dataset%nmax = nmax
    
                if (dataset%id == dataset_user) then
                    !< @brief user-defined dataset
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                    if (i + 6 > num_args) then
                        write(*, *) ' FAILURE! > incorrect arguments for [user] dataset'
                        ierr = 1        ! signal ERROR
                        return
                    end if
    
                    !< filename
                    call get_command_argument(i + 2, dataset%filename)
Evgeny Mortikov's avatar
Evgeny Mortikov committed

                    !< surface type
                    call get_command_argument(i + 3, arg)
                    dataset%surface = get_surface_id(arg)
                    if (dataset%surface == -1) then
                        write(*, *) ' FAILURE! > unknown surface [key]: ', trim(arg)
                        ierr = 1        ! signal ERROR
                        return
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                    end if
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                    call get_command_argument(i + 4, arg)
                    call str2real(dataset%h, arg, status)
                    if (status /= 0) then
                        write(*, *) ' FAILURE! > expecting real h [value]'
                        ierr = 1        ! signal ERROR
                        return 
                    end if
                    if (dataset%h <= 0) then
                        write(*, *) ' FAILURE! > h [value] should be positive'
                        ierr = 1        ! signal ERROR
                        return
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                    call get_command_argument(i + 5, arg)
                    call str2real(dataset%z0_m, arg, status)
                    if (status /= 0) then
                        write(*, *) ' FAILURE! > expecting real z0(m) [value]'
                        ierr = 1        ! signal ERROR
                        return
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                    call get_command_argument(i + 6, arg)
                    call str2real(dataset%z0_h, arg, status)
                    if (status /= 0) then
                        write(*, *) ' FAILURE! > expecting real z0(h) [value]'
                        ierr = 1        ! signal ERROR
                        return
                    end if
                    
                end if 
            else if (trim(arg) == trim(arg_key_output)) then
                if (i == num_args) then
                    write(*, *) ' FAILURE! > missing output [key] argument'
                    ierr = 1        ! signal ERROR
                    return 
                end if
                call get_command_argument(i + 1, arg)
                filename_out = trim(arg)
                is_output_set = 1
            else if (trim(arg) == trim(arg_key_nmax)) then
                if (i == num_args) then
                    write(*, *) ' FAILURE! > missing nmax [key] argument'
                    ierr = 1        ! signal ERROR
                    return
                end if
                call get_command_argument(i + 1, arg)
                call str2int(dataset%nmax, arg, status)
                if (status /= 0) then
                    write(*, *) ' FAILURE! > expecting int nmax [value]'
                    ierr = 1        ! signal ERROR
                    return
                end if
                if (dataset%nmax <= 0) then
                    write(*, *) ' FAILURE! > nmax [value] should be positive'
                    ierr = 1        ! signal ERROR
                    return
                end if
            else if (trim(arg) == trim(arg_key_config)) then
                if (i == num_args) then
                    write(*, *) ' FAILURE! > missing configuration file [key] argument'
                    ierr = 1        ! signal ERROR
                    return
                end if
                
                call get_command_argument(i + 1, arg)
#ifdef USE_CONFIG_PARSER
                call c_config_run(trim(arg)//C_NULL_CHAR, status)
                if (status == 0) then
                    write(*, *) ' FAILURE! > unable to parse configuration file: ', trim(arg)
                    ierr = 1        ! signal ERROR
                    return
                end if 
    
                call c_config_is_varname("model.id"//C_NULL_CHAR, status)
                if (status /= 0) then
                    call c_config_get_string("model.id"//C_NULL_CHAR, config_field, status)
                    if (status == 0) then
                        ierr = 1        ! signal ERROR
                        return
                    end if
                    model = get_model_id(char_array2str(config_field))
                    if (model == -1) then
                        write(*, *) ' FAILURE! > unknown model [key]: ', trim(char_array2str(config_field))
                        ierr = 1        ! signal ERROR
                        return
                    end if
                end if
    
                call c_config_is_varname("dataset.id"//C_NULL_CHAR, status)            
                if (status /= 0) then
                    call c_config_get_string("dataset.id"//C_NULL_CHAR, config_field, status)
                    if (status == 0) then
                        ierr = 1        ! signal ERROR
                        return
                    end if
                    id = get_dataset_id(char_array2str(config_field))
                    if (id == -1) then
                        write(*, *) ' FAILURE! > unknown dataset [key]: ', trim(char_array2str(config_field))
                        ierr = 1        ! signal ERROR
                        return
                    end if
                    !< save nmax if previously set
                    nmax = dataset%nmax
                    call set_dataset(dataset, id)
                    dataset%nmax = nmax
    
                    call c_config_is_varname("dataset.filename"//C_NULL_CHAR, status)
                    if ((status /= 0).or.(dataset%id == dataset_user)) then
                        !< mandatory in user dataset
                        call c_config_get_string("dataset.filename"//C_NULL_CHAR, config_field, status)
                        if (status == 0) then
                            ierr = 1        ! signal ERROR
                            return
                        end if
                        dataset%filename = char_array2str(config_field)
                    end if
Evgeny Mortikov's avatar
Evgeny Mortikov committed

                    call c_config_is_varname("dataset.surface"//C_NULL_CHAR, status)
                    if ((status /= 0).or.(dataset%id == dataset_user)) then
                        !< mandatory in user dataset
                        call c_config_get_string("dataset.surface"//C_NULL_CHAR, config_field, status)
                        if (status == 0) then
                            ierr = 1        ! signal ERROR
                            return
                        end if
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                        dataset%surface = get_surface_id(char_array2str(config_field))
                        if (dataset%surface == -1) then
                            write(*, *) ' FAILURE! > unknown surface [key]: ', trim(char_array2str(config_field))
                            ierr = 1        ! signal ERROR
                            return
Evgeny Mortikov's avatar
Evgeny Mortikov committed
                        end if 
                    endif 
    
                    call c_config_is_varname("dataset.h"//C_NULL_CHAR, status)
                    if ((status /= 0).or.(dataset%id == dataset_user)) then
                        !< mandatory in user dataset
                        call c_config_get_float("dataset.h"//C_NULL_CHAR, dataset%h, status)
                        if (status == 0) then
                            ierr = 1        ! signal ERROR
                            return
                        end if
                    end if
    
                    call c_config_is_varname("dataset.z0_m"//C_NULL_CHAR, status)
                    if ((status /= 0).or.(dataset%id == dataset_user)) then
                        !< mandatory in user dataset
                        call c_config_get_float("dataset.z0_m"//C_NULL_CHAR, dataset%z0_m, status)
                        if (status == 0) then
                            ierr = 1        ! signal ERROR
                            return
                        end if
                    end if
    
                    call c_config_is_varname("dataset.z0_h"//C_NULL_CHAR, status)
                    if ((status /= 0).or.(dataset%id == dataset_user)) then
                        !< mandatory in user dataset
                        call c_config_get_float("dataset.z0_h"//C_NULL_CHAR, dataset%z0_h, status)
                        if (status == 0) then
                            ierr = 1        ! signal ERROR
                            return
                        end if
                    end if
                end if
    
                call c_config_is_varname("dataset.nmax"//C_NULL_CHAR, status)
                if (status /= 0) then
                    call c_config_get_int("dataset.nmax"//C_NULL_CHAR, dataset%nmax, status)
                    if (status == 0) then
                        ierr = 1        ! signal ERROR
                        return
                    end if
                end if 
            
                call c_config_is_varname("output.filename"//C_NULL_CHAR, status)
                if (status /= 0) then
                    call c_config_get_string("output.filename"//C_NULL_CHAR, config_field, status)
                    if (status == 0) then
                        ierr = 1        ! signal ERROR
                        return
                    end if
                    filename_out = char_array2str(config_field)
                
                    is_output_set = 1
                end if
    
                if (allocated(config_field)) deallocate(config_field)
#endif
            end if        
        end do
    
        !< set auto-defined output filename 
        if (is_output_set == 0) then 
            filename_out = 'output-' // trim(get_dataset_tag(dataset%id)) // '.txt'
        end if
    
    end subroutine
    

end module sfx_run