diff --git a/.idea/misc.xml b/.idea/misc.xml
index ca8c3e84ff9323893952ecd05e24828cf3b455ac..473d3b4a54f672281993f4993c671fb8d9c52083 100644
--- a/.idea/misc.xml
+++ b/.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
diff --git a/.idea/vcs.xml b/.idea/vcs.xml
index 35eb1ddfbbc029bcab630581847471d7f238ec53..f9e8820ca83979a0fa3891e7e7bfdac8c44d5784 100644
--- a/.idea/vcs.xml
+++ b/.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
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7ef0b8a71a241b21e8bbe7dec1937987b134ad9f..5df81955e23cc0dd653e531e852d5aee9a188ce1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,13 +1,17 @@
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 )
\ No newline at end of file
+set(lib_files
+ 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
diff --git a/src/config-utils.f90 b/src/config-utils.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b93137308298905c83857ac4105ccd0e06edabfe
--- /dev/null
+++ b/src/config-utils.f90
@@ -0,0 +1,160 @@
+! 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
diff --git a/src/diag_pbl.f90 b/src/diag_pbl.f90
new file mode 100644
index 0000000000000000000000000000000000000000..da216f1f39d8bb3baeec2e00b8864d913f483401
--- /dev/null
+++ b/src/diag_pbl.f90
@@ -0,0 +1,220 @@
+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
diff --git a/src/external_forcing.f90 b/src/external_forcing.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a3d72bea0ea491d119a63b87a26bede5d8a101c8
--- /dev/null
+++ b/src/external_forcing.f90
@@ -0,0 +1,114 @@
+! 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
diff --git a/src/parkinds.f90 b/src/parkinds.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ac95e5f40e38f505834e8068579708063d0825a8
--- /dev/null
+++ b/src/parkinds.f90
@@ -0,0 +1,22 @@
+! 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
diff --git a/src/pbl_grid.f90 b/src/pbl_grid.f90
new file mode 100644
index 0000000000000000000000000000000000000000..3624fed417ed5fc3638e9ac4aa00202064330dab
--- /dev/null
+++ b/src/pbl_grid.f90
@@ -0,0 +1,100 @@
+! 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
diff --git a/src/pbl_turb_data.f90 b/src/pbl_turb_data.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9e55a8af34844f7a2d37ca1d16e092a0e8e48202
--- /dev/null
+++ b/src/pbl_turb_data.f90
@@ -0,0 +1,62 @@
+! 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
diff --git a/src/phys_fluid.f90 b/src/phys_fluid.f90
new file mode 100644
index 0000000000000000000000000000000000000000..9cb9d2d767cbdc838b7526c11a64a92917edfbe5
--- /dev/null
+++ b/src/phys_fluid.f90
@@ -0,0 +1,50 @@
+! 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
diff --git a/src/rrtm_interface.f90 b/src/rrtm_interface.f90
new file mode 100644
index 0000000000000000000000000000000000000000..e2556f41156258d56651e9c482efecf4ad7bd7ef
--- /dev/null
+++ b/src/rrtm_interface.f90
@@ -0,0 +1,242 @@
+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
diff --git a/src/scm_io_default.f90 b/src/scm_io_default.f90
new file mode 100644
index 0000000000000000000000000000000000000000..b8eb908b771da936571382c9fed605a4a6e0df1f
--- /dev/null
+++ b/src/scm_io_default.f90
@@ -0,0 +1,146 @@
+
+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
diff --git a/src/scm_sfx_data.f90 b/src/scm_sfx_data.f90
new file mode 100644
index 0000000000000000000000000000000000000000..98df36546cd15401d8abc7f31c6ec55504c68f05
--- /dev/null
+++ b/src/scm_sfx_data.f90
@@ -0,0 +1,15 @@
+! 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
diff --git a/src/scm_state_data.f90 b/src/scm_state_data.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ce08dc9ad7ffae5c7c214ffd012b5c67d7b3a6bd
--- /dev/null
+++ b/src/scm_state_data.f90
@@ -0,0 +1,87 @@
+! 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
diff --git a/src/state_utils.f90 b/src/state_utils.f90
new file mode 100644
index 0000000000000000000000000000000000000000..311241cc30eb23d4c4cba32bcc00c3605cba2f2b
--- /dev/null
+++ b/src/state_utils.f90
@@ -0,0 +1,80 @@
+! 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
diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90
new file mode 100644
index 0000000000000000000000000000000000000000..91478bce2a6c83b0676d29679f4967e6d81c7c8a
--- /dev/null
+++ b/src/tests/gabls1.f90
@@ -0,0 +1,535 @@
+! 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