gabls1 experiment added

.idea/misc.xml

@@ -3,6 +3,10 @@
   <component name="Black">
     <option name="sdkName" value="Python 3.8 (untitled)" />
   </component>
+  <component name="CMakePythonSetting">
+    <option name="pythonIntegrationState" value="NO" />
+  </component>
+  <component name="CMakeWorkspace" PROJECT_DIR="$PROJECT_DIR$" />
   <component name="JsBuildToolPackageJson" sorting="DEFINITION_ORDER" />
   <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (untitled)" project-jdk-type="Python SDK" />
 </project>
\ No newline at end of file

.idea/vcs.xml

@@ -2,5 +2,6 @@
 <project version="4">
   <component name="VcsDirectoryMappings">
     <mapping directory="" vcs="Git" />
+    <mapping directory="$PROJECT_DIR$/cmake-build-debug/_deps/inmcm60_sfx-src" vcs="Git" />
   </component>
 </project>
\ No newline at end of file

CMakeLists.txt

 cmake_minimum_required(VERSION 3.23)
 option(BUILD_DOC    "Build documentation"    OFF)
-option(USE_CONFIG_PARSER "Build config parser"    ON)
-option(USE_NETCDF "Enable netcdf library" ON)
+option(USE_NETCDF "Enable netcdf library for tests output" OFF)
 option(WITH_RRTMG "with(out) RRTMG radiation transfer" OFF)
+option(WITH_TESTS "Generate and compile test cases" OFF)
 
+#project cmake includes
+include(FetchContent)
+# Project main definitions
 project(scm_abl)
 enable_language(Fortran)
-enable_language(C CXX)
+enable_language(C)
 set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/modules)
 
 if(BUILD_DOC)
     find_package(Doxygen)
@@ -24,4 +28,92 @@ if(BUILD_DOC)
 endif(BUILD_DOC)
 
 
