LCOV - code coverage report
Current view: | top level - chemistry/utils - prescribed_strataero.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 58 | 218 | 26.6 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 4 | 7 | 57.1 %

              Line data    Source code
           1             : !-------------------------------------------------------------------
           2             : ! manages reading and interpolation of prescribed stratospheric aerosols
           3             : ! Created by: Francis Vitt
           4             : !-------------------------------------------------------------------
           5             : module prescribed_strataero
           6             :
           7             :   use shr_kind_mod,     only : r8 => shr_kind_r8
           8             :   use cam_abortutils,   only : endrun
           9             :   use spmd_utils,       only : masterproc
          10             :   use tracer_data,      only : trfld, trfile
          11             :   use cam_logfile,      only : iulog
          12             :
          13             :   implicit none
          14             :   private
          15             :   save
          16             :
          17             :   type(trfld), pointer :: fields(:)
          18             :   type(trfile)         :: file
          19             :
          20             :   public :: prescribed_strataero_readnl
          21             :   public :: prescribed_strataero_register
          22             :   public :: prescribed_strataero_init
          23             :   public :: prescribed_strataero_adv
          24             :   public :: write_prescribed_strataero_restart
          25             :   public :: read_prescribed_strataero_restart
          26             :   public :: has_prescribed_strataero
          27             :   public :: init_prescribed_strataero_restart
          28             :
          29             :   logical :: has_prescribed_strataero = .false.
          30             :   character(len=16), parameter :: mmr_name = 'VOLC_MMR'
          31             :   character(len=16), parameter :: rad_name = 'VOLC_RAD_GEOM'
          32             :   character(len=16), parameter :: sad_name = 'VOLC_SAD'
          33             :   character(len=16), parameter :: mass_name = 'VOLC_MASS'
          34             :   character(len=16), parameter :: mass_column_name = 'VOLC_MASS_C'
          35             :   character(len=16), parameter :: dens_name = 'VOLC_DENS'
          36             :
          37             :   character(len=16), parameter :: mmr_name1 = 'VOLC_MMR1'
          38             :   character(len=16), parameter :: mmr_name2 = 'VOLC_MMR2'
          39             :   character(len=16), parameter :: mmr_name3 = 'VOLC_MMR3'
          40             :   character(len=16), parameter :: rad_name1 = 'VOLC_RAD_GEOM1'
          41             :   character(len=16), parameter :: rad_name2 = 'VOLC_RAD_GEOM2'
          42             :   character(len=16), parameter :: rad_name3 = 'VOLC_RAD_GEOM3'
          43             :   character(len=16), parameter :: mass_name1 = 'VOLC_MASS1'
          44             :   character(len=16), parameter :: mass_name2 = 'VOLC_MASS2'
          45             :   character(len=16), parameter :: mass_name3 = 'VOLC_MASS3'
          46             :   character(len=16), parameter :: mass_column_name1 = 'VOLC_MASS_C1'
          47             :   character(len=16), parameter :: mass_column_name2 = 'VOLC_MASS_C2'
          48             :   character(len=16), parameter :: mass_column_name3 = 'VOLC_MASS_C3'
          49             :   character(len=16), parameter :: dens_name1 = 'VOLC_DENS1'
          50             :   character(len=16), parameter :: dens_name2 = 'VOLC_DENS2'
          51             :   character(len=16), parameter :: dens_name3 = 'VOLC_DENS3'
          52             :
          53             :   ! These variables are settable via the namelist (with longer names)
          54             :   character(len=32)  :: specifier(7) = ' '
          55             :   character(len=256) :: filename = 'NONE'
          56             :   character(len=256) :: filelist = ''
          57             :   character(len=256) :: datapath = ''
          58             :   character(len=32)  :: data_type = 'SERIAL'
          59             :   logical            :: rmv_file = .false.
          60             :   integer            :: cycle_yr  = 0
          61             :   integer            :: fixed_ymd = 0
          62             :   integer            :: fixed_tod = 0
          63             :   integer            :: rad_ndx1 = -1
          64             :   integer            :: rad_ndx2 = -1
          65             :   integer            :: rad_ndx3 = -1
          66             :   integer            :: sad_ndx = -1
          67             :   integer            :: mmr_ndx1 = -1
          68             :   integer            :: mmr_ndx2 = -1
          69             :   integer            :: mmr_ndx3 = -1
          70             :
          71             :   logical            :: prescribed_strataero_use_chemtrop = .false.
          72             :   logical            :: three_mode = .true.
          73             :   integer :: rad_fld_no=-1, sad_fld_no=-1
          74             :
          75             : contains
          76             :
          77             : !-------------------------------------------------------------------
          78             : !-------------------------------------------------------------------
          79        1536 : subroutine prescribed_strataero_readnl(nlfile)
          80             :
          81             :    use namelist_utils,  only: find_group_name
          82             :    use units,           only: getunit, freeunit
          83             :    use mpishorthand
          84             :
          85             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          86             :
          87             :    ! Local variables
          88             :    integer :: unitn, ierr
          89             :    character(len=*), parameter :: subname = 'prescribed_strataero_readnl'
          90             :
          91             :    character(len=32)  :: prescribed_strataero_specifier(7)
          92             :    character(len=256) :: prescribed_strataero_file
          93             :    character(len=256) :: prescribed_strataero_filelist
          94             :    character(len=256) :: prescribed_strataero_datapath
          95             :    character(len=32)  :: prescribed_strataero_type
          96             :    logical            :: prescribed_strataero_rmfile
          97             :    integer            :: prescribed_strataero_cycle_yr
          98             :    integer            :: prescribed_strataero_fixed_ymd
          99             :    integer            :: prescribed_strataero_fixed_tod
         100             :
         101             :    namelist /prescribed_strataero_nl/ &
         102             :       prescribed_strataero_specifier, &
         103             :       prescribed_strataero_file,      &
         104             :       prescribed_strataero_filelist,  &
         105             :       prescribed_strataero_datapath,  &
         106             :       prescribed_strataero_type,      &
         107             :       prescribed_strataero_rmfile,    &
         108             :       prescribed_strataero_cycle_yr,  &
         109             :       prescribed_strataero_fixed_ymd, &
         110             :       prescribed_strataero_fixed_tod, &
         111             :       prescribed_strataero_use_chemtrop
         112             :    !-----------------------------------------------------------------------------
         113             :
         114             :    ! Initialize namelist variables from local module variables.
         115       12288 :    prescribed_strataero_specifier= specifier
         116        1536 :    prescribed_strataero_file     = filename
         117        1536 :    prescribed_strataero_filelist = filelist
         118        1536 :    prescribed_strataero_datapath = datapath
         119        1536 :    prescribed_strataero_type     = data_type
         120        1536 :    prescribed_strataero_rmfile   = rmv_file
         121        1536 :    prescribed_strataero_cycle_yr = cycle_yr
         122        1536 :    prescribed_strataero_fixed_ymd= fixed_ymd
         123        1536 :    prescribed_strataero_fixed_tod= fixed_tod
         124             :
         125             :    ! Read namelist
         126        1536 :    if (masterproc) then
         127           2 :       unitn = getunit()
         128           2 :       open( unitn, file=trim(nlfile), status='old' )
         129           2 :       call find_group_name(unitn, 'prescribed_strataero_nl', status=ierr)
         130           2 :       if (ierr == 0) then
         131           0 :          read(unitn, prescribed_strataero_nl, iostat=ierr)
         132           0 :          if (ierr /= 0) then
         133           0 :             call endrun(subname // ':: ERROR reading namelist')
         134             :          end if
         135             :       end if
         136           2 :       close(unitn)
         137           2 :       call freeunit(unitn)
         138             :    end if
         139             :
         140             : #ifdef SPMD
         141             :    ! Broadcast namelist variables
         142        1536 :    call mpibcast(prescribed_strataero_specifier,len(prescribed_strataero_specifier)*7, mpichar, 0, mpicom)
         143        1536 :    call mpibcast(prescribed_strataero_file,     len(prescribed_strataero_file),        mpichar, 0, mpicom)
         144        1536 :    call mpibcast(prescribed_strataero_filelist, len(prescribed_strataero_filelist),    mpichar, 0, mpicom)
         145        1536 :    call mpibcast(prescribed_strataero_datapath, len(prescribed_strataero_datapath),    mpichar, 0, mpicom)
         146        1536 :    call mpibcast(prescribed_strataero_type,     len(prescribed_strataero_type),        mpichar, 0, mpicom)
         147        1536 :    call mpibcast(prescribed_strataero_rmfile,   1, mpilog,  0, mpicom)
         148        1536 :    call mpibcast(prescribed_strataero_cycle_yr, 1, mpiint,  0, mpicom)
         149        1536 :    call mpibcast(prescribed_strataero_fixed_ymd,1, mpiint,  0, mpicom)
         150        1536 :    call mpibcast(prescribed_strataero_fixed_tod,1, mpiint,  0, mpicom)
         151        1536 :    call mpibcast(prescribed_strataero_use_chemtrop, 1, mpilog,  0, mpicom)
         152             : #endif
         153             :
         154             :    ! Update module variables with user settings.
         155       12288 :    specifier(:) = prescribed_strataero_specifier(:)
         156        1536 :    filename   = prescribed_strataero_file
         157        1536 :    filelist   = prescribed_strataero_filelist
         158        1536 :    datapath   = prescribed_strataero_datapath
         159        1536 :    data_type  = prescribed_strataero_type
         160        1536 :    rmv_file   = prescribed_strataero_rmfile
         161        1536 :    cycle_yr   = prescribed_strataero_cycle_yr
         162        1536 :    fixed_ymd  = prescribed_strataero_fixed_ymd
         163        1536 :    fixed_tod  = prescribed_strataero_fixed_tod
         164             :
         165             :    ! Turn on prescribed volcanics if user has specified an input dataset.
         166        1536 :    if (len_trim(filename) > 0 .and. filename.ne.'NONE') has_prescribed_strataero = .true.
         167             :
         168        1536 : end subroutine prescribed_strataero_readnl
         169             :
         170             : !-------------------------------------------------------------------
         171             : !-------------------------------------------------------------------
         172        1536 :   subroutine prescribed_strataero_register()
         173             :     use ppgrid,         only: pver,pcols
         174             :     use physics_buffer, only: pbuf_add_field, dtype_r8
         175             :     use pio,            only: var_desc_t, file_desc_t, pio_closefile, pio_inq_varid, pio_seterrorhandling, &
         176             :                               PIO_INTERNAL_ERROR, PIO_BCAST_ERROR, PIO_NOERR
         177             :     use cam_pio_utils,  only: cam_pio_openfile
         178             :     use ioFileMod, only : getfil
         179             :
         180             :     type(var_desc_t)  :: varid
         181             :     type(file_desc_t) :: file_handle
         182             :     character(len=256) :: filepath, filen
         183             :     integer :: ierr
         184             :
         185        1536 :     if (has_prescribed_strataero) then
         186             :
         187           0 :        filepath = trim(datapath)//'/'//trim(filename)
         188             :
         189           0 :        call getfil( filepath, filen, 0 )
         190           0 :        call cam_pio_openfile( file_handle, filen, 0 )
         191             :
         192           0 :        call pio_seterrorhandling(file_handle, PIO_BCAST_ERROR)
         193             :
         194           0 :        ierr = pio_inq_varid( file_handle, 'so4mass_a1', varid )
         195           0 :        three_mode = three_mode .and. (ierr.eq.PIO_NOERR)
         196           0 :        ierr = pio_inq_varid( file_handle, 'so4mass_a2', varid )
         197           0 :        three_mode = three_mode .and. (ierr.eq.PIO_NOERR)
         198           0 :        ierr = pio_inq_varid( file_handle, 'so4mass_a3', varid )
         199           0 :        three_mode = three_mode .and. (ierr.eq.PIO_NOERR)
         200           0 :        ierr = pio_inq_varid( file_handle, 'diamwet_a1', varid )
         201           0 :        three_mode = three_mode .and. (ierr.eq.PIO_NOERR)
         202           0 :        ierr = pio_inq_varid( file_handle, 'diamwet_a2', varid )
         203           0 :        three_mode = three_mode .and. (ierr.eq.PIO_NOERR)
         204           0 :        ierr = pio_inq_varid( file_handle, 'diamwet_a3', varid )
         205           0 :        three_mode = three_mode .and. (ierr.eq.PIO_NOERR)
         206             :
         207           0 :        call pio_seterrorhandling(file_handle, PIO_INTERNAL_ERROR)
         208             :
         209           0 :        call pio_closefile( file_handle )
         210             :
         211           0 :        if (three_mode) then
         212           0 :           call pbuf_add_field(mmr_name1, 'physpkg', dtype_r8,(/pcols,pver/), mmr_ndx1)
         213           0 :           call pbuf_add_field(mmr_name2, 'physpkg', dtype_r8,(/pcols,pver/), mmr_ndx2)
         214           0 :           call pbuf_add_field(mmr_name3, 'physpkg', dtype_r8,(/pcols,pver/), mmr_ndx3)
         215           0 :           call pbuf_add_field(rad_name1, 'physpkg', dtype_r8,(/pcols,pver/), rad_ndx1)
         216           0 :           call pbuf_add_field(rad_name2, 'physpkg', dtype_r8,(/pcols,pver/), rad_ndx2)
         217           0 :           call pbuf_add_field(rad_name3, 'physpkg', dtype_r8,(/pcols,pver/), rad_ndx3)
         218           0 :           call pbuf_add_field(sad_name, 'physpkg', dtype_r8,(/pcols,pver/), sad_ndx)
         219             :           specifier(1:7) = (/'VOLC_MMR1:so4mass_a1            ', &
         220             :                              'VOLC_MMR2:so4mass_a2            ', &
         221             :                              'VOLC_MMR3:so4mass_a3            ', &
         222             :                              'VOLC_RAD_GEOM1:diamwet_a1       ', &
         223             :                              'VOLC_RAD_GEOM2:diamwet_a2       ', &
         224             :                              'VOLC_RAD_GEOM3:diamwet_a3       ', &
         225           0 :                              'VOLC_SAD:SAD_AERO               ' /)
         226           0 :           rad_fld_no = 4
         227           0 :           sad_fld_no = 7
         228             :        else
         229           0 :           if (masterproc) then
         230           0 :              write(iulog, *) ' pbuf add mmr_name = '//trim(mmr_name)
         231             :           end if
         232           0 :           call pbuf_add_field(mmr_name, 'physpkg', dtype_r8,(/pcols,pver/), mmr_ndx1)
         233           0 :           call pbuf_add_field(rad_name, 'physpkg', dtype_r8,(/pcols,pver/), rad_ndx1)
         234           0 :           call pbuf_add_field(sad_name, 'physpkg', dtype_r8,(/pcols,pver/), sad_ndx)
         235             :           specifier(1:3) = (/'VOLC_MMR:H2SO4_mass             ', &
         236             :                              'VOLC_RAD_GEOM:rmode             ', &
         237           0 :                              'VOLC_SAD:sad                    ' /)
         238           0 :           rad_fld_no = 2
         239           0 :           sad_fld_no = 3
         240             :        endif
         241             :     endif
         242             :
         243        1536 :   endsubroutine prescribed_strataero_register
         244             :
         245             : !-------------------------------------------------------------------
         246             : !-------------------------------------------------------------------
         247        1536 :   subroutine prescribed_strataero_init()
         248             :
         249        1536 :     use tracer_data, only : trcdata_init
         250             :     use cam_history, only : addfld, horiz_only
         251             :     use error_messages, only: handle_err
         252             :
         253        1536 :     if ( has_prescribed_strataero ) then
         254           0 :        if ( masterproc ) then
         255           0 :           write(iulog,*) 'stratospheric aerosol is prescribed in :'//trim(filename)
         256             :        endif
         257             :     else
         258             :        return
         259             :     endif
         260             :
         261           0 :     allocate(file%in_pbuf(size(specifier)))
         262           0 :     file%in_pbuf(:) = .true.
         263           0 :     file%geop_alt = .true.
         264             :
         265             :     call trcdata_init( specifier, filename, filelist, datapath, fields, file, &
         266           0 :                        rmv_file, cycle_yr, fixed_ymd, fixed_tod, data_type)
         267             :
         268           0 :     if (three_mode) then
         269           0 :        call addfld(dens_name1, (/ 'lev' /), 'I','molecules/cm3', 'prescribed volcanic aerosol number density in Mode 1' )
         270           0 :        call addfld(dens_name2, (/ 'lev' /), 'I','molecules/cm3', 'prescribed volcanic aerosol number density in Mode 2' )
         271           0 :        call addfld(dens_name3, (/ 'lev' /), 'I','molecules/cm3', 'prescribed volcanic aerosol number density in Mode 3' )
         272           0 :        call addfld(mmr_name1, (/ 'lev' /), 'I','kg/kg', 'prescribed volcanic aerosol dry mass mixing ratio in Mode 1' )
         273           0 :        call addfld(mmr_name2, (/ 'lev' /), 'I','kg/kg', 'prescribed volcanic aerosol dry mass mixing ratio in Mode 2' )
         274           0 :        call addfld(mmr_name3, (/ 'lev' /), 'I','kg/kg', 'prescribed volcanic aerosol dry mass mixing ratio in Mode 3' )
         275           0 :        call addfld(rad_name1, (/ 'lev' /), 'I','m', 'volcanic aerosol geometric-mode radius in Mode 1' )
         276           0 :        call addfld(rad_name2, (/ 'lev' /), 'I','m', 'volcanic aerosol geometric-mode radius in Mode 2' )
         277           0 :        call addfld(rad_name3, (/ 'lev' /), 'I','m', 'volcanic aerosol geometric-mode radius in Mode 3' )
         278           0 :        call addfld(mass_name1, (/ 'lev' /), 'I','kg/m^2', 'volcanic aerosol vertical mass path in layer in Mode 1' )
         279           0 :        call addfld(mass_name2, (/ 'lev' /), 'I','kg/m^2', 'volcanic aerosol vertical mass path in layer in Mode 2' )
         280           0 :        call addfld(mass_name3, (/ 'lev' /), 'I','kg/m^2', 'volcanic aerosol vertical mass path in layer in Mode 3' )
         281           0 :        call addfld(mass_column_name1, horiz_only, 'I','kg/m^2', 'volcanic aerosol column mass in Mode 1' )
         282           0 :        call addfld(mass_column_name2, horiz_only, 'I','kg/m^2', 'volcanic aerosol column mass in Mode 2' )
         283           0 :        call addfld(mass_column_name3, horiz_only, 'I','kg/m^2', 'volcanic aerosol column mass IN Mode 3' )
         284             :     else
         285           0 :        call addfld(dens_name, (/ 'lev' /), 'I','molecules/cm3', 'prescribed volcanic aerosol number density' )
         286           0 :        call addfld(mmr_name,  (/ 'lev' /), 'I','kg/kg', 'prescribed volcanic aerosol dry mass mixing ratio' )
         287           0 :        call addfld(rad_name,  (/ 'lev' /), 'I','m', 'volcanic aerosol geometric-mode radius' )
         288           0 :        call addfld(mass_name, (/ 'lev' /), 'I','kg/m^2', 'volcanic aerosol vertical mass path in layer' )
         289           0 :        call addfld(mass_column_name, horiz_only, 'I','kg/m^2', 'volcanic aerosol column mass' )
         290             :     endif
         291           0 :     call addfld(sad_name, (/ 'lev' /), 'I','cm2/cm3', 'stratospheric aerosol surface area density' )
         292             :
         293        1536 :   end subroutine prescribed_strataero_init
         294             :
         295             : !-------------------------------------------------------------------
         296             : !-------------------------------------------------------------------
         297      741888 :   subroutine prescribed_strataero_adv( state, pbuf2d)
         298             :
         299        1536 :     use tracer_data,  only : advance_trcdata
         300             :     use physics_types,only : physics_state
         301             :     use ppgrid,       only : begchunk, endchunk
         302             :     use ppgrid,       only : pcols, pver
         303             :     use string_utils, only : to_lower, GLC
         304             :     use cam_history,  only : outfld
         305             :     use physconst,    only : mwdry                ! molecular weight dry air ~ kg/kmole
         306             :     use physconst,    only : boltz, gravit        ! J/K/molecule
         307             :     use tropopause,   only : tropopause_findChemTrop
         308             :
         309             :     use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk
         310             :     use physconst,      only : pi
         311             :
         312             :     type(physics_state), intent(in)    :: state(begchunk:endchunk)
         313             :
         314             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         315             :
         316      370944 :     type(physics_buffer_desc), pointer :: pbuf_chnk(:)
         317             :
         318             :     integer :: c,ncol,i,k
         319             :     real(r8) :: to_mmr(pcols,pver)
         320             :     real(r8), parameter :: molmass = 4._r8/3._r8*98.0_r8 !convert dry mass to wet mass of h2so4
         321             :     real(r8) :: volcmass1(pcols,pver)
         322             :     real(r8) :: volcmass2(pcols,pver)
         323             :     real(r8) :: volcmass3(pcols,pver)
         324             :     real(r8) :: columnmass1(pcols)
         325             :     real(r8) :: columnmass2(pcols)
         326             :     real(r8) :: columnmass3(pcols)
         327             :     integer  :: tropLev(pcols)
         328             :     real(r8) :: area_fact, radius_fact
         329             :
         330      370944 :     real(r8), pointer :: mass1(:,:)
         331      370944 :     real(r8), pointer :: mass2(:,:)
         332      370944 :     real(r8), pointer :: mass3(:,:)
         333      370944 :     real(r8), pointer :: area(:,:)
         334      370944 :     real(r8), pointer :: radius1(:,:)
         335      370944 :     real(r8), pointer :: radius2(:,:)
         336      370944 :     real(r8), pointer :: radius3(:,:)
         337             :
         338             :     !WACCM-derived relation between mass concentration and wet aerosol radius in meters
         339             :     real(r8),parameter :: radius_conversion = 1.9e-4_r8
         340             :
         341             :     logical :: zero_aerosols
         342             :     real(r8), parameter :: rad2deg = 180._r8/pi                ! radians to degrees conversion factor
         343             :
         344      370944 :     if( .not. has_prescribed_strataero ) return
         345             :
         346           0 :     call advance_trcdata( fields, file, state, pbuf2d )
         347             :
         348             :     ! copy prescribed tracer fields into state svariable with the correct units
         349           0 :     do c = begchunk,endchunk
         350             :
         351           0 :        pbuf_chnk => pbuf_get_chunk(pbuf2d, c)
         352             :
         353           0 :        ncol = state(c)%ncol
         354             :
         355           0 :        select case ( to_lower(trim(fields(1)%units(:GLC(fields(1)%units)))) )
         356             :        case ("molecules/cm3air", "molec/cm3","/cm3","molecules/cm3","cm^-3","cm**-3")
         357           0 :           to_mmr(:ncol,:) = (molmass*1.e6_r8*boltz*state(c)%t(:ncol,:))/(mwdry*state(c)%pmiddry(:ncol,:))
         358             :        case ('kg/kg','mmr','kg kg-1')
         359           0 :           to_mmr(:ncol,:) = 1._r8 ! input file must have converted to wet sulfate mass (=4/3*dry mass)
         360             :        case ('mol/mol','mole/mole','vmr','fraction')
         361           0 :           to_mmr(:ncol,:) = molmass/mwdry
         362             :        case default
         363           0 :           write(iulog,*) 'prescribed_strataero_adv: mass units = ',trim(fields(1)%units) ,' are not recognized'
         364           0 :           call endrun('prescribed_strataero_adv: mass units are not recognized')
         365             :        end select
         366             :
         367           0 :        if (mmr_ndx1>0) call pbuf_get_field(pbuf_chnk, mmr_ndx1, mass1)
         368           0 :        if (mmr_ndx2>0) call pbuf_get_field(pbuf_chnk, mmr_ndx2, mass2)
         369           0 :        if (mmr_ndx3>0) call pbuf_get_field(pbuf_chnk, mmr_ndx3, mass3)
         370             :
         371           0 :        if (three_mode) then
         372           0 :           call outfld( dens_name1, mass1(:,:), pcols, state(c)%lchnk)
         373           0 :           call outfld( dens_name2, mass2(:,:), pcols, state(c)%lchnk)
         374           0 :           call outfld( dens_name3, mass3(:,:), pcols, state(c)%lchnk)
         375             :        else
         376           0 :           call outfld( dens_name, mass1(:,:), pcols, state(c)%lchnk)
         377             :        endif
         378             :
         379           0 :        if (mmr_ndx1>0) mass1(:ncol,:) = to_mmr(:ncol,:) * mass1(:ncol,:) ! mmr
         380           0 :        if (mmr_ndx2>0) mass2(:ncol,:) = to_mmr(:ncol,:) * mass2(:ncol,:) ! mmr
         381           0 :        if (mmr_ndx3>0) mass3(:ncol,:) = to_mmr(:ncol,:) * mass3(:ncol,:) ! mmr
         382             :
         383           0 :        if (rad_ndx1>0) call pbuf_get_field(pbuf_chnk, rad_ndx1, radius1)
         384           0 :        if (rad_ndx2>0) call pbuf_get_field(pbuf_chnk, rad_ndx2, radius2)
         385           0 :        if (rad_ndx3>0) call pbuf_get_field(pbuf_chnk, rad_ndx3, radius3)
         386             :
         387           0 :        select case ( to_lower(trim(fields(rad_fld_no)%units(:GLC(fields(rad_fld_no)%units)))) )
         388             :        case ("m","meters")
         389           0 :           radius_fact = 1._r8
         390             :        case ("cm","centimeters")
         391           0 :           radius_fact = 1.e-2_r8
         392             :        case default
         393           0 :           write(iulog,*) 'prescribed_strataero_adv: radius units = ',trim(fields(rad_fld_no)%units) ,' are not recognized'
         394           0 :           call endrun('prescribed_strataero_adv: radius units are not recognized')
         395             :        end select
         396             :
         397             :        !MAM output is diamter so we need to half the value
         398           0 :        if (three_mode) then
         399           0 :           radius1(:ncol,:) = radius_fact*radius1(:ncol,:)*0.5_r8
         400           0 :           radius2(:ncol,:) = radius_fact*radius2(:ncol,:)*0.5_r8
         401           0 :           radius3(:ncol,:) = radius_fact*radius3(:ncol,:)*0.5_r8
         402             :        else
         403           0 :           radius1(:ncol,:) = radius_fact*radius1(:ncol,:)
         404             :        endif
         405             :
         406           0 :        call pbuf_get_field(pbuf_chnk, sad_ndx, area)
         407             :
         408           0 :        select case ( to_lower(trim(fields(sad_fld_no)%units(:7))) )
         409             :        case ("um2/cm3")
         410           0 :           area_fact = 1.e-8_r8
         411             :        case ("cm2/cm3")
         412           0 :           area_fact = 1._r8
         413             :        case default
         414           0 :           write(iulog,*) 'prescribed_strataero_adv: surface area density units = ',&
         415           0 :                          trim(fields(rad_fld_no)%units) ,' are not recognized'
         416           0 :           call endrun('prescribed_strataero_adv: surface area density units are not recognized')
         417             :        end select
         418           0 :        area(:ncol,:) = area_fact*area(:ncol,:)
         419             :
         420             :        ! this definition of tropopause is consistent with what is used in chemistry
         421             :        !REMOVECAM - no longer need this when CAM is retired and pcols no longer exists
         422           0 :        tropLev = 0
         423             :        !REMOVECAM_END
         424           0 :        call tropopause_findChemTrop(state(c), tropLev)
         425             :
         426           0 :        do i = 1,ncol
         427           0 :           do k = 1,pver
         428           0 :              zero_aerosols = k >= tropLev(i)
         429           0 :              if ( .not.prescribed_strataero_use_chemtrop .and. abs( state(c)%lat(i)*rad2deg ) > 50._r8 ) then
         430           0 :                 zero_aerosols = state(c)%pmid(i,k) >= 30000._r8
         431             :              endif
         432             :              ! set to zero at and below tropopause
         433           0 :              if ( zero_aerosols ) then
         434           0 :                 if (mmr_ndx1>0) mass1(i,k) = 0._r8
         435           0 :                 if (mmr_ndx2>0) mass2(i,k) = 0._r8
         436           0 :                 if (mmr_ndx3>0) mass3(i,k) = 0._r8
         437           0 :                 if (rad_ndx1>0) radius1(i,k) = 0._r8
         438           0 :                 if (rad_ndx2>0) radius2(i,k) = 0._r8
         439           0 :                 if (rad_ndx3>0) radius3(i,k) = 0._r8
         440           0 :                 area(i,k) = 0._r8
         441             :              endif
         442             :           enddo
         443             :        enddo
         444             :
         445           0 :        volcmass1(:ncol,:) = mass1(:ncol,:)*state(c)%pdel(:ncol,:)/gravit
         446           0 :        columnmass1(:ncol) = sum(volcmass1(:ncol,:), 2)
         447             :
         448           0 :        if (three_mode) then
         449           0 :           volcmass2(:ncol,:) = mass2(:ncol,:)*state(c)%pdel(:ncol,:)/gravit
         450           0 :           volcmass3(:ncol,:) = mass3(:ncol,:)*state(c)%pdel(:ncol,:)/gravit
         451           0 :           columnmass2(:ncol) = sum(volcmass2(:ncol,:), 2)
         452           0 :           columnmass3(:ncol) = sum(volcmass3(:ncol,:), 2)
         453           0 :           call outfld( mmr_name1,         mass1(:,:),     pcols, state(c)%lchnk)
         454           0 :           call outfld( mmr_name2,         mass2(:,:),     pcols, state(c)%lchnk)
         455           0 :           call outfld( mmr_name3,         mass3(:,:),     pcols, state(c)%lchnk)
         456           0 :           call outfld( mass_name1,        volcmass1(:,:), pcols, state(c)%lchnk)
         457           0 :           call outfld( mass_name2,        volcmass2(:,:), pcols, state(c)%lchnk)
         458           0 :           call outfld( mass_name3,        volcmass3(:,:), pcols, state(c)%lchnk)
         459           0 :           call outfld( mass_column_name1, columnmass1(:), pcols, state(c)%lchnk)
         460           0 :           call outfld( mass_column_name2, columnmass2(:), pcols, state(c)%lchnk)
         461           0 :           call outfld( mass_column_name3, columnmass3(:), pcols, state(c)%lchnk)
         462           0 :           call outfld( rad_name1,         radius1(:,:),   pcols, state(c)%lchnk)
         463           0 :           call outfld( rad_name2,         radius2(:,:),   pcols, state(c)%lchnk)
         464           0 :           call outfld( rad_name3,         radius3(:,:),   pcols, state(c)%lchnk)
         465             :        else
         466           0 :           call outfld( mmr_name,         mass1(:,:),     pcols, state(c)%lchnk)
         467           0 :           call outfld( mass_name,        volcmass1(:,:), pcols, state(c)%lchnk)
         468           0 :           call outfld( mass_column_name, columnmass1(:), pcols, state(c)%lchnk)
         469           0 :           call outfld( rad_name,         radius1(:,:),   pcols, state(c)%lchnk)
         470             :        endif
         471             :
         472           0 :        call outfld( sad_name,         area(:,:),     pcols, state(c)%lchnk)
         473             :
         474             :     enddo
         475             :
         476      370944 :   end subroutine prescribed_strataero_adv
         477             :
         478             : !-------------------------------------------------------------------
         479           0 :   subroutine init_prescribed_strataero_restart( piofile )
         480      370944 :     use pio, only : file_desc_t
         481             :     use tracer_data, only : init_trc_restart
         482             :
         483             :     type(file_desc_t),intent(inout) :: pioFile     ! pio File pointer
         484             :
         485           0 :     call init_trc_restart( 'prescribed_strataero', piofile, file )
         486             :
         487           0 :   end subroutine init_prescribed_strataero_restart
         488             : !-------------------------------------------------------------------
         489           0 :   subroutine write_prescribed_strataero_restart( piofile )
         490           0 :     use tracer_data, only : write_trc_restart
         491             :     use pio, only : file_desc_t
         492             :
         493             :     type(file_desc_t) :: piofile
         494             :
         495           0 :     call write_trc_restart( piofile, file )
         496             :
         497           0 :   end subroutine write_prescribed_strataero_restart
         498             : !-------------------------------------------------------------------
         499             : !-------------------------------------------------------------------
         500           0 :   subroutine read_prescribed_strataero_restart( pioFile )
         501           0 :     use tracer_data, only : read_trc_restart
         502             :     use pio, only : file_desc_t
         503             :
         504             :     type(file_desc_t) :: piofile
         505             :
         506           0 :     call read_trc_restart( 'prescribed_strataero', piofile, file )
         507             :
         508           0 :   end subroutine read_prescribed_strataero_restart
         509             :
         510             : end module prescribed_strataero

Generated by: LCOV version 1.14