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