! Created by Andrey Debolskiy on 12.10.2024.
program gabls2
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 pbl_turb_common
use pbl_dry_contrgradient
use pbl_solver
use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z, &
get_sigma_from_z_theta, get_theta_v, get_density_theta_vector
use diag_pbl
use pbl_grid
use sfx_data, only: meteoDataType, sfxDataType
use sfx_esm, only: get_surface_fluxes_most => get_surface_fluxes, &
numericsType_most => numericsType
use config_parser
use config_utils
implicit none
!local functions
real, external:: get_ts_gabls2
!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-gabls2-ex.txt'
real, allocatable :: ttemp(:,:), utemp(:,:),vtemp(:,:), ttime(:)
real, allocatable:: w_sub(:)
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(pblContrGradDataType):: cbl;
type(turbBLDataType):: turb
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 = 3.00000
vg = -9.0
dt = 10.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_gabls2.dat'
fname_nc = 'gabls2.nc'
series_f%fname = 'phys2.dat'
call set_file(series_f, series_f%fname)
scan_f%fname='time_scan2.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
write(*, *) ' FAILURE! > could not initialize config file: ', fname_config
ierr = 1 ! signal ERROR
stop
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
write(*, *) ' FAILURE! > could not set: time.begin '
ierr = 1 ! signal ERROR
stop
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
write(*, *) ' FAILURE! > could not set: time.end '
ierr = 1 ! signal ERROR
stop
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
write(*, *) ' FAILURE! > could not set: time.dt '
ierr = 1 ! signal ERROR
stop
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("output.dt"//C_NULL_CHAR, output_dt, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: output.dt '
ierr = 1 ! signal ERROR
stop
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
write(*, *) ' FAILURE! > could not set: Ug '
ierr = 1 ! signal ERROR
stop
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
write(*, *) ' FAILURE! > could not set: Vg '
ierr = 1 ! signal ERROR
stop
end if
end if
call c_config_is_varname("fcoriolis"//C_NULL_CHAR, status)
if ((status /= 0)) then
!< mandatory in user dataset
call c_config_get_float("fcoriolis"//C_NULL_CHAR, f_cor, status)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: fcoriolis '
ierr = 1 ! signal ERROR
stop
end if
end if
call set_fluid_default(fluid_params)
call get_fluid_params(fluid_params, status, ierr)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: Fluid params '
ierr = 1 ! signal ERROR
stop
end if
!z coordinate
call get_grid_params(grid, status, ierr)
if (status == 0) then
write(*, *) ' FAILURE! > could not set: grid '
ierr = 1 ! signal ERROR
stop
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)
bl_old%cor= f_cor
bl%cor = f_cor
bl%land_type = 0.0
bl%p0=fluid_params%p0
!setup surface
z0 = 0.03
zh = z0 / 10.0
! init blData states
call scm_data_allocate(bl, kmax)
call scm_data_allocate(bl_old, kmax)
call cbl_allocate(cbl, kmax)
call init_state(bl, ug, vg, fluid_params%tref)
call init_state(bl_old, ug, vg, fluid_params%tref)
call init_theta_gabls2(bl%theta,grid%z_cell, grid%kmax)
call scm_data_copy(bl_old,bl)
allocate(w_sub(0:kmax),source=0.0)
!finish updating grid
call get_sigma_from_z_theta( grid, bl, fluid_params%p0 * 100.0, fluid_params)
call theta2ta(bl_old%temp, bl_old%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
call theta2ta(bl%temp, bl%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl%kmax)
call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
call scm_data_copy(bl_old,bl)
!Initialize turbulent closure
call turb_data_allocate(turb, kmax)
! getclosurefromconfig
call get_closure_params(turb, status, ierr)
!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 begin', nt
call to_file_1d_2var('temperature2.dat', grid%z_cell, bl%theta, grid%kmax)
do while (time_current < time_end)
call put_tscan( time_current/3600.0, grid%z_cell,bl_old, bl_old%kmax, scan_f)
tsurf = get_ts_gabls2(time_current)
meteo_cell%U = sqrt(bl_old%u(kmax)**2 + bl_old%v(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
!meteo_cell%z0_t= z0 / 10.0
call get_surface_fluxes_most(sfx_cell, meteo_cell, numerics_sfx)
call get_density_theta_vector( bl, fluid_params, grid)
bl%surf%cm2u = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U
taux = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U * bl_old%u(kmax)
tauy = bl%rho(kmax)* sfx_cell%Cm**2 * meteo_cell%U * bl_old%v(kmax)
hflx = - bl%rho(kmax) &
* sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U * meteo_cell%dT
bl%surf%hs = hflx
bl%surf%es = sfx_cell%Cm * sfx_cell%Ct * meteo_cell%U *meteo_cell%dQ
turb%ustr = sfx_cell%Cm * meteo_cell%U
umst = turb%ustr
!write(*,*) 'surf', bl%surf%hs, turb%ustr, tsurf, bl%theta(kmax)
call get_density_theta_vector( bl, fluid_params, grid)
call get_theta_v(bl%theta_v, bl%theta, bl%qv, fluid_params, grid%kmax)
call get_coeffs_general(turb, bl, fluid_params, 2, grid)
!write(*,*) 'bl%kpbl', bl%kpbl, grid%z_cell(bl%kpbl), sfx_cell%Cm, sfx_cell%Ct, tsurf
call to_file_1d_3var('coeffs1.dat',grid%z_cell, turb%km, turb%km, bl_old%kmax)
!do k=1,grid%kmax
! write(*,*) 'kdiffs', k, bl%vdctq(k), bl%vdcuv(k)
!end do
!call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
call solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
call new_cntrg(cbl, bl, fluid_params, grid, dt)
!write(*,*) 'dttheta,', bl%theta - bl_old%theta
call scm_data_copy(bl_old,bl)
call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
!write(*,*) 'aftercoriolis,'
call get_wsub(w_sub, time_current, grid)
call apply_subsidence(bl, w_sub, fluid_params, grid, dt)
!call to_file_1d_2var('temperature.dat', grid%z_cell, bl_old%U - bl%U, bl_old%kmax)
call scm_data_copy(bl_old,bl)
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 theta2ta(bl_old%temp, bl_old%theta, &
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) = turb%ustr
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('temperature2.dat', grid%z_cell, bl%theta, bl_old%kmax)
call close_file(series_f)
call close_file(scan_f)
call deallocate_pbl_grid(grid)
call scm_data_deallocate(bl)
call scm_data_deallocate(bl_old)
call pbl_data_deallocate(turb)
call cbl_deallocate(cbl)
end program gabls2
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(inout) :: bl
integer kmax
kmax= bl%kmax
write(*,*) 'blkmax=', kmax
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 init_theta_gabls2(thold, z, kmax)
implicit none
integer,intent(in):: kmax
real, dimension(kmax):: thold, z
integer k
real T1, T2, Z1, Z2
!write(*,*) 'theta_init'
do k = kmax, 1, -1
if (z(k) <= 200.0 ) then
z1 = 0.0
z2 = 200.0
T1 = 288.0
T2 = 286.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 200.0 .and. z(k)<= 850.0 ) then
z1 = 200.0
z2 = 850.0
T1 = 286.0
T2 = 286.0
thold(k) = t2
elseif (z(k) > 850.0 .and. z(k)<= 900.0) then
z1 = 850.0
z2 = 900.0
T1 = 286.0
T2 = 288.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 900.0 .and. z(k)<= 1000.0) then
z1 = 900.0
z2 = 1000.0
T1 = 288.0
T2 = 292.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 1000.0 .and. z(k)<= 2000.0) then
z1 = 1000.0
z2 = 2000.0
T1 = 292.0
T2 = 300.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 2000.0 .and. z(k)<= 3500.0) then
z1 = 2000.0
z2 = 3500.0
T1 = 300.0
T2 = 310.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 3500.0 .and. z(k)<= 4000.0) then
z1 = 3500.0
z2 = 4000.0
T1 = 310.0
T2 = 312.0
thold(k) = (t2 * (z(k) - z1) + t1 * (z2 - z(k))) / (z2 - z1)
elseif (z(k) > 4000.0) then
thold(k) = 312.0
end if
!write(*,*) z(k), thold(k)
end do
end subroutine init_theta_gabls2
real function get_ts_gabls2(curtime)
implicit none
real curtime
real hrs, ts
hrs = curtime / 3600.0 + 16.0
ts = 4.4
if (hrs <= 17.4) then
ts = -10.0 - 25.0 * cos(0.22 * hrs + 0.2)
end if
if ((hrs > 17.4) .and. (hrs <= 30.0)) then
ts = -0.54 * hrs + 15.2
end if
if ((hrs > 30.0) .and. (hrs <= 41.9)) then
ts = -7.0 - 25.0 * cos(0.21 * hrs + 1.8)
end if
if ((hrs > 41.9) .and. (hrs <= 53.3)) then
ts = -0.37 * hrs + 18.0;
end if
if ((hrs > 53.3) .and. (hrs <= 65.6)) then
ts = -4.0 - 25.0 * cos(0.22 * hrs + 2.5)
end if
ts = ts + 273.15
get_ts_gabls2 = ts
end function get_ts_gabls2
subroutine add_coriolis(bl, bl_old, 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(inout):: bl
type(stateBLDataType), intent(in):: bl_old
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 get_wsub(w, curtime, grid)
use pbl_grid, only: pblgridDataType
implicit none
type(pblgridDataType), intent(in):: grid
real, intent(in):: curtime
real, dimension(*), intent(out):: w
integer k
if(curtime >= 24.0* 3600.0) then
do k= grid%kmax, 1, -1
if (grid%z_edge(k) <= 1000.0) then
w(k) = - 0.005 * grid%z_edge(k) / 1000.0
else
w(k) = - 0.005
end if
write(*,*) 'ws', k, w(k)
end do
else
w(1:grid%kmax) = 0.0
end if
end subroutine get_wsub
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
stop
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