diff --git a/config-examples/config-gabls1-ex.txt b/config-examples/config-gabls1-ex.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5f885834cbdf2e3f3d943aca4226f69cd21ee8a6
--- /dev/null
+++ b/config-examples/config-gabls1-ex.txt
@@ -0,0 +1,29 @@
+
+Ug = 8.0;
+Vg = 0.0;
+time{
+ begin = 0.0;
+ end = 9.0 * 3600.0;
+ dt = 20.0;
+ out_dt = 360.0;
+ }
+
+fluid {
+ pref= 1013.4; #hPa
+ tref = 265.0; #K
+ g = 9.81; # m/s2
+ beta = g /tref;
+ kappa = 0.4; # von Karman constant
+ }
+
+surface_model {
+
+ id = "esm"; # optional: "esm" (default), "sheba", "most", "log"
+ type = "land";
+ z0_m = 0.01; # aerodynamic roughness [m]
+ z0_h = -1; # no prescribed value # -> using scheme assigned by surface type
+ }
+
+output{
+ dt = 0.2 *3600.0;
+}
\ No newline at end of file
diff --git a/src/config-utils.f90 b/src/config-utils.f90
index a917496a327c0ffadb1f72a762ef1d4a7bee5bfe..d3026263401b6b2e6d450cb8539a0352b10ee787 100644
--- a/src/config-utils.f90
+++ b/src/config-utils.f90
@@ -8,7 +8,7 @@ module config_utils
implicit none
integer, public, save:: is_config_initialized = 0
public init_config, get_fluid_params, get_grid_params
- public get_geo_forcing, get_heat_forcing
+ !public get_geo_forcing, get_heat_forcing
contains
diff --git a/src/pbl_fo_turb.f90 b/src/pbl_fo_turb.f90
index fe140e683046c6692f05c15e51f2641a0f4e8e3d..959dbccb60bed4a3e57a342506f1b31bed6a659c 100644
--- a/src/pbl_fo_turb.f90
+++ b/src/pbl_fo_turb.f90
@@ -5,7 +5,44 @@ module pbl_fo_turb
public get_fo_coeffs
public get_turb_length
- public
+ public get_fo_diff
+ public get_rig
contains
+ subroutine get_fo_diff(turb, bl, fluid, grid)
+ use pbl_turb_data, only: turbBLDataType
+ use scm_state_data, only: stateBLDataType
+ use pbl_grid, only: pblgridDataType
+ use phys_fluid, only: fluidParamsDataType
+ end subroutine get_fo_diff
+ subroutine get_fo_coeffs(turb, bl, fluid, grid)
+ use pbl_turb_data, only: turbBLDataType
+ use scm_state_data, only: stateBLDataType
+ use pbl_grid, only: pblgridDataType
+ use phys_fluid, only: fluidParamsDataType
+ end subroutine get_fo_coeffs
+
+ subroutine get_turb_length(turb, bl, fluid, grid, hbl_option)
+ use pbl_turb_data, only: turbBLDataType
+ use scm_state_data, only: stateBLDataType
+ use pbl_grid, only: pblgridDataType
+ use phys_fluid, only: fluidParamsDataType
+ use diag_pbl
+ implicit none
+ type(fluidParamsDataType), intent(in):: fluid
+ type(stateBLDataType), intent(in):: bl
+ type(pblgridDataType), intent(in)::grid
+ type(turbBLDataTyoe), intent(inout):: turb
+
+ call get_rig(turb, bl, fluid, grid)
+
+
+ end subroutine get_turb_length
+
+ subroutine get_rig(turb, bl, fluid, grid)
+ use pbl_turb_data, only: turbBLDataType
+ use scm_state_data, only: stateBLDataType
+ use pbl_grid, only: pblgridDataType
+ use phys_fluid, only: fluidParamsDataType
+ end subroutine get_rig
end module pbl_fo_turb
\ No newline at end of file
diff --git a/src/pbl_grid.f90 b/src/pbl_grid.f90
index da0333f030e16b9a171f097f509adf4c73ec3de8..90a0d0a992106fec0a76e8b05e1efad16cc7c153 100644
--- a/src/pbl_grid.f90
+++ b/src/pbl_grid.f90
@@ -10,7 +10,9 @@ module pbl_grid
real(kind=rf), allocatable :: dzc(:)
real(kind=rf), allocatable :: dze(:)
real(kind=rf), allocatable :: sigma(:)
- real(kind=rf), allocatable :: dsigma(:)
+ real(kind=rf), allocatable :: sigma_edge(:)
+ real(kind=rf), allocatable :: dsigmac(:)
+ real(kind=rf), allocatable :: dsigmae(:)
integer(kind=im) :: kmax
end type pblgridDataType
!grid types
@@ -23,11 +25,8 @@ module pbl_grid
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
- public :: set_inmcm_21
-
- public :: get_sig_from_z
- public :: get_z_from_sig
+ !public :: set_from_file_z
+ !public :: set_inmcm_21
contains
subroutine allocate_pbl_grid(grid, nk)
@@ -40,7 +39,9 @@ module pbl_grid
allocate(grid%dzc(nk))
allocate(grid%dze(nk))
allocate(grid%sigma(nk))
- allocate(grid%dsigma(nk))
+ allocate(grid%dsigmac(nk))
+ allocate(grid%sigma_edge(0:nk))
+ allocate(grid%dsigmae(nk))
end subroutine allocate_pbl_grid
subroutine deallocate_pbl_grid(grid)
@@ -51,7 +52,9 @@ module pbl_grid
deallocate(grid%dzc)
deallocate(grid%dze)
deallocate(grid%sigma)
- deallocate(grid%dsigma)
+ deallocate(grid%dsigmac)
+ deallocate(grid%sigma_edge)
+ deallocate(grid%dsigmae)
end subroutine deallocate_pbl_grid
subroutine set_pbl_grid_via_edge(grid,zpos_cell,zpos_edge,nk)
@@ -72,11 +75,11 @@ module pbl_grid
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)
+ grid%dze(k+1) = zpos_edge(k+1) - zpos_edge(k)
+ grid%dzc(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%dze(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
@@ -97,8 +100,11 @@ module pbl_grid
grid%z_cell(k) = 0.5 * (grid%z_edge(k+1) + grid%z_edge(k))
end do
grid%z_edge(0) = h_top
-
+ grid%dzc(1:nk) = dz
+ grid%dze(1:nk) = dz
end subroutine set_pbl_grid_uniform
+
+
end module pbl_grid
\ No newline at end of file
diff --git a/src/state_utils.f90 b/src/state_utils.f90
index 2d0eebcfe5e874660856f6a1d28e3b6fc7abf9c4..d4c3af35bbe0efa953110dbe6c1fba8e194b815f 100644
--- a/src/state_utils.f90
+++ b/src/state_utils.f90
@@ -3,11 +3,41 @@
module state_utils
public geo
+ public get_geopotential_moist
public theta2ta
public ta2theta
+ public get_exner
+ public get_sigma_from_z
+ public get_z_from_sigma
+ public get_density_theta
+ public get_density
public dqsdt_sc, qmax_sc, qmax
public DQSDT2, DQSDT, CLDT
contains
+ subroutine get_geopotential_moist(F,t, qv, sig, fs, rd, eps, kmax)
+ implicit none
+ integer, intent(in):: kmax
+ real, intent(in):: t(kmax)
+ real, intent(in):: qv(kmax)
+ real, intent(in):: sig(kmax)
+ real, intent(in):: rd
+ real, intent(in):: eps
+ real, intent(in):: fs
+ real, intent(out)::f(kmax)
+
+ integer ka, k
+ real thalf, qhalf
+
+ F(kmax) = fs - log(sig(kmax)) * rd * t(kmax) * (1.0 + eps * qv(kmax)) / (1.0 + qv(kmax))
+
+ do k = kmax-1, 1, -1
+ ka = kmax+1
+ thalf = 0.5 * (t(ka) + t(k))
+ qhalf = 0.5 * (qv(ka) + qv(k))
+ F(k) = F(ka) - (log(sig(ka)) -log(sig(k))) &
+ * rd * thalf * (1.0 + eps * qhalf) / (1.0 + qhalf)
+ end do
+ end subroutine get_geopotential_moist
SUBROUTINE GEO(T,FS,F,SIG, R, LL)
implicit none
@@ -136,4 +166,113 @@ module state_utils
RETURN
END FUNCTION DQSDT2
+ real function get_density( p, t, qv, fluid)
+ use phys_fluid, only: fluidParamsDataType
+ implicit none
+ type(fluidParamsDataType),intent(in):: fluid
+ real, intent(in):: p, t, qv
+
+ get_density = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
+ end function get_density
+
+ real function get_density_theta( p, th, qv, fluid)
+ use phys_fluid, only: fluidParamsDataType
+ implicit none
+ type(fluidParamsDataType),intent(in):: fluid
+ real, intent(in):: p, th, qv
+
+ real t
+ t = (p / fluid%p0)**(fluid%R_d / fluid%cp)
+ get_density_theta = (p / (fluid%R_d * t)) * (1.0 + qv) / (1 + qv * fluid%eps_v)
+ end function get_density_theta
+
+ subroutine get_exner( exnr, p, fluid, kmax)
+ use phys_fluid, only: fluidParamsDataType
+ implicit none
+ integer,intent(in):: kmax
+ type(fluidParamsDataType), intent(in):: fluid
+ real,dimension(kmax), intent(in):: p
+ real,dimension(kmax), intent(out):: exnr
+
+ integer k
+
+ do k = kmax,1,-1
+ exnr(k) = (p(k)/fluid%p0)**(fluid%R_d / fluid%cp)
+ end do
+ end subroutine get_exner
+
+ subroutine get_exner_sig( exnr, ps, fluid, grid, kmax)
+ use pbl_grid, only : pblgridDataType
+ use phys_fluid, only: fluidParamsDataType
+ implicit none
+ integer,intent(in):: kmax
+ type(fluidParamsDataType), intent(in):: fluid
+ type(pblgridDataType), intent(in)::grid
+ real, intent(in):: ps
+ real,dimension(kmax), intent(out):: exnr
+
+ integer k
+ do k = kmax,1,-1
+ exnr(k) = (ps * grid%sigma(k)/fluid%p0)**(fluid%R_d / fluid%cp)
+ end do
+ end subroutine get_exner_sig
+
+ subroutine get_sigma_from_z( grid, bl, ps, fluid)
+ use pbl_grid, only : pblgridDataType
+ use phys_fluid, only: fluidParamsDataType
+ use scm_state_data, only: stateBLDataType
+ implicit none
+ type(fluidParamsDataType), intent(in):: fluid
+ type(stateBLDataType), intent(in):: bl
+ real, intent(in):: ps
+ type(pblgridDataType), intent(inout)::grid
+
+
+ real dp, rho, p
+ integer k, kmax
+ kmax = grid%kmax
+ grid%sigma_edge(kmax) = 1.0
+ rho = get_density(ps, bl%temp(kmax), bl%qv(kmax), fluid)
+ dp = - rho * fluid%g * grid%dzc(kmax)
+ grid%sigma(kmax) = (ps + dp/2) / ps
+ grid%sigma_edge(kmax) = (ps + dp) / ps
+ p = ps + dp
+ do k = kmax-1,1,-1
+ rho = get_density(p, bl%temp(k), bl%qv(k), fluid)
+ dp = - rho * fluid%g * grid%dzc(k)
+ grid%sigma_edge(k) = grid%sigma_edge(k+1) + dp / ps
+ rho = get_density(p, 0.5*(bl%temp(k+1) + bl%temp(k)), &
+ 0.5*(bl%qv(k+1) + bl%qv(k)), fluid)
+ grid%sigma(k) = grid%sigma(k+1) &
+ - rho * fluid%g * grid%dze(k)
+ p = grid%sigma_edge(k) * ps
+ end do
+ grid%dsigmac(2:kmax) = grid%sigma(2:kmax) - grid%sigma(1:kmax-1)
+ grid%sigma_edge(0) = 0.0
+ grid%dsigmae(1:kmax) = grid%sigma_edge(1:kmax) - grid%sigma_edge(1:kmax-1)
+ end subroutine get_sigma_from_z
+
+ subroutine get_z_from_sigma( grid, bl, zs, fluid, kmax)
+ use pbl_grid, only : pblgridDataType
+ use phys_fluid, only: fluidParamsDataType
+ use scm_state_data, only: stateBLDataType
+ implicit none
+ integer,intent(in):: kmax
+ type(fluidParamsDataType), intent(in):: fluid
+ type(stateBLDataType), intent(in):: bl
+ real, intent(in):: zs
+ type(pblgridDataType), intent(inout)::grid
+
+ real fs
+ real, dimension(kmax):: f
+ integer k
+ call get_geopotential_moist(f,bl%temp, bl%qv, grid%sigma, &
+ fluid%g * zs, fluid%R_d, fluid%eps_v, kmax)
+ grid%z_cell(1:kmax) = f(1:kmax) / fluid%g
+ grid%z_edge(kmax) = zs
+ do k = kmax - 1, 0, -1
+ grid%z_edge(k) = grid%z_cell(kmax + 1) + 0.5 *( grid%z_cell(kmax + 1) - grid%z_edge(kmax + 1))
+ end do
+ grid%dzc(2:kmax) = grid%z_cell(2:kmax) - grid%z_cell(1:kmax-1)
+ end subroutine get_z_from_sigma
end module state_utils
\ No newline at end of file
diff --git a/src/tests/gabls1.f90 b/src/tests/gabls1.f90
index 8bd572b69dec5a21ddfd4604cd6ffa9ae1cc4319..0ddb027467bd967ddd9ce3419f6338765abae95e 100644
--- a/src/tests/gabls1.f90
+++ b/src/tests/gabls1.f90
@@ -11,7 +11,7 @@ program gabls1
use scm_state_data
use pbl_turb_data
use pbl_solver
- use state_utils, only : geo, theta2ta, ta2theta
+ use state_utils, only : geo, theta2ta, ta2theta, get_sigma_from_z
use diag_pbl
use pbl_grid
use sfx_data, only: meteoDataType, sfxDataType
@@ -193,6 +193,8 @@ program gabls1
END DO
call theta2ta(bl_old%temp, bl_old%theta, fluid_params%pref / 1000.0, grid%sigma, 287.5 / 1003.5, bl_old%kmax)
+ !finish updating grid
+ call get_sigma_from_z( grid, bl, fluid_params%p0 * 100.0, fluid_params)
!io setup
@@ -302,45 +304,6 @@ subroutine init_state(bl, ug_,vg_,tref)
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