From 33db9380b3eb8e0034a8efab44dffe053e37d673 Mon Sep 17 00:00:00 2001 From: Andrey Debolskiy <and.debol@gmail.com> Date: Fri, 29 Nov 2024 11:25:01 +0300 Subject: [PATCH] gabls1 experiment added --- .idea/misc.xml | 4 + .idea/vcs.xml | 1 + CMakeLists.txt | 100 +++++++- src/config-utils.f90 | 160 ++++++++++++ src/diag_pbl.f90 | 220 ++++++++++++++++ src/external_forcing.f90 | 114 +++++++++ src/parkinds.f90 | 22 ++ src/pbl_grid.f90 | 100 ++++++++ src/pbl_turb_data.f90 | 62 +++++ src/phys_fluid.f90 | 50 ++++ src/rrtm_interface.f90 | 242 ++++++++++++++++++ src/scm_io_default.f90 | 146 +++++++++++ src/scm_sfx_data.f90 | 15 ++ src/scm_state_data.f90 | 87 +++++++ src/state_utils.f90 | 80 ++++++ src/tests/gabls1.f90 | 535 +++++++++++++++++++++++++++++++++++++++ 16 files changed, 1934 insertions(+), 4 deletions(-) create mode 100644 src/config-utils.f90 create mode 100644 src/diag_pbl.f90 create mode 100644 src/external_forcing.f90 create mode 100644 src/parkinds.f90 create mode 100644 src/pbl_grid.f90 create mode 100644 src/pbl_turb_data.f90 create mode 100644 src/phys_fluid.f90 create mode 100644 src/rrtm_interface.f90 create mode 100644 src/scm_io_default.f90 create mode 100644 src/scm_sfx_data.f90 create mode 100644 src/scm_state_data.f90 create mode 100644 src/state_utils.f90 create mode 100644 src/tests/gabls1.f90 diff --git a/.idea/misc.xml b/.idea/misc.xml index ca8c3e8..473d3b4 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 35eb1dd..f9e8820 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 7ef0b8a..5df8195 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 0000000..b931373 --- /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 0000000..da216f1 --- /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 0000000..a3d72be --- /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 0000000..ac95e5f --- /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 0000000..3624fed --- /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 0000000..9e55a8a --- /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 0000000..9cb9d2d --- /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 0000000..e2556f4 --- /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 0000000..b8eb908 --- /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 0000000..98df365 --- /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 0000000..ce08dc9 --- /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 0000000..311241c --- /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 0000000..91478bc --- /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 -- GitLab