uEMEP_read_weekly_shipping_asi_data Subroutine

public subroutine uEMEP_read_weekly_shipping_asi_data()

Uses

  • proc~~uemep_read_weekly_shipping_asi_data~~UsesGraph proc~uemep_read_weekly_shipping_asi_data uEMEP_read_weekly_shipping_asi_data module~uemep_definitions uEMEP_definitions proc~uemep_read_weekly_shipping_asi_data->module~uemep_definitions

Arguments

None

Calls

proc~~uemep_read_weekly_shipping_asi_data~~CallsGraph proc~uemep_read_weekly_shipping_asi_data uEMEP_read_weekly_shipping_asi_data proc~date_to_number date_to_number proc~uemep_read_weekly_shipping_asi_data->proc~date_to_number proc~datestr_to_date datestr_to_date proc~uemep_read_weekly_shipping_asi_data->proc~datestr_to_date proc~lb2lambert2_uemep lb2lambert2_uEMEP proc~uemep_read_weekly_shipping_asi_data->proc~lb2lambert2_uemep proc~ll2ps_spherical LL2PS_spherical proc~uemep_read_weekly_shipping_asi_data->proc~ll2ps_spherical proc~nxtdat nxtdat proc~uemep_read_weekly_shipping_asi_data->proc~nxtdat proc~proj2ll PROJ2LL proc~uemep_read_weekly_shipping_asi_data->proc~proj2ll sngl sngl proc~date_to_number->sngl proc~laea2ll LAEA2LL proc~proj2ll->proc~laea2ll proc~lambert2lb2_uemep lambert2lb2_uEMEP proc~proj2ll->proc~lambert2lb2_uemep proc~ltm2ll ltm2ll proc~proj2ll->proc~ltm2ll proc~ps2ll_spherical PS2LL_spherical proc~proj2ll->proc~ps2ll_spherical proc~rdm2ll RDM2LL proc~proj2ll->proc~rdm2ll proc~utm2ll utm2ll proc~proj2ll->proc~utm2ll dcos dcos proc~ltm2ll->dcos dsin dsin proc~ltm2ll->dsin dtan dtan proc~ltm2ll->dtan proc~utm2ll->dcos proc~utm2ll->dsin proc~utm2ll->dtan

Called by

proc~~uemep_read_weekly_shipping_asi_data~~CalledByGraph proc~uemep_read_weekly_shipping_asi_data uEMEP_read_weekly_shipping_asi_data proc~uemep_calculate_emissions_for_emep uEMEP_calculate_emissions_for_EMEP proc~uemep_calculate_emissions_for_emep->proc~uemep_read_weekly_shipping_asi_data program~uemep uEMEP program~uemep->proc~uemep_read_weekly_shipping_asi_data program~uemep->proc~uemep_calculate_emissions_for_emep

