! 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