Skip to content
Snippets Groups Projects
gabls2.f90 22.5 KiB
Newer Older
  • Learn to ignore specific revisions
  • ! 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_most, 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)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            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
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            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 = bl%rho(grid%kmax) * 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
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            !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, 1,  grid)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            !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 solve_diffusion(bl, bl_old, turb, fluid_params, grid, dt)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            !call cbl_apply_cntrg(cbl, bl, fluid_params, grid, dt)
            !write(*,*) 'dttheta,', bl%theta - bl_old%theta
    
            call add_coriolis(bl, bl_old, ug, vg, f_cor, dt, grid)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            !write(*,*) 'aftercoriolis,'
    
            call    get_wsub(w_sub, time_current, grid)
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
            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
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                    w(k) = - 0.005 * grid%z_edge(k) / 1000.0
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                    w(k) = - 0.005
    
    Debolskiy Andrey's avatar
    Debolskiy Andrey committed
                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