LCOV - code coverage report
Current view: | top level - physics/cam - tracers.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 19 | 200 | 9.5 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 3 | 8 | 37.5 %

              Line data    Source code
           1             : !======================================================================
           2             : ! Module implements passive test tracers.
           3             : !
           4             : ! Two options:
           5             : !
           6             : ! 1) Specify only the number of tracers desired by setting the -nadv_tt option to configure.
           7             : !    This results in setting up the desired number of tracers making use of the tracers_suite
           8             : !    module to generate tracer names and initialize mixing ratios.
           9             : !
          10             : ! 2) Specify both the number of tracers using configure's -nadv_tt option, and specify
          11             : !    the same number of tracer names using the test_tracer_names namelist variable.  This
          12             : !    specifies a set of passive tracers that are either initialized from analytic expressions
          13             : !    or by reading values from the IC file.  The tracers for which analytic expressions are
          14             : !    available are the following:
          15             : !
          16             : !    test_tracer_names     description
          17             : !    -----------------     -----------
          18             : !    TT_SLOT               Non-smooth scalar field (slotted cylinder)
          19             : !    TT_SLOT1              1st slotted cylinder for 3 tracer test (Lauritzen and Thuburn, 2012, QJRMS)
          20             : !    TT_SLOT2              2nd slotted cylinder for 3 tracer test (Lauritzen and Thuburn, 2012, QJRMS)
          21             : !    TT_SLOT3              3rd slotted cylinder for 3 tracer test: TT_SLOT3=1-(TT_SLOT1+TT_SLOT2)
          22             : !    TT_GBALL              Smooth Gaussian "ball"
          23             : !    TT_TANH               Zonally constant, tanh function of latitude.
          24             : !    TT_EM8                Constant field of size 1.e-8.
          25             : !    TT_Y2_2               Approximately Y^2_2 spherical harmonic
          26             : !    TT_Y32_16             Approximately Y^32_16 spherical harmonic
          27             : !    TT_LATP2              Zonally constant, latitude + 2.
          28             : !    TT_LONP2              Meridionally constant, cos(longitude) + 2.
          29             : !    TT_COSB               Cosine bell (Kent et al., 2012, MWR)
          30             : !    TT_CCOSB              Non-linearly correlated Cosine bell (Lauritzen and Thuburn,2012, QJRMS)
          31             : !    TT_lCCOSB             Linearly correlated Cosine bell
          32             : !
          33             : !======================================================================
          34             :
          35             : module tracers
          36             :
          37             : use shr_kind_mod,    only: r8 => shr_kind_r8
          38             : use shr_sys_mod,     only: shr_sys_flush
          39             : use spmd_utils,      only: masterproc
          40             : use ppgrid,          only: pver
          41             : use physconst,       only: mwdry, cpair
          42             : use constituents,    only: cnst_add, cnst_name, cnst_longname
          43             : use tracers_suite,   only: get_tracer_name, init_cnst_tr
          44             : use cam_history,     only: addfld, add_default
          45             : use cam_logfile,     only: iulog
          46             : use cam_abortutils,  only: endrun
          47             :
          48             : implicit none
          49             : private
          50             : save
          51             :
          52             : public :: &
          53             :    tracers_readnl,            &! read namelist
          54             :    tracers_register,          &! register constituent
          55             :    tracers_implements_cnst,   &! true if named constituent is implemented by this package
          56             :    tracers_init_cnst,         &! initialize constituent field
          57             :    tracers_init               ! initialize history fields, datasets
          58             :
          59             : integer, parameter :: num_names_max = 30
          60             : integer, parameter :: num_analytic  = 14
          61             :
          62             : ! Data from namelist variables
          63             : integer           :: test_tracer_num = 0
          64             : character(len=16) :: test_tracer_names(num_names_max)
          65             :
          66             : logical :: tracers_flag       = .false.  ! true => turn on test tracer code
          67             : logical :: tracers_suite_flag = .false.  ! true => test tracers provided by tracers_suite module
          68             :
          69             : integer :: ixtrct=-999                   ! index of 1st constituent
          70             :
          71             : character(len=16), parameter :: analytic_names(num_analytic) =     &
          72             :      (/'TT_SLOT1        ', 'TT_SLOT2        ', 'TT_SLOT3        ', &
          73             :        'TT_SLOT         ', 'TT_GBALL        ', 'TT_TANH         ', &
          74             :        'TT_EM8          ', 'TT_Y2_2         ', 'TT_Y32_16       ', &
          75             :        'TT_LATP2        ', 'TT_LONP2        ', 'TT_COSB         ', &
          76             :        'TT_CCOSB        ', 'TT_lCCOSB       '/)
          77             :
          78             :
          79             : logical :: analytic_tracer(num_names_max)
          80             :
          81             : !======================================================================
          82             : contains
          83             : !======================================================================
          84             :
          85        1536 : subroutine tracers_readnl(nlfile)
          86             :
          87             :    use namelist_utils,  only: find_group_name
          88             :    use units,           only: getunit, freeunit
          89             :    use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_character
          90             :
          91             :    ! args
          92             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          93             :
          94             :    ! Local variables
          95             :    integer :: unitn, ierr
          96             :    integer :: i, j
          97             :    integer :: num_names
          98             :    character(len=*), parameter :: subname = 'tracers_readnl'
          99             :
         100             :    namelist /test_tracers_nl/ test_tracer_num, test_tracer_names
         101             :    !-----------------------------------------------------------------------------
         102             :
         103       47616 :    test_tracer_names = (/ (' ', i=1,num_names_max) /)
         104             :
         105        1536 :    if (masterproc) then
         106           2 :       unitn = getunit()
         107           2 :       open( unitn, file=trim(nlfile), status='old' )
         108           2 :       call find_group_name(unitn, 'test_tracers_nl', status=ierr)
         109           2 :       if (ierr == 0) then
         110           0 :          read(unitn, test_tracers_nl, iostat=ierr)
         111           0 :          if (ierr /= 0) then
         112           0 :             call endrun(subname // ':: ERROR reading namelist')
         113             :          end if
         114             :       end if
         115           2 :       close(unitn)
         116           2 :       call freeunit(unitn)
         117             :    end if
         118             :
         119             :    call mpi_bcast(test_tracer_names, len(test_tracer_names)*num_names_max, mpi_character, &
         120        1536 :                   mstrid, mpicom, ierr)
         121        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: test_tracer_names")
         122             :
         123        1536 :    call mpi_bcast(test_tracer_num,   1, mpi_integer, mstrid, mpicom, ierr)
         124        1536 :    if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: iradsw")
         125             :
         126             :    ! If any tracers have been specified then turn on the tracers module
         127        1536 :    if (test_tracer_num > 0) then
         128           0 :       tracers_flag = .true.
         129             :    else
         130        1536 :       return
         131             :    end if
         132             :
         133             :    ! Determine the number of tracer names supplied:
         134           0 :    num_names = 0
         135           0 :    analytic_tracer = .false.
         136           0 :    do i = 1, num_names_max
         137           0 :       if (len_trim(test_tracer_names(i)) > 0) then
         138           0 :          num_names = num_names + 1
         139             :
         140             :          ! Does the tracer have an analytic IC?
         141           0 :          do j = 1, num_analytic
         142           0 :             if (trim(test_tracer_names(i)) == trim(analytic_names(j))) then
         143           0 :                analytic_tracer(i) = .true.
         144           0 :                exit
         145             :             end if
         146             :          end do
         147             :       else
         148             :          exit
         149             :       end if
         150             :    end do
         151             :
         152           0 :    if (num_names > 0) then
         153             :       ! If test_tracer_names have been specified, the test_tracer_num should
         154             :       ! equal the number of names supplied.
         155           0 :       if (num_names /= test_tracer_num) then
         156           0 :          write(iulog, *) subname//' number of names, number of tracers: ', num_names, test_tracer_num
         157           0 :          call endrun(subname // ':: number of names does not match number of tracers')
         158             :       end if
         159             :    else
         160             :       ! If no names have been supplied then
         161             :       ! the tracers will be provided by the tracers_suite module.
         162           0 :       tracers_suite_flag = .true.
         163             :    end if
         164             :
         165             :    ! Print summary to log file
         166           0 :    if (masterproc) then
         167             :
         168           0 :       write(iulog, *) 'Test Tracers Module'
         169           0 :       write(iulog, *) '  Number of Test Tracers:', test_tracer_num
         170           0 :       if (tracers_suite_flag) then
         171           0 :          write(iulog, *) '  Tracers will be provided by tracers_suite module.'
         172             :       else
         173           0 :          do i = 1, num_names
         174           0 :             if (analytic_tracer(i)) then
         175             :                write(iulog, *) '  '//trim(test_tracer_names(i))//&
         176           0 :                   ' will be initialized from an analytic expression'
         177             :             else
         178             :                write(iulog, *) '  '//trim(test_tracer_names(i))//&
         179           0 :                   ' will be initialized from the IC file'
         180             :             end if
         181             :          end do
         182             :       end if
         183             :    end if
         184             :
         185             : end subroutine  tracers_readnl
         186             :
         187             : !======================================================================
         188             :
         189        1536 : subroutine tracers_register()
         190             :
         191             :    ! Register advected tracers.
         192             :
         193             :    ! Local variables
         194             :    integer  :: m, mm
         195             :    logical  :: read_from_file
         196             :    real(r8) :: minc
         197             :    character(len=16) :: name
         198             :    !-----------------------------------------------------------------------
         199             :
         200        1536 :    if (.not. tracers_flag) return
         201             :
         202           0 :    minc = -1.e36_r8   ! min mixing ratio (disable qneg3)
         203             :
         204           0 :    do m = 1, test_tracer_num
         205             :
         206           0 :       read_from_file = .true.
         207           0 :       if (tracers_suite_flag) then
         208           0 :          name = get_tracer_name(m)  ! get name from suite file
         209           0 :          read_from_file = .false.
         210             :       else
         211           0 :          name = test_tracer_names(m)
         212           0 :          if (analytic_tracer(m)) read_from_file = .false.
         213             :       end if
         214             :
         215             :       ! add constituent name to list of advected, save index number ixtrct
         216             :       call cnst_add(name, mwdry, cpair, minc, mm, &
         217           0 :                     readiv=read_from_file, mixtype='dry')
         218           0 :       if (m == 1) ixtrct = mm  ! save index number of first tracer
         219             :    end do
         220             :
         221             : end subroutine tracers_register
         222             :
         223             : !======================================================================
         224             :
         225           0 : function tracers_implements_cnst(name)
         226             :
         227             :    ! return true if specified constituent is implemented by this package
         228             :
         229             :    ! Arguments
         230             :    character(len=*), intent(in) :: name   ! constituent name
         231             :    logical :: tracers_implements_cnst        ! return value
         232             :
         233             :    ! Local variables
         234             :    integer :: m
         235             :    character(len=16) :: trc_name
         236             :    !-----------------------------------------------------------------------
         237             :
         238           0 :    tracers_implements_cnst = .false.
         239           0 :    if (.not. tracers_flag) return
         240             :
         241           0 :    do m = 1, test_tracer_num
         242             :
         243           0 :       if (tracers_suite_flag) then
         244           0 :          trc_name = get_tracer_name(m)
         245             :       else
         246           0 :          trc_name = test_tracer_names(m)
         247             :       end if
         248             :
         249           0 :       if (name == trc_name) then
         250           0 :          tracers_implements_cnst = .true.
         251             :          return
         252             :       end if
         253             :    end do
         254             :
         255             : end function tracers_implements_cnst
         256             :
         257             : !===============================================================================
         258             :
         259           0 : subroutine tracers_init_cnst(name, latvals, lonvals, mask, q, z)
         260             :
         261             :    ! Initialize test tracer mixing ratio
         262             :
         263             :    character(len=*), intent(in)  :: name
         264             :    real(r8),         intent(in)  :: latvals(:) ! lat in radians (ncol)
         265             :    real(r8),         intent(in)  :: lonvals(:) ! lon in radians (ncol)
         266             :    logical,          intent(in)  :: mask(:)    ! Only initialize where .true.
         267             :    real(r8),         intent(out) :: q(:,:)     ! kg tracer/kg dry air (gcol, plev
         268             :    real(r8),optional,intent(in)  :: z(:,:)     ! height of full level pressure
         269             :
         270             :    ! Local
         271             :    integer :: m
         272             :    logical :: found
         273             :    character(len=*), parameter :: subname = 'tracers_init_cnst'
         274             :    !----------------------------------------------------------------------------
         275             :
         276           0 :    if (.not. tracers_flag) return
         277             :
         278           0 :    found = .false.
         279           0 :    if (tracers_suite_flag) then
         280             :
         281           0 :      do m = 1, test_tracer_num
         282           0 :        if (name ==  get_tracer_name(m))  then
         283           0 :          call init_cnst_tr(m, latvals, lonvals, mask, q)
         284           0 :          found = .true.
         285           0 :          exit
         286             :        endif
         287             :      end do
         288             :
         289             :    else
         290             :
         291           0 :      do m = 1, test_tracer_num
         292           0 :        if (name == test_tracer_names(m)) then
         293           0 :          if (analytic_tracer(m)) then
         294           0 :            if (present(z)) then
         295           0 :              call test_func_set(name, latvals, lonvals, mask, q, z=z)
         296             :            else
         297           0 :              call test_func_set(name, latvals, lonvals, mask, q)
         298             :            end if
         299             :            found = .true.
         300             :            exit
         301             :          else
         302             :            ! The initial values were supposed to be read from the IC file.  This
         303             :            ! call should not have been made in that case, so it appears that a requested
         304             :            ! tracer is not on the IC file.
         305           0 :            write(iulog, *) subname//': ERROR: tracer ', trim(name), ' should be on IC file'
         306           0 :            call endrun(subname//': ERROR: tracer missing from IC file')
         307             :          end if
         308             :        end if
         309             :      end do
         310             :    end if
         311             :
         312             :    if (.not. found) then
         313             :      ! unrecognized tracer name
         314           0 :      write(iulog, *) subname//': ERROR: ', trim(name), ' not recognized'
         315           0 :      call endrun(subname//': ERROR: tracer name not recognized')
         316             :    end if
         317             :
         318             : end subroutine tracers_init_cnst
         319             :
         320             : !===============================================================================
         321             :
         322        1536 : subroutine tracers_init()
         323             :
         324             :    ! Add tracers to history output
         325             :
         326             :    ! Local
         327             :    integer m, mm
         328             :    character(len=16) :: name   ! constituent name
         329             :
         330        1536 :    if (.not. tracers_flag ) return
         331             :
         332           0 :    do m = 1,test_tracer_num
         333           0 :       mm = ixtrct + m - 1
         334           0 :       call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm))
         335           0 :       call add_default(cnst_name(mm), 1, ' ')
         336             :    end do
         337             : end subroutine tracers_init
         338             :
         339             : !=========================================================================================
         340             :
         341           0 : subroutine test_func_set(name, latvals, lonvals, mask, q, z)
         342             :
         343             :    ! use test_func code to set array q
         344             :
         345             :    character(len=*), intent(in)  :: name
         346             :    real(r8),         intent(in)  :: latvals(:) ! lat in radians (ncol)
         347             :    real(r8),         intent(in)  :: lonvals(:) ! lon in radians (ncol)
         348             :    logical,          intent(in)  :: mask(:)    ! Only initialize where .true.
         349             :    real(r8),         intent(out) :: q(:,:)     ! kg tracer/kg dry air (gcol, plev
         350             :    real(r8),optional,intent(in)  :: z(:,:)     ! height of full level pressure
         351             :
         352             :    ! local variables
         353             :    integer :: i, k
         354             :    !----------------------------------------------------------------------------
         355             :
         356           0 :    do i = 1, size(mask)
         357           0 :      if (mask(i)) then
         358           0 :        if (present(z)) then
         359           0 :          do k = 1, size(q, 2)
         360           0 :            q(i,k) = test_func(name, latvals(i), lonvals(i), k, z=z(i,k))
         361             :          end do
         362             :        else
         363           0 :          do k = 1, size(q, 2)
         364           0 :            q(i,k) = test_func(name, latvals(i), lonvals(i), k)
         365             :          end do
         366             :        end if
         367             :      end if
         368             :    end do
         369             :
         370           0 : end subroutine test_func_set
         371             :
         372             : !=========================================================================================
         373             :
         374           0 : function test_func(name, lat, lon, k, z) result(fout)
         375             :
         376             :    ! Analytic test functions.
         377             :
         378             :    use physconst, only: pi
         379             :    use hycoef,    only: hyam, hybm, ps0
         380             :
         381             :    character(len=*), intent(in) :: name
         382             :    real(r8),         intent(in) :: lon     ! radians
         383             :    real(r8),         intent(in) :: lat     ! radians
         384             :    integer,          intent(in) :: k       ! vertical index for layer mid-points
         385             :    real(r8),optional,intent(in) :: z       ! height
         386             :
         387             :    real(r8) :: fout
         388             :
         389             :    real(r8), parameter :: psurf_moist = 100000.0_r8 ! moist surface pressure
         390             :    real(r8), parameter :: deg2rad = pi/180._r8
         391             :    real(r8), parameter :: z0 = 4500.0_r8 !Eq 9 in Kent et al. (2012)
         392             :    real(r8), parameter :: zz = 1000.0_r8  !Eq 9 in Kent et al. (2012)
         393             :
         394             :    real(r8) :: lon1, lat1, R0, Rg1
         395             :    real(r8) :: eta, eta_c, d1, f1, f2, c
         396             :    real(r8) :: cos_tmp, sin_tmp
         397             :
         398             :    !----------------------------------------------------------------------------
         399             :
         400           0 :    select case(name)
         401             :    case('TT_SLOT')
         402             :      !
         403             :      !   Non-smooth scalar field (slotted cylinder): min=0.1 and max=1
         404             :      !
         405           0 :      R0 = 0.25_r8
         406           0 :      lon1  = pi/9.0_r8
         407           0 :      lat1  = 2.0_r8*pi/9.0_r8
         408           0 :      Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
         409           0 :      c = 1.0_r8
         410             :
         411           0 :      if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
         412             :        fout = c
         413             :      elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
         414           0 :           .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
         415             :        fout = c
         416             :      else
         417           0 :        fout = 0.1_r8
         418             :      endif
         419             :    case('TT_SLOT1')
         420             :      !
         421             :      !   Slotted cylinder for 3 tracer test
         422             :      !
         423           0 :      R0 = 0.25_r8
         424           0 :      lon1  = pi/9.0_r8
         425           0 :      lat1  = 2.0_r8*pi/9.0_r8
         426           0 :      Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
         427           0 :      c = 1.0_r8/3.0_r8
         428             :
         429           0 :      if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
         430             :        fout = c
         431             :      elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
         432           0 :           .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
         433             :        fout = c
         434             :      else
         435           0 :        fout = 0.1_r8
         436             :      endif
         437             :    case('TT_SLOT2')
         438             :      !
         439             :      !   Slotted cylinder for 3 tracer test
         440             :      !
         441           0 :      R0 = 0.25_r8
         442           0 :      lon1  = pi/9.0_r8
         443           0 :      lat1  = 2.0_r8*pi/9.0_r8+pi/18.0_r8
         444           0 :      Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
         445           0 :      c = 2.0_r8/3.0_r8
         446             :
         447           0 :      if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
         448             :        fout = c
         449             :      elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
         450           0 :           .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
         451             :        fout = c
         452             :      else
         453           0 :        fout = 0.1_r8
         454             :      endif
         455             :    case('TT_SLOT3')
         456             :      !
         457             :      !   Slotted cylinder for 3 tracer test
         458             :      !   Residual of two slotted cylinders so that TT_SLOT1+TT_SLOT2+TT_SLOT3=1
         459             :      !   Lauritzen and Thuburn (2012, QJRMS)
         460             :      !
         461           0 :      R0 = 0.25_r8
         462           0 :      lon1  = pi/9.0_r8
         463           0 :      lat1  = 2.0_r8*pi/9.0_r8
         464           0 :      Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
         465           0 :      c = 1.0_r8/3.0_r8
         466             :
         467           0 :      if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
         468           0 :        f1 = c
         469             :      elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
         470           0 :           .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
         471           0 :        f1 = c
         472             :      else
         473           0 :        f1 = 0.1_r8
         474             :      endif
         475             :
         476           0 :      R0 = 0.25_r8
         477           0 :      lon1  = pi/9.0_r8
         478           0 :      lat1  = 2.0_r8*pi/9.0_r8+pi/18.0_r8
         479           0 :      Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
         480           0 :      c = 2.0_r8/3.0_r8
         481           0 :      if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
         482             :        f2 = c
         483             :      elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
         484           0 :           .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
         485             :        f2 = c
         486             :      else
         487           0 :        f2 = 0.1_r8
         488             :      endif
         489             :
         490           0 :      fout = 1.0_r8-(f1+f2)
         491             :    case('TT_GBALL')
         492             :      !
         493             :      ! Smooth Gaussian "ball"
         494             :      !
         495           0 :      R0    = 10.0_r8           ! radius of the perturbation
         496           0 :      lon1  = 20.0_r8*deg2rad   ! longitudinal position, 20E
         497           0 :      lat1  = 40.0_r8*deg2rad  ! latitudinal position, 40N
         498           0 :      eta_c = 0.6_r8
         499           0 :      sin_tmp = SIN(lat1)*SIN(lat)
         500           0 :      cos_tmp = COS(lat1)*COS(lat)
         501             :
         502           0 :      Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) )    ! great circle distance
         503           0 :      eta =  (hyam(k)*ps0 + hybm(k)*psurf_moist)/psurf_moist
         504           0 :      fout = EXP(- ((Rg1*R0)**2 + ((eta-eta_c)/0.1_r8)**2))
         505           0 :      IF (ABS(fout) < 1.0E-8_r8) fout = 0.0_r8
         506             :
         507             :    case('TT_TANH')
         508             :      !
         509             :      !
         510             :      !
         511           0 :      fout = 0.5_r8 * ( tanh( 3.0_r8*abs(lat)-pi ) + 1.0_r8)
         512             :
         513             :    case('TT_EM8')
         514           0 :      fout = 1.0e-8_r8
         515             :
         516             :    case('TT_Y2_2')
         517             :      !
         518             :      ! approximately Y^2_2 spherical harmonic
         519             :      !
         520           0 :      fout = 0.5_r8 + 0.5_r8*(cos(lat)*cos(lat)*cos(2.0_r8*lon))
         521             :
         522             :    case('TT_Y32_16')
         523             :      !
         524             :      ! approximately Y32_16 spherical harmonic
         525             :      !
         526           0 :      fout = 0.5_r8 + 0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
         527             :
         528             :    case('TT_LATP2')
         529           0 :      fout = 2.0_r8 + lat
         530             :
         531             :    case('TT_LONP2')
         532           0 :      fout = 2.0_r8 + cos(lon)
         533             :    case('TT_COSB')
         534             :      !
         535             :      ! Cosine bell inspired by Kent et al., 2012, MWR; only one bell and location changed
         536             :      ! https://journals.ametsoc.org/doi/pdf/10.1175/MWR-D-11-00150.1
         537             :      !
         538           0 :      R0    = 0.5_r8           ! radius of the perturbation
         539           0 :      lon1  = pi/9.0_r8
         540           0 :      lat1  = 2.0_r8*pi/9.0_r8
         541             :
         542           0 :      sin_tmp = SIN(lat1)*SIN(lat)
         543           0 :      cos_tmp = COS(lat1)*COS(lat)
         544           0 :      Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) )    ! great circle distance
         545             :
         546           0 :      if (present(z)) THEN
         547           0 :        d1 = MIN(1.0_r8,(Rg1/R0)**2+((z-z0)/zz)**2)
         548             :      else
         549           0 :        d1 = MIN(1.0_r8,(Rg1/R0)**2)
         550             :      end if
         551             :
         552           0 :      if (Rg1 < R0) then
         553           0 :        fout = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1))
         554             :      else
         555             :        fout = 0.1_r8
         556             :      end if
         557             :    case('TT_CCOSB')
         558             :      !
         559             :      ! Correlated cosine bell
         560             :      !
         561           0 :      R0    = 0.5_r8           ! radius of the perturbation
         562           0 :      lon1  = pi/9.0_r8
         563           0 :      lat1  = 2.0_r8*pi/9.0_r8
         564             :
         565           0 :      sin_tmp = SIN(lat1)*SIN(lat)
         566           0 :      cos_tmp = COS(lat1)*COS(lat)
         567           0 :      Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) )    ! great circle distance
         568             :
         569           0 :      if (present(z)) THEN
         570           0 :        d1 = MIN(1.0_r8,(Rg1/R0)**2+((z-z0)/zz)**2)
         571             :      else
         572           0 :        d1 = MIN(1.0_r8,(Rg1/R0)**2)
         573             :      end if
         574             :
         575           0 :      if (Rg1 < R0) then
         576           0 :        f1 = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1))
         577             :      else
         578           0 :        f1 = 0.1_r8
         579             :      end if
         580             :
         581           0 :      fout=corr_fct(f1)
         582             :    case('TT_lCCOSB')
         583             :      !
         584             :      ! Correlated cosine bell
         585             :      !
         586           0 :      R0    = 0.5_r8           ! radius of the perturbation
         587           0 :      lon1  = pi/9.0_r8
         588           0 :      lat1  = 2.0_r8*pi/9.0_r8
         589             :
         590           0 :      sin_tmp = SIN(lat1)*SIN(lat)
         591           0 :      cos_tmp = COS(lat1)*COS(lat)
         592           0 :      Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) )    ! great circle distance
         593             :
         594           0 :      if (present(z)) THEN
         595           0 :        d1 = MIN(1.0_r8,(Rg1/R0)**2+((z-z0)/zz)**2)
         596             :      else
         597           0 :        d1 = MIN(1.0_r8,(Rg1/R0)**2)
         598             :      end if
         599             :
         600           0 :      if (Rg1 < R0) then
         601           0 :        f1 = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1))
         602             :      else
         603           0 :        f1 = 0.1_r8
         604             :      end if
         605             :
         606           0 :      fout=1.0_r8-f1
         607             :
         608             :    case default
         609           0 :      call endrun("test_func: ERROR: name not recognized")
         610             :    end select
         611           0 : end function test_func
         612             :
         613             : !
         614             : ! correlation function (Lauritzen and Thuburn, 2012, QJRMS)
         615             : !
         616           0 : function corr_fct(x) result(out)
         617             :   real(r8), intent(in) :: x
         618             :   real(r8)             :: out
         619           0 :   out = -0.8_r8*x**2+0.9_r8
         620           0 : end function corr_fct
         621             : end module tracers

Generated by: LCOV version 1.14