LCOV - code coverage report
Current view: | top level - physics/cam - phys_debug_util.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 19 | 84 | 22.6 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 2 | 3 | 66.7 %

              Line data    Source code
           1             : module phys_debug_util
           2             :
           3             : !----------------------------------------------------------------------------------------
           4             :
           5             : ! Module to facilitate debugging of physics parameterizations.
           6             : !
           7             : ! The user requests a location for debugging in lat/lon coordinates
           8             : ! (degrees).  The initialization routine does a global search to find the
           9             : ! column in the physics grid closest to the requested location.  The local
          10             : ! indices of that column in the physics decomposition are stored as module
          11             : ! data.  The user code then passes the local chunk index of the chunked
          12             : ! data into the subroutine that will write diagnostic information for the
          13             : ! column.  The function phys_debug_col returns the local column index if
          14             : ! the column of interest is contained in the chunk, and zero otherwise.
          15             : ! Printing is done only if a column index >0 is returned.
          16             : !
          17             : ! Phil Rasch, B. Eaton, Feb 2008
          18             : !----------------------------------------------------------------------------------------
          19             :
          20             : use shr_kind_mod,   only: r8 => shr_kind_r8
          21             : use spmd_utils,     only: masterproc, iam, mpicom, npes
          22             : use cam_logfile,    only: iulog
          23             : use cam_abortutils, only: endrun
          24             :
          25             : implicit none
          26             : private
          27             : save
          28             :
          29             : real(r8), parameter :: uninitr8 = huge(1._r8)
          30             :
          31             : ! Public methods
          32             : public phys_debug_readnl  ! read namelist input
          33             : public phys_debug_init    ! initialize the method to a chunk and column
          34             : public phys_debug_col     ! return local column index in debug chunk
          35             :
          36             : ! Namelist variables
          37             : real(r8) :: phys_debug_lat = uninitr8 ! latitude of requested debug column location in degrees
          38             : real(r8) :: phys_debug_lon = uninitr8 ! longitude of requested debug column location in degrees
          39             :
          40             :
          41             : integer :: debchunk = -999            ! local index of the chuck we will debug
          42             : integer :: debcol   = -999            ! the column within the chunk we will debug
          43             :
          44             : !================================================================================
          45             : contains
          46             : !================================================================================
          47             :
          48        1536 : subroutine phys_debug_readnl(nlfile)
          49             :
          50             :    use namelist_utils,  only: find_group_name
          51             :    use units,           only: getunit, freeunit
          52             :    use mpishorthand
          53             :
          54             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          55             :
          56             :    ! Local variables
          57             :    integer :: unitn, ierr
          58             :    character(len=*), parameter :: subname = 'phys_debug_readnl'
          59             :
          60             :    namelist /phys_debug_nl/ phys_debug_lat, phys_debug_lon
          61             :    !-----------------------------------------------------------------------------
          62             :
          63        1536 :    if (masterproc) then
          64           2 :       unitn = getunit()
          65           2 :       open( unitn, file=trim(nlfile), status='old' )
          66           2 :       call find_group_name(unitn, 'phys_debug_nl', status=ierr)
          67           2 :       if (ierr == 0) then
          68           0 :          read(unitn, phys_debug_nl, iostat=ierr)
          69           0 :          if (ierr /= 0) then
          70           0 :             call endrun(subname // ':: ERROR reading namelist')
          71             :          end if
          72             :       end if
          73           2 :       close(unitn)
          74           2 :       call freeunit(unitn)
          75             :       ! Check inputs
          76           2 :       if (phys_debug_lat /= uninitr8) then
          77           0 :          if (abs(phys_debug_lat) > 90.0_r8) then
          78           0 :             write(iulog, *) subname, ': phys_debug_lat out of range [-90., 90.]'
          79           0 :             call endrun(subname//': phys_debug_lat out of range [-90., 90.]')
          80             :          end if
          81             :       else
          82           2 :          write(iulog, *) subname, ': phys_debug_lat = ', phys_debug_lat
          83             :       end if
          84           2 :       if (phys_debug_lon /= uninitr8) then
          85           0 :          if ((phys_debug_lon < 0.0_r8) .or. (phys_debug_lon > 360.0_r8)) then
          86           0 :             write(iulog, *) subname, ': phys_debug_lon out of range [0., 360.]'
          87           0 :             call endrun(subname//': phys_debug_lon out of range [0., 360.]')
          88             :          end if
          89             :       else
          90           2 :          write(iulog, *) subname, ': phys_debug_lon = ', phys_debug_lon
          91             :       end if
          92             :    end if
          93             :
          94             : #ifdef SPMD
          95             :    ! Broadcast namelist variables
          96        1536 :    call mpibcast(phys_debug_lat, 1, mpir8, 0, mpicom)
          97        1536 :    call mpibcast(phys_debug_lon, 1, mpir8, 0, mpicom)
          98             : #endif
          99             :
         100        1536 : end subroutine phys_debug_readnl
         101             :
         102             : !==============================================================================
         103             :
         104        1536 : subroutine phys_debug_init()
         105             :    use mpi,       only: mpi_real8, mpi_integer, mpi_min, mpi_max
         106             :    use physconst, only: pi
         107             :    use ppgrid,    only: begchunk, endchunk
         108             :    use phys_grid, only: get_ncols_p, get_rlat_p, get_rlon_p
         109             :
         110             :    integer             :: owner, lchunk, icol, ncol
         111             :    integer             :: lchunk_min, icol_min, minlondist
         112             :    real(r8)            :: deblat, deblon
         113             :    real(r8)            :: latmin, lonmin
         114             :    real(r8)            :: lat, lon, dist, temp1, temp2
         115             :    real(r8)            :: mindist
         116             :    real(r8), parameter :: maxangle = pi / 4.0_r8
         117             :    real(r8), parameter :: rad2deg = 180.0_r8 / pi
         118             :    real(r8), parameter :: deg2rad = pi / 180.0_r8
         119             :    real(r8), parameter :: maxtol = 0.99999_r8 ! max cos value
         120             :    real(r8), parameter :: maxlat = pi * maxtol / 2.0_r8
         121             :    !---------------------------------------------------------------------------
         122             :
         123             :    ! If no debug column specified then do nothing
         124        1536 :    if ((phys_debug_lat == uninitr8) .or. (phys_debug_lon == uninitr8)) then
         125             :       return
         126             :    end if
         127             :
         128             :    ! User has specified a column location for debugging.  Find the closest
         129             :    ! column in the physics grid.
         130           0 :    mindist = 2.0_r8 * pi
         131           0 :    deblat = pi
         132           0 :    deblon = 3.0_r8 * pi
         133           0 :    latmin = phys_debug_lat * deg2rad
         134           0 :    lonmin = phys_debug_lon * deg2rad
         135           0 :    lchunk_min = -1
         136           0 :    icol_min = -1
         137           0 :    do lchunk = begchunk, endchunk
         138           0 :       ncol = get_ncols_p(lchunk)
         139           0 :       do icol = 1, ncol
         140           0 :          lat = get_rlat_p(lchunk, icol)
         141           0 :          lon = get_rlon_p(lchunk, icol)
         142           0 :          if ( (abs(lat - latmin) <= maxangle) .and.                           &
         143           0 :               (abs(lon - lonmin) <= maxangle)) then
         144             :             ! maxangle could be pi but why waste all those trig functions?
         145           0 :             if ((lat == latmin) .and. (lon == lonmin)) then
         146           0 :                dist = 0.0_r8
         147             :             else
         148             :                temp1 = (sin(latmin) * sin(lat)) +                             &
         149           0 :                     (cos(latmin) * cos(lat) * cos(lon - lonmin))
         150           0 :                if (temp1 > maxtol) then
         151             :                   ! Use haversine formula
         152           0 :                   temp1 = sin((latmin - lat) / 2.0_r8)
         153           0 :                   temp2 = sin((lonmin - lon) / 2.0_r8)
         154             :                   dist = (temp1 * temp1) +                                    &
         155           0 :                        (cos(latmin)* cos(lat) * temp2 * temp2)
         156           0 :                   dist = 2.0_r8 * asin(sqrt(dist))
         157             :                else
         158           0 :                   dist = acos(temp1)
         159             :                end if
         160             :             end if
         161           0 :             if ( (dist < mindist) .or.                                        &
         162             :                  ((dist == mindist) .and.                                     &
         163             :                  (abs(lon - lonmin) < abs(deblon - lonmin)))) then
         164           0 :                lchunk_min = lchunk
         165           0 :                icol_min = icol
         166           0 :                mindist = dist
         167           0 :                deblon = lon
         168           0 :                deblat = lat
         169           0 :                if (dist == 0.0_r8) then
         170             :                   exit
         171             :                end if
         172             :             end if
         173             :          end if
         174             :       end do
         175             :    end do
         176             :    ! We need to find the minimum mindist and use only that value
         177           0 :    dist = mindist
         178           0 :    call MPI_allreduce(dist, mindist, 1, mpi_real8, mpi_min, mpicom, icol)
         179             :    ! Special case for pole points
         180           0 :    if (deblat > pi / 2.0_r8) then
         181           0 :       temp1 = 0.0_r8
         182             :    else
         183           0 :       temp1 = abs(deblat)
         184             :    end if
         185           0 :    call MPI_allreduce(temp1, lat, 1, mpi_real8, mpi_max, mpicom, icol)
         186           0 :    if ((abs(latmin) > maxlat) .or. (lat > maxlat)) then
         187           0 :       if (dist == mindist) then
         188             :          ! Only distance winners can compete
         189           0 :          lon = abs(deblon - lonmin)
         190             :       else
         191           0 :          lon = 3.0_r8 * pi
         192             :       end if
         193           0 :       call MPI_allreduce(lon, minlondist, 1, mpi_real8, mpi_min, mpicom, icol)
         194             :       ! Kill the losers
         195           0 :       if (lon /= minlondist) then
         196           0 :          dist = dist + 1.0_r8
         197             :       end if
         198             :    end if
         199             :    ! Now, only task(s) which have real minimum distance should be owner
         200           0 :    if (dist == mindist) then
         201           0 :       lchunk = iam
         202             :    else
         203           0 :       lchunk = npes + 2
         204             :    end if
         205           0 :    call MPI_allreduce(lchunk, owner, 1, mpi_integer, mpi_min, mpicom, icol)
         206             :    ! If the column is owned by this process then save its local indices
         207           0 :    if (iam == owner) then
         208           0 :       debchunk         = lchunk_min
         209           0 :       debcol           = icol_min
         210           0 :       deblat           = get_rlat_p(lchunk_min, icol_min) * rad2deg
         211           0 :       deblon           = get_rlon_p(lchunk_min, icol_min) * rad2deg
         212           0 :       write(iulog,*) 'phys_debug_init: debugging column at lat=', deblat, '  lon=', deblon
         213             :    end if
         214             :
         215        1536 : end subroutine phys_debug_init
         216             :
         217             : !================================================================================
         218             :
         219           0 : integer function phys_debug_col(chunk)
         220             :
         221             :    integer,  intent(in) :: chunk
         222             :    !-----------------------------------------------------------------------------
         223             :
         224           0 :    if (chunk == debchunk) then
         225           0 :       phys_debug_col = debcol
         226             :    else
         227             :       phys_debug_col = 0
         228             :    endif
         229             :
         230        1536 : end function phys_debug_col
         231             :
         232             : !================================================================================
         233             :
         234             : end module phys_debug_util

Generated by: LCOV version 1.14