LCOV - code coverage report
Current view: | top level - physics/cam - qneg_module.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 149 | 167 | 89.2 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 10 | 10 | 100.0 %

              Line data    Source code
           1             : module qneg_module
           2             :
           3             :   use shr_kind_mod,        only: r8 => shr_kind_r8, CS => SHR_KIND_CS
           4             :   use perf_mod,            only: t_startf, t_stopf
           5             :   use cam_logfile,         only: iulog
           6             :   use cam_abortutils,      only: endrun
           7             :   use shr_sys_mod,         only: shr_sys_flush
           8             :   use cam_history_support, only: max_fieldname_len
           9             :   use ppgrid,              only: pcols
          10             :   use constituents,        only: pcnst, cnst_name
          11             :
          12             :   implicit none
          13             :   private
          14             :   save
          15             :
          16             :   ! Public interface.
          17             :
          18             :   public :: qneg_readnl
          19             :   public :: qneg_init
          20             :   public :: qneg3
          21             :   public :: qneg4
          22             :   public :: qneg_print_summary
          23             :
          24             :   ! Private module variables
          25             :   character(len=8) :: print_qneg_warn
          26             :   logical          :: log_warnings = .false.
          27             :   logical          :: collect_stats = .false.
          28             :   logical          :: timestep_reset = .false.
          29             :
          30             :   real(r8), parameter :: tol = 1.e-12_r8
          31             :   real(r8), parameter :: worst_reset = 1.e35_r8
          32             :
          33             :   ! Diagnostic field names
          34             :   integer, parameter               :: num_diag_fields = (2 * pcnst) + 1
          35             :   character(len=max_fieldname_len) :: diag_names(num_diag_fields)
          36             :   logical             :: cnst_out_calc = .false.
          37             :   logical             :: cnst_outfld(num_diag_fields) = .false.
          38             :
          39             :   ! Summary buffers
          40             :   integer, parameter  :: num3_bins = 24
          41             :   integer, parameter  :: num4_bins = 4
          42             :   character(len=CS)   :: qneg3_warn_labels(num3_bins) = ''
          43             :   character(len=CS)   :: qneg4_warn_labels(num4_bins) = ''
          44             :   integer             :: qneg3_warn_num(pcnst, num3_bins)   = 0
          45             :   integer             :: qneg4_warn_num(num4_bins)          = 0
          46             :   real(r8)            :: qneg3_warn_worst(pcnst, num3_bins) = worst_reset
          47             :   real(r8)            :: qneg4_warn_worst(num4_bins)        = worst_reset
          48             :
          49             :   private             :: calc_cnst_out
          50             :   private             :: find_index3
          51             :   private             :: find_index4
          52             :   interface reset_stats
          53             :      module procedure reset_stats_scalar
          54             :      module procedure reset_stats_array
          55             :   end interface reset_stats
          56             :
          57             : contains
          58             :
          59        1536 :   subroutine qneg_readnl(nlfile)
          60             :     use namelist_utils,  only: find_group_name
          61             :     use units,           only: getunit, freeunit
          62             :     use spmd_utils,      only: mpicom, mstrid=>masterprocid, mpi_character, masterproc
          63             :     ! File containing namelist input.
          64             :     character(len=*), intent(in) :: nlfile
          65             :
          66             :     ! Local variables
          67             :     integer :: unitn, ierr
          68             :     character(len=*), parameter :: sub = 'qneg_readnl'
          69             :
          70             :     namelist /qneg_nl/ print_qneg_warn
          71             :
          72        1536 :     print_qneg_warn = ''
          73             :
          74        1536 :     if (masterproc) then
          75           2 :        unitn = getunit()
          76           2 :        open( unitn, file=trim(nlfile), status='old' )
          77           2 :        call find_group_name(unitn, 'qneg_nl', status=ierr)
          78           2 :        if (ierr == 0) then
          79           2 :           read(unitn, qneg_nl, iostat=ierr)
          80           2 :           if (ierr /= 0) then
          81           0 :              call endrun(sub // ':: ERROR reading namelist qneg_nl')
          82             :           end if
          83             :        end if
          84           2 :        close(unitn)
          85           2 :        call freeunit(unitn)
          86             :     end if
          87             :
          88        1536 :     call mpi_bcast(print_qneg_warn, len(print_qneg_warn), mpi_character, mstrid, mpicom, ierr)
          89        1536 :     if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: print_qneg_warn")
          90             :
          91        3072 :     select case(trim(print_qneg_warn))
          92             :     case('summary')
          93        1536 :        collect_stats = .true.
          94        1536 :        timestep_reset = .false.
          95             :     case('timestep')
          96           0 :        collect_stats = .true.
          97           0 :        timestep_reset = .true.
          98             :     case('off')
          99           0 :        collect_stats = .false.
         100           0 :        timestep_reset = .false.
         101             :     case default
         102        3072 :        call endrun(sub//" FATAL: '"//trim(print_qneg_warn)//"' is not a valid value for print_qneg_warn")
         103             :     end select
         104             :
         105        1536 :     if (masterproc) then
         106           2 :        if (collect_stats) then
         107           2 :           if (timestep_reset) then
         108           0 :              write(iulog, *) sub, ": QNEG statistics will be collected and printed for each timestep"
         109             :           else
         110           2 :              write(iulog, *) sub, ": QNEG statistics will be collected and printed at the end of the run"
         111             :           end if
         112             :        else
         113           0 :           write(iulog, *) sub, ": QNEG statistics will not be collected"
         114             :        end if
         115             :     end if
         116             :
         117        1536 :   end subroutine qneg_readnl
         118             :
         119        1536 :   subroutine qneg_init()
         120             :     use cam_history,    only: addfld, horiz_only
         121             :     use constituents,   only: cnst_longname
         122             :
         123             :     integer :: index
         124             :
         125       64512 :     do index = 1, pcnst
         126       62976 :        diag_names(index) = trim(cnst_name(index))//'_qneg3'
         127             :        call addfld(diag_names(index), (/ 'lev' /), 'I', 'kg/kg',              &
         128      125952 :             trim(cnst_longname(index))//' QNEG3 error (cell)')
         129       62976 :        diag_names(pcnst+index) = trim(cnst_name(index))//'_qneg3_col'
         130             :        call addfld(diag_names(pcnst+index), horiz_only, 'I', 'kg/kg',         &
         131       64512 :             trim(cnst_longname(index))//' QNEG3 error (column)')
         132             :     end do
         133        1536 :     diag_names((2*pcnst) + 1) = 'qflux_exceeded'
         134             :     call addfld(diag_names((2*pcnst) + 1), horiz_only, 'I', 'kg/m^2/s',     &
         135        1536 :          'qflux excess (QNEG4)')
         136             :
         137        1536 :   end subroutine qneg_init
         138             :
         139        1536 :   subroutine calc_cnst_out()
         140        1536 :     use cam_history, only: hist_fld_active, history_initialized
         141             :     integer :: index
         142             :
         143        1536 :     if (history_initialized()) then
         144             :        ! to protect against routines that call qneg3 too early
         145      129024 :        do index = 1, num_diag_fields
         146      129024 :           cnst_outfld(index) = hist_fld_active(trim(diag_names(index)))
         147             :        end do
         148        1536 :        cnst_out_calc = .true.
         149             :     end if
         150             :
         151        1536 :   end subroutine calc_cnst_out
         152             :
         153   743173128 :   integer function find_index3(nam) result(index)
         154             :     ! Find a valid or new index for 'nam' entries
         155             :     character(len=*),  intent(in) :: nam
         156             :
         157             :     integer                      :: i
         158             :
         159   743173128 :     index = -1
         160  9419391144 :     do i = 1, num3_bins
         161  9419391144 :        if (trim(nam) == trim(qneg3_warn_labels(i))) then
         162             :           ! We found this entry, return its index
         163             :           index = i
         164   743173128 :           exit
         165  8676248736 :        else if (len_trim(qneg3_warn_labels(i)) == 0) then
         166             :           ! We have run out of known entries, use a new one and reset its stats
         167       30720 :           qneg3_warn_labels(i) = nam
         168       30720 :           index = i
         169       30720 :           call reset_stats(qneg3_warn_num(:, index), qneg3_warn_worst(:,index))
         170       30720 :           exit
         171             :        end if
         172             :     end do
         173   743174664 :   end function find_index3
         174             :
         175     1489176 :   integer function find_index4(nam) result(index)
         176             :     ! Find a valid or new index for 'nam' entries
         177             :     character(len=*),  intent(in) :: nam
         178             :
         179             :     integer                      :: i
         180             :
         181     1489176 :     index = -1
         182     1489176 :     do i = 1, num4_bins
         183     1489176 :        if (trim(nam) == trim(qneg4_warn_labels(i))) then
         184             :           ! We found this entry, return its index
         185             :           index = i
         186     1489176 :           exit
         187        1536 :        else if (len_trim(qneg4_warn_labels(i)) == 0) then
         188             :           ! We have run out of known entries, use a new one and reset its stats
         189        1536 :           qneg4_warn_labels(i) = nam
         190        1536 :           index = i
         191             :           call reset_stats(qneg4_warn_num(index), qneg4_warn_worst(index))
         192        1536 :           exit
         193             :        end if
         194             :     end do
         195     1489176 :   end function find_index4
         196             :
         197   743173128 :   subroutine qneg3 (subnam, idx, ncol, ncold, lver, lconst_beg, &
         198   743173128 :        lconst_end, qmin, q)
         199             :     !-----------------------------------------------------------------------
         200             :     !
         201             :     ! Purpose:
         202             :     ! Check moisture and tracers for minimum value, reset any below
         203             :     ! minimum value to minimum value and return information to allow
         204             :     ! warning message to be printed. The global average is NOT preserved.
         205             :     !
         206             :     ! Method:
         207             :     ! <Describe the algorithm(s) used in the routine.>
         208             :     ! <Also include any applicable external references.>
         209             :     !
         210             :     ! Author: J. Rosinski
         211             :     !
         212             :     !-----------------------------------------------------------------------
         213             :     use cam_history, only: outfld
         214             :
         215             :     !------------------------------Arguments--------------------------------
         216             :     !
         217             :     ! Input arguments
         218             :     !
         219             :     character(len=*), intent(in) :: subnam ! name of calling routine
         220             :
         221             :     integer, intent(in) :: idx          ! chunk/latitude index
         222             :     integer, intent(in) :: ncol         ! number of atmospheric columns
         223             :     integer, intent(in) :: ncold        ! declared number of atmospheric columns
         224             :     integer, intent(in) :: lver         ! number of vertical levels in column
         225             :     integer, intent(in) :: lconst_beg   ! beginning constituent
         226             :     integer, intent(in) :: lconst_end   ! ending    constituent
         227             :
         228             :     real(r8), intent(in) :: qmin(lconst_beg:lconst_end)      ! Global minimum constituent concentration
         229             :
         230             :     !
         231             :     ! Input/Output arguments
         232             :     !
         233             :     real(r8), intent(inout) :: q(ncold,lver,lconst_beg:lconst_end) ! moisture/tracer field
         234             :     !
         235             :     !---------------------------Local workspace-----------------------------
         236             :     !
         237             :     integer  :: nvals            ! number of values found < qmin
         238             :     integer  :: i, k             ! longitude, level indices
         239             :     integer  :: index            ! For storing stats
         240             :     integer  :: m                ! constituent index
         241             :     integer  :: iw,kw            ! i,k indices of worst violator
         242             :
         243             :     logical  :: found            ! true => at least 1 minimum violator found
         244             :
         245  1486346256 :     real(r8) :: badvals(ncold, lver) ! Collector for outfld calls
         246  1486346256 :     real(r8) :: badcols(ncold)  ! Column sum for outfld
         247             :     real(r8) :: worst           ! biggest violator
         248             :     !
         249             :     !-----------------------------------------------------------------------
         250             :     !
         251             :
         252   743173128 :     call t_startf ('qneg3')
         253             :     ! The first time we call this, we need to determine whether to call outfld
         254   743173128 :     if (.not. cnst_out_calc) then
         255        1536 :        call calc_cnst_out()
         256             :     end if
         257             :
         258   743173128 :     if (collect_stats) then
         259   743173128 :        index = find_index3(trim(subnam))
         260             :     else
         261             :        index = -1
         262             :     end if
         263             :
         264  1605975696 :     do m = lconst_beg, lconst_end
         265             :        nvals = 0
         266 81103441392 :        found = .false.
         267 81103441392 :        worst = worst_reset
         268 >13649*10^8 :        badvals(:,:) = 0.0_r8
         269 81103441392 :        iw = -1
         270 81103441392 :        kw = -1
         271             :        !
         272             :        ! Test all field values for being less than minimum value. Set q = qmin
         273             :        ! for all such points. Trace offenders and identify worst one.
         274             :        !
         275 81103441392 :        do k = 1, lver
         276 >13406*10^8 :           do i = 1, ncol
         277 >13398*10^8 :              if (q(i,k,m) < qmin(m)) then
         278  7550970901 :                 found = .true.
         279             :                 nvals = nvals + 1
         280  7550970901 :                 badvals(i, k) = q(i, k, m)
         281  7550970901 :                 if (index > 0) then
         282  7550970901 :                    qneg3_warn_num(m, index) = qneg3_warn_num(m, index) + 1
         283             :                 end if
         284  7550970901 :                 if (q(i,k,m) < worst) then
         285   486955362 :                    worst = q(i,k,m)
         286   486955362 :                    iw = i
         287   486955362 :                    kw = k
         288   486955362 :                    if (index > 0) then
         289   486955362 :                       qneg3_warn_worst(m, index) = worst
         290             :                    end if
         291             :                 end if
         292  7550970901 :                 q(i,k,m) = qmin(m)
         293             :              end if
         294             :           end do
         295             :        end do
         296             :        ! Maybe output bad values
         297   862802568 :        if ((cnst_outfld(m)) .and. (worst < worst_reset)) then
         298           0 :           call outfld(trim(diag_names(m)), badvals, pcols, idx)
         299             :        end if
         300  1605975696 :        if ((cnst_outfld(pcnst+m)) .and. (worst < worst_reset)) then
         301           0 :           do i = 1, pcols
         302           0 :              badcols(i) = SUM(badvals(i,:))
         303             :           end do
         304           0 :           call outfld(trim(diag_names(pcnst+m)), badcols, pcols, idx)
         305             :        end if
         306             :     end do
         307   743173128 :     call t_stopf ('qneg3')
         308             :
         309   743173128 :   end subroutine qneg3
         310             :
         311     1489176 :   subroutine qneg4 (subnam, lchnk, ncol, ztodt,                            &
         312     1489176 :        qbot, srfrpdel, shflx, lhflx, qflx)
         313             :     !-----------------------------------------------------------------------
         314             :     !
         315             :     ! Purpose:
         316             :     ! Check if moisture flux into the ground is exceeding the total
         317             :     ! moisture content of the lowest model layer (creating negative moisture
         318             :     ! values).  If so, then subtract the excess from the moisture and
         319             :     ! latent heat fluxes and add it to the sensible heat flux.
         320             :     !
         321             :     ! Method:
         322             :     ! <Describe the algorithm(s) used in the routine.>
         323             :     ! <Also include any applicable external references.>
         324             :     !
         325             :     ! Author: J. Olson
         326             :     !
         327             :     !-----------------------------------------------------------------------
         328   743173128 :     use physconst,    only: gravit, latvap
         329             :     use constituents, only: qmin
         330             :     use cam_history,  only: outfld
         331             :
         332             :     !
         333             :     ! Input arguments
         334             :     !
         335             :     character(len=*), intent(in) :: subnam   ! name of calling routine
         336             :     !
         337             :     integer, intent(in) :: lchnk             ! chunk index
         338             :     integer, intent(in) :: ncol              ! number of atmospheric columns
         339             :     !
         340             :     real(r8), intent(in) :: ztodt            ! two times model timestep (2 delta-t)
         341             :     real(r8), intent(in) :: qbot(ncol,pcnst) ! moisture at lowest model level
         342             :     real(r8), intent(in) :: srfrpdel(ncol)   ! 1./(pint(K+1)-pint(K))
         343             :     !
         344             :     ! Input/Output arguments
         345             :     !
         346             :     real(r8), intent(inout) :: shflx(ncol)   ! Surface sensible heat flux (J/m2/s)
         347             :     real(r8), intent(inout) :: lhflx(ncol)   ! Surface latent   heat flux (J/m2/s)
         348             :     real(r8), intent(inout) :: qflx (ncol,pcnst) ! surface water flux (kg/m^2/s)
         349             :     !
         350             :     !---------------------------Local workspace-----------------------------
         351             :     !
         352             :     integer :: i                 ! column index
         353             :     integer :: iw                ! i index of worst violator
         354             :     integer :: index             ! caller bin index
         355             :     !
         356             :     real(r8):: worst             ! biggest violator
         357     2978352 :     real(r8):: excess(ncol)     ! Excess downward sfc latent heat flux
         358             :     !
         359             :     !-----------------------------------------------------------------------
         360             :
         361     1489176 :     call t_startf ('qneg4')
         362             :     ! The first time we call this, we need to determine whether to call outfld
         363     1489176 :     if (.not. cnst_out_calc) then
         364           0 :        call calc_cnst_out()
         365             :     end if
         366             :
         367     1489176 :     if (collect_stats) then
         368     1489176 :        index = find_index4(trim(subnam))
         369             :     else
         370             :        index = -1
         371             :     end if
         372             :
         373             :     !
         374             :     ! Compute excess downward (negative) q flux compared to a theoretical
         375             :     ! maximum downward q flux.  The theoretical max is based upon the
         376             :     ! given moisture content of lowest level of the model atmosphere.
         377             :     !
         378     1489176 :     worst = worst_reset
         379    24865776 :     do i = 1, ncol
         380    23376600 :        excess(i) = qflx(i,1) - (qmin(1) - qbot(i,1))/(ztodt*gravit*srfrpdel(i))
         381             :        !
         382             :        ! If there is an excess downward (negative) q flux, then subtract
         383             :        ! excess from "qflx" and "lhflx" and add to "shflx".
         384             :        !
         385    24865776 :        if (excess(i) < 0._r8) then
         386           2 :           if (excess(i) < worst) then
         387           2 :              iw = i
         388           2 :              worst = excess(i)
         389             :           end if
         390           2 :           qflx (i,1) = qflx (i,1) - excess(i)
         391           2 :           lhflx(i) = lhflx(i) - excess(i)*latvap
         392           2 :           shflx(i) = shflx(i) + excess(i)*latvap
         393           2 :           if (index > 0) then
         394           2 :              qneg4_warn_num(index) = qneg4_warn_num(index) + 1
         395             :           end if
         396             :        end if
         397             :     end do
         398             :     ! Maybe output bad values
         399     1489176 :     if ((cnst_outfld((2*pcnst)+1)) .and. (worst < worst_reset)) then
         400           0 :        do i = 1, ncol
         401           0 :           if (excess(i) > 0.0_r8) then
         402           0 :              excess(i) = 0.0_r8
         403             :           end if
         404             :        end do
         405           0 :        call outfld(trim(diag_names((2*pcnst)+1)), excess(1:ncol), ncol, lchnk)
         406             :     end if
         407     1489176 :     call t_stopf ('qneg4')
         408             :
         409     1489176 :   end subroutine qneg4
         410             :
         411      369408 :   subroutine qneg_print_summary(end_of_run)
         412     1489176 :     use spmd_utils, only: mpicom, masterprocid, masterproc
         413             :     use spmd_utils, only: MPI_MIN, MPI_SUM, MPI_INTEGER, MPI_REAL8
         414             :
         415             :     logical, intent(in) :: end_of_run
         416             :
         417             :     integer             :: global_warn_num(pcnst)
         418             :     real(r8)            :: global_warn_worst(pcnst)
         419             :     integer             :: index, m
         420             :     integer             :: ierr
         421             :
         422      369408 :     if (collect_stats) then
         423      369408 :        if (timestep_reset .or. end_of_run) then
         424       38400 :           do index = 1, num3_bins
         425             :              ! QNEG3
         426       36864 :              call reset_stats(global_warn_num(:), global_warn_worst(:))
         427             :              call MPI_REDUCE(qneg3_warn_num(:, index), global_warn_num(:),    &
         428       36864 :                   pcnst, MPI_INTEGER, MPI_SUM, masterprocid, mpicom, ierr)
         429             :              call MPI_REDUCE(qneg3_warn_worst(:, index), global_warn_worst(:),&
         430       36864 :                   pcnst, MPI_REAL8, MPI_MIN, masterprocid, mpicom, ierr)
         431       36864 :              if (masterproc) then
         432        2016 :                 do m = 1, pcnst
         433        1968 :                    if ( (global_warn_num(m) > 0) .and.                        &
         434             :                         (abs(global_warn_worst(m)) > tol)) then
         435          67 :                       write(iulog, 9100) trim(qneg3_warn_labels(index)),      &
         436          67 :                            trim(cnst_name(m)), global_warn_num(m),            &
         437         134 :                            global_warn_worst(m)
         438             :                    end if
         439        2016 :                    call shr_sys_flush(iulog)
         440             :                 end do
         441             :              end if
         442       38400 :              call reset_stats(qneg3_warn_num(:,index), qneg3_warn_worst(:,index))
         443             :           end do
         444        7680 :           do index = 1, num4_bins
         445             :              ! QNEG4
         446        6144 :              call reset_stats(qneg4_warn_num(index), qneg4_warn_worst(index))
         447        6144 :              call reset_stats(global_warn_num(1), global_warn_worst(1))
         448             :              call MPI_REDUCE(qneg4_warn_num(index), global_warn_num(1),       &
         449        6144 :                   1, MPI_INTEGER, MPI_SUM, masterprocid, mpicom, ierr)
         450             :              call MPI_REDUCE(qneg4_warn_worst(index), global_warn_worst(1),   &
         451        6144 :                   1, MPI_REAL8, MPI_MIN, masterprocid, mpicom, ierr)
         452        6144 :              if (masterproc) then
         453           8 :                 if ( (global_warn_num(1) > 0) .and.                           &
         454             :                      (abs(global_warn_worst(1)) > tol)) then
         455           0 :                    write(iulog, 9101) trim(qneg4_warn_labels(index)),          &
         456           0 :                         global_warn_num(1), global_warn_worst(1)
         457             :                 end if
         458           8 :                 call shr_sys_flush(iulog)
         459             :              end if
         460       13824 :              call reset_stats(qneg4_warn_num(index), qneg4_warn_worst(index))
         461             :           end do
         462             :        end if
         463             :     end if
         464             :
         465      369408 :     return
         466             : 9100 format(' QNEG3 from ', a, ':', a, &
         467             :          ' Min. mixing ratio violated at ', i9, ' points. Worst = ', e10.1)
         468             : 9101 format(' QNEG4 from ',a,': moisture flux exceeded at', &
         469             :           i9, ' points. Worst = ', e10.1)
         470             :   end subroutine qneg_print_summary
         471             :
         472      104448 :   subroutine reset_stats_array(num_array, worst_array)
         473             :     ! Private routine to reset statistics
         474             :     integer,  intent(inout) :: num_array(:)
         475             :     real(r8), intent(inout) :: worst_array(:)
         476             :
         477     4386816 :     num_array(:)    = 0
         478     4386816 :     worst_array(:)  = worst_reset
         479      104448 :   end subroutine reset_stats_array
         480             :
         481       19968 :   subroutine reset_stats_scalar(num, worst)
         482             :     ! Private routine to reset statistics
         483             :     integer,  intent(inout) :: num
         484             :     real(r8), intent(inout) :: worst
         485             :
         486       19968 :     num    = 0
         487       19968 :     worst  = worst_reset
         488       19968 :   end subroutine reset_stats_scalar
         489             :
         490             : end module qneg_module

Generated by: LCOV version 1.14