Source Code

    subroutine uEMEP_read_weekly_shipping_asi_data

        use uEMEP_definitions

        implicit none

        integer i
        character(256) temp_str1
        integer unit_in
        logical :: exists
        real totalnoxemission,totalparticulatematteremission
        real y_ship,x_ship
        integer i_ship_index,j_ship_index
        integer source_index,subsource_index
        integer t
        integer, allocatable :: count_subgrid(:,:,:)

        integer i_pollutant

        integer a(6)
        character(256) format_temp,week_of_year_str
        double precision date_num,date_num_start
        integer week_of_year,ship_week,ship_counts
        logical nxtdat_flag
        real ship_delta_x,ship_delta_y
        real lat_ship,lon_ship

        write(unit_logfile,'(A)') ''
        write(unit_logfile,'(A)') '================================================================'
        write(unit_logfile,'(A)') 'Reading weekly shipping asi data  (uEMEP_read_weekly_shipping_asi_data)'
        write(unit_logfile,'(A)') '================================================================'

        source_index=shipping_index
        n_subsource(source_index)=1
        proxy_emission_subgrid(:,:,source_index,:)=0.
        t=1

        allocate (count_subgrid(emission_subgrid_dim(x_dim_index,source_index),emission_subgrid_dim(y_dim_index,source_index),n_pollutant_loop))
        count_subgrid=0

        !Determine week of year approximately using date string. Needs proper function
        format_temp='yyyymmdd'
        call datestr_to_date(config_date_str,format_temp,a)
        date_num=date_to_number(a,ref_year_meteo)
        a(2)=1;a(3)=1;a(4)=1;
        date_num_start=date_to_number(a,ref_year_meteo)
        week_of_year=1+int((date_num-date_num_start)/7.)
        week_of_year=max(min(week_of_year,52),1)
        !write(*,*) week_of_year
        write(week_of_year_str,'(i2)') week_of_year
        write(unit_logfile,'(a)') 'Week of year: '//trim(ADJUSTL(week_of_year_str))

        pathfilename_ship(1)=trim(pathname_ship(1))//trim(filename_ship(1))//'_week_'//trim(ADJUSTL(week_of_year_str))//'.txt'

        !Test existence of the shipping filename. If does not exist then use default
        inquire(file=trim(pathfilename_ship(1)),exist=exists)
        if (.not.exists) then
            write(unit_logfile,'(A,A)') ' ERROR: Shipping file does not exist: ', trim(pathfilename_ship(1))
            stop
        endif

        !Open the file for reading
        unit_in=20
        open(unit_in,file=pathfilename_ship(1),access='sequential',status='old',readonly)
        write(unit_logfile,'(a)') ' Opening shipping file '//trim(pathfilename_ship(1))

        rewind(unit_in)

        subsource_index=1

        !Skip over lines starting with #
        call nxtdat(unit_in,nxtdat_flag)

        !Read data
        read(unit_in,*) temp_str1,ship_delta_x
        read(unit_in,*) temp_str1,ship_delta_y
        read(unit_in,*) temp_str1,ship_week
        read(unit_in,*) temp_str1,ship_counts

        !Skip header
        read(unit_in,*) temp_str1
        do i=1,ship_counts
            read(unit_in,*) x_ship,y_ship,totalnoxemission,totalparticulatematteremission

            !Special case when saving emissions, convert to either latlon or lambert
            if (save_emissions_for_EMEP(shipping_index)) then
                call PROJ2LL(x_ship,y_ship,lon_ship,lat_ship,projection_attributes,projection_type)
                ! call utm2ll_modern(1, utm_zone,y_ship,x_ship,lat_ship,lon_ship)
                if (EMEP_projection_type.eq.LL_projection_index) then
                    x_ship=lon_ship
                    y_ship=lat_ship
                elseif (EMEP_projection_type.eq.LCC_projection_index) then
                    call lb2lambert2_uEMEP(x_ship,y_ship,lon_ship,lat_ship,EMEP_projection_attributes)
                elseif (EMEP_projection_type.eq.PS_projection_index) then
                    call LL2PS_spherical(x_ship,y_ship,lon_ship,lat_ship,EMEP_projection_attributes)
                endif
            endif

            i_ship_index=1+floor((x_ship-emission_subgrid_min(x_dim_index,source_index))/emission_subgrid_delta(x_dim_index,source_index))
            j_ship_index=1+floor((y_ship-emission_subgrid_min(y_dim_index,source_index))/emission_subgrid_delta(y_dim_index,source_index))

            !Add to subgrid
            if (i_ship_index.ge.1.and.i_ship_index.le.emission_subgrid_dim(x_dim_index,source_index) &
                .and.j_ship_index.ge.1.and.j_ship_index.le.emission_subgrid_dim(y_dim_index,source_index)) then
                do i_pollutant=1,n_pollutant_loop
                    if  (totalnoxemission.gt.0.and.pollutant_loop_index(i_pollutant).eq.nox_nc_index) then
                        proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)=proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)+totalnoxemission
                        count_subgrid(i_ship_index,j_ship_index,i_pollutant)=count_subgrid(i_ship_index,j_ship_index,i_pollutant)+1
                    elseif (totalparticulatematteremission.gt.0.and.(pollutant_loop_index(i_pollutant).eq.pm25_nc_index.or.pollutant_loop_index(i_pollutant).eq.pm10_nc_index)) then
                        proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)=proxy_emission_subgrid(i_ship_index,j_ship_index,source_index,i_pollutant)+totalparticulatematteremission
                        count_subgrid(i_ship_index,j_ship_index,i_pollutant)=count_subgrid(i_ship_index,j_ship_index,i_pollutant)+1
                    endif
                enddo
            endif

        enddo

        write(unit_logfile,'(A,I)') 'Shipping counts = ',ship_counts
        do i_pollutant=1,n_pollutant_loop
            write(unit_logfile,'(A,es12.3)') 'Total emission (g/hr) '//trim(pollutant_file_str(pollutant_loop_index(i_pollutant)))//' = ',sum(proxy_emission_subgrid(1:emission_subgrid_dim(x_dim_index,source_index),1:emission_subgrid_dim(y_dim_index,source_index),source_index,i_pollutant))
        enddo

        close(unit_in)

        deallocate (count_subgrid)

    end subroutine uEMEP_read_weekly_shipping_asi_data