doc/app/iofiles/pids/aerosol: create_basic_salsa_dynamic_driver.py

#------------------------------------------------------------------------------#
# This file is part of the PALM model system.
#
# PALM is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# PALM is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# PALM. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright 2019-2019 University of Helsinki
#------------------------------------------------------------------------------#

import math
from   netCDF4 import Dataset
import numpy as np
import datetime


class SalsaDynamicDriver:
    """ This is an example script to generate dynamic drivers that include aerosol 
    backgound concentrations for PALM. 

    You can use it as a starting point for creating your setup specific
    driver.  
    """
    
    
    def __init__(self):
        """ Open the salsa dynamic driver as NetCDF4 file. Here, you have to
        give the full path to the static driver that shall be created. Existing
        file with same name is deleted.
        """
        print('Opening file...')
        self.nc_file = Dataset('salsa_dynamic_driver.nc', 'w', format='NETCDF4')


    def write_global_attributes(self):
        """ Write global attributes to static driver. """
        print("Writing global attributes...")

        # optional global attributes
        self.nc_file.origin_lon = 55.0 # used to initialize coriolis parameter
        self.nc_file.origin_lat = 0.0 # (overwrite initialization_parameters)
        self.nc_file.origin_time = '2000-06-21 12:00:00 +00'
        self.nc_file.origin_x = 308124
        self.nc_file.origin_y = 6098908
        self.nc_file.origin_z = 0.0
        self.nc_file.rotation_angle = 0.0
        self.nc_file.author = 'Your Name'
        self.nc_file.comment = 'Miscellaneous information about the data ' \
                               'or methods to produce it.'
        self.nc_file.creation_date = str(datetime.datetime.now())
        self.nc_file.institution = 'INAR/Physics, University of Helsinki'
        self.nc_file.history = ''
        self.nc_file.palm_revision = ''
        self.nc_file.title = 'Salsa dynamic driver for some arbitrary PALM setup'


    def define_dimensions(self):
        """ Set dimensions on which variables are defined. """
        print("Writing dimensions...")

        # General grid parameters:

        self.nx = 19
        self.ny = 19
        self.nz = 20
        dx = 2
        dy = 2
        dz = 2

        # time steps in the background concentration data:
        self.times = np.arange( 0.0, 7201.0, 3600.0 )


        # Aerosol composition and size distribution:

        # number of chemical components
        self.ncomposition_index = 7

        # number of aerosol size bins in the subrange 1 and 2
        nbin   = [1, 7]

        # subrange diameter limit (e.g. subrange 1: 3-10 nm, subrange 2: 10 nm - 2.5 um
        reglim = [3.0e-9, 10.0e-9, 2.5e-6]


        # Mandatory dimensions
        self.nmax_string_length = 25

        self.ntime = len( self.times )

        self.nc_file.createDimension('x' ,self.nx+1)
        self.x = self.nc_file.createVariable('x', 'f8', ('x',))
        self.x.units = 'm'
        self.x.standard_name = 'x coordinate of cell centers'
        self.x[:] = np.arange( 0.5*dx, ( self.nx+1 ) * dx, dx )

        self.nc_file.createDimension('y', self.ny+1)
        self.y = self.nc_file.createVariable('y', 'f8', ('y',))
        self.y.units = 'm'
        self.y.standard_name = 'y coordinate of cell centers'
        self.y[:] = np.arange( 0.5*dy, ( self.ny+1 ) * dy, dy )

        self.nc_file.createDimension('z', self.nz)
        self.z = self.nc_file.createVariable('z', 'f8', ('z',))
        self.z.units = 'm'
        self.z.standard_name = 'z coordinate of cell centers'
        self.z[:] = np.arange( 0.5*dz, ( self.nz ) * dz, dz )  # NO STRETCHING!

        self.nc_file.createDimension('max_string_length', self.nmax_string_length )
        self.max_string_length = self.nc_file.createVariable('max_string_length', 'i4',
                                                             ('max_string_length',))
        self.max_string_length.units = ''
        self.max_string_length[:] = np.linspace( 1, self.nmax_string_length, self.nmax_string_length )

        self.nc_file.createDimension('composition_index', self.ncomposition_index )
        self.composition_index = self.nc_file.createVariable('composition_index', 'i4', 
                                                             ('composition_index',))
        self.composition_index.units = ''
        self.composition_index.long_name = 'composition index'
        self.composition_index[:] = np.linspace( 1, self.ncomposition_index, self.ncomposition_index )

        self.nc_file.createDimension('time' , self.ntime )
        self.time = self.nc_file.createVariable('time', 'f8', ('time',))
        self.time.units = 's'
        self.time.standard_name = 'time since utc init'
        self.time[:] = self.times

        self.ndmid, self.bin_limits = define_bins( nbin, reglim )
        self.nc_file.createDimension('Dmid' , np.sum( nbin ) )
        self.Dmid = self.nc_file.createVariable('Dmid', 'f8', ('Dmid',))
        self.Dmid.units = 'm'
        self.Dmid.standard_name = 'mean diamater per size bin'
        self.Dmid[:] = self.ndmid

    def add_variables(self):
        """ Uncomment variables below as you like. 

        Be aware that some variables depend on others. For a description of 
        each variable, please have a look at the documentation at
        https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/pids/aerosol

        An example of how you modify the variables is given below:

        building_2d_array = np.ones((self.ny+1,self.nx+1)) * -9999.0
        south_wall, north_wall, left_wall, right_wall = 20, 25, 20, 25
        building_2d_array[south_wall:north_wall,left_wall:right_wall] = 50
        nc_buildings_2d = self.nc_file.createVariable(
            'buildings_2d', 'f4', ('y','x'),fill_value=-9999.0)
        nc_buildings_2d.lod = 1
        nc_buildings_2d[:,:] = building_2d_array
        """
        print("Writing variables...")

        composition_name_list = ['H2SO4                    ','OC                       ',
                                 'BC                       ','DU                       ',
                                 'SS                       ','HNO3                     ', 
                                 'NH3                      ']

        lsf_mf_a = np.ones( ( self.ntime, self.nz, self.ncomposition_index ) )
        lsf_mf_b = np.ones( ( self.ntime, self.nz, self.ncomposition_index ) )
        lsf_psd  = np.ones( ( self.ntime, self.nz, len( self.ndmid ) ) )

        dpg       = np.array([20.3, 72.0])
        n_lognorm = np.array([18960.0, 13750.0])
        sigmag    = np.array([1.7, 1.6])

        # sectional size distribution
        nsect = np.zeros( len( self.ndmid ), dtype=float )
        for l in range( 1, len( self.ndmid )+1 ): 
          d1 = self.bin_limits[l-1]
          d2 = self.bin_limits[l]
          delta_d = ( d2 - d1 ) / 10.0
          for ib in range( 1, len( self.ndmid )+1 ):
            d1 = self.bin_limits[l-1] + ( ib - 1 ) * delta_d
            d2 = d1 + delta_d
            dmidi = ( d1 + d2 ) / 2.0
            deltadp = np.log10( d2 / d1 )
            nsect[l-1] = np.sum( n_lognorm * 1.0E6 * deltadp / ( np.sqrt( 2.0 * np.pi ) * 
                                                                 np.log10( sigmag ) ) *
                         np.exp( -np.log10( dmidi / ( 1.0E-9 * dpg ) )**2.0 / 
                           ( 2.0 * np.log10( sigmag ) ** 2.0 ) ) )

        # set a constant profile:
        for t in range( self.ntime ):
          for k in range( self.nz ):
            lsf_mf_a[t,k,:] = np.array( [0.001, 0.699, 0.1, 0.0, 0.0, 0.1, 0.1] )
            lsf_mf_b[t,k,:] = np.array( [0.0,   0.0,   0.0, 0.0, 0.0, 0.0, 0.0] )
            lsf_psd[t,k,:]  = nsect


        # Save into the file
        nc_composition_name = self.nc_file.createVariable( 'composition_name', 'S1',
                                                         ('composition_index','max_string_length',) )
        nc_composition_name[:] = list( map( lambda x : list(x), composition_name_list ) )
        nc_composition_name.long_name = 'aerosol composition name'
        nc_composition_name.standard_name = 'composition_name'


        nc_init_mf_a = self.nc_file.createVariable( 'init_atmosphere_mass_fracs_a', 'f4',
                                                    ('z','composition_index',), fill_value=-9999.0 )
        nc_init_mf_a[:] = lsf_mf_a[0,:,:]
        nc_init_mf_a.long_name = "initial mass fraction profile: a bins"
        nc_init_mf_a.standard_name = 'init_atmosphere_mass_fracs_a'


        nc_init_mf_b = self.nc_file.createVariable( 'init_atmosphere_mass_fracs_b', 'f4',
                                                    ('z','composition_index',), fill_value=-9999.0 )
        nc_init_mf_b[:] = lsf_mf_b[0,:,:]
        nc_init_mf_b.long_name = "initial mass fraction profile: b bins"
        nc_init_mf_b.standard_name = 'init_atmosphere_mass_fracs_b'


        nc_init_psd = self.nc_file.createVariable( 'init_atmosphere_aerosol', 'f4', ('z','Dmid',),
                                                   fill_value=-9999.0 )
        nc_init_psd[:] = lsf_psd[0,:,:]
        nc_init_psd.unit = '#/m3'
        nc_init_psd.long_name = 'initial vertical profile of aerosol concentration'
        nc_init_psd.standard_name = 'init_atmosphere_aerosol'
        nc_init_psd.lod = 1

        ls_names = ['left','right','north','south']
        ls_dims = [ self.y, self.y, self.x, self.x ]
        ls_dimsname = ['y','y','x','x']

        nvi = 0
        for nv in ls_names:
          dummy = np.zeros( [ self.ntime, self.nz, len( ls_dims[nvi] ), self.ncomposition_index ] )

          for i in range( len( ls_dims[nvi] ) ):
            dummy[:,:,i,:] = lsf_mf_a
          namev = 'ls_forcing_{}_mass_fracs_a'.format( nv )
          nc_lsf_mf_a = self.nc_file.createVariable( namev, 'f4', ('time','z', ls_dimsname[nvi],
                                                                   'composition_index',),
                                                     fill_value=-9999.0 )
          nc_lsf_mf_a[:] = dummy
          nc_lsf_mf_a.long_name = "boundary conditions of mass fraction profile: a bins"
          nc_lsf_mf_a.standard_name = namev

          for i in range( len( ls_dims[nvi] ) ):
            dummy[:,:,i,:] = lsf_mf_b
          namev = 'ls_forcing_{}_mass_fracs_b'.format( nv )
          nc_lsf_mf_b = self.nc_file.createVariable( namev, 'f4', ('time','z',ls_dimsname[nvi],
                                                                   'composition_index',),
                                                     fill_value=-9999.0 )
          nc_lsf_mf_b[:] = dummy
          nc_lsf_mf_b.long_name = "boundary conditions of mass fraction profile: b bins"
          nc_lsf_mf_b.standard_name = namev

          nvi += 1

        dummy = np.zeros( [ self.ntime, self.ny+1, self.nx+1, self.ncomposition_index ] )
        for j in range( self.ny ):
          for i in range( self.nx ):
            dummy[:,j,i,:] = lsf_mf_a[:,-1,:]
        namev = 'ls_forcing_top_mass_fracs_a'
        nc_lsf_mf_a = self.nc_file.createVariable( namev, 'f4', ('time','y','x','composition_index',), 
                                                   fill_value=-9999.0 )
        nc_lsf_mf_a[:] = dummy
        nc_lsf_mf_a.long_name = "boundary conditions of mass fraction profile: a bins"
        nc_lsf_mf_a.standard_name = namev

        for j in range( self.ny ):
          for i in range( self.nx ):
            dummy[:,j,i,:] = lsf_mf_b[:,-1,:]
        namev = 'ls_forcing_top_mass_fracs_b'
        nc_lsf_mf_b = self.nc_file.createVariable( namev, 'f4', ('time','y','x','composition_index',), 
                                                   fill_value=-9999.0 )
        nc_lsf_mf_b[:] = dummy
        nc_lsf_mf_b.long_name = "boundary conditions of mass fraction profile: b bins"
        nc_lsf_mf_b.standard_name = namev


        nvi = 0
        for nv in ls_names:
          dummy = np.zeros( [ self.ntime, self.nz, len( ls_dims[nvi] ), len( self.ndmid ) ] )
          for i in range( len( ls_dims[nvi] ) ):
            dummy[:,:,i,:] = lsf_psd
          namev = 'ls_forcing_{}_aerosol'.format( nv )
          nc_lsf_psd = self.nc_file.createVariable( namev, 'f4', ('time','z',ls_dimsname[nvi],'Dmid',), 
                                                    fill_value=-9999.0 )
          nc_lsf_psd[:] = dummy
          nc_lsf_psd.long_name = 'boundary condition of aerosol concentration'
          nc_lsf_psd.standard_name = namev
          nc_lsf_psd.unit = '#/m3'
          nc_lsf_psd.lod = 1
          nvi += 1

        dummy =  np.zeros( [ self.ntime, self.ny+1, self.nx+1, len( self.ndmid ) ], dtype=float ) 
        for j in range( self.ny ):
          for i in range( self.nx ):
            dummy[:,j,i,:] = lsf_psd[:,-1,:]
        namev = 'ls_forcing_top_aerosol'
        nc_lsf_psd = self.nc_file.createVariable( 'ls_forcing_top_aerosol', 'f4', 
                                                 ('time','y','x','Dmid',), fill_value=-9999.0 )
        nc_lsf_psd[:] = dummy
        nc_lsf_psd.long_name = 'boundary condition of aerosol concentration'
        nc_lsf_psd.standard_name = 'ls_forcing_top_aerosol'
        nc_lsf_psd.unit = '#/m3'
        nc_lsf_psd.lod = 1
        nvi += 1


    def finalize(self):
        """ Close file """
        print("Closing file...")
        
        self.nc_file.close()

