LCOV - code coverage report
Current view: | top level - chemistry/utils - prescribed_aero.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 114 | 218 | 52.3 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 7 | 13 | 53.8 %

              Line data    Source code
           1             : !-------------------------------------------------------------------
           2             : ! Manages reading and interpolation of prescribed aerosol concentrations.
           3             : ! This places the concentration fields in the physics buffer so that
           4             : ! radiation package can access them.
           5             : !
           6             : ! This has been generalized so that the field names in the data files
           7             : ! and the field names in the physics buffer can be arbitrary.
           8             : !
           9             : ! The prescribed_aero_specifier namelist variable specifies a list of
          10             : ! variable names of the concentration fields in the netCDF dataset (ncdf_fld_name)
          11             : ! and the corresponding names used in the physics buffer:
          12             : !
          13             : ! prescribed_aero_specifier = 'pbuf_name1:ncdf_fld_name1','pbuf_name2:ncdf_fld_name2', ...
          14             : !
          15             : ! If there is no ":" then the specified name is used as both the
          16             : ! pbuf_name and ncdf_fld_name
          17             : !
          18             : ! Created by: Francis Vitt
          19             : !-------------------------------------------------------------------
          20             : module prescribed_aero
          21             :
          22             :   use shr_kind_mod,     only : r8 => shr_kind_r8
          23             :   use cam_abortutils,   only : endrun
          24             :   use spmd_utils,       only : masterproc
          25             :   use tracer_data,      only : trfld, trfile
          26             :   use cam_logfile,      only : iulog
          27             :   use pio,              only : var_desc_t
          28             :
          29             :   implicit none
          30             :   private
          31             :   save
          32             :
          33             :   type(trfld), pointer :: fields(:)
          34             :   type(trfile)         :: file
          35             :
          36             :   public :: prescribed_aero_init
          37             :   public :: prescribed_aero_adv
          38             :   public :: write_prescribed_aero_restart
          39             :   public :: read_prescribed_aero_restart
          40             :   public :: has_prescribed_aero
          41             :   public :: prescribed_aero_register
          42             :   public :: init_prescribed_aero_restart
          43             :   public :: prescribed_aero_readnl
          44             :
          45             :   logical :: has_prescribed_aero = .false.
          46             :
          47             :   ! Decides if its a modal aerosol simulation or not
          48             :   logical :: clim_modal_aero     = .false.
          49             :
          50             :   integer, parameter, public :: N_AERO = 50
          51             :
          52             :   integer :: number_flds
          53             :
          54             :   character(len=256) :: filename = 'NONE'
          55             :   character(len=256) :: filelist = ' '
          56             :   character(len=256) :: datapath = ' '
          57             :   character(len=32)  :: datatype = 'SERIAL'
          58             :   logical            :: rmv_file = .false.
          59             :   integer            :: cycle_yr = 0
          60             :   integer            :: fixed_ymd = 0
          61             :   integer            :: fixed_tod = 0
          62             :   character(len=32)  :: specifier(N_AERO) = ''
          63             :
          64             :   ! prescribed aerosol names
          65             :   character(len=16), allocatable :: pbuf_names(:)
          66             :
          67             :   integer :: aero_cnt
          68             :   integer :: aero_cnt_c = 0 ! clound borne species count for modal aerosols (zero otherwise)
          69             :
          70             :   ! Normal random number which persists from one tiemstep to the next
          71             :   real(r8) :: randn_persists
          72             :
          73             :   ! Following definitions are added to
          74             :   ! allow randn_persists to persist during restart runs
          75             :   type(var_desc_t) :: randn_persists_desc
          76             :   character(len=16), parameter :: randn_persists_name = 'prescraero_randn'
          77             :
          78             : contains
          79             :
          80             : !-------------------------------------------------------------------
          81             : ! registers aerosol fields to the phys buffer
          82             : !-------------------------------------------------------------------
          83        2304 :   subroutine prescribed_aero_register()
          84             :
          85             :     use ppgrid,         only: pver,pcols
          86             :
          87             :     use physics_buffer, only : pbuf_add_field, dtype_r8
          88             :     integer :: i,idx
          89             :
          90        1536 :     if (has_prescribed_aero) then
          91       21504 :        do i = 1,aero_cnt
          92       21504 :           call pbuf_add_field(pbuf_names(i),'physpkg',dtype_r8,(/pcols,pver/),idx)
          93             :        enddo
          94             :     endif
          95             :
          96        1536 :   endsubroutine prescribed_aero_register
          97             :
          98             : !-------------------------------------------------------------------
          99             : ! parses prescribed_aero_specifier namelist option
         100             : !-------------------------------------------------------------------
         101        1536 :   subroutine prescribed_aero_init()
         102             :
         103        1536 :     use tracer_data, only : trcdata_init
         104             :     use cam_history, only : addfld
         105             :     use ppgrid,      only : pver
         106             :     use error_messages, only: handle_err
         107             :     use ppgrid,         only: pcols, pver, begchunk, endchunk
         108             :     use physics_buffer, only : physics_buffer_desc
         109             :
         110             :     implicit none
         111             :
         112             :     ! local vars
         113             :     character(len=32)  :: spec_a
         114             :     integer :: ndx, istat, i, i_c, j
         115             :
         116        1536 :     if ( has_prescribed_aero ) then
         117        1536 :        if ( masterproc ) then
         118           2 :           write(iulog,*) 'aero is prescribed in :'//trim(filename)
         119             :        endif
         120             :     else
         121             :        return
         122             :     endif
         123             :
         124             :
         125        1536 :     allocate (file%in_pbuf(size(specifier)))
         126        1536 :     if (clim_modal_aero) then
         127           0 :       file%in_pbuf(:) = .false.
         128           0 :       do i = 1,N_AERO
         129           0 :           do j=1,aero_cnt
         130           0 :               if(specifier(i) .eq. pbuf_names(j)) then
         131           0 :                   file%in_pbuf(i) = .true.
         132           0 :                   exit
         133             :               endif
         134             :           enddo
         135             :       enddo
         136             :     else
         137       78336 :       file%in_pbuf(:) = .true.
         138             :     endif
         139             :     call trcdata_init( specifier, filename, filelist, datapath, fields, file, &
         140        1536 :                        rmv_file, cycle_yr, fixed_ymd, fixed_tod, datatype)
         141             :
         142        1536 :     number_flds = 0
         143        1536 :     if (associated(fields)) number_flds = size( fields )
         144             :
         145        1536 :     if( number_flds < 1 ) then
         146           0 :        if ( masterproc ) then
         147           0 :           write(iulog,*) 'There are no prescribed aerosols'
         148           0 :           write(iulog,*) ' '
         149             :        endif
         150           0 :        return
         151             :     end if
         152             :
         153             :     ! Following loop add fields for output. For modal aerosols, first the cld borne aersols
         154             :     ! are added and then their interstitial counterparts are added. The loop exits once all the cld borne
         155             :     ! aerosols (and their interstitial counterparts) are added.  Note that the units(fields(i)%units),
         156             :     ! will be same for both interstitial and cloud borne species.
         157             :     ! All other aerosol treatments(bulk) are left unchanged
         158             :     i_c = 0
         159       21504 :     fldloop:do i = 1,number_flds
         160       21504 :        if(clim_modal_aero .and. index(trim(fields(i)%fldnam),'_c') > 1) then ! Only cloud borne species
         161           0 :           call addfld(trim(fields(i)%fldnam)//'_D', (/ 'lev' /), 'A',trim(fields(i)%units), 'prescribed aero' )
         162           0 :           call spec_c_to_a(trim(fields(i)%fldnam),spec_a)
         163           0 :           call addfld(trim(spec_a)//'_D', (/ 'lev' /), 'A',trim(fields(i)%units), 'prescribed aero' )
         164           0 :           i_c = i_c + 1
         165           0 :           if(i_c >= aero_cnt_c) exit fldloop
         166             :        else
         167       39936 :           call addfld(trim(fields(i)%fldnam)//'_D', (/ 'lev' /), 'A',trim(fields(i)%units), 'prescribed aero' )
         168             :        endif
         169             :     enddo fldloop
         170             :
         171        1536 :   end subroutine prescribed_aero_init
         172             :
         173             : !-------------------------------------------------------------------
         174             : ! reads namelist options
         175             : !-------------------------------------------------------------------
         176        1536 : subroutine prescribed_aero_readnl(nlfile)
         177             :
         178        1536 :    use namelist_utils,  only: find_group_name
         179             :    use units,           only: getunit, freeunit
         180             :    use mpishorthand
         181             :    use rad_constituents, only: rad_cnst_get_info ! Added to query if it is a modal aero sim or not
         182             :
         183             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
         184             :
         185             :    ! Local variables
         186             :    integer :: unitn, ierr, nmodes, aero_loop_end
         187             :    logical :: skip_spec
         188             :    character(len=*), parameter :: subname = 'prescribed_aero_readnl'
         189             :
         190             :    character(len=32)  :: prescribed_aero_specifier(N_AERO)
         191             :    character(len=256) :: prescribed_aero_file
         192             :    character(len=256) :: prescribed_aero_filelist
         193             :    character(len=256) :: prescribed_aero_datapath
         194             :    character(len=32)  :: prescribed_aero_type
         195             :    logical            :: prescribed_aero_rmfile
         196             :    integer            :: prescribed_aero_cycle_yr
         197             :    integer            :: prescribed_aero_fixed_ymd
         198             :    integer            :: prescribed_aero_fixed_tod
         199             :    integer :: i,k
         200             :
         201             :    namelist /prescribed_aero_nl/ &
         202             :       prescribed_aero_specifier, &
         203             :       prescribed_aero_file,      &
         204             :       prescribed_aero_filelist,  &
         205             :       prescribed_aero_datapath,  &
         206             :       prescribed_aero_type,      &
         207             :       prescribed_aero_rmfile,    &
         208             :       prescribed_aero_cycle_yr,  &
         209             :       prescribed_aero_fixed_ymd, &
         210             :       prescribed_aero_fixed_tod
         211             :    !-----------------------------------------------------------------------------
         212             :
         213             :    ! Initialize namelist variables from local module variables.
         214       78336 :    prescribed_aero_specifier= specifier
         215        1536 :    prescribed_aero_file     = filename
         216        1536 :    prescribed_aero_filelist = filelist
         217        1536 :    prescribed_aero_datapath = datapath
         218        1536 :    prescribed_aero_type     = datatype
         219        1536 :    prescribed_aero_rmfile   = rmv_file
         220        1536 :    prescribed_aero_cycle_yr = cycle_yr
         221        1536 :    prescribed_aero_fixed_ymd= fixed_ymd
         222        1536 :    prescribed_aero_fixed_tod= fixed_tod
         223             :
         224             :    ! Read namelist
         225        1536 :    if (masterproc) then
         226           2 :       unitn = getunit()
         227           2 :       open( unitn, file=trim(nlfile), status='old' )
         228           2 :       call find_group_name(unitn, 'prescribed_aero_nl', status=ierr)
         229           2 :       if (ierr == 0) then
         230           2 :          read(unitn, prescribed_aero_nl, iostat=ierr)
         231           2 :          if (ierr /= 0) then
         232           0 :             call endrun(subname // ':: ERROR reading namelist')
         233             :          end if
         234             :       end if
         235           2 :       close(unitn)
         236           2 :       call freeunit(unitn)
         237             :    end if
         238             :
         239             : #ifdef SPMD
         240             :    ! Broadcast namelist variables
         241        1536 :    call mpibcast(prescribed_aero_specifier,len(prescribed_aero_specifier(1))*N_AERO,     mpichar, 0, mpicom)
         242        1536 :    call mpibcast(prescribed_aero_file,     len(prescribed_aero_file),     mpichar, 0, mpicom)
         243        1536 :    call mpibcast(prescribed_aero_filelist, len(prescribed_aero_filelist), mpichar, 0, mpicom)
         244        1536 :    call mpibcast(prescribed_aero_datapath, len(prescribed_aero_datapath), mpichar, 0, mpicom)
         245        1536 :    call mpibcast(prescribed_aero_type,     len(prescribed_aero_type),     mpichar, 0, mpicom)
         246        1536 :    call mpibcast(prescribed_aero_rmfile,   1, mpilog,  0, mpicom)
         247        1536 :    call mpibcast(prescribed_aero_cycle_yr, 1, mpiint,  0, mpicom)
         248        1536 :    call mpibcast(prescribed_aero_fixed_ymd,1, mpiint,  0, mpicom)
         249        1536 :    call mpibcast(prescribed_aero_fixed_tod,1, mpiint,  0, mpicom)
         250             : #endif
         251             :
         252             :    ! Update module variables with user settings.
         253       78336 :    specifier  = prescribed_aero_specifier
         254        1536 :    filename   = prescribed_aero_file
         255        1536 :    filelist   = prescribed_aero_filelist
         256        1536 :    datapath   = prescribed_aero_datapath
         257        1536 :    datatype   = prescribed_aero_type
         258        1536 :    rmv_file   = prescribed_aero_rmfile
         259        1536 :    cycle_yr   = prescribed_aero_cycle_yr
         260        1536 :    fixed_ymd  = prescribed_aero_fixed_ymd
         261        1536 :    fixed_tod  = prescribed_aero_fixed_tod
         262             :
         263             :    ! Turn on prescribed aerosols if user has specified an input dataset.
         264        1536 :    has_prescribed_aero = len_trim(filename) > 0 .and. filename.ne.'NONE'
         265             :
         266        1536 :    if ( .not. has_prescribed_aero) return
         267             :
         268             :    ! Determine whether its a 'modal' aerosol simulation  or not
         269        1536 :    call rad_cnst_get_info(0, nmodes=nmodes)
         270        1536 :    clim_modal_aero = (nmodes > 0)
         271             :
         272             :    ! For modal aerosols, interstitial species(*_a) are diagnosed from
         273             :    ! their *_logm and *_logv counterparts (e.g. soa_a1 is diagnosed from
         274             :    ! soa_a1_logm and soa_a1_logv). Therefore, only *_logm and *_logv and cloud
         275             :    ! borne (*_c) species are specified in the build-namelist. In the following
         276             :    ! cnt_loop, we will count the cloud borne species and *_logm species(in lieu
         277             :    ! of *_a species). We will skip *_logv species. This will ensure that aero_cnt
         278             :    ! variable is the sum of cloud borne and interstitial species (which we will
         279             :    ! manually add in pbuf_names array later). We are also counting cloud borne
         280             :    ! (*_c) species which will help adding the same number of interstitial species
         281             :    ! in pbuf_names array
         282             :
         283             :    ! determine which prescibed aerosols are specified
         284        1536 :    aero_cnt   = 0
         285        1536 :    aero_cnt_c = 0 ! cloud borne species count
         286       21504 :    cnt_loop: do i = 1,N_AERO
         287       21504 :       if ( len_trim(specifier(i))==0 ) exit cnt_loop
         288       19968 :       skip_spec = .FALSE.
         289       19968 :       if(clim_modal_aero) then
         290             :          ! For modal aerosols, skip counting species ending with *_logv
         291           0 :          if(index(specifier(i),'_c') >= 1) aero_cnt_c = aero_cnt_c + 1
         292           0 :          if(index(specifier(i),'_logv') >= 1) skip_spec = .TRUE.
         293             :       endif
         294       21504 :       if(.NOT.skip_spec)aero_cnt = aero_cnt+1
         295             :    end do cnt_loop
         296             :
         297        1536 :    has_prescribed_aero = aero_cnt>0
         298        1536 :    if ( .not. has_prescribed_aero) return
         299             :
         300        4608 :    allocate(pbuf_names(aero_cnt))
         301             :
         302             :    !  For modal aerosols, the following loop will add the cloud borne
         303             :    ! species directly and interstitial species through the "add_interstitial_spec"
         304             :    ! call. Interstitial species are added at the end of the cloud borne species in
         305             :    ! pbuf_names array.
         306             :    ! Note that aero_cnt_c should be zero for all other aerosol trearments
         307             :    ! except the modal aerosols (e.g. bulk)
         308             :
         309        1536 :    if(.NOT. clim_modal_aero) aero_cnt_c = 0
         310        1536 :    aero_loop_end = aero_cnt - aero_cnt_c
         311             :
         312       21504 :    do i = 1,aero_loop_end
         313       19968 :       k = scan( specifier(i),':')
         314       21504 :       if (k>1) then
         315       19968 :          pbuf_names(i) = trim(specifier(i)(:k-1))
         316       19968 :          if(clim_modal_aero)call add_interstitial_spec(aero_loop_end,i)
         317             :       else
         318           0 :          pbuf_names(i) = trim(specifier(i))
         319           0 :          if(clim_modal_aero)call add_interstitial_spec(aero_loop_end,i)
         320             :       endif
         321             :    enddo
         322             :
         323        1536 : end subroutine prescribed_aero_readnl
         324             :
         325             : !-------------------------------------------------------------------
         326             : ! Add interstitial aerosols in pbuf_names array for modal aerosols
         327             : !-------------------------------------------------------------------
         328           0 : subroutine add_interstitial_spec(aero_loop_end,i_in)
         329             :   implicit none
         330             :
         331             :   !Arguments
         332             :   integer, intent(in) :: i_in, aero_loop_end
         333             :
         334             :   !Local
         335             :   character(len=32) :: spec_a
         336             :
         337             :   !  Replace 'c' with 'a' in species name
         338           0 :   call spec_c_to_a(pbuf_names(i_in), spec_a)
         339           0 :   pbuf_names(aero_loop_end+i_in) = spec_a
         340        1536 : end subroutine add_interstitial_spec
         341             :
         342             : !-------------------------------------------------------------------
         343             : ! A generic subroutine which replaces 'c' in the cloud borne
         344             : ! species name with 'a' to make it interstital species
         345             : !-------------------------------------------------------------------
         346           0 : subroutine spec_c_to_a (spec_c_in,spec_a_out)
         347             :   implicit none
         348             :
         349             :   !Arguments
         350             :   character(len=*), intent(in)  :: spec_c_in
         351             :   character(len=*), intent(out) :: spec_a_out
         352             :
         353             :   !Local
         354             :   character(len=32)   :: name
         355             :   character(len=1000) :: errMsg
         356             :   integer             :: k_c, k_cp1
         357             :
         358           0 :   k_c = index(trim(adjustl(spec_c_in)),'_c')
         359           0 :   if(k_c >= 1) then
         360           0 :      k_cp1             = k_c + 1
         361           0 :      name              = trim(adjustl(spec_c_in))
         362           0 :      name(k_cp1:k_cp1) = 'a'
         363           0 :      spec_a_out        = trim(adjustl(name))
         364             :   else
         365           0 :      write(errMsg,*) "Input string (",trim(spec_c_in)," is not a cld borne aerosol,", &
         366           0 :           "cannot form interstitial species name"
         367           0 :      call endrun(trim(errMsg))
         368             :   endif
         369           0 : end subroutine spec_c_to_a
         370             :
         371             : !-------------------------------------------------------------------
         372             : ! advances the aerosol fields to the current time step
         373             : !-------------------------------------------------------------------
         374      741888 :   subroutine prescribed_aero_adv( state, pbuf2d )
         375             :
         376             :     use tracer_data,  only : advance_trcdata
         377             :     use physics_types,only : physics_state
         378             :     use ppgrid,       only : begchunk, endchunk
         379             :     use ppgrid,       only : pcols, pver
         380             :     use string_utils, only : to_lower, GLC
         381             :     use cam_history,  only : outfld
         382             :
         383             :     use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk, &
         384             :          pbuf_get_index
         385             :
         386             :     implicit none
         387             :
         388             :     type(physics_state), intent(in)    :: state(begchunk:endchunk)
         389             :
         390             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         391             :
         392      370944 :     type(physics_buffer_desc), pointer :: pbuf_chnk(:)
         393             :
         394             :     character(len=32) :: spec_a
         395             :     integer :: c,ncol, i, i_c, spec_a_ndx, errcode
         396      370944 :     real(r8),pointer :: outdata(:,:)
         397             :     logical :: cld_borne_aero = .FALSE.
         398             :
         399      370944 :     if( .not. has_prescribed_aero ) return
         400             :
         401      370944 :     call advance_trcdata( fields, file, state, pbuf2d )
         402             :
         403             :     ! Diagnose interstital species using random sampling
         404      370944 :     if ( clim_modal_aero ) then
         405           0 :       call rand_sample_prescribed_aero(state, pbuf2d )
         406             :     endif
         407             :
         408             :     ! set the tracer fields with the correct units
         409      370944 :     i_c = 0
         410     5193216 :     fldloop : do i = 1,number_flds
         411             :
         412             : !$OMP PARALLEL DO PRIVATE (C, NCOL, OUTDATA,PBUF_CHNK)
         413    24262056 :        do c = begchunk,endchunk
         414    19439784 :           ncol = state(c)%ncol
         415    19439784 :           pbuf_chnk => pbuf_get_chunk(pbuf2d, c)
         416    24262056 :           if(clim_modal_aero .and. index(trim(fields(i)%fldnam),'_c') > 1) then ! Only cloud borne species
         417           0 :              call pbuf_get_field(pbuf_chnk, fields(i)%pbuf_ndx, outdata)
         418           0 :              call outfld( trim(fields(i)%fldnam)//'_D', outdata(:ncol,:), ncol, state(c)%lchnk )
         419             :
         420           0 :              call spec_c_to_a(trim(fields(i)%fldnam),spec_a)
         421           0 :              spec_a_ndx = pbuf_get_index(trim(fields(i)%fldnam),errcode)
         422           0 :              call pbuf_get_field(pbuf_chnk, spec_a_ndx, outdata)
         423           0 :              call outfld( trim(spec_a)//'_D', outdata(:ncol,:), ncol, state(c)%lchnk )
         424           0 :              cld_borne_aero = .TRUE.
         425             :           else
         426    19439784 :              call pbuf_get_field(pbuf_chnk, fields(i)%pbuf_ndx, outdata)
         427    19439784 :              call outfld( trim(fields(i)%fldnam)//'_D', outdata(:ncol,:), ncol, state(c)%lchnk )
         428             :           endif
         429             :        enddo
         430     5193216 :        if(cld_borne_aero)then
         431           0 :           i_c = i_c + 1
         432           0 :           cld_borne_aero = .FALSE.
         433           0 :           if(i_c >= aero_cnt_c) exit fldloop
         434             :        endif
         435             :     enddo fldloop
         436             :
         437      741888 :   end subroutine prescribed_aero_adv
         438             :
         439             : !-------------------------------------------------------------------
         440           0 :   subroutine rand_sample_prescribed_aero(state, pbuf2d)
         441             :
         442             :     !Purpose: Generates log normal distribution for the interstitial species.
         443             :     !Note that only the interstitial aerosols are diagnosed here
         444             :     !
         445             :     !Written by Balwinder Singh (12/14/2012)
         446             :     !Adapted from Jin-Ho Yoon
         447             :     !
         448             :     !Update log:
         449             :
         450      370944 :     use physics_types,  only : physics_state
         451             :     use ppgrid,         only : begchunk, endchunk, pver
         452             :     use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk, &
         453             :          pbuf_get_index
         454             :     use ppgrid,       only : pcols, pver
         455             :
         456             :     !Arguments
         457             :     type(physics_state), intent(in)    :: state(begchunk:endchunk)
         458             :     type(physics_buffer_desc), pointer :: pbuf2d(:,:)
         459             :
         460             :     !Local
         461             :     real(r8), parameter :: mean_max_val = 5.0_r8
         462             :     real(r8), parameter :: std_max_val  = 3.0_r8
         463             :
         464             :     integer  :: c, ncol, i, kc, klog, logv_ndx, logm_ndx, a_ndx, icol, kpver
         465             :     real(r8) :: logm2, variance, std, mean_max, std_max, randn
         466           0 :     type(physics_buffer_desc), pointer :: pbuf_chnk(:)
         467           0 :     real(r8), pointer :: data(:,:)
         468             :     real(r8) :: logm_data(pcols,pver),logv_data(pcols,pver)
         469             :
         470           0 :     randn = randn_prescribed_aero()
         471           0 :     do i = 1, aero_cnt
         472             :        !Species with '_a' are updated using random sampling.
         473             :        !Cloud borne species (ends with _c1,_c2, etc.) have to be skipped
         474           0 :        kc   = index(pbuf_names(i),'_c')
         475             :        ! if kc>1, species is cloud borne species
         476           0 :        if( kc > 1) cycle
         477             :
         478             :        !If species is ending with _a1, a2 etc., then find indicies of their
         479             :        !logv and lom conterparts
         480           0 :        logv_ndx = logvm_get_index(pbuf_names(i),'v')
         481           0 :        logm_ndx = logvm_get_index(pbuf_names(i),'m')
         482           0 :        a_ndx    = pbuf_get_index(trim(adjustl(pbuf_names(i))))
         483             :
         484           0 :        do c = begchunk, endchunk
         485           0 :           ncol = state(c)%ncol
         486           0 :           pbuf_chnk => pbuf_get_chunk(pbuf2d, c)
         487             :
         488             :           !Get the species mixing ratio nad its logv and logm values
         489           0 :           call pbuf_get_field(pbuf_chnk, a_ndx, data)
         490           0 :           logv_data = fields(logv_ndx)%data(:,:,c)
         491           0 :           logm_data = fields(logm_ndx)%data(:,:,c)
         492             :
         493           0 :           do icol = 1, ncol
         494           0 :              do kpver = 1, pver
         495           0 :                 logm2    = logm_data(icol,kpver) * logm_data(icol,kpver)
         496             :
         497             :                 !Compute variance
         498           0 :                 variance = max( 0.0_r8, (logv_data(icol,kpver) - logm2) )
         499             :
         500             :                 !Standard deviation
         501           0 :                 std      = sqrt(variance)
         502             :
         503             :                 !Bounds to keep mixing ratios from going unphysical
         504           0 :                 mean_max = exp(logm_data(icol,kpver)) * mean_max_val
         505           0 :                 std_max  = exp(logm_data(icol,kpver)  + std_max_val * std )
         506             :
         507           0 :                 data(icol,kpver) = min(exp(logm_data(icol,kpver)+randn*std),mean_max,std_max)
         508             :
         509             :              enddo !pver
         510             :           enddo !col
         511             :        enddo !chunk
         512             :     enddo !flds
         513             :
         514           0 :   end subroutine rand_sample_prescribed_aero
         515             : !-------------------------------------------------------------------
         516           0 :   function randn_prescribed_aero()
         517             :     !Pupose: This function generates a new normally distributed random
         518             :     ! number at end end of each day. This random number then stays the same
         519             :     ! for the whole day
         520             :     !
         521             :     !Written by Balwinder Singh (12/14/2012)
         522             :     !Adapted from Jin-Ho Yoon
         523             :     !
         524             :     !Update log:
         525             :
         526           0 :     use time_manager,   only : is_end_curr_day, is_first_step, get_nstep
         527             :
         528             :     integer, parameter :: rconst1_1 = 5000000
         529             :     integer, parameter :: rconst1_2 = 50
         530             :     integer, parameter :: rconst2_1 = 10000000
         531             :     integer, parameter :: rconst2_2 = 10
         532             :
         533             :     integer :: i, seed_size, nstep
         534           0 :     integer, allocatable :: seed(:)
         535             :
         536             :     real(r8) :: randn_prescribed_aero
         537             :     real(r8) :: randu1, randu2
         538             :
         539             :     !Use same random number for the entire day and generate a new normally
         540             :     !distributed random number at the start of the new day
         541           0 :     if(is_first_step() .or. is_end_curr_day()) then
         542             :        !Generate two uniformly distributed random numbers (between 0 and 1)
         543           0 :        call random_seed(size=seed_size)
         544           0 :        allocate(seed(seed_size))
         545             :
         546             :        ! Using nstep as a seed to generate same sequence
         547           0 :        nstep = get_nstep()
         548           0 :        do i = 1, seed_size
         549           0 :           seed(i) = rconst1_1*nstep + rconst1_2*(i-1)
         550             :        end do
         551           0 :        call random_seed(put=seed)
         552           0 :        call random_number (randu1)
         553             :
         554           0 :        do i = 1, seed_size
         555           0 :           seed(i) = rconst2_1*nstep + rconst2_2*(i-1)
         556             :        end do
         557           0 :        call random_seed(put=seed)
         558           0 :        call random_number (randu2)
         559           0 :        deallocate(seed)
         560             :
         561             :        !Normal distribution (Mean = 0, standard dev = 1)
         562           0 :        randn_prescribed_aero = boxMuller(randu1,randu2)
         563           0 :        randn_persists = randn_prescribed_aero
         564             :     else
         565             :        !Use the previously generated random number
         566           0 :        randn_prescribed_aero = randn_persists
         567             :     endif
         568           0 :   end function randn_prescribed_aero
         569             : !-------------------------------------------------------------------
         570           0 :   function logvm_get_index(name,type) result(index)
         571             :     implicit none
         572             :
         573             :     !Args
         574             :     character(len=*), intent(in) :: name, type
         575             :
         576             :     !Loc
         577             :     character(len=64)   :: tmp_name
         578             :     character(len=1000) :: msgStr
         579             :     integer :: index, i
         580             :
         581             :
         582           0 :     index = -1
         583           0 :     tmp_name = trim(adjustl(name))//'_log'//trim(adjustl(type))
         584             :
         585           0 :     fldloop: do i = 1, number_flds
         586           0 :        if(fields(i)%fldnam == tmp_name) then
         587             :           index = i
         588             :           exit fldloop
         589             :        endif
         590             :     enddo fldloop
         591             :
         592           0 :     if(index == -1) then
         593           0 :        write(msgStr,*) "Prescribed_aero.F90: ",tmp_name," doesn't exist in the fields%fldnam"
         594           0 :        call endrun(msgStr)
         595             :     endif
         596             :
         597           0 :   end function logvm_get_index
         598             : !-------------------------------------------------------------------
         599           0 :   function boxMuller(ru1,ru2) result(rn)
         600             :     use physconst,      only : pi
         601             :     implicit none
         602             :
         603             :     !Args
         604             :     real(r8), intent(in)  :: ru1, ru2
         605             :
         606             :     !Loc
         607             :     real(r8), parameter :: pi2 = 2._r8 * pi
         608             :     real(r8) :: ur, theta, rn
         609             :
         610             :     !Based on Box Muller Method, generate normally distributed random numbers
         611           0 :     ur    = sqrt(-2._r8 * log(ru1))
         612           0 :     theta = pi2 * ru2
         613             :
         614             :     !Normal distribution (Mean = 0, standard dev = 1)
         615           0 :     rn = ur * cos(theta)
         616             :
         617           0 :   end function boxMuller
         618             :
         619             : !-------------------------------------------------------------------
         620             : !-------------------------------------------------------------------
         621        1536 :   subroutine init_prescribed_aero_restart( piofile )
         622             :     use pio, only : file_desc_t, pio_def_var, pio_double
         623             :     use tracer_data, only : init_trc_restart
         624             :     implicit none
         625             :     type(file_desc_t),intent(inout) :: pioFile     ! pio File pointer
         626             :     integer :: ierr
         627             :
         628             :     ! For allowing randn_persists to persist during reststarts
         629        1536 :     ierr = pio_def_var(piofile, randn_persists_name, pio_double, randn_persists_desc)
         630             :
         631        1536 :     call init_trc_restart( 'prescribed_aero', piofile, file )
         632             :
         633        1536 :   end subroutine init_prescribed_aero_restart
         634             : !-------------------------------------------------------------------
         635        1536 :   subroutine write_prescribed_aero_restart( piofile )
         636        1536 :     use tracer_data, only : write_trc_restart
         637             :     use pio, only : file_desc_t, pio_put_var
         638             :     implicit none
         639             :
         640             :     type(file_desc_t) :: piofile
         641             :     integer :: ierr
         642             :
         643             :     !  For allowing randn_persists to persist during reststarts
         644        3072 :     ierr = pio_put_var(piofile, randn_persists_desc, (/randn_persists/))
         645             :
         646        1536 :     call write_trc_restart( piofile, file )
         647             :
         648        1536 :   end subroutine write_prescribed_aero_restart
         649             :
         650             : !-------------------------------------------------------------------
         651             : !-------------------------------------------------------------------
         652         768 :   subroutine read_prescribed_aero_restart( pioFile )
         653        1536 :     use tracer_data, only : read_trc_restart
         654             :     use pio, only : file_desc_t, pio_inq_varid, pio_get_var
         655             :     implicit none
         656             :
         657             :     type(file_desc_t) :: piofile
         658             :     integer :: ierr
         659             :
         660             :     ! For allowing randn_persists to persist during reststarts
         661         768 :     ierr = pio_inq_varid(pioFile, randn_persists_name, randn_persists_desc)
         662         768 :     ierr = pio_get_var(pioFile, randn_persists_desc, randn_persists)
         663             :
         664         768 :     call read_trc_restart( 'prescribed_aero', piofile, file )
         665             :
         666         768 :   end subroutine read_prescribed_aero_restart
         667             :
         668             : end module prescribed_aero

Generated by: LCOV version 1.14