LCOV - code coverage report
Current view: | top level - physics/cam - aoa_tracers.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 29 | 110 | 26.4 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 5 | 8 | 62.5 %

              Line data    Source code
           1             : !===============================================================================
           2             : ! Age of air test tracers
           3             : ! provides dissipation rate and surface fluxes for diagnostic constituents
           4             : !===============================================================================
           5             :
           6             : module aoa_tracers
           7             :
           8             :   use shr_kind_mod, only: r8 => shr_kind_r8
           9             :   use spmd_utils,   only: masterproc
          10             :   use ppgrid,       only: pcols, pver
          11             :   use constituents, only: pcnst, cnst_add, cnst_name, cnst_longname
          12             :   use cam_logfile,  only: iulog
          13             :   use ref_pres,     only: pref_mid_norm
          14             :   use time_manager, only: get_curr_date, get_start_date
          15             :   use time_manager, only: is_leapyear, timemgr_get_calendar_cf, get_calday
          16             :
          17             :   implicit none
          18             :   private
          19             :
          20             :   ! Public interfaces
          21             :   public :: aoa_tracers_register         ! register constituents
          22             :   public :: aoa_tracers_implements_cnst  ! true if named constituent is implemented by this package
          23             :   public :: aoa_tracers_init_cnst        ! initialize constituent field
          24             :   public :: aoa_tracers_init             ! initialize history fields, datasets
          25             :   public :: aoa_tracers_timestep_init    ! place to perform per timestep initialization
          26             :   public :: aoa_tracers_timestep_tend    ! calculate tendencies
          27             :   public :: aoa_tracers_readnl           ! read namelist options
          28             :
          29             :   ! Private module data
          30             :
          31             :   integer, parameter :: ncnst=3  ! number of constituents implemented by this module
          32             :
          33             :   ! constituent names
          34             :   character(len=6), parameter :: c_names(ncnst) = (/'AOAMF ', 'HORZ  ', 'VERT  '/)
          35             :
          36             :   ! constituent source/sink names
          37             :   character(len=8), parameter :: src_names(ncnst) = (/'AOAMFSRC', 'HORZSRC ', 'VERTSRC '/)
          38             :
          39             :   integer :: ifirst = -1 ! global index of first constituent
          40             :   integer :: ixaoa  = -1 ! global index for AOAMFSRC tracer
          41             :   integer :: ixht   = -1 ! global index for HORZ tracer
          42             :   integer :: ixvt   = -1 ! global index for VERT tracer
          43             :
          44             :   ! Data from namelist variables
          45             :   logical :: aoa_tracers_flag  = .false.    ! true => turn on test tracer code, namelist variable
          46             :   logical :: aoa_read_from_ic_file = .true. ! true => tracers initialized from IC file
          47             :
          48             :   real(r8),  parameter ::  treldays = 15._r8
          49             :   real(r8),  parameter ::  vert_offset = 10._r8
          50             :
          51             :   ! 15-days used for diagnostic of transport circulation and K-tensors
          52             :   ! relaxation (in the original papers PM-1987 and YSGD-2000) => Zonal Mean
          53             :   ! to evaluate eddy-fluxes for 2D-diagnostics, here relaxation to the GLOBAL MEAN  IC
          54             :   ! it may help to keep gradients but will rule-out 2D-transport diagnostics
          55             :   ! in km  to avoid negative values of  vertical tracers
          56             :   ! VERT(k) = -7._r8*alog(hyam(k)+hybm(k)) + vert_offset
          57             :
          58             :   ! PM-1987:
          59             :   ! Plumb, R. A., and J. D. Mahlman (1987), The zonally averaged transport
          60             :   ! characteristics of the GFDL general circulation/transport model,
          61             :   ! J. Atmos.Sci.,44, 298-327
          62             :
          63             :   ! YSGD-2000:
          64             :   ! Yudin, Valery A., Sergey P. Smyshlyaev, Marvin A. Geller, Victor L. Dvortsov, 2000:
          65             :   ! Transport Diagnostics of GCMs and Implications for 2D Chemistry-Transport Model of
          66             :   ! Troposphere and Stratosphere. J. Atmos. Sci., 57, 673-699.
          67             :   ! doi: http://dx.doi.org/10.1175/1520-0469(2000)057<0673:TDOGAI>2.0.CO;2
          68             :
          69             :   real(r8) :: qrel_vert(pver) = -huge(1._r8)  ! = -7._r8*log(pref_mid_norm(k)) + vert_offset
          70             :
          71             :   integer :: yr0 = -huge(1)
          72             :   real(r8) :: calday0 = -huge(1._r8)
          73             :   real(r8) :: years = -huge(1._r8)
          74             :
          75             : !===============================================================================
          76             : contains
          77             : !===============================================================================
          78             :
          79             : !================================================================================
          80        1536 :   subroutine aoa_tracers_readnl(nlfile)
          81             :
          82             :     use namelist_utils, only: find_group_name
          83             :     use cam_abortutils, only: endrun
          84             :     use spmd_utils,     only: mpicom, masterprocid, mpi_logical, mpi_success
          85             :
          86             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          87             :
          88             :     ! Local variables
          89             :     integer :: unitn, ierr
          90             :     character(len=*), parameter :: subname = 'aoa_tracers_readnl'
          91             :
          92             :     namelist /aoa_tracers_nl/ aoa_tracers_flag, aoa_read_from_ic_file
          93             :
          94             :     !-----------------------------------------------------------------------------
          95             :
          96        1536 :     if (masterproc) then
          97           2 :        open( newunit=unitn, file=trim(nlfile), status='old' )
          98           2 :        call find_group_name(unitn, 'aoa_tracers_nl', status=ierr)
          99           2 :        if (ierr == 0) then
         100           0 :           read(unitn, aoa_tracers_nl, iostat=ierr)
         101           0 :           if (ierr /= 0) then
         102           0 :              call endrun(subname // ':: ERROR reading namelist')
         103             :           end if
         104             :        end if
         105           2 :        close(unitn)
         106             :     end if
         107             :
         108        1536 :     call mpi_bcast(aoa_tracers_flag, 1, mpi_logical, masterprocid, mpicom, ierr)
         109        1536 :     if (ierr/=mpi_success) then
         110           0 :        call endrun(subname//': MPI_BCAST ERROR: aoa_tracers_flag')
         111             :     end if
         112        1536 :     call mpi_bcast(aoa_read_from_ic_file, 1, mpi_logical, masterprocid, mpicom, ierr)
         113        1536 :     if (ierr/=mpi_success) then
         114           0 :        call endrun(subname//': MPI_BCAST ERROR: aoa_read_from_ic_file')
         115             :     end if
         116             :
         117        1536 :   endsubroutine aoa_tracers_readnl
         118             :
         119             : !================================================================================
         120             :
         121        1536 :   subroutine aoa_tracers_register
         122             :     !-----------------------------------------------------------------------
         123             :     !
         124             :     ! Purpose: register advected constituents
         125             :     !
         126             :     !-----------------------------------------------------------------------
         127             :     use physconst,  only: cpair, mwdry
         128             :     !-----------------------------------------------------------------------
         129             :
         130             :     integer :: k
         131             :
         132        1536 :     if (.not. aoa_tracers_flag) return
         133             :
         134             :     call cnst_add(c_names(1), mwdry, cpair, 0._r8, ixaoa, readiv=aoa_read_from_ic_file, &
         135           0 :                   longname='mixing ratio LB tracer')
         136             :
         137             :     call cnst_add(c_names(2), mwdry, cpair, 1._r8, ixht,   readiv=aoa_read_from_ic_file, &
         138           0 :                   longname='horizontal tracer')
         139             :     call cnst_add(c_names(3), mwdry, cpair, 0._r8, ixvt,   readiv=aoa_read_from_ic_file, &
         140           0 :                   longname='vertical tracer')
         141             :
         142           0 :     ifirst = ixaoa
         143             :
         144           0 :     do k = 1,pver
         145           0 :        qrel_vert(k) = -7._r8*log(pref_mid_norm(k)) + vert_offset
         146             :     enddo
         147             :
         148             :   end subroutine aoa_tracers_register
         149             :
         150             : !===============================================================================
         151             :
         152           0 :   function aoa_tracers_implements_cnst(name)
         153             :     !-----------------------------------------------------------------------
         154             :     !
         155             :     ! Purpose: return true if specified constituent is implemented by this package
         156             :     !
         157             :     !-----------------------------------------------------------------------
         158             :
         159             :     character(len=*), intent(in) :: name   ! constituent name
         160             :     logical :: aoa_tracers_implements_cnst        ! return value
         161             :
         162             :     !---------------------------Local workspace-----------------------------
         163             :     integer :: m
         164             :     !-----------------------------------------------------------------------
         165             :
         166           0 :     aoa_tracers_implements_cnst = .false.
         167             :
         168           0 :     if (.not. aoa_tracers_flag) return
         169             :
         170           0 :     do m = 1, ncnst
         171           0 :        if (name == c_names(m)) then
         172           0 :           aoa_tracers_implements_cnst = .true.
         173             :           return
         174             :        end if
         175             :     end do
         176             :
         177             :   end function aoa_tracers_implements_cnst
         178             :
         179             : !===============================================================================
         180             :
         181           0 :   subroutine aoa_tracers_init_cnst(name, latvals, lonvals, mask, q)
         182             :
         183             :     !-----------------------------------------------------------------------
         184             :     !
         185             :     ! Purpose: initialize test tracers mixing ratio fields
         186             :     !  This subroutine is called at the beginning of an initial run ONLY
         187             :     !
         188             :     !-----------------------------------------------------------------------
         189             :
         190             :     character(len=*), intent(in)  :: name
         191             :     real(r8),         intent(in)  :: latvals(:) ! lat in degrees (ncol)
         192             :     real(r8),         intent(in)  :: lonvals(:) ! lon in degrees (ncol)
         193             :     logical,          intent(in)  :: mask(:)    ! Only initialize where .true.
         194             :     real(r8),         intent(out) :: q(:,:)   ! kg tracer/kg dry air (gcol, plev)
         195             :
         196             :     integer :: m
         197             :     !-----------------------------------------------------------------------
         198             :
         199           0 :     if (.not. aoa_tracers_flag) return
         200             :
         201           0 :     do m = 1, ncnst
         202           0 :        if (name ==  c_names(m))  then
         203             :           ! pass global constituent index
         204           0 :           call init_cnst_3d(ifirst+m-1, latvals, lonvals, mask, q)
         205             :        endif
         206             :     end do
         207             :
         208             :   end subroutine aoa_tracers_init_cnst
         209             :
         210             : !===============================================================================
         211             :
         212        1536 :   subroutine aoa_tracers_init
         213             :
         214             :     !-----------------------------------------------------------------------
         215             :     !
         216             :     ! Purpose: initialize age of air constituents
         217             :     !          (declare history variables)
         218             :     !-----------------------------------------------------------------------
         219             :
         220             :     use cam_history,    only: addfld, add_default
         221             :
         222             :     integer :: m, mm
         223             :     integer :: yr, mon, day, sec, ymd
         224             :
         225             :     !-----------------------------------------------------------------------
         226             :
         227        1536 :     if (.not. aoa_tracers_flag) return
         228             :
         229             :     ! Set names of tendencies and declare them as history variables
         230             :
         231           0 :     do m = 1, ncnst
         232           0 :        mm = ifirst+m-1
         233           0 :        call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm))
         234           0 :        call addfld(src_names(m),  (/ 'lev' /), 'A', 'kg/kg/s', trim(cnst_name(mm))//' source/sink')
         235             :
         236           0 :        call add_default (cnst_name(mm), 1, ' ')
         237           0 :        call add_default (src_names(m),  1, ' ')
         238             :     end do
         239             :
         240           0 :     call get_start_date(yr, mon, day, sec)
         241             :
         242           0 :     ymd = yr*10000 + mon*100 + day
         243             :
         244           0 :     yr0 = yr
         245           0 :     calday0 = get_calday(ymd, sec)
         246             :
         247        1536 :   end subroutine aoa_tracers_init
         248             :
         249             : !===============================================================================
         250             :
         251      370944 :   subroutine aoa_tracers_timestep_init( phys_state )
         252             :     !-----------------------------------------------------------------------
         253             :     ! Provides a place to reinitialize diagnostic constituents HORZ and VERT
         254             :     !-----------------------------------------------------------------------
         255             :
         256        1536 :     use ppgrid,         only: begchunk, endchunk
         257             :     use physics_types,  only: physics_state
         258             :
         259             :     type(physics_state), intent(inout), dimension(begchunk:endchunk), optional :: phys_state
         260             :
         261             :     integer c, i, k, ncol
         262             :     integer yr, mon, day, tod,  ymd
         263             :     real(r8) :: calday, dpy
         264             :     !--------------------------------------------------------------------------
         265             :
         266      370944 :     if (.not. aoa_tracers_flag) return
         267             :
         268           0 :     call get_curr_date (yr,mon,day,tod)
         269             :
         270           0 :     if ( day == 1 .and. tod == 0) then
         271           0 :        if (masterproc) then
         272           0 :          write(iulog,*) 'AGE_OF_AIR_CONSTITUENTS: RE-INITIALIZING HORZ/VERT CONSTITUENTS'
         273             :        endif
         274             :
         275           0 :        do c = begchunk, endchunk
         276           0 :           ncol = phys_state(c)%ncol
         277           0 :           do k = 1, pver
         278           0 :              do i = 1, ncol
         279           0 :                 phys_state(c)%q(i,k,ixht) = 2._r8 + sin(phys_state(c)%lat(i))
         280           0 :                 phys_state(c)%q(i,k,ixvt) = qrel_vert(k)
         281             :              end do
         282             :           end do
         283             :        end do
         284             :
         285             :     end if
         286             :
         287           0 :     ymd = yr*10000 + mon*100 + day
         288           0 :     calday = get_calday(ymd, tod)
         289             :
         290           0 :     dpy = 365._r8
         291           0 :     if (timemgr_get_calendar_cf() == 'gregorian' .and. is_leapyear(yr)) then
         292           0 :        dpy = 366._r8
         293             :     end if
         294           0 :     years = (yr-yr0) + (calday-calday0)/dpy
         295             :
         296      741888 :   end subroutine aoa_tracers_timestep_init
         297             :
         298             : !===============================================================================
         299             :
         300     5956704 :   subroutine aoa_tracers_timestep_tend(state, ptend, dt)
         301             :
         302      370944 :     use physics_types, only: physics_state, physics_ptend, physics_ptend_init
         303             :     use cam_history,   only: outfld
         304             :
         305             :     ! Arguments
         306             :     type(physics_state), intent(in)    :: state              ! state variables
         307             :     type(physics_ptend), intent(out)   :: ptend              ! package tendencies
         308             :     real(r8),            intent(in)    :: dt                 ! timestep size (sec)
         309             :
         310             :     !----------------- Local workspace-------------------------------
         311             :
         312             :     integer :: i, k
         313             :     integer :: lchnk                          ! chunk identifier
         314             :     integer :: ncol                           ! no. of column in chunk
         315             :     real(r8) :: qrel                          ! value to be relaxed to
         316             :     real(r8) :: xhorz                         ! updated value of HORZ
         317             :     real(r8) :: xvert                         ! updated value of VERT
         318             :     logical  :: lq(pcnst)
         319             :     real(r8) :: teul                          ! relaxation in  1/sec*dt/2 = k*dt/2
         320             :     real(r8) :: wimp                          !     1./(1.+ k*dt/2)
         321             :     real(r8) :: wsrc                          !  teul*wimp
         322             :
         323             :     real(r8) :: xmmr
         324             :     real(r8), parameter :: mmr0 = 1.0e-6_r8   ! initial lower boundary mmr
         325             :     real(r8), parameter :: per_yr = 0.02_r8   ! fractional increase per year
         326             :
         327             :     !------------------------------------------------------------------
         328             :
         329     1489176 :     teul = .5_r8*dt/(86400._r8 * treldays)   ! 1/2 for the semi-implicit scheme if dt=time step
         330     1489176 :     wimp = 1._r8/(1._r8 +teul)
         331     1489176 :     wsrc = teul*wimp
         332             :
         333     1489176 :     if (.not. aoa_tracers_flag) then
         334     1489176 :        call physics_ptend_init(ptend,state%psetcols,'none') !Initialize an empty ptend for use with physics_update
         335             :        return
         336             :     end if
         337             :
         338           0 :     lq(:)     = .FALSE.
         339           0 :     lq(ixaoa) = .TRUE.
         340           0 :     lq(ixht)  = .TRUE.
         341           0 :     lq(ixvt)  = .TRUE.
         342             :
         343           0 :     call physics_ptend_init(ptend,state%psetcols, 'aoa_tracers', lq=lq)
         344             :
         345           0 :     lchnk = state%lchnk
         346           0 :     ncol  = state%ncol
         347             :
         348             :     ! AOAMF
         349           0 :     xmmr = mmr0*(1._r8 + per_yr*years)
         350           0 :     ptend%q(1:ncol,pver,ixaoa) = (xmmr - state%q(1:ncol,pver,ixaoa)) / dt
         351             :
         352           0 :     do k = 1, pver
         353           0 :        do i = 1, ncol
         354             :
         355             :           ! HORZ
         356           0 :           qrel              = 2._r8 + sin(state%lat(i))          ! qrel  should zonal mean
         357           0 :           xhorz             = state%q(i,k,ixht)*wimp + wsrc*qrel ! Xnew = weight*3D-tracer + (1.-weight)*1D-tracer
         358           0 :           ptend%q(i,k,ixht) = (xhorz - state%q(i,k,ixht)) / dt   ! Xnew = weight*3D-tracer + (1.-weight)*2D-tracer  zonal mean
         359             :                                                                  !  Can be still used .... to diagnose fluxes OT-tracers
         360             :           ! VERT
         361           0 :           qrel              = qrel_vert(k)                       ! qrel  should zonal mean
         362           0 :           xvert             = wimp*state%q(i,k,ixvt) + wsrc*qrel
         363           0 :           ptend%q(i,k,ixvt) = (xvert - state%q(i,k,ixvt)) / dt
         364             :
         365             :        end do
         366             :
         367             :     end do
         368             :
         369             :     ! record tendencies on history files
         370           0 :     call outfld (src_names(1), ptend%q(:,:,ixaoa),  pcols, lchnk)
         371           0 :     call outfld (src_names(2), ptend%q(:,:,ixht),   pcols, lchnk)
         372           0 :     call outfld (src_names(3), ptend%q(:,:,ixvt),   pcols, lchnk)
         373             :
         374     1489176 :   end subroutine aoa_tracers_timestep_tend
         375             :
         376             : !===========================================================================
         377             :
         378           0 :   subroutine init_cnst_3d(m, latvals, lonvals, mask, q)
         379             :
         380             :     integer,  intent(in)  :: m          ! global constituent index
         381             :     real(r8), intent(in)  :: latvals(:) ! lat in degrees (ncol)
         382             :     real(r8), intent(in)  :: lonvals(:) ! lon in degrees (ncol)
         383             :     logical,  intent(in)  :: mask(:)    ! Only initialize where .true.
         384             :     real(r8), intent(out) :: q(:,:)     ! kg tracer/kg dry air (gcol,plev)
         385             :
         386             :     integer :: j, k, gsize
         387             :     !-----------------------------------------------------------------------
         388             :
         389           0 :     if (masterproc) then
         390           0 :        write(iulog,*) 'AGE-OF-AIR CONSTITUENTS: INITIALIZING ',cnst_name(m),m
         391             :     end if
         392             :
         393           0 :     if (m == ixaoa) then
         394             :
         395             :        ! AOAMF
         396           0 :        q(:,:) = 0.0_r8
         397             :
         398           0 :     else if (m == ixht) then
         399             :
         400             :        ! HORZ
         401           0 :        gsize = size(q, 1)
         402           0 :        do j = 1, gsize
         403           0 :           q(j,:) = 2._r8 + sin(latvals(j))
         404             :        end do
         405             :
         406           0 :     else if (m == ixvt) then
         407             :
         408             :        ! VERT
         409           0 :        do k = 1, pver
         410           0 :           do j = 1, size(q,1)
         411           0 :              q(j,k) = qrel_vert(k)
         412             :           end do
         413             :        end do
         414             :
         415             :     end if
         416             :
         417     1489176 :   end subroutine init_cnst_3d
         418             :
         419             : !=====================================================================
         420             :
         421             : end module aoa_tracers

Generated by: LCOV version 1.14