-set(lib_files )
+set(lib_files
\ No newline at end of file
+    src/parkinds.f90
+    src/pbl_grid.f90
+        src/phys_fluid.f90
+    src/state_utils.f90
+    src/scm_state_data.f90
+    src/pbl_turb_data.f90
+    src/diag_pbl.f90)
+
+add_library(abl_lib ${lib_files})
+
+
+target_compile_options(abl_lib
+        PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -cpp>)
+if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug")
+    if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+        target_compile_options(abl_lib
+                PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
+    elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+        target_compile_options(abl_lib
+                PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -O0 -fbacktrace -ffpe-trap=zero,overflow,underflow >)
+    endif()
+endif()
+
+if ("${CMAKE_BUILD_TYPE}" STREQUAL "Release")
+    if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+        target_compile_options(abl_lib
+                PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
+    elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+        target_compile_options(abl_lib
+                PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -fbacktrace -ffpe-trap=zero,overflow,underflow >)
+    endif()
+endif()
+
+if (WITH_TESTS)
+    # RRTMG option
+    if (WITH_RRTMG)
+        add_definitions(-DUSE_RRTMG)
+        add_subdirectory(radiation)
+    endif()
+    # adding sfx and config-parser stuff
+    #fetch sfx model
+    FetchContent_Declare(
+            inmcm60_sfx
+            GIT_REPOSITORY http://tesla.parallel.ru/inmcm-mirror/sfx.git
+            GIT_TAG        origin/main)
+    FetchContent_MakeAvailable(inmcm60_sfx)
+    FetchContent_GetProperties(inmcm60_sfx)
+
+    if(NOT inmcm60_sfx_POPULATED)
+        FetchContent_Populate(sfx_lib)
+        add_subdirectory(${inmcm60_sfx_SOURCE_DIR} ${inmcm60_sfx_BINARY_DIR}  EXCLUDE_FROM_ALL)
+    endif()
+
+    #Config parser from sfx_lib is needed
+    if (WITH_RRTMG)
+        list(APPEND test_files  src/rrtm_interface.f90)
+    endif()
+    #gabls1 experiment
+    add_executable(gabls1 src/tests/gabls1.f90 src/config-utils.f90 src/scm_io_default.f90 src/scm_sfx_data.f90)
+    target_include_directories(gabls1 PUBLIC ${CMAKE_BINARY_DIR}/modules/)
+    target_link_libraries( gabls1 PRIVATE abl_lib )
+    target_link_libraries(gabls1 PRIVATE sfx_lib)
+    target_link_libraries(gabls1 PRIVATE config_parser_F config_parser_CXX)
+    if(WITH_RRTMG)
+        target_link_libraries(gabls1 PRIVATE rrtm)
+    endif()
+    target_compile_options(gabls1
+            PRIVATE $<$<COMPILE_LANGUAGE:Fortran>: -cpp>)
+    if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug")
+        if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+            target_compile_options(gabls1
+                    PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
+        elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+            target_compile_options(gabls1
+                    PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -O0 -fbacktrace -ffpe-trap=zero,overflow,underflow >)
+        endif()
+    endif()
+
+    if ("${CMAKE_BUILD_TYPE}" STREQUAL "Release")
+        if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+            target_compile_options(gabls1
+                    PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -traceback -check all -ftrapuv -debug all >)
+        elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+            target_compile_options(gabls1
+                    PUBLIC $<$<COMPILE_LANGUAGE:Fortran>: -g -fbacktrace -ffpe-trap=zero,overflow,underflow >)
+        endif()
+    endif()
+endif()
\ No newline at end of file

src/config-utils.f90

+! Created by Andrey Debolskiy on 18.10.2024.
+
+module config_utils
+    use config_parser
+    use ISO_C_BINDING, only: C_NULL_CHAR
+    use parkinds, only: rf=>kind_rf, im=>kind_im
+
+    implicit none
+    integer, public, save:: is_config_initialized = 0
+    public init_config, get_fluid_params, get_grid_params
+
+
+    contains
+    subroutine init_config(fname,status, ierr)
+
+        implicit none
+        character(len = *), intent(in) ::fname
+        integer,intent(out):: status, ierr
+        call c_config_run(trim(fname)//C_NULL_CHAR, status)
+        if (status == 0) then
+            write(*, *) ' FAILURE! > unable to parse configuration file: ', trim(fname)
+            return
+        end if
+        is_config_initialized = 1
+    end subroutine init_config
+
+    subroutine get_fluid_params(fluid_params, status, ierr)
+        use phys_fluid, only: fluidParamsDataType
+
+        type(fluidParamsDataType), intent(inout):: fluid_params
+        integer,intent(out):: status, ierr
+        ierr = 0
+
+        if( is_config_initialized /= 0 ) then
+            ! Fluid params
+            call c_config_is_varname("fluid.tref"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                call c_config_get_float("fluid.tref"//C_NULL_CHAR, fluid_params%tref, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+            call c_config_is_varname("fluid.pref"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                !< mandatory in user dataset
+                call c_config_get_float("fluid.pref"//C_NULL_CHAR, fluid_params%pref, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+            call c_config_is_varname("fluid.beta"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                !< mandatory in user dataset
+                call c_config_get_float("fluid.beta"//C_NULL_CHAR, fluid_params%beta, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+            call c_config_is_varname("fluid.g"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                !< mandatory in user dataset
+                call c_config_get_float("fluid.g"//C_NULL_CHAR, fluid_params%g, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+            call c_config_is_varname("fluid.kappa"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                call c_config_get_float("fluid.kappa"//C_NULL_CHAR, fluid_params%kappa, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+
+        end if
+    end subroutine get_fluid_params
+
+    subroutine get_grid_params(grid, status, ierr)
+    use pbl_grid, only: pblgridDataType, &
+                grid_inmcm21_tag, grid_inmcm73_tag, &
+            grid_streached_tag, grid_uniform_tag, &
+            set_pbl_grid_uniform
+
+    implicit none
+
+    type(pblgridDataType), intent(inout):: grid
+    integer,intent(out):: status, ierr
+    character(len=50) :: tag
+    character, allocatable :: config_field(:)
+    real(kind=rf):: h_bot, h_top
+    integer(kind=im):: nz
+
+
+    ierr = 0
+
+    if( is_config_initialized /= 0 ) then
+        ! grid type
+        call c_config_is_varname("grid.type"//C_NULL_CHAR, status)
+        if ((status /= 0)) then
+            call c_config_get_string("grid.type"//C_NULL_CHAR, config_field, status)
+            if (status == 0) then
+                ierr = 1        ! signal ERROR
+                return
+            end if
+        end if
+        if (trim(tag) == trim(grid_uniform_tag)) then
+            call c_config_is_varname("grid.nz"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                call c_config_get_int("grid.nz"//C_NULL_CHAR, nz, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+
+            call c_config_is_varname("grid.h_bot"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                call c_config_get_float("grid.h_bot"//C_NULL_CHAR, h_bot, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+            call c_config_is_varname("grid.h_top"//C_NULL_CHAR, status)
+            if ((status /= 0)) then
+                call c_config_get_float("grid.h_top"//C_NULL_CHAR, h_top, status)
+                if (status == 0) then
+                    ierr = 1        ! signal ERROR
+                    return
+                end if
+            end if
+            call set_pbl_grid_uniform(grid, h_bot, h_top, nz)
+        end if
+    else
+        status  = 0
+        ierr = 2
+    end if
+    end subroutine get_grid_params
+
+    !> @brief character array to string conversion
+    function char_array2str(char_array) result(str)
+        ! ----------------------------------------------------------------------------
+        implicit none
+        character, intent(in) :: char_array(:)
+        character(len=:), allocatable :: str
+        integer :: i
+        ! ----------------------------------------------------------------------------
+
+        str = ""
+        do i = 1, size(char_array)
+            str = str(:) // char_array(i)
+        end do
+
+    end function
+end module config_utils
\ No newline at end of file

src/diag_pbl.f90

+module diag_pbl
+
+contains
+    subroutine get_hpbl_old(hpbla, kpbl, theta_v, z, z_surf, kl, cor, ustar)
+        implicit none
+        real, intent(in):: cor, ustar, z_surf
+        integer, intent(in):: kl
+        real, intent(in), dimension(KL):: theta_v, z
+        real, intent(out):: hpbla
+        integer, intent(out):: kpbl
+
+        !local variabls
+        real :: dz_low, hdyn, dz_hdyn, dz_conv
+        integer:: k, kpbld, kpblc
+
+
+
+        dz_low = z(kl) - z_surf
+        hdyn = AMIN1(dz_low , 0.5E0 * ustar/cor)
+
+        kpblc = kl
+        kpbld = kl
+
+        do k = kl-1,1,-1
+            dz_hdyn =  dz_low - HDYN
+            dz_conv = theta_v(k) - theta_v(kl)
+            if(kpbld.EQ.kl .AND. dz_hdyn.GE.0.E0) kpbld = k
+            if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k
+        enddo
+        kpbl = MIN0(kpblc, kpbld, KL-2)
+        hpbla = z(kpbl) - z_surf
+    end subroutine get_hpbl_old
+
+    subroutine get_hpbl(hpbla, kpbl, theta_v, z, z_surf, kl, cor, ustar)
+        implicit none
+        real, intent(in):: cor, ustar, z_surf
+        integer, intent(in):: kl
+        real, intent(in), dimension(KL):: theta_v, z
+        real, intent(out):: hpbla
+        integer, intent(out):: kpbl
+
+        !local variabls
+        real :: dz_low, hdyn, dz_hdyn, dz_conv
+        integer:: k, kpbld, kpblc
+
+
+
+        dz_low = z(kl) - z_surf
+        hdyn = AMIN1(dz_low , 0.5E0 * ustar/cor)
+
+        kpblc = kl
+        kpbld = kl
+
+        do k = kl-1,1,-1
+            dz_hdyn =  dz_low - HDYN
+            dz_conv = theta_v(k) - theta_v(kl)
+            if(kpbld.EQ.kl .AND. dz_hdyn.GE.0.E0) kpbld = k
+            if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k
+        enddo
+        kpbl = MIN0(kpblc, kpbld, KL-2)
+        hpbla = z(kpbl) - z_surf
+    end subroutine get_hpbl
+
+    subroutine diag_pblh_inmcm(z,thetav,lat,zs,ustar,kl,hpbl)
+    ! This subroutine calculates height of the PBL above ground level according to INMCM algorithm
+    ! input:
+    ! z - array with heights above sea level
+    ! zs - height of surface above sea level
+    ! thetav - array with virtual potential temperatures
+    ! lat - latitude
+    ! ustar - ustar at the lowest model level
+    ! kl - number of vertical levels
+    !!! IMPORTANT:  In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
+    ! output:
+    ! hpbl - height of PBL
+
+    implicit none
+    integer,intent(in)::kl
+    real,intent(in),dimension(kl)::z,thetav
+    real,intent(in)::lat,zs,ustar
+    real,intent(out)::hpbl
+    real,parameter::omega=7.2921/(10**5)
+    real,parameter::cormin=5.0/(10**5)
+    real cor,yzkl,hdyn,ydelz,ydeltv
+    integer KPBL,KPBLC,KPBLD,k
+
+    cor = 2 * omega * sin (lat * 3.14 / 180)
+    cor = max(cormin,abs(cor))
+
+    yzkl = z(kl)-zs
+    hdyn = min(yzkl,0.5 * ustar/cor)
+
+    ydelz = yzkl - hdyn
+
+    if(ydelz.ge.0)then
+    KPBLD = kl-1
+    else
+    KPBLD = kl
+    endif
+
+    KPBLC = kl
+    do k=kl-1,1,-1
+    ydeltv = thetav(k) - thetav(kl)
+    if(KPBLC.eq.kl.and.ydeltv.gt.0) KPBLC = k
+    enddo
+
+    KPBL = min(KPBLC,KPBLD,KL-2)
+
+    call get_hpbl(hpbl, kpbl, thetav, z, zs, kl, cor, ustar)
+    !hpbl = z(KPBL) - zs
+    return
+    end subroutine diag_pblh_inmcm
+
+
+
+    subroutine diag_pblh_rib(theta,z,u,kl,zs,hpbl,nkpbl)
+    ! This subroutine calculates PBL depth according to (Troen and Mahrt 1986) as the lowest level where Rib>Ric
+    ! Ric varies between 0.15 and 0.5
+    ! Rib = (g/theta0)*(theta(z)-theta(s))*z/u(z)**2
+    !input:
+    ! theta - array with (virtual) potential temperature
+    ! z - array with heights above sea level
+    ! u - array with wind speed
+    ! zs - height of surface above sea level
+    ! kl - number of vertical levels
+    ! output:
+    ! hpbl - PBL height above ground level
+    !!! IMPORTANT:  In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
+    implicit none
+    integer,intent(in)::kl
+    real,intent(in),dimension(kl)::z,theta,u
+    real,intent(in)::zs
+    real,intent(out)::hpbl
+    integer,intent(out)::nkpbl
+    real,parameter:: g = 9.81
+    real,parameter:: Ric = 0.25
+    real,parameter:: theta0 = 300.0
+    real,parameter:: upper_bound = 6000 !upper boundary of PBLH
+    real,dimension(kl)::Rib
+    real du
+    integer k,KPBL
+
+    KPBL=kl
+
+    do k=kl-1,1,-1
+    if(z(k).lt.upper_bound)then
+
+    du = u(k)-u(kl)
+    if(du.eq.0) du = 0.1
+
+    Rib(k) = (g / theta0) * (theta(k) - theta(kl)) * (z(k) - z(kl)) / (du**2)
+
+    if(Rib(k).ge.Ric.and.KPBL.eq.kl)then
+     KPBL = k
+     if(KPBL.le.kl-2)then
+     hpbl = z(k) - (z(k)-z(k+1))*(Rib(k)-Ric)/(Rib(k)-Rib(k+1))
+     else
+     hpbl = z(k)
+     endif
+    endif
+
+    endif
+    enddo
+
+    if(KPBL.eq.kl) hpbl=z(kl)
+
+    hpbl = hpbl - zs
+    nkpbl = kl - KPBL + 1
+
+    end subroutine diag_pblh_rib
+
+    subroutine diag_pblh_rh(t,z,q,p,kl,zs,hpbl,nkpbl)
+        ! This subroutine calculates PBL height as height where minimal gradient of relative humidity is found
+        !input:
+        ! t - air temperature (not potential)
+        ! z - heights above sea level
+        ! q - specific humidity
+        ! p - pressure
+        ! u - wind speed
+        ! zs - height of surface above sea level
+        ! kl - number of vertical levels; index of the lowest level
+        ! output:
+        ! hpbl - PBL height above ground level
+        ! nkpbl - number of levels inside PBL
+        !!! IMPORTANT:  In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
+    implicit none
+        integer,intent(in)::kl
+        real,intent(in),dimension(kl)::t,z,q,p
+        real,intent(in)::zs
+        real,intent(out)::hpbl
+        integer,intent(out)::nkpbl
+        real,dimension(kl):: Es,rh
+        real,parameter::E0 = 611.0
+        real,parameter:: upper_bound = 6000 !maximum PBLH
+        integer i,j,k,KPBL
+        real mindrh,drhdz
+
+          Es(:) = E0*10**(7.45*(t(:)-273.15)/(235+t(:)-273.15))
+          rh(:) = 100*q(:)*p(:)/(0.622*Es(:))
+
+        KPBL=kl
+
+        mindrh = 100.0
+        do k=kl-1,1,-1
+            if(z(k).lt.upper_bound)then
+
+                drhdz = rh(k) - rh(k+1)
+
+                if(drhdz.lt.mindrh)then
+                    KPBL = k
+                    mindrh = drhdz
+                endif
+            endif
+        enddo
+
+        hpbl = z(KPBL) - zs
+        nkpbl = kl - KPBL + 1
+    end subroutine diag_pblh_rh
+
+end module diag_pbl
\ No newline at end of file

src/external_forcing.f90

+! Created by Andrey Debolskiy on 15.08.2024.
+
+module external_forcing
+    implicit none
+    !type declarations
+    type forcing_1d
+        real, allocatable :: x(:) ! 0 - closed, 1 - opened
+        real, allocatable :: val(:)
+        integer :: xsize
+    end type forcing_1d
+    type forcing_2d
+        real, allocatable :: x(:) ! 0 - closed, 1 - opened
+        real, allocatable :: y(:) ! 0 - closed, 1 - opened
+        real, allocatable :: val(:,:)
+        integer :: xsize, ysize
+    end type forcing_2d
+
+    public :: ext_forc_dealloc_1d, ext_forc_dealloc_2d
+    public :: get_1d_interpolated, get_1d_to_1d_interpolated, &
+                read_2d_focing
+
+
+    contains
+        real function get_1d_interpolated(xq, forcing1d) result(val)
+            real, intent(in) :: xq
+            type (forcing_1d), intent(in) :: forcing1d
+            !local
+            integer i
+            real x1,x2, val1, val2
+
+            if (xq < forcing1d%x(1)) then
+                val = forcing1d%val(1)
+                return
+            end if
+
+            do i = 1, forcing1d%xsize
+                if (forcing1d%x(i) >= xq) exit
+            end do
+
+            if (i == forcing1d%xsize) then
+                val = forcing1d%val(i)
+                return
+            else
+                x1 = forcing1d%x(i-1)
+                x2 = forcing1d%x(i)
+                val1  = forcing1d%val(i-1)
+                val2  = forcing1d%val(i)
+                val = val1 + (val2 - val1) * (xq - x1) / (x2 -x1)
+                return
+            end if
+        end function get_1d_interpolated
+
+        subroutine get_1d_to_1d_interpolated(valq, xq, nq, forcing1d)
+            integer, intent(in):: nq
+            real, dimension(nq), intent(in) :: xq
+            real, dimension(nq), intent(out) :: valq
+            type (forcing_1d), intent(in) :: forcing1d
+            !local
+            integer k
+            do k = 1, nq
+                valq(k) = get_1d_interpolated(xq(k), forcing1d)
+            enddo
+        end subroutine get_1d_to_1d_interpolated
+
+        subroutine read_2d_focing(out_forc, filename, filename_x, filename_y)
+            implicit none
+            type (forcing_2d), intent(out):: out_forc
+            character(*), intent(in):: filename, filename_x, filename_y
+
+            !local vars
+            integer nx,ny, i, j
+            real x,y
+            open(10, FILE=trim(filename))
+            read(10,*) nx
+            read(10,*) ny
+
+            out_forc%xsize=nx
+            out_forc%ysize=ny
+            allocate(out_forc%val(nx,ny))
+            do i  = 1, nx
+                READ(10,*) (out_forc%val(i,j),j=1,ny)
+            end do
+            close(10)
+            ! read coords
+            allocate(out_forc%x(nx))
+            open(10, FILE=trim(filename_x))
+            do i=1,nx
+                read(10,*) out_forc%x(i)
+            end do
+            close(10)
+
+            allocate(out_forc%y(ny))
+            open(10, FILE=trim(filename_y))
+            do i=1,ny
+                read(10,*) out_forc%y(i)
+            end do
+            close(10)
+        end subroutine read_2d_focing
+
+        subroutine ext_forc_dealloc_1d(forc_1d)
+            implicit none
+            type (forcing_1d), intent(inout) :: forc_1d
+                deallocate(forc_1d%val)
+                deallocate(forc_1d%x)
+        end subroutine ext_forc_dealloc_1d
+
+        subroutine ext_forc_dealloc_2d(forc_2d)
+            implicit none
+            type (forcing_2d), intent(inout) :: forc_2d
+            deallocate(forc_2d%val)
+            deallocate(forc_2d%x)
+            deallocate(forc_2d%y)
+        end subroutine ext_forc_dealloc_2d
+end module external_forcing
\ No newline at end of file

src/parkinds.f90

+! Created by Andrey Debolskiy on 03.07.2024.
+! this is the module to keep track of c/c++ default size type compliance
+module parkinds
+    use, intrinsic :: iso_c_binding
+    implicit none
+    save
+
+    !     integer kinds
+    !     ------------
+    integer, parameter :: kind_ib = c_long  ! 8 byte integer
+    integer, parameter :: kind_im = c_int ! 4 byte integer
+
+    !
+    !     real kinds
+    !     ----------
+    !
+    integer, parameter :: kind_rd = c_double ! 8 byte real
+    integer, parameter :: kind_rf = c_float  ! 4 byte real
+
+    integer, parameter :: kind_bool = c_bool
+
+end module parkinds
\ No newline at end of file

src/pbl_grid.f90

+! Created by Andrey Debolskiy on 03.07.2024.
+
+module pbl_grid
+    use parkinds, only: rf=>kind_rf, im=>kind_im
+    implicit none
+
+    type, public :: pblgridDataType
+        real(kind=rf), allocatable :: z_edge(:)
+        real(kind=rf), allocatable :: z_cell(:)
+        real(kind=rf), allocatable :: dzc(:)
+        real(kind=rf), allocatable :: dze(:)
+        real(kind=rf), allocatable :: sigma(:)
+        real(kind=rf), allocatable :: dsigma(:)
+        integer(kind=im) :: kmax
+    end type pblgridDataType
+    !grid types
+    character(len = 16), parameter :: grid_uniform_tag = 'uniform'
+    character(len = 16), parameter :: grid_streached_tag = 'streached'
+    character(len = 16), parameter :: grid_inmcm21_tag = 'inmcm21'
+    character(len = 16), parameter :: grid_inmcm73_tag = 'inmcm73'
+    ! public interface
+    ! --------------------------------------------------------------------------------
+    public :: allocate_pbl_grid, deallocate_pbl_grid
+    public :: set_pbl_grid_via_edge  !(grid,zpos,nk)
+    public :: set_pbl_grid_uniform
+    !public :: set_from_file()
+
+    contains
+        subroutine allocate_pbl_grid(grid, nk)
+            implicit none
+            integer(kind=im), intent(in) :: nk
+            type(pblgridDataType), intent(inout) :: grid
+            grid%kmax = nk
+            allocate(grid%z_edge(0:nk)) ! to be consistent with other places in INMCM
+            allocate(grid%z_cell(nk))
+            allocate(grid%dzc(nk))
+            allocate(grid%dze(nk))
+            allocate(grid%sigma(nk))
+            allocate(grid%dsigma(nk))
+        end subroutine allocate_pbl_grid
+
+        subroutine deallocate_pbl_grid(grid)
+            implicit none
+            type(pblgridDataType), intent(inout) :: grid
+            deallocate(grid%z_cell)
+            deallocate(grid%z_edge)
+            deallocate(grid%dzc)
+            deallocate(grid%dze)
+            deallocate(grid%sigma)
+            deallocate(grid%dsigma)
+        end subroutine deallocate_pbl_grid
+
+        subroutine set_pbl_grid_via_edge(grid,zpos_cell,zpos_edge,nk)
+            implicit none
+            type(pblgridDataType), intent(inout) :: grid
+            integer(kind=im), intent(in) :: nk
+            real(kind=rf), intent(in) :: zpos_edge(0:nk)
+            real(kind=rf), intent(in) :: zpos_cell(0:nk)
+            !local  variables
+            integer(kind=im):: k
+
+            if(grid%kmax /= nk) then
+                write(*,*) 'Grid size does not match in set_pbl_grid! grid%kmax =',grid%kmax,' nk = ',nk
+                stop ! not sure how to do error handling in fortran
+            end if
+            grid%z_edge(grid%kmax) = zpos_edge(nk)
+            grid%z_cell(grid%kmax) = zpos_cell(nk)
+            do k = nk-1, 1, -1
+                grid%z_edge(k) = zpos_edge(k)
+                grid%z_cell(k) = zpos_cell(k)
+                grid%dzc(k+1) = zpos_edge(k+1) - zpos_edge(k)
+                grid%dze(k+1) = zpos_cell(k+1) - zpos_cell(k)
+            end do
+            grid%z_edge(0) = zpos_edge(0)
+            grid%dzc(1) =zpos_edge(1) - zpos_edge(0)
+            grid%dzc(1) = 2.0 * (zpos_cell(1) - zpos_edge(0))
+        end subroutine set_pbl_grid_via_edge
+
+        subroutine set_pbl_grid_uniform(grid, h_bottom, h_top, nk)
+            implicit none
+            type(pblgridDataType), intent(out):: grid
+            integer(kind=im), intent(in) :: nk
+            real(kind=rf), intent(in) :: h_bottom, h_top
+
+            !local
+            real(kind=rf):: dz, stuff
+            integer(kind=im) :: k
+            call allocate_pbl_grid(grid, nk)
+            dz = (h_top - h_bottom) / real(nk, kind = kind(dz))
+            grid%z_edge(nk) = h_bottom
+            do k = nk-1, 1, -1
+                grid%z_edge(k) = grid%z_edge(k+1) + dz
+                grid%z_cell(k) = 0.5 * (grid%z_edge(k+1) + grid%z_edge(k))
+            end do
+            grid%z_edge(0) = h_top
+
+        end subroutine set_pbl_grid_uniform
+
+
+end module pbl_grid
\ No newline at end of file

src/pbl_turb_data.f90

+! Created by Andrey Debolskiy on 26.11.2024.
+
+module pbl_turb_data
+    type turbBLDataType
+        real, allocatable:: uw(:), vw(:), thetaw(:), thetavw(:), qvw(:)
+        real, allocatable:: s2(:),n2(:)
+        real, allocatable :: qcw(:), cloud_frac_w(:)
+        real, allocatable:: tke(:), eps(:), l_turb(:)
+        real, allocatable:: rig(:), prt(:), rif(:)
+        real, allocatable:: km(:), kh(:)
+        integer :: kmax
+    end type turbBLDataType
+    public turb_data_allocate, turb_data_deallocate
+    contains
+    subroutine turb_data_allocate(turb_data, kmax)
+    implicit none
+        type(turbBLDataType),intent(inout):: turb_data
+        integer:: kmax
+
+        turb_data%kmax = kmax
+        allocate(turb_data%uw(kmax))
+        allocate(turb_data%vw(kmax))
+        allocate(turb_data%thetaw(kmax))
+        allocate(turb_data%thetavw(kmax))
+        allocate(turb_data%qvw(kmax))
+        allocate(turb_data%qcw(kmax))
+        allocate(turb_data%cloud_frac_w(kmax))
+        allocate(turb_data%s2(kmax))
+        allocate(turb_data%n2(kmax))
+        allocate(turb_data%tke(kmax))
+        allocate(turb_data%eps(kmax))
+        allocate(turb_data%l_turb(kmax))
+        allocate(turb_data%rig(kmax))
+        allocate(turb_data%rif(kmax))
+        allocate(turb_data%prt(kmax))
+        allocate(turb_data%km(kmax))
+        allocate(turb_data%kh(kmax))
+    end subroutine turb_data_allocate
+
+    subroutine scm_data_deallocate(turb_data)
+        implicit none
+        type(turbBLDataType), intent(inout):: turb_data
+        deallocate(turb_data%uw)
+        deallocate(turb_data%vw)
+        deallocate(turb_data%thetaw)
+        deallocate(turb_data%thetavw)
+        deallocate(turb_data%qvw)
+        deallocate(turb_data%qcw)
+        deallocate(turb_data%cloud_frac_w)
+        deallocate(turb_data%s2)
+        deallocate(turb_data%n2)
+        deallocate(turb_data%km)
+        deallocate(turb_data%kh)
+        deallocate(turb_data%tke)
+        deallocate(turb_data%eps)
+        deallocate(turb_data%l_turb)
+        deallocate(turb_data%rig)
+        deallocate(turb_data%rif)
+        deallocate(turb_data%prt)
+    end subroutine scm_data_deallocate
+
+end module pbl_turb_data
\ No newline at end of file

src/phys_fluid.f90

+! Created by Andrey Debolskiy on 03.07.2024.
+
+module phys_fluid
+    use parkinds, only:  rf => kind_rf
+implicit none
+    real(kind=rf), parameter :: pref_default = 1000.00 !< reference pressure hPa
+    real(kind=rf), parameter :: tref_default = 287.15 !< reference temperature
+    real(kind=rf), parameter  :: R_d_default = 287.052874  !< gas constant for dry air
+    real(kind=rf), parameter  :: R_v_default = 461.501 !< gas constant for water vapor
+    real(kind=rf), parameter  :: eps_v_default = 1.6077177 !< epilon R_d / R_v
+    real(kind=rf), parameter  :: cp_default = 1.005e3 !< heat capacity of air in isobaric proccess
+    real(kind=rf), parameter  :: g_default = 9.806  !< gravity acceleraion constant
+    real(kind=rf), parameter  :: omega_default = 7.292115e-5  !< gravity acceleraion constant
+    real(kind=rf), parameter  :: kappa_default = 0.4  !< von Karman constant
+    real(kind=rf), parameter  :: appa_default = 0.28607164179  !< dry adiapatic parameter R_d/cp
+
+    type, public :: fluidParamsDataType
+        real(kind=rf) :: pref
+        real(kind=rf) :: p0
+        real(kind=rf) :: tref
+        real(kind=rf) :: R_d
+        real(kind=rf) :: R_v
+        real(kind=rf) :: eps_v
+        real(kind=rf) :: cp
+        real(kind=rf) :: g
+        real(kind=rf) :: omega
+        real(kind=rf) :: kappa
+        real(kind=rf) :: appa
+        real(kind=rf) :: beta
+    end type fluidParamsDataType
+    public :: set_fluid_default
+contains
+    subroutine set_fluid_default(fp)
+        implicit none
+        type(fluidParamsDataType), intent(inout):: fp
+        fp%pref = pref_default
+        fp%p0 = fp%pref
+        fp%tref = tref_default
+        fp%R_d = R_d_default
+        fp%R_v = R_v_default
+        fp%eps_v = eps_v_default
+        fp%cp = cp_default
+        fp%g = g_default
+        fp%omega = omega_default
+        fp%kappa = kappa_default
+        fp%appa = fp%R_d / fp%cp
+        fp%beta = fp%g / fp%tref
+    end subroutine set_fluid_default
+
+end module phys_fluid
\ No newline at end of file

src/rrtm_interface.f90

+MODULE INM_rrtm_interface
+
+  use parkinds, only: rf=>kind_rf, im=>kind_im
+  use rrtmg_lw_rad, only: rrtmg_lw
+  use rrtmg_sw_rad, only: rrtmg_sw
+  use parrrtm, only : nbndlw, ngptlw
+  use parrrsw, only : nbndsw, ngptsw  
+  
+  IMPLICIT NONE
+
+  PRIVATE
+  
+  PUBLIC :: rrtm_init
+
+  type radiationDataType
+        real, allocatable :: uflx(:,:) !Total sky upward flux (W/m2)
+        real, allocatable :: dflx(:,:) !Total sky downward flux (W/m2)
+        real, allocatable :: uflxc(:,:) !Clear sky upward flux (W/m2)
+        real, allocatable :: dflxc(:,:) !Clear sky downward flux (W/m2)
+        real, allocatable :: hr(:,:) !Total sky radiative heating rate (K/d)    
+        real, allocatable :: hrc(:,:) !Clear sky radiative heating rate (K/d)
+  end type radiationDataType
+ 
+  !real,dimension(1,34) :: uflx,dflx,uflxc,dflxc       ! Total sky and Clear sky longwave upward,downward flux (W/m2)
+  !real,dimension(1,34) :: swuflx,swdflx,swuflxc,swdflxc       ! Total sky and Clear sky shortwave upward,downward flux (W/m2)
+  !real,dimension(1,33) :: hr,hrc           ! Total sky and Clear sky longwave radiative heating rate (K/d)
+  !real,dimension(1,33) :: swhr,swhrc           ! Total sky and Clear sky longwave radiative heating rate (K/d)
+	
+
+CONTAINS
+
+SUBROUTINE rrtm_init (t_init,t_lay,p_init,p_lay,q_lay,o3_lay,o2_lay,tsfc,nlay,swuflx,swdflx,swuflxc,swdflxc,swhr,swhrc,uflx,dflx,uflxc,dflxc,hr,hrc)
+
+!INPUT
+
+      type (forcing_1d), intent(in) :: t_init,t_lay,p_init,p_lay,q_init,o3_lay,o2_lay
+      integer, intent(in) :: nlay
+      real(kind=rf), intent(in) :: tsfc ! Surface temperature (K)
+
+!OUTUP
+      real,dimension(ncol,nlay+1), intent(out) :: uflx,dflx,uflxc,dflxc       ! Total sky and Clear sky longwave upward,downward flux (W/m2)
+      real,dimension(ncol,nlay+1), intent(out) :: swuflx,swdflx,swuflxc,swdflxc       ! Total sky and Clear sky shortwave upward,downward flux (W/m2)
+      real,dimension(ncol,nlay), intent(out) :: hr,hrc           ! Total sky and Clear sky longwave radiative heating rate (K/d)
+      real,dimension(ncol,nlay), intent(out) :: swhr,swhrc           ! Total sky and Clear sky longwave radiative heating rate (K/d)
+
+
+!Local
+!----- Input for rrtm-----
+! Note: All volume mixing ratios are in dimensionless units of mole fraction obtained
+! by scaling mass mixing ratio (g/g) with the appropriate molecular weights (g/mol) 
+   
+    !real, allocatable, intent(inout) :: uflx(:,:),dflx(:,:),uflxc(:,:),dflxc(:,:),swuflx(:,:),swdflx(:,:),swuflxc(:,:),swdflxc(:,:)
+    !uflx,dflx,uflxc,dflxc         - Total sky and Clear sky longwave upward,downward flux (W/m2)
+    !swuflx,swdflx,swuflxc,swdflxc - Total sky and Clear sky shortwave upward,downward flux (W/m2)
+    !hr,hrc                        - Total sky and Clear sky longwave radiative heating rate (K/d)
+    !swhr,swhrc                    - Total sky and Clear sky longwave radiative heating rate (K/d)
+
+
+    integer,parameter :: ncol=1             
+	  integer,parameter:: dyofyr=150          ! Day of the year (used to get Earth/Sun)
+	  integer, parameter :: isolvar = 0         ! Flag for solar variability method
+	  real,parameter :: scon =  1360.85 ! Solar constant (W/m2)
+	  integer,parameter :: iaer =10         ! Aerosol option flag
+                                                      !    0: No aerosol
+                                                      !    6: ECMWF method
+                                                      !    10:Input aerosol optical 
+                                                      !       properties
+      integer,parameter :: icld=0        ! Cloud overlap method
+                                                      !    0: Clear only
+                                                      !    1: Random
+                                                      !    2: Maximum/random
+                                                      !    3: Maximum
+	  integer,parameter:: idrv=0           ! Flag for calculation of dFdT, the change
+                                                      !    in upward flux as a function of 
+                                                      !    surface temperature [0=off, 1=on]
+                                                      !    0: Normal forward calculation
+                                                      !    1: Normal forward calculation with
+                                                      !       duflx_dt and duflxc_dt output												  
+													  
+
+      real:: play(ncol,nlay),plev(ncol,nlay+1),tlay(ncol,nlay),tlev(ncol,nlay+1)             ! Layer,Interface pressures and temperatures (hPa, mb)  (K)        
+      real,dimension(ncol,nlay):: h2ovmr,o3vmr,co2vmr,ch4vmr,n2ovmr,o2vmr,cfc11vmr,cfc12vmr,cfc22vmr,ccl4vmr       ! H2O,o3,co2,ch4,n20,o2 volume mixing ratio
+       real:: emis(ncol,nbndlw)          ! Surface emissivity
+      real:: asdir,aldir,asdif,aldif        ! UV/vis,Near-IR surface albedo direct,diffuse rad
+      real :: adjes              ! Flux adjustment for Earth/Sun distance
+      real :: coszen(ncol)          ! Cosine of solar zenith angle
+
+      integer :: inflgsw         ! Flag for cloud optical properties
+      integer :: iceflgsw        ! Flag for ice particle specification
+      integer :: liqflgsw        ! Flag for liquid droplet specification
+      integer :: inflglw         ! Flag for cloud optical properties
+      integer :: iceflglw        ! Flag for ice particle specification
+      integer :: liqflglw        ! Flag for liquid droplet specification
+
+    real,dimension(ngptlw,ncol,nlay) :: cldfmcl_lw,taucmcl_lw,ssacmcl_lw,asmcmcl_lw,fsfcmcl_lw,ciwpmcl_lw,clwpmcl_lw
+	  real,dimension(ncol,nlay)  :: reicmcl_lw,relqmcl_lw
+
+      real :: cldfmcl_sw(ngptsw,ncol,nlay)     ! Cloud fraction
+      real :: taucmcl_sw (ngptsw,ncol,nlay)     ! In-cloud optical depth
+      real :: ssacmcl_sw (ngptsw,ncol,nlay)     ! In-cloud single scattering albedo
+                                                      !    Dimensions: (ngptsw,ncol,nlay)
+      real :: asmcmcl_sw (ngptsw,ncol,nlay)     ! In-cloud asymmetry parameter
+                                                      !    Dimensions: (ngptsw,ncol,nlay)
+      real :: fsfcmcl_sw (ngptsw,ncol,nlay)     ! In-cloud forward scattering fraction
+                                                      !    Dimensions: (ngptsw,ncol,nlay)
+      real :: ciwpmcl_sw (ngptsw,ncol,nlay)     ! In-cloud ice water path (g/m2)
+                                                      !    Dimensions: (ngptsw,ncol,nlay)
+      real :: clwpmcl_sw (ngptsw,ncol,nlay)     ! In-cloud liquid water path (g/m2)
+                                                      !    Dimensions: (ngptsw,ncol,nlay)
+      real :: reicmcl_sw (ncol,nlay)       ! Cloud ice effective radius (microns)
+                                                      !    Dimensions: (ncol,nlay)
+                                                      ! specific definition of reicmcl depends on setting of iceflgsw:
+                                                      ! iceflgsw = 0: (inactive)
+                                                      ! 
+                                                      ! iceflgsw = 1: ice effective radius, r_ec, (Ebert and Curry, 1992),
+                                                      !               r_ec range is limited to 13.0 to 130.0 microns
+                                                      ! iceflgsw = 2: ice effective radius, r_k, (Key, Streamer Ref. Manual, 1996)
+                                                      !               r_k range is limited to 5.0 to 131.0 microns
+                                                      ! iceflgsw = 3: generalized effective size, dge, (Fu, 1996),
+                                                      !               dge range is limited to 5.0 to 140.0 microns
+                                                      !               [dge = 1.0315 * r_ec]													  												  
+      real :: relqmcl_sw (ncol,nlay)       ! Cloud water drop effective radius (microns)
+	  
+      real:: tauaer_lw (ncol,nlay,nbndlw)     ! Aerosol longwave
+      real:: tauaer_sw (ncol,nlay,nbndsw), ssaaer_sw (ncol,nlay,nbndsw),asmaer_sw (ncol,nlay,nbndsw),ecaer_sw (ncol,nlay,nbndsw)     ! Aerosol  
+  
+
+
+! ----- Local -----
+integer :: i,j,k
+
+
+  
+  do i=1,nlay
+    tlay(1,i)=t_lay%val(i)
+	  play(1,i)=p_lay%val(i)
+    h2ovmr(1,i)=q_lay%val(i)
+    o3vmr(1,i)=o3_lay%val(i)
+    o2vmr(1,i)=o2_lay%val(i)
+	enddo
+  do i=1,nlay+1
+    tlev(1,i)=p_init%val(i)
+	  plev(1,i)=p_init%val(i)
+	enddo
+!for longwave
+ !flags
+ inflglw    = 0
+ iceflglw   = 0
+ liqflglw   = 0
+ !variables
+  do i=1,ncol
+     do j=1,nlay
+        do k=1,nbndlw
+          cldfmcl_lw(k,i,j) = 0
+          taucmcl_lw(k,i,j) = 0
+          ssacmcl_lw(k,i,j) = 0
+          asmcmcl_lw(k,i,j) = 0
+          fsfcmcl_lw(k,i,j) = 0
+          ciwpmcl_lw(k,i,j) = 0
+          clwpmcl_lw(k,i,j) = 0
+          tauaer_lw(i,j,k)  = 0
+ 
+        enddo
+        co2vmr(i,j)     = 0
+        ch4vmr(i,j)     = 0
+        n2ovmr(i,j)     = 0
+        cfc11vmr(i,j)   = 0
+        cfc12vmr(i,j)   = 0
+        cfc22vmr(i,j)   = 0
+        ccl4vmr(i,j)    = 0
+        reicmcl_lw(i,j) = 0
+        relqmcl_lw(i,j) = 0
+      enddo
+  enddo
+	call rrtmg_lw &
+            (ncol    ,nlay    ,icld    ,idrv    , &
+             play    ,plev    ,tlay    ,tlev    ,tsfc    , & 
+             h2ovmr  ,o3vmr   ,co2vmr  ,ch4vmr  ,n2ovmr  ,o2vmr , &
+             cfc11vmr,cfc12vmr,cfc22vmr,ccl4vmr ,emis    , &
+             inflglw ,iceflglw,liqflglw,cldfmcl_lw , &
+             taucmcl_lw ,ciwpmcl_lw ,clwpmcl_lw ,reicmcl_lw ,relqmcl_lw , &
+             tauaer_lw  , &
+             uflx    ,dflx    ,hr      ,uflxc   ,dflxc,  hrc)
+
+  !for shortwave
+  !flags
+  inflgsw    = 0
+  iceflgsw   = 0
+  liqflgsw   = 0
+  !variables
+  asdir = 0.2
+  asdif = 0.2
+  aldir = 0.2
+  aldif = 0.2
+  coszen = 0.6
+  adjes = 1
+  !arrays
+    do i=1,ncol
+     do j=1,nlay
+        do k=1,nbndsw
+          cldfmcl_sw(k,i,j) = 0._rb
+          taucmcl_sw(k,i,j) = 0._rb
+          ssacmcl_sw(k,i,j) = 0._rb
+          asmcmcl_sw(k,i,j) = 0._rb
+          fsfcmcl_sw(k,i,j) = 0._rb
+          ciwpmcl(k,i,j)    = 0._rb
+          clwpmcl(k,i,j)    = 0._rb
+          tauaer_sw(i,j,k)     = 0._rb
+          ssaaer_sw(i,j,k)     = 0._rb
+          asmaer_sw(i,j,k)     = 0._rb
+          ecaer_sw(i,j,k)     = 0._rb
+        enddo
+        co2vmr(i,j)     = 0._rb
+        ch4vmr(i,j)     = 0._rb
+        n2ovmr(i,j)     = 0._rb
+        cfc11vmr(i,j)   = 0._rb
+        cfc12vmr(i,j)   = 0._rb
+        cfc22vmr(i,j)   = 0._rb
+        ccl4vmr(i,j)    = 0._rb
+        reicmcl(i,j) = 0._rb
+        relqmcl(i,j) = 0._rb
+      enddo
+  enddo
+	call rrtmg_sw &
+            (ncol    ,nlay    ,icld    ,iaer    , &
+             play    ,plev    ,tlay    ,tlev    ,tsfc   , &
+             h2ovmr , o3vmr   ,co2vmr  ,ch4vmr  ,n2ovmr ,o2vmr , &
+             asdir   ,asdif   ,aldir   ,aldif   , &
+             coszen  ,adjes   ,dyofyr  ,scon    ,isolvar, &
+             inflgsw ,iceflgsw,liqflgsw,cldfmcl_sw , &
+             taucmcl_sw  ,ssacmcl_sw  ,asmcmcl_sw  ,fsfcmcl_sw  , &
+             ciwpmcl ,clwpmcl ,reicmcl ,relqmcl , &
+             tauaer_sw  ,ssaaer_sw  ,asmaer_sw  ,ecaer_sw   , &
+             swuflx  ,swdflx  ,swhr    ,swuflxc ,swdflxc ,swhrc)	 
+
+   
+  END SUBROUTINE rrtm_init
+  
+  END MODULE INM_rrtm_interface
+  
+  
+  
\ No newline at end of file

src/scm_io_default.f90

+
+module scm_io_default
+    implicit none
+    !type declaration
+    type io_struct
+        character(len = 160) :: fname
+        integer :: status ! 0 - closed, 1 - opened
+        integer :: unit_id
+    end type io_struct
+
+    character(len = 160) tmp_str
+    integer, public, parameter :: nunits_max = 70
+    public
+
+    contains
+        subroutine to_file_1d_2var(fname, var1, var2, n)
+            implicit none
+            character(*), intent(in):: fname
+            integer, intent(in):: n
+            real, dimension(n), intent(in):: var1, var2
+            integer i
+            integer istat
+            character(len = 7) sta
+
+            open(10, FILE=trim(fname))
+            do i =1, n
+                write(10,*) var1(i), var2(i)
+            end do
+            close(10,iostat = istat,STATUS = sta)
+        end subroutine to_file_1d_2var
+
+        subroutine to_file_1d_3var(fname, var1, var2, var3, n)
+            implicit none
+            character(*), intent(in):: fname
+            integer, intent(in):: n
+            real, dimension(n), intent(in):: var1, var2, var3
+            integer i
+            integer istat
+            character(len = 7) sta
+
+            open(10, FILE=trim(fname))
+            do i =1, n
+                write(10,*) var1(i), var2(i), var3(i)
+            end do
+            close(10,iostat = istat,STATUS = sta)
+        end subroutine to_file_1d_3var
+
+        subroutine set_file( f, fname )
+            implicit none
+            type(io_struct), intent(inout):: f
+            character(*), intent(in):: fname
+
+                f%unit_id = get_file_unit()
+                open(f%unit_id, FILE=trim(fname))
+                f%status = 1
+                f%fname = trim(fname)
+
+            write(*,*) 'file opened ', f%fname
+            write(*,*) ' unit: ',f%unit_id
+            write(*,*) 'max_units: ', nunits_max
+        end subroutine set_file
+
+        subroutine write_series(stamp, nlength, f)
+            implicit none
+            type(io_struct), intent(inout):: f
+            integer, intent(in):: nlength
+            real, intent(in), dimension(nlength)::stamp
+            write(f%unit_id,*) stamp(:)
+        end subroutine write_series
+
+        subroutine write_timescan(stamp,nz, nlength, f)
+            implicit none
+            type(io_struct), intent(in):: f
+            integer, intent(in):: nlength, nz
+            real, intent(in), dimension(nlength, nz)::stamp
+            integer k
+            do k=1,nz
+            write(f%unit_id,*) stamp(:, k)
+            end do
+        end subroutine  write_timescan
+
+        subroutine close_file(f)
+            implicit none
+            type(io_struct), intent(inout):: f
+
+            close(f%unit_id)
+            write(*,*) 'file closed ', f%fname
+            f%status = 0
+        end subroutine close_file
+
+        !   get_file_unit returns a unit number that is not in use
+        integer function get_file_unit ()
+            integer lu, iostat
+            integer, save:: m
+            logical, save:: initialized = .true.
+            logical   opened
+
+            if (initialized) then
+                m = nunits_max
+                initialized = .false.
+            end if
+
+            if (m < 8 ) then
+                m = 2 * nunits_max
+            end if
+
+            do lu = m,7,-1
+                inquire (unit=lu, opened=opened, iostat=iostat)
+                if (iostat.ne.0) cycle
+                if (.not.opened) exit
+            end do
+            !
+            get_file_unit = lu
+            return
+        end function get_file_unit
+
+    subroutine read_1d_plain(x, val, nrows, fname)
+        implicit none
+        real, allocatable, intent(inout) :: x(:), val(:)
+        integer, intent(inout):: nrows
+        character(*), intent(in):: fname
+
+        integer ::  i, io
+
+        ! get number of rows
+        nrows = 0
+        open (1, file = trim(fname))
+        do
+            read(1,*,iostat=io)
+            if (io/=0) exit
+            nrows = nrows + 1
+        end do
+        close (1)
+        !check if arrays are already allocated
+        if (allocated(x)) deallocate(x)
+        if (allocated(val)) deallocate(val)
+        !allocate arrays
+        allocate(x(nrows))
+        allocate(val(nrows))
+        ! reopen file and read data
+        open (1, file = trim(fname))
+        do i = 1,nrows
+            read(1,*) x(i), val(i)
+        end do
+    end subroutine read_1d_plain
+end module scm_io_default

src/scm_sfx_data.f90

+! Created by Andrey Debolskiy on 20.11.2024.
+
+module scm_sfx_data
+
+    implicit none
+
+    real taux, tauy
+    real cu, ct
+    real hflx, qflx
+    real ust, vst, umst
+    real u10, v10
+    real ROLIM, WAREA, WWET, ATGR
+
+
+end module scm_sfx_data
\ No newline at end of file

src/scm_state_data.f90

+! Created by Andrey Debolskiy on 20.11.2024.
+
+module scm_state_data
+    implicit none
+    type stateBLDataType
+        real, allocatable:: u(:), v(:), temp(:), theta(:) ,qv(:)
+        real, allocatable:: rho(:)
+        real, allocatable :: lwp(:), cloud_frac(:), s_e(:)
+        real, allocatable:: km(:), kh(:)
+        real, allocatable:: vdctq(:)
+        real, allocatable:: vdcuv(:)
+        real :: HPBLA_diag
+        real :: p0
+        integer :: ktvdm=10
+        integer :: ktvd
+        integer :: kpbl
+        integer :: ktop
+        integer :: kmax
+    end type stateBLDataType
+    public scm_data_allocate
+    public scm_data_deallocate
+    public scm_data_copy
+    contains
+        subroutine scm_data_allocate(bl_data, kmax)
+            implicit none
+            type(stateBLDataType),intent(inout):: bl_data
+            integer:: kmax
+            bl_data%kmax = kmax
+            allocate(bl_data%u(kmax))
+            allocate(bl_data%v(kmax))
+            allocate(bl_data%temp(kmax))
+            allocate(bl_data%theta(kmax))
+            allocate(bl_data%qv(kmax))
+            allocate(bl_data%lwp(kmax))
+            allocate(bl_data%cloud_frac(kmax))
+            allocate(bl_data%s_e(kmax))
+            allocate(bl_data%km(kmax))
+            allocate(bl_data%kh(kmax))
+            allocate(bl_data%vdcuv(kmax))
+            allocate(bl_data%vdctq(kmax))
+        end subroutine scm_data_allocate
+
+        subroutine scm_data_deallocate(bl_data)
+            implicit none
+            type(stateBLDataType), intent(inout):: bl_data
+            deallocate(bl_data%u)
+            deallocate(bl_data%v)
+            deallocate(bl_data%temp)
+            deallocate(bl_data%theta)
+            deallocate(bl_data%qv)
+            deallocate(bl_data%lwp)
+            deallocate(bl_data%cloud_frac)
+            deallocate(bl_data%s_e)
+            deallocate(bl_data%km)
+            deallocate(bl_data%kh)
+            deallocate(bl_data%vdcuv)
+            deallocate(bl_data%vdctq)
+        end subroutine scm_data_deallocate
+
+        subroutine scm_data_copy(bl, bl_old)
+            implicit none
+            type(stateBLDataType),intent(in):: bl_old
+            type(stateBLDataType),intent(out):: bl
+            if (bl%kmax /= bl_old%kmax) then
+                write(*,*) 'error in copy BLData: size is not compatible'
+            else
+                bl%u  = bl_old%u
+                bl%v = bl_old%v
+                bl%temp = bl_old%temp
+                bl%theta = bl_old%theta
+                bl%qv = bl_old%qv
+                bl%lwp = bl_old%lwp
+                bl%cloud_frac= bl_old%cloud_frac
+                bl%s_e = bl_old%s_e
+                bl%km = bl_old%km
+                bl%kh = bl_old%kh
+                bl%vdcuv = bl_old%vdcuv
+                bl%vdctq = bl_old%vdctq
+                bl%hpbla_diag = bl_old%hpbla_diag
+                bl%p0 = bl_old%p0
+                bl%ktvdm = bl_old%ktvdm
+                bl%kpbl = bl_old%kpbl
+                bl%ktop = bl_old%ktop
+            end if
+
+        end subroutine scm_data_copy
+end module scm_state_data
\ No newline at end of file

src/state_utils.f90

+! Created by Andrey Debolskiy on 17.12.2023.
+
+module state_utils
+
+
+    contains
+
+        SUBROUTINE GEO(T,FS,F,SIG, R, LL)
+        implicit none
+        integer, intent(in):: LL
+        real, intent(in):: T(LL)
+        real, intent(in):: SIG(LL)
+        real, intent(in):: FS
+        real, intent(out):: F(LL)
+        real, intent(in):: R
+
+        integer k, ka
+
+        F(LL)=FS-ALOG(SIG(LL))*R*T(LL)
+
+        do  K=LL-1,1,-1
+
+            KA=K+1
+            F(K)=F(KA)+0.5E0*R*(ALOG(SIG(KA))-ALOG(SIG(K)))* &
+             (T(KA)+T(K))
+
+        end do
+
+        END SUBROUTINE GEO
+
+
+
+
+
+            SUBROUTINE NANOLOVKA(X,INAN)
+            implicit none
+            real, intent(in):: x
+            integer, intent(out):: inan
+                inan=1
+
+            IF(X.LE.0.0) THEN
+                INAN=0
+            END IF
+
+            IF(X.GE.0.0) THEN
+                INAN=0
+            END IF
+
+            END SUBROUTINE NANOLOVKA
+
+        subroutine theta2ta(ta, theta, p0, sig, appa, kl)
+            implicit none
+            real, dimension(kl), intent(in):: theta
+            real, dimension(kl), intent(in)::sig
+            real, intent(in):: p0, appa
+            integer, intent(in):: kl
+            real, dimension(kl), intent(out):: ta
+
+            integer k
+
+            do k=1,kl
+                ta(k) = theta(k) * (p0 * sig(k))**appa
+            end do
+        end subroutine theta2ta
+
+        subroutine ta2theta(theta, ta, p0, sig, appa, kl)
+            implicit none
+            real, dimension(kl), intent(in):: ta
+            real, dimension(kl), intent(in)::sig
+            real, intent(in):: p0, appa
+            integer, intent(in):: kl
+            real, dimension(kl), intent(out):: theta
+
+            integer k
+
+            do k=1,kl
+                theta(k) = ta(k) * (p0 * sig(k))**(-appa)
+            end do
+        end subroutine ta2theta
+end module state_utils
\ No newline at end of file

src/tests/gabls1.f90

+! Created by Andrey Debolskiy on 12.10.2024.
+program gabls1
+    use ISO_C_BINDING, only: C_NULL_CHAR
+    use parkinds, only: rf=>kind_rf, im=>kind_im
+    use phys_fluid
+    use scm_io_default
+#ifdef USE_NETCDF
+    use io, only: variable_metadata, write_2d_real_nc, write_1d_real_nc
+#endif
+    use scm_sfx_data, only: taux, tauy, umst, hflx, qflx, cu
+    use scm_state_data
+    use pbl_turb_data
+    use state_utils, only : geo, theta2ta, ta2theta
+    use diag_pbl
+    use pbl_grid
+    use sfx_data, only: meteoDataType, sfxDataType
+    use sfx_most, only: get_surface_fluxes_most => get_surface_fluxes, &
+            numericsType_most => numericsType
+    use config_parser
+    use config_utils
+
+
+    implicit none
+    !local varaibles
+    real time_begin, time_end, time_current
+    real dt, out_dt, output_dt
+    integer status, ierr
+    logical bstatus
+    character(len = 160) ::fname
+    character(len = 128) :: fname_nc
+    character(len = 128) ::fname_config = '../config-examples/config-gabls1-ex.txt'
+    real, allocatable :: ttemp(:,:), utemp(:,:),vtemp(:,:), ttime(:)
+    real(kind=rf):: ug, vg
+    real(kind=rf):: z0, zh, f_cor
+    real(kind=rf):: tsurf
+    real(kind=rf):: shir
+    type(fluidParamsDataType) fluid_params
+    type(pblgridDataType) grid
+
+    type(stateBLDataType):: bl, bl_old;
+
+    type(meteoDataType) :: meteo_cell
+    type(sfxDataType) :: sfx_cell
+    type(numericsType_most) :: numerics_sfx
+
+
+    !diagnostic variables
+    real hpbl
+    integer nkpbl
+
+    !io variables
+    real, dimension (5):: series_val
+    type (io_struct) :: series_f, scan_f
+#ifdef USE_NETCDF
+    type(variable_metadata) :: meta_z, meta_t
+    type(variable_metadata) :: meta_temperature, meta_uwind, meta_vwind
+#endif
+    integer k, nt, kmax
+    ! init experiment params
+    numerics_sfx%maxiters_charnock = 10
+
+    ug = 8.00000
+    vg = 0.0
+    dt = 1.0;
+!    call set_fluid_default(fluid_params)
+!    fluid_params%tref = 265.0
+!    fluid_params%pref = 1013.2500
+!    fluid_params%beta = fluid_params%g / fluid_params%tref
+!    fluid_params%kappa = 0.4
+!    fluid_params%p0= fluid_params%pref
+
+
+
+    !set output filenames
+    fname = 'test.dat'
+    fname_nc = 'gabls1.nc'
+    series_f%fname = 'phys.dat'
+    call set_file(series_f, series_f%fname)
+
+    scan_f%fname='time_scan.dat'
+    call set_file(scan_f, scan_f%fname)
+#ifdef USE_NETCDF
+    call set_meta_nc(meta_z, meta_t, meta_temperature, meta_uwind, meta_vwind)
+#endif
+
+
+
+    call init_config(fname_config,status, ierr)
+    if (status == 0) then
+        ierr = 1        ! signal ERROR
+        return
+    end if
+
+    call c_config_is_varname("time.begin"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("time.begin"//C_NULL_CHAR, time_begin, status)
+        if (status == 0) then
+            ierr = 1        ! signal ERROR
+            return
+        end if
+    end if
+    call c_config_is_varname("time.end"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("time.end"//C_NULL_CHAR, time_end, status)
+        if (status == 0) then
+            ierr = 1        ! signal ERROR
+            return
+        end if
+    end if
+    call c_config_is_varname("time.dt"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("time.dt"//C_NULL_CHAR, dt, status)
+        if (status == 0) then
+            ierr = 1        ! signal ERROR
+            return
+        end if
+    end if
+
+    call c_config_is_varname("time.out_dt"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("time.out_dt"//C_NULL_CHAR, output_dt, status)
+        if (status == 0) then
+            ierr = 1        ! signal ERROR
+            return
+        end if
+    end if
+
+    call c_config_is_varname("Ug"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("Ug"//C_NULL_CHAR, ug, status)
+        if (status == 0) then
+            ierr = 1        ! signal ERROR
+            return
+        end if
+    end if
+        call c_config_is_varname("Vg"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_float("Vg"//C_NULL_CHAR, vg, status)
+        if (status == 0) then
+            ierr = 1        ! signal ERROR
+            return
+        end if
+    end if
+    call set_fluid_default(fluid_params)
+    call get_fluid_params(fluid_params, status, ierr)
+    if (status == 0) then
+        ierr = 1        ! signal ERROR
+        return
+    end if
+    !z coordinate
+    call get_grid_params(grid, status, ierr)
+    if (status == 0) then
+        ierr = 1        ! signal ERROR
+        return
+    end if
+    kmax = grid%kmax
+    time_current = time_begin
+    ! hack to insure shir array is good
+    shir = 3.141592653589793238 * 72.3798 / 180.0
+    f_cor = 2.0 * fluid_params%omega * sin(shir)
+    write(*,*) 'F coriolis: =  ', f_cor
+    !fs(1, 1) = 0.0
+    !setup surface
+    z0 = 0.01
+    zh = z0 / 10.0
+    ! init blData states
+    call scm_data_allocate(bl, kmax)
+    call scm_data_allocate(bl_old, kmax)
+    call init_state(bl, ug, vg, fluid_params%tref)
+    call scm_data_copy(bl_old,bl)
+
+
+
+    !TODO why zet is not initialized in beirit properly?????
+    !CALL GEO(, FS(1, 1), grid%z_cell, grid%kmax)  !change for clima
+    DO K = 1, kmax
+        grid%z_cell(k) = grid%z_cell(k) / fluid_params%g
+
+        if (grid%z_cell(k) > 100.0) then
+            bl_old%theta(k) = fluid_params%tref + 0.01* (grid%z_cell(k) - 100.0)
+        else
+            bl_old%theta(k) = fluid_params%tref
+        end if
+        !write(*,*) k, sig(k)* fluid_params%pref* 100.0, grid%z_cell(k)
+    END DO
+    call theta2ta(bl_old%temp, bl_old%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
+
+    !io setup
+    nt = floor((time_end - time_begin)/output_dt)
+    allocate(ttime(nt))
+    allocate(ttemp(kmax,nt),utemp(kmax,nt),vtemp(kmax,nt))
+#ifdef USE_NETCDF
+    call write_1d_real_nc(grid%z_cell,fname_nc,meta_z)
+#endif
+    out_dt = 0.0
+    nt=0
+    write(*,*)'nt=0', nt
+    do while (time_current < time_end)
+        tsurf = fluid_params%tref - 0.25 * (time_current - time_begin) / 3600.0;
+        meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%u(kmax)**2)
+        meteo_cell%dT = -tsurf + bl_old%theta(kmax)
+        meteo_cell%Tsemi = 0.5 * (tsurf + bl_old%theta(kmax))
+        meteo_cell%dQ  = 0.0
+        meteo_cell%h = grid%z_cell(kmax)
+        meteo_cell%z0_m = z0
+
+        call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
+        cu = 1000.0 * sfx_cell%Cm**2
+        taux = 1.1 * sfx_cell%Cm**2 * meteo_cell%U *  bl_old%u(kmax)
+        tauy = 1.1 * sfx_cell%Cm**2 * meteo_cell%U *  bl_old%v(kmax)
+        hflx =  sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
+        umst = sfx_cell%Cm * meteo_cell%U
+
+        !call vdiff_scm()
+        call add_coriolis(bl, ug, vg, f_cor, dt, grid)
+
+        time_current = time_current + dt
+        if (time_current >=out_dt) then
+            nt = nt+1
+            ttime(nt) = time_current/3600.0
+            out_dt = out_dt + output_dt
+            write(*,*)'nt= ', nt
+
+
+            call ta2theta( bl_old%theta, bl_old%temp, &
+                    fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
+            call diag_pblh_inmcm(grid%z_cell,bl_old%theta,shir,0.0,umst,bl_old%kmax,hpbl)
+            !create output
+            series_val(1) = time_current/3600.0
+            series_val(2) = hflx
+            series_val(3) = umst
+            series_val(4) = hpbl
+            series_val(5) = 0.0
+
+            call write_series(series_val, 5, series_f)
+
+            call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
+#ifdef USE_NETCDF
+            do k = 1, kl
+                ttemp(k,nt) = thold(k)
+                utemp(k,nt) = uold(k)
+                vtemp(k,nt) = vold(k)
+            end do
+#endif
+            !write(*,*) ccld * 10000.0
+        end if
+    end do
+#ifdef USE_NETCDF
+    write(*,*)'nt= ', nt
+    call write_1d_real_nc(ttime,fname_nc,meta_t)
+    write(*,*)'here'
+    call  write_2d_real_nc(ttemp,fname_nc,meta_temperature)
+    call  write_2d_real_nc(utemp,fname_nc,meta_uwind)
+    call  write_2d_real_nc(vtemp,fname_nc,meta_vwind)
+#endif
+
+    call ta2theta(bl_old%theta, bl_old%temp, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
+    call to_file_1d_3var(fname, grid%z_cell, bl_old%u, bl_old%v, bl_old%kmax)
+    call to_file_1d_3var('coeffs.dat',grid%z_cell, bl_old%km, bl_old%kh, bl_old%kmax)
+    call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%theta, bl_old%kmax)
+    call close_file(series_f)
+    call close_file(scan_f)
+
+    call deallocate_pbl_grid(grid)
+    !deallocate(ttemp)
+    !deallocate(ttime)
+    !deallocate(utemp)
+    !deallocate(vtemp)
+end program gabls1
+
+
+subroutine init_state(bl, ug_,vg_,tref)
+    use parkinds, only: rf=>kind_rf, im=>kind_im
+    use scm_state_data, only:stateBLDataType
+
+    implicit none
+    real(kind=rf):: ug_,vg_,tref
+    type(stateBLDataType), intent(out) :: bl
+    bl%u(1:bl%kmax) = ug_
+    bl%v(1:bl%kmax) = vg_
+    bl%temp(1:bl%kmax) = tref
+    bl%theta(1:bl%kmax) = tref
+    bl%qv(1:bl%kmax) = 0.0
+    bl%lwp(1:bl%kmax) = 0.0
+    bl%cloud_frac(1:bl%kmax) = 0.0
+    bl%s_e(1:bl%kmax) = 0.0
+    bl%km(1:bl%kmax) = 0.0
+    bl%kh(1:bl%kmax) = 0.0
+    bl%vdcuv(1:bl%kmax) = 0.0
+    bl%vdctq(1:bl%kmax) = 0.0
+end subroutine init_state
+
+
+!subroutine update_sfx_fluxes(tsurf, z0, zh, z, tair, uair, vair, beta, kappa)
+!    use phys
+!    use sfx_scm, only : taux, tauy, hflx, umst, cu
+!    implicit none
+!    real, intent(in) :: tsurf, z0, zh, z, tair, uair, vair, beta, kappa
+!    real umod, ustar, tstar
+!    real zeta, rib, cm
+!    real psim,psih,z_d,b
+!    integer i
+!    cm= 5.0
+!    z_d = z /z0
+!    b = z0/zh
+!    umod = sqrt(uair**2 + vair**2)
+!    rib = (9.8 * 0.5  / (tair + tsurf))* (tair - tsurf) * z / (Umod * Umod)
+!    ustar = kappa * umod / log(z / z0)
+!    umst =ustar
+!    tstar = kappa * (tair - tsurf) / log(z / zh)
+!    zeta = 0.0
+!    if (abs(rib) < 0.001) then
+!        zeta = 0.0
+!        psim =  log(z_d)
+!        psih =  log(z/zh)
+!    else
+!        do i = 1,1
+!            psih =  log(z_d) +  cm * zeta * (z_d * b - 1.0)  / (z_d * b)
+!            psim =  log(z_d) +  cm * zeta * (z_d - 1.0)  / z_d
+!            zeta = (Rib * psim * psim) / (kappa* psih)
+!        end do
+!    end if
+!    ustar = kappa * umod  / psim
+!    tstar = kappa * (tair - tsurf) / psih
+!    taux = (ustar / umod) * ustar * uair
+!    tauy = (ustar / umod) * ustar * vair
+!    cu =  ustar **2 / umod
+!    hflx = 1.0 * tstar * ustar
+!    hsn(1, 1) = hflx
+!    hlt(1, 1) = 0.0!
+
+!end subroutine update_sfx_fluxes
+
+subroutine add_coriolis(bl, ug, vg, f, dt, grid)
+    use scm_state_data
+    use pbl_grid, only: pblgridDataType
+    implicit none
+    real, intent(in) :: ug, vg, f, dt
+    type(stateBLDataType), intent(out):: bl
+    type(pblgridDataType), intent(in):: grid
+
+    integer k
+        do k = 1, grid%kmax
+        bl%v(k) = bl%v(k)  - dt * f * (bl%u(k) - ug)
+        bl%u(k) = bl%u(k) + dt * f * ( bl%v(k) - vg)
+    end do
+
+end subroutine add_coriolis
+
+subroutine put_tscan( time, z, bl, nl, f)
+    use parkinds, only : rf=>kind_rf, im => kind_im
+    use scm_io_default
+    use scm_state_data, only:stateBLDataType
+    implicit none
+    type(stateBLDataType), intent(in):: bl
+    real(kind=rf), intent(in), dimension(nl):: z
+    real(kind=rf),intent(in) :: time
+    type (io_struct),intent(in) :: f
+    integer(kind=im), intent(in) :: nl
+
+
+    !local
+    real(kind=rf), dimension(5,nl)::stamp
+    integer k
+
+    ! copy to stamp
+    do k=1,nl
+        stamp(1,k) = time
+        stamp(2,k) = z(k)
+        stamp(3,k) = bl%theta(k)
+        stamp(4,k) = bl%u(k)
+        stamp(5,k) = bl%v(k)
+    end do
+
+    ! call to write timestamp
+    call write_timescan(stamp,nl, 5, f)
+
+end subroutine put_tscan
+
+subroutine get_surface_from_config(model, surface, z0,zh )
+    use iso_c_binding, only: C_NULL_CHAR
+    use config_parser
+    use sfx_common, only: char_array2str
+    use sfx_config
+    use sfx_surface, only: get_surface_id
+    use parkinds, only: rf=>kind_rf
+    implicit none
+    character, allocatable :: config_field(:)
+
+    integer,intent(out) :: model, surface
+    real(kind=rf), intent(out) :: z0,zh
+
+    integer status, ierr
+
+
+    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("surface.type"//C_NULL_CHAR, status)
+    if ((status /= 0)) then
+        !< mandatory in user dataset
+        call c_config_get_string("surface.type"//C_NULL_CHAR, config_field, status)
+        if (status == 0) then
+            ierr = 1        ! signal ERROR
+            return
+        end if
+        surface = get_surface_id(char_array2str(config_field))
+        if (surface == -1) then
+            write(*, *) ' FAILURE! > unknown surface [key]: ', trim(char_array2str(config_field))
+            ierr = 1        ! signal ERROR
+            return
+        end if
+    endif
+
+
+end subroutine get_surface_from_config
+#ifdef USE_NETCDF
+subroutine set_meta_nc(z_meta, t_meta, theta_meta, uwind_meta, vwind_meta)
+    use io, only: variable_metadata
+    type(variable_metadata), intent(inout):: z_meta, t_meta, theta_meta, uwind_meta, vwind_meta
+    z_meta = variable_metadata( &
+            name = 'Z', &
+            dim_names = [character(len=32) :: 'Z', '', '', ''], &
+            char_name = [character(len=32) :: &
+                    'long_name', &
+                    'units', &
+                    '', '', '', '', '', '', '', ''], &
+            char_value = [character(len=32) :: &
+                    'Height', &
+                    'm', &
+                    '', '', '', '', '', '', '', ''], &
+            real_name = [character(len=32) :: &
+                    'missing_value', &
+                    'scale_factor', &
+                    '', '', '', '', '', '', '', ''], &
+            real_value = [&
+                    -9999.9, &
+                            1.0, &
+                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+    t_meta = variable_metadata( &
+            name = 'Time', &
+            dim_names = [character(len=32) :: 'Time', '', '', ''], &
+            char_name = [character(len=32) :: &
+                    'long_name', &
+                    'units', &
+                    '', '', '', '', '', '', '', ''], &
+            char_value = [character(len=32) :: &
+                    'Time', &
+                    'hours since 00-00-00', &
+                    '', '', '', '', '', '', '', ''], &
+            real_name = [character(len=32) :: &
+                    'missing_value', &
+                    'scale_factor', &
+                    '', '', '', '', '', '', '', ''], &
+            real_value = [&
+                    -9999.9, &
+                            1.0, &
+                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+    theta_meta = variable_metadata( &
+            name = 'th', &
+            dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
+            char_name = [character(len=32) :: &
+                    'long_name', &
+                    'units', &
+                    '', '', '', '', '', '', '', ''], &
+            char_value = [character(len=32) :: &
+                    'Potential Temperature', &
+                    'K', &
+                    '', '', '', '', '', '', '', ''], &
+            real_name = [character(len=32) :: &
+                    'missing_value', &
+                    'scale_factor', &
+                    '', '', '', '', '', '', '', ''], &
+            real_value = [&
+                    -9999.9, &
+                            1.0, &
+                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+    uwind_meta = variable_metadata( &
+            name = 'ua', &
+            dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
+            char_name = [character(len=32) :: &
+                    'long_name', &
+                    'units', &
+                    '', '', '', '', '', '', '', ''], &
+            char_value = [character(len=32) :: &
+                    'Longitude wind component', &
+                    'm/s', &
+                    '', '', '', '', '', '', '', ''], &
+            real_name = [character(len=32) :: &
+                    'missing_value', &
+                    'scale_factor', &
+                    '', '', '', '', '', '', '', ''], &
+            real_value = [&
+                    -9999.9, &
+                            1.0, &
+                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+    vwind_meta = variable_metadata( &
+            name = 'va', &
+            dim_names = [character(len=32) :: 'Z', 'Time', '', ''], &
+            char_name = [character(len=32) :: &
+                    'long_name', &
+                    'units', &
+                    '', '', '', '', '', '', '', ''], &
+            char_value = [character(len=32) :: &
+                    'Latitude wind component', &
+                    'm/s', &
+                    '', '', '', '', '', '', '', ''], &
+            real_name = [character(len=32) :: &
+                    'missing_value', &
+                    'scale_factor', &
+                    '', '', '', '', '', '', '', ''], &
+            real_value = [&
+                    -9999.9, &
+                            1.0, &
+                            0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+end subroutine set_meta_nc
+#endif
\ No newline at end of file