LCOV - code coverage report
Current view: | top level - chemistry/utils - aircraft_emit.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 35 | 161 | 21.7 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 4 | 8 | 50.0 %

              Line data    Source code
           1             : module aircraft_emit
           2             : !-----------------------------------------------------------------------
           3             : !
           4             : ! Purpose:
           5             : ! Manages reading and interpolation of aircraft aerosols
           6             : !
           7             : ! Authors: Chih-Chieh (Jack) Chen and Cheryl Craig -- February 2010
           8             : !
           9             : !-----------------------------------------------------------------------
          10             :   use perf_mod,     only : t_startf, t_stopf
          11             :
          12             :   use shr_kind_mod,     only : r8 => shr_kind_r8
          13             :   use cam_abortutils,   only : endrun
          14             :   use spmd_utils,       only : masterproc
          15             :   use tracer_data,      only : trfld, trfile
          16             :   use cam_logfile,      only : iulog
          17             :
          18             :   implicit none
          19             :   private
          20             :   save
          21             :
          22             :   public :: aircraft_emit_init
          23             :   public :: aircraft_emit_adv
          24             :   public :: aircraft_emit_register
          25             :   public :: aircraft_emit_readnl
          26             :   public :: get_aircraft
          27             :
          28             :   type :: forcing_air
          29             :      real(r8)              :: mw
          30             :      character(len=256) :: filelist
          31             :      character(len=256) :: filename
          32             :      real(r8), pointer     :: times(:)
          33             :      real(r8), pointer     :: levi(:)
          34             :      character(len=11)  :: species
          35             :      character(len=8)  :: units
          36             :      integer                   :: nsectors
          37             :      character(len=32),pointer :: sectors(:)
          38             :      type(trfld),pointer       :: fields(:)
          39             :      type(trfile)              :: file
          40             :   end type forcing_air
          41             :
          42             :   type(forcing_air),allocatable :: forcings_air(:)
          43             :
          44             :   integer, parameter :: N_AERO = 13
          45             :   character(len=13)    :: aero_names(N_AERO) = (/'ac_SLANT_DIST','ac_TRACK_DIST','ac_HC        ','ac_NOX       ','ac_PMNV      ',&
          46             :                                                  'ac_PMSO      ','ac_PMFO      ','ac_FUELBURN  ','ac_CO2       ','ac_H2O       ',&
          47             :                                                  'ac_SOX       ','ac_CO        ','ac_BC        '/)
          48             :
          49             :   real(r8), parameter :: molmass(N_AERO) = 1._r8
          50             :
          51             :   logical :: advective_tracer(N_AERO) = .false.
          52             :   character(len=3) :: mixtype(N_AERO) = 'wet'
          53             :
          54             :   real(r8) :: cptmp = 666.0_r8
          55             :   real(r8) :: qmin = 0.0_r8
          56             :   logical :: cam_outfld = .false.
          57             :
          58             :   integer            :: index_map(N_AERO)
          59             :   character(len=256) :: air_specifier(N_AERO)=''
          60             :   character(len=256) :: air_datapath=''
          61             :   character(len=24)  :: air_type = 'CYCLICAL_LIST' ! 'CYCLICAL_LIST'
          62             :
          63             :   logical            :: rmv_file = .false.
          64             :
          65             :   integer :: number_flds
          66             :
          67             :   integer :: aircraft_cnt = 0
          68             :   character(len=16) :: spc_name_list(N_AERO)
          69             :   character(len=256) :: spc_flist(N_AERO),spc_fname(N_AERO)
          70             :   logical :: dist(N_AERO)
          71             :
          72             : contains
          73             :
          74           0 :   subroutine get_aircraft(cnt, spc_name_list_out)
          75             :   integer, intent(out) :: cnt
          76             :   character(len=16), optional, intent(out) :: spc_name_list_out(N_AERO)
          77             :   integer :: i
          78             :
          79           0 :   spc_name_list_out = ''
          80             :
          81           0 :   cnt = aircraft_cnt
          82           0 :   if( cnt>0 ) then
          83           0 :     do i=1,cnt
          84           0 :        spc_name_list_out(i) = spc_name_list(i)
          85             :     end do
          86             :   end if
          87             :
          88           0 :   end subroutine get_aircraft
          89             :
          90        1536 :   subroutine aircraft_emit_register()
          91             :
          92             : !------------------------------------------------------------------
          93             : ! **** Add the aircraft aerosol data to the physics buffer ****
          94             : !------------------------------------------------------------------
          95             :     use ppgrid,         only: pver,pcols
          96             :     use physics_buffer, only : pbuf_add_field, dtype_r8
          97             :     use tracer_data,    only: incr_filename
          98             :     use constituents,   only: cnst_add
          99             :
         100             :     integer :: i,idx, mm, ind, n
         101             :     character(len=16) :: spc_name
         102             :     character(len=256) :: filelist, curr_filename
         103             :     character(len=128) :: long_name
         104             :     logical            :: has_fixed_ubc=.false.
         105             :     logical            :: read_iv=.false.
         106             :
         107             :     !------------------------------------------------------------------
         108             :     ! Return if air_specifier is blank (no aircraft data to process)
         109             :     !------------------------------------------------------------------
         110        1536 :     dist(:) = .false.
         111        1536 :     aircraft_cnt = 0
         112        1536 :     if (air_specifier(1) == "") return
         113             :
         114             : ! count aircraft emission species used in the simulation
         115           0 :     count_emis: do n=1,N_AERO
         116             :
         117           0 :         if( len_trim(air_specifier(n) ) == 0 ) then
         118             :             exit count_emis
         119             :         endif
         120             :
         121           0 :         i = scan(air_specifier(n),'->')
         122           0 :         spc_name = trim(adjustl(air_specifier(n)(:i-1)))
         123           0 :         filelist = trim(adjustl(air_specifier(n)(i+2:)))
         124             :
         125           0 :         mm = get_aircraft_ndx(spc_name)
         126           0 :         if( mm < 1 ) then
         127           0 :          call endrun('aircraft_emit_register: '//trim(spc_name)//' is not in the aircraft emission dataset')
         128             :         endif
         129             :
         130           0 :         if (trim(spc_name) == 'ac_SLANT_DIST'.or. trim(spc_name) == 'ac_TRACK_DIST') dist(n) = .true.
         131             :
         132           0 :         aircraft_cnt = aircraft_cnt + 1
         133           0 :         call pbuf_add_field(aero_names(mm),'physpkg',dtype_r8,(/pcols,pver/),idx)
         134             :
         135           0 :         spc_flist(aircraft_cnt) = filelist
         136           0 :         spc_name_list(aircraft_cnt) = spc_name
         137           0 :         index_map(aircraft_cnt) = mm
         138             :
         139           0 :         curr_filename=''
         140             :         spc_fname(aircraft_cnt) = incr_filename( curr_filename, filenames_list=spc_flist(aircraft_cnt), &
         141           0 :                                                  datapath=air_datapath)
         142             :
         143           0 :         if( advective_tracer(mm) ) then
         144           0 :           long_name = 'aircraft_'//trim(spc_name)
         145             :           call cnst_add(aero_names(mm),molmass(mm),cptmp,qmin,ind,longname=long_name,readiv=read_iv, &
         146           0 :                         mixtype=mixtype(mm),cam_outfld=cam_outfld,fixed_ubc=has_fixed_ubc)
         147             :         endif
         148             :
         149             :     enddo count_emis
         150             : ! count aircraft emission species used in the simulation
         151             :
         152        1536 :   endsubroutine aircraft_emit_register
         153             :
         154        1536 :   subroutine aircraft_emit_init()
         155             : !-------------------------------------------------------------------
         156             : ! **** Initialize the aircraft aerosol data handling ****
         157             : !-------------------------------------------------------------------
         158        1536 :     use cam_history,    only: addfld, add_default
         159             :     use tracer_data,    only: trcdata_init
         160             :     use phys_control,   only: phys_getopts
         161             :
         162             :     implicit none
         163             :
         164             :     character(len=16)  :: spc_name
         165             :
         166             :     integer :: astat, m
         167             :
         168             :     logical :: history_chemistry
         169             :
         170        1536 :     call phys_getopts(history_chemistry_out=history_chemistry)
         171             :
         172             :     !------------------------------------------------------------------
         173             :     ! Return if aircraft_cnt is zero (no aircraft data to process)
         174             :     !------------------------------------------------------------------
         175        1536 :     if (aircraft_cnt == 0 ) return
         176             :
         177           0 :     if (masterproc) write(iulog,*) ' '
         178             :
         179             :     !-----------------------------------------------------------------------
         180             :     !       allocate forcings type array
         181             :     !-----------------------------------------------------------------------
         182           0 :     allocate( forcings_air(aircraft_cnt), stat=astat )
         183           0 :     if( astat/= 0 ) then
         184           0 :        write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air array; error = ',astat
         185           0 :        call endrun
         186             :     end if
         187             :
         188             :     !-----------------------------------------------------------------------
         189             :     !       setup the forcings_air type array
         190             :     !-----------------------------------------------------------------------
         191           0 :     species_loop : do m = 1,aircraft_cnt
         192             :
         193           0 :           allocate( forcings_air(m)%sectors(1), stat=astat )
         194           0 :           if( astat/= 0 ) then
         195           0 :              write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air%sectors array; error = ',astat
         196           0 :              call endrun
         197             :           end if
         198             :
         199           0 :           allocate( forcings_air(m)%fields(1), stat=astat )
         200           0 :           if( astat/= 0 ) then
         201           0 :              write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air%fields array; error = ',astat
         202           0 :              call endrun
         203             :           end if
         204             :
         205           0 :           spc_name = spc_name_list(m)
         206             :           !-----------------------------------------------------------------------
         207             :           !         default settings
         208             :           !-----------------------------------------------------------------------
         209           0 :           forcings_air(m)%file%stepTime    = .true.  ! Aircraft data is not to be interpolated in time
         210           0 :           forcings_air(m)%file%cyclical_list    = .true.  ! Aircraft data cycles over the filename list
         211           0 :           forcings_air(m)%file%weight_by_lat     = .true.  ! Aircraft data -  interpolated with latitude weighting
         212           0 :           forcings_air(m)%file%conserve_column = .true. ! Aircraft data - vertically interpolated to conserve the total column
         213           0 :           forcings_air(m)%file%dist = dist(m)
         214           0 :           forcings_air(m)%species          = spc_name
         215           0 :           forcings_air(m)%sectors          = spc_name ! Only one species per file for aircraft data
         216           0 :           forcings_air(m)%nsectors         = 1
         217           0 :           forcings_air(m)%filelist         = spc_flist(m)
         218             : !         forcings_air(m)%file%curr_filename    = spc_fname(m)
         219           0 :           forcings_air(m)%filename         = spc_fname(m)
         220             :     end do species_loop
         221             :
         222           0 :     if (masterproc) then
         223             :        !-----------------------------------------------------------------------
         224             :        !            diagnostics
         225             :        !-----------------------------------------------------------------------
         226           0 :        write(iulog,*) ' '
         227           0 :        write(iulog,*) 'aircraft_emit_init: diagnostics'
         228           0 :        write(iulog,*) ' '
         229           0 :        write(iulog,*) 'aircraft_emit timing specs'
         230           0 :        write(iulog,*) 'type = ',air_type
         231           0 :        write(iulog,*) ' '
         232           0 :        write(iulog,*) 'there are ',aircraft_cnt,' species of aircraft emission'
         233           0 :        do m = 1,aircraft_cnt
         234           0 :           write(iulog,*) ' '
         235           0 :           write(iulog,*) 'forcing type ',m
         236           0 :           write(iulog,*) 'species = ',trim(forcings_air(m)%species)
         237           0 :           write(iulog,*) 'filelist= ',trim(forcings_air(m)%filelist)
         238             :        end do
         239           0 :        write(iulog,*) ' '
         240             :     endif
         241             :
         242             :     !------------------------------------------------------------------
         243             :     !       Initialize the aircraft file processing
         244             :     !------------------------------------------------------------------
         245           0 :     do m=1,aircraft_cnt
         246             :
         247           0 :        allocate (forcings_air(m)%file%in_pbuf(size(forcings_air(m)%sectors)))
         248           0 :        forcings_air(m)%file%in_pbuf(:) = .true.
         249           0 :        call trcdata_init( forcings_air(m)%sectors, forcings_air(m)%filename, forcings_air(m)%filelist, air_datapath, &
         250           0 :                           forcings_air(m)%fields, forcings_air(m)%file, rmv_file, 0, 0, 0, air_type)
         251             :
         252             :
         253           0 :        number_flds = 0
         254           0 :        if (associated(forcings_air(m)%fields)) number_flds = size( forcings_air(m)%fields )
         255             :
         256           0 :        if( number_flds < 1 ) then
         257           0 :           if ( masterproc ) then
         258           0 :              write(iulog,*) 'There are no aircraft aerosols'
         259           0 :              write(iulog,*) ' '
         260           0 :              call endrun
         261             :           endif
         262             :        end if
         263             :
         264           0 :        spc_name = spc_name_list(m)
         265           0 :        call addfld( trim(spc_name), (/ 'lev' /), 'A', forcings_air(m)%fields(1)%units, &
         266           0 :                     'aircraft emission '//trim(spc_name) )
         267           0 :        if (history_chemistry) then
         268           0 :           call add_default( trim(spc_name), 1, ' ' )
         269             :        end if
         270             :    end do
         271             :
         272             :
         273        3072 :   end subroutine aircraft_emit_init
         274             :
         275             :
         276             :
         277      741888 :   subroutine aircraft_emit_adv( state, pbuf2d)
         278             : !-------------------------------------------------------------------
         279             : ! **** Advance to the next aircraft data ****
         280             : !-------------------------------------------------------------------
         281             :
         282        1536 :     use tracer_data,  only : advance_trcdata
         283             :     use physics_types,only : physics_state
         284             :     use ppgrid,       only : begchunk, endchunk
         285             :     use ppgrid,       only : pcols, pver
         286             :     use string_utils, only : to_lower, GLC
         287             :     use cam_history,  only : outfld
         288             :     use physconst,    only : mwdry       ! molecular weight dry air ~ kg/kmole
         289             :     use physconst,    only : boltz                ! J/K/molecule
         290             : ! C.-C. Chen
         291             : !    use physconst,    only : gravit, rearth
         292             :     use phys_grid,    only : get_wght_all_p
         293             :
         294             :     use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk
         295             :
         296             :     implicit none
         297             :
         298             :     type(physics_state), intent(in)    :: state(begchunk:endchunk)
         299             :
         300             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         301             :
         302      370944 :     type(physics_buffer_desc), pointer :: pbuf_chnk(:)
         303             :     integer :: ind,c,ncol,i,caseid,m,n
         304             :     real(r8) :: to_mmr(pcols,pver)
         305      370944 :     real(r8),pointer :: tmpptr(:,:)
         306             :
         307             : ! C.-C. Chen
         308             :     real(r8) :: wght(pcols)
         309             :
         310             :    !------------------------------------------------------------------
         311             :    ! Return if aircraft_cnt is zero (no aircraft data to process)
         312             :    !------------------------------------------------------------------
         313      370944 :     if (aircraft_cnt == 0 ) return
         314           0 :     call t_startf('All_aircraft_emit_adv')
         315             :
         316             :    !-------------------------------------------------------------------
         317             :    !    For each field, read more data if needed and interpolate it to the current model time
         318             :    !-------------------------------------------------------------------
         319           0 :     do m = 1, aircraft_cnt
         320           0 :        call advance_trcdata( forcings_air(m)%fields, forcings_air(m)%file, state, pbuf2d)
         321             :
         322             :    !-------------------------------------------------------------------
         323             :    !    set the tracer fields with the correct units
         324             :    !-------------------------------------------------------------------
         325           0 :        do i = 1,number_flds
         326             :
         327             : ! C.-C. Chen, adding case 4  for kg/sec
         328           0 :           select case ( to_lower(trim(forcings_air(m)%fields(i)%units(:GLC(forcings_air(m)%fields(i)%units)))) )
         329             :           case ("molec/cm3","/cm3","molecules/cm3","cm^-3","cm**-3")
         330           0 :              caseid = 1
         331             :           case ('kg/kg','mmr')
         332           0 :              caseid = 2
         333             :           case ('mol/mol','mole/mole','vmr','fraction')
         334           0 :              caseid = 3
         335             :           case ('kg/kg/sec')
         336           0 :               caseid = 4
         337             :           case ('kg m-2 s-1')
         338           0 :               caseid = 5
         339             :           case ('m/sec' )
         340           0 :               caseid = 6
         341             :           case default
         342           0 :              print*, 'aircraft_emit_adv: units = ',trim(forcings_air(m)%fields(i)%units) ,' are not recognized'
         343           0 :              call endrun('aircraft_emit_adv: units are not recognized')
         344             :           end select
         345             :
         346           0 :           ind = index_map(i)
         347             :
         348             : !$OMP PARALLEL DO PRIVATE (C, NCOL, TO_MMR, tmpptr, pbuf_chnk)
         349           0 :           do c = begchunk,endchunk
         350           0 :              ncol = state(c)%ncol
         351             :
         352             : ! C.-C. Chen, turning emission data to mixing ratio
         353           0 :              call get_wght_all_p(c,ncol,wght(:ncol))
         354             :
         355           0 :              if (caseid == 1) then
         356           0 :                 to_mmr(:ncol,:) = (molmass(ind)*1.e6_r8*boltz*state(c)%t(:ncol,:))/(mwdry*state(c)%pmiddry(:ncol,:))
         357           0 :              elseif (caseid == 2) then
         358           0 :                 to_mmr(:ncol,:) = 1._r8
         359           0 :              elseif (caseid == 4) then
         360             : !                do n=1,ncol
         361             : !                  to_mmr(n,:) = 1.0_r8/(rearth*rearth*wght(n)*state(c)%pdel(n,:)/gravit)
         362             : !                end do
         363           0 :                 to_mmr(:ncol,:) = 1.0_r8
         364           0 :              elseif (caseid == 5) then
         365           0 :                 to_mmr(:ncol,:) = 1.0_r8
         366           0 :              elseif (caseid == 6) then
         367           0 :                 to_mmr(:ncol,:) = 1.0_r8
         368             :              else
         369           0 :                 to_mmr(:ncol,:) = molmass(ind)/mwdry
         370             :              endif
         371           0 :              pbuf_chnk => pbuf_get_chunk(pbuf2d, c)
         372           0 :              call pbuf_get_field(pbuf_chnk, forcings_air(m)%fields(i)%pbuf_ndx, tmpptr )
         373             :
         374           0 :              tmpptr(:ncol,:) = tmpptr(:ncol,:)*to_mmr(:ncol,:)
         375             :
         376           0 :              call outfld( forcings_air(m)%fields(i)%fldnam, &
         377           0 :                   tmpptr(:ncol,:), ncol, state(c)%lchnk )
         378             :
         379             :           enddo
         380             :        enddo
         381             :     enddo
         382             :
         383           0 :     call t_stopf('All_aircraft_emit_adv')
         384      370944 :   end subroutine aircraft_emit_adv
         385             :
         386        1536 :   subroutine aircraft_emit_readnl(nlfile)
         387             : !-------------------------------------------------------------------
         388             : ! **** Read in the aircraft_emit namelist *****
         389             : !-------------------------------------------------------------------
         390      370944 :    use namelist_utils,  only: find_group_name
         391             :    use units,           only: getunit, freeunit
         392             :    use mpishorthand
         393             :
         394             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         395             :
         396             :    ! Local variables
         397             :    integer :: unitn, ierr
         398             :    character(len=*), parameter :: subname = 'aircraft_emit_readnl'
         399             :
         400             :    character(len=256) :: aircraft_specifier(N_AERO)
         401             :    character(len=256) :: aircraft_datapath
         402             :    character(len=24)  :: aircraft_type
         403             :
         404             :    namelist /aircraft_emit_nl/  aircraft_specifier, aircraft_datapath, aircraft_type
         405             :    !-----------------------------------------------------------------------------
         406             :
         407             :    ! Initialize namelist variables from local module variables.
         408       21504 :    aircraft_specifier= air_specifier
         409        1536 :    aircraft_datapath = air_datapath
         410        1536 :    aircraft_type     = air_type
         411             :
         412             :    ! Read namelist
         413        1536 :    if (masterproc) then
         414           2 :       unitn = getunit()
         415           2 :       open( unitn, file=trim(nlfile), status='old' )
         416           2 :       call find_group_name(unitn, 'aircraft_emit_nl', status=ierr)
         417           2 :       if (ierr == 0) then
         418           0 :          read(unitn, aircraft_emit_nl, iostat=ierr)
         419           0 :          if (ierr /= 0) then
         420           0 :             call endrun(subname // ':: ERROR reading namelist')
         421             :          end if
         422             :       end if
         423           2 :       close(unitn)
         424           2 :       call freeunit(unitn)
         425             :    end if
         426             :
         427             : #ifdef SPMD
         428             :    ! Broadcast namelist variables
         429        1536 :    call mpibcast(aircraft_specifier,len(aircraft_specifier(1))*N_AERO,     mpichar, 0, mpicom)
         430        1536 :    call mpibcast(aircraft_datapath, len(aircraft_datapath),                mpichar, 0, mpicom)
         431        1536 :    call mpibcast(aircraft_type,     len(aircraft_type),                    mpichar, 0, mpicom)
         432             : #endif
         433             :
         434             :    ! Update module variables with user settings.
         435       21504 :    air_specifier  = aircraft_specifier
         436        1536 :    air_datapath   = aircraft_datapath
         437        1536 :    air_type       = aircraft_type
         438             :
         439        1536 :  end subroutine aircraft_emit_readnl
         440             :
         441           0 :  integer function get_aircraft_ndx( name )
         442             :
         443             :     implicit none
         444             :     character(len=*), intent(in) :: name
         445             :
         446             :     integer :: i
         447             :
         448           0 :     get_aircraft_ndx = 0
         449           0 :     do i = 1,N_AERO
         450           0 :       if ( trim(name) == trim(aero_names(i)) ) then
         451           0 :         get_aircraft_ndx = i
         452           0 :         return
         453             :       endif
         454             :     enddo
         455             :
         456             :   end function get_aircraft_ndx
         457             :
         458           0 : end module aircraft_emit

Generated by: LCOV version 1.14