Line data    Source code

       1             : ! Include shortname defintions, so that the F77 code does not have to be modified to
       2             : ! reference the CARMA structure.
       3             : #include "carma_globaer.h"
       4             : 
       5             : !! This routine checks the column for errors where the sum of the core masses
       6             : !! is larger than the concentration element. This implies a negative mass
       7             : !! of the concentration element, which is not physical.
       8             : !!
       9             : !! This routine attempts to conserve mass by using mass from the positive
      10             : !! values to offset the negative values.
      11             : !!
      12             : !! @author Charles Bardeen
      13             : !! version Feb-2023
      14     7354368 : subroutine fixcorecol(carma, cstate, rc)
      15             : 
      16             :   ! types
      17             :   use carma_precision_mod
      18             :   use carma_enums_mod
      19             :   use carma_constants_mod
      20             :   use carma_types_mod
      21             :   use carmastate_mod
      22             :   use carma_mod
      23             : 
      24             :   implicit none
      25             : 
      26             :   type(carma_type), intent(in)         :: carma   !! the carma object
      27             :   type(carmastate_type), intent(inout) :: cstate  !! the carma state object
      28             :   integer, intent(inout)               :: rc      !! return code, negative indicates failure
      29             : 
      30             :   ! Declare local variables
      31             :   integer        :: ibin
      32             :   integer        :: igroup
      33             :   integer        :: iepart
      34             :   integer        :: iz
      35             :   integer        :: icore
      36             :   integer        :: i
      37    14708736 :   real(kind=f)   :: total_core(NZ)
      38    14708736 :   real(kind=f)   :: concgas_md(NZ)
      39             :   real(kind=f)   :: total_mass
      40             :   real(kind=f)   :: missing_mass
      41             :   real(kind=f)   :: factor
      42             : 
      43     7354368 :   rc = RC_OK
      44             : 
      45             :   ! Advection can cause error in tracer/tracer relationship that can cause
      46             :   ! negative values for the concentration element. Use a mass conserving
      47             :   ! fixer to make sure there are no negative values.
      48             : 
      49             :   ! Find out the total amount of the concentration mass and the amount needed to
      50             :   ! fill in any negative values.
      51             :   !
      52    22063104 :   do igroup =  1, NGROUP
      53             : 
      54             :     ! We only need to do this where there are core masses.
      55    22063104 :     if (ncore(igroup) .gt. 0) then
      56             : 
      57     7354368 :       iepart = ienconc(igroup)
      58             : 
      59   154441728 :       do ibin = 1, NBIN
      60             : 
      61             :         ! Find the total of the core masses
      62  4853882880 :         total_core(:) = 0._f
      63             : 
      64  1470873600 :         do i = 1, ncore(igroup)
      65  1323786240 :           icore = icorelem(i,igroup)
      66 43832033280 :           total_core(:) = total_core(:) + pc(:, ibin, icore)
      67             :         end do
      68             : 
      69             :         ! Are there any places where the sum of the core masses is
      70             :         ! greater than the concentration element?
      71  4514076885 :         if (any(total_core(:) > pc(:, ibin, iepart)*rmass(ibin, igroup))) then
      72             : 
      73             :           ! Determine the total mass and the missing mass.
      74    13593248 :           total_mass = 0._f
      75    13593248 :           missing_mass = 0._f
      76             : 
      77   448577184 :           do iz = 1, NZ
      78   434983936 :             concgas_md(iz) = pc(iz, ibin, iepart)*rmass(ibin, igroup) - total_core(iz)
      79   448577184 :             if (concgas_md(iz) > 0._f) then
      80   389533736 :               total_mass = total_mass + concgas_md(iz) * dz(iz)
      81    45450200 :             else if (concgas_md(iz) < 0._f) then
      82    38487097 :               missing_mass = missing_mass - concgas_md(iz) * dz(iz)
      83             :             end if
      84             :           end do
      85             : 
      86             :           ! Is there enough mass to fill in the holes?
      87    13593248 :           if (total_mass >= missing_mass) then
      88             : 
      89             :             ! Do we need to add any fudge factor to this for
      90             :             ! roundoff, ...?
      91    13483936 :             factor = (total_mass - missing_mass) / total_mass
      92             : 
      93   444969888 :             do iz = 1, NZ
      94             : 
      95             :               ! Scale the positive concentration mass by enough to fill
      96             :               ! the gaps.
      97   444969888 :               if (concgas_md(iz) > 0._f) then
      98             : 
      99   387709393 :                 pc(iz, ibin, iepart) = (total_core(iz) + factor * concgas_md(iz)) / rmass(ibin, igroup)
     100             : 
     101             :               ! Fill the negative values with zero.
     102    43776559 :               else if (concgas_md(iz) < 0._f) then
     103             : 
     104    37717572 :                 pc(iz, ibin, iepart) = total_core(iz) / rmass(ibin, igroup)
     105             :               end if
     106             :             end do
     107             :           else
     108             : 
     109             :             ! Since there isn't enough mass in the column to fix the negative values, just
     110             :             ! zero out the column.
     111     3607296 :             pc(:, ibin, iepart) = total_core(:) / rmass(ibin, igroup)
     112             :           end if
     113             :         end if
     114             :       end do
     115             :     end if
     116             :   end do
     117             : 
     118     7354368 :   return
     119     7354368 : end subroutine fixcorecol