reworked some functions

config-examples/config-gabls1-ex.txt

+
+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

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

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

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

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

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