def define_bins( nbin, reglim ):
    """ This function defines the sectional bin limits based on number of bins
        (nbin) and diameter limits (reglim) per subrange
    """

    nbins = np.sum( nbin ) # = subrange 1 + subrange 2

    # Log-normal to sectional

    vlolim = np.zeros( nbins )
    vhilim = np.zeros( nbins )
    dmid   = np.zeros( nbins )
    bin_limits = np.zeros( nbins )

    # Sectional bin limits
    ratio_d = reglim[1] / reglim[0]
    for b in range( nbin[0] ):
      vlolim[b] = np.pi / 6.0 * ( reglim[0] * ratio_d **( float(b) / nbin[0] ) )**3
      vhilim[b] = np.pi / 6.0 * ( reglim[0] * ratio_d **( float(b+1) / nbin[0] ) )**3
      dmid[b] = np.sqrt( ( 6.0 * vhilim[b] / np.pi )**0.33333333 * ( 6.0 * vlolim[b] / np.pi 
                                                                    )**0.33333333 )

    ratio_d = reglim[2] / reglim[1]
    for b in np.arange( nbin[0], np.sum( nbin ),1 ):
      c = b-nbin[0]
      vlolim[b] = np.pi / 6.0 * ( reglim[1] * ratio_d ** ( float(c) / nbin[1] ) )**3
      vhilim[b] = np.pi / 6.0 * ( reglim[1] * ratio_d ** ( float(c+1) / nbin[1] ) ) ** 3
      dmid[b] = np.sqrt( ( 6.0 * vhilim[b] / np.pi )**0.33333333 * ( 6.0 * vlolim[b] / np.pi 
                                                                    )**0.33333333 )

    bin_limits = ( 6.0 * vlolim / np.pi )**0.33333333
    bin_limits = np.append( bin_limits, reglim[-1] )

    return dmid, bin_limits


if __name__ == '__main__':
    driver = SalsaDynamicDriver()
    driver.write_global_attributes()
    driver.define_dimensions()
    driver.add_variables()
    driver.finalize()