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 calculates new gas concentrations.
       6             : !!
       7             : !! @author Andy Ackerman, Bill McKie, Chuck Bardeen
       8             : !! @version Dec-1995, Sep-1997, Nov-2009
       9   213271518 : subroutine gsolve(carma, cstate, iz, previous_ice, previous_liquid, scale_threshold, rc)
      10             : 
      11             :   ! types
      12             :   use carma_precision_mod
      13             :   use carma_enums_mod
      14             :   use carma_constants_mod
      15             :   use carma_types_mod
      16             :   use carmastate_mod
      17             :   use carma_mod
      18             : 
      19             :   implicit none
      20             : 
      21             :   type(carma_type), intent(in)         :: carma   !! the carma object
      22             :   type(carmastate_type), intent(inout) :: cstate  !! the carma state object
      23             :   integer, intent(in)                  :: iz      !! z index
      24             :   real(kind=f), intent(in)             :: previous_ice(NGAS)      !! total ice at the start of substep
      25             :   real(kind=f), intent(in)             :: previous_liquid(NGAS)   !! total liquid at the start of substep
      26             :   real(kind=f)                         :: scale_threshold !! Scaling factor for convergence thresholds
      27             :   integer, intent(inout)               :: rc      !! return code, negative indicates failure
      28             : 
      29             :   ! Local Variables
      30             :   integer                              :: igas    !! gas index
      31             :   real(kind=f)                         :: gc_cgs
      32             :   real(kind=f)                         :: rvap
      33   426543036 :   real(kind=f)                         :: total_ice(NGAS)      ! total ice
      34   426543036 :   real(kind=f)                         :: total_liquid(NGAS)   ! total liquid
      35             :   real(kind=f)                         :: threshold            ! convergence threshold
      36             : 
      37             : 
      38             :   1 format(/,'gsolve::ERROR - negative gas concentration for ',a,' : iz=',i4,',xc=', &
      39             :               f7.2,',yc=',f7.2,',gc=',e10.3,',gasprod=',e10.3,',supsati=',e10.3, &
      40             :               ',supsatl=',e10.3,',t=',f6.2)
      41             :   2 format('gsolve::ERROR - conditions at beginning of the step : gc=',e10.3,',supsati=',e17.10, &
      42             :               ',supsatl=',e17.10,',t=',f6.2,',d_gc=',e10.3,',d_t=',f6.2)
      43             :   3 format(/,'microfast::WARNING - gas concentration change exceeds threshold: ',a,' : iz=',i4,',xc=', &
      44             :               f7.2,',yc=',f7.2, ', (gc-gcl)/gcl=', e10.3)
      45             : 
      46             : 
      47             :   ! Determine the total amount of condensate for each gas.
      48   213271518 :   call totalcondensate(carma, cstate, iz, total_ice, total_liquid, rc)
      49             : 
      50   639814554 :   do igas = 1,NGAS
      51             : 
      52             :     ! We do not seem to be conserving mass and energy, so rather than relying upon gasprod
      53             :     ! and rlheat, recalculate the total change in condensate to determine the change
      54             :     ! in gas and energy.
      55             :     !
      56             :     ! This is because in the old scheme, the particles were solved for implicitly, but the
      57             :     ! gas and latent heat were solved for explicitly using the same rates.
      58   426543036 :     gasprod(igas) = ((previous_ice(igas) - total_ice(igas)) + (previous_liquid(igas) - total_liquid(igas))) / dtime
      59           0 :     rlprod        = rlprod - ((previous_ice(igas) - total_ice(igas)) * (rlhe(iz,igas) + rlhm(iz,igas)) + &
      60   426543036 :                        (previous_liquid(igas) - total_liquid(igas)) * (rlhe(iz,igas))) / (CP * rhoa(iz) * dtime)
      61             : 
      62             :     ! Don't let the gas concentration go negative.
      63   426543036 :     gc(iz,igas) = gc(iz,igas) + dtime * gasprod(igas)
      64             : 
      65   426543036 :     if (gc(iz,igas) < 0.0_f) then
      66    20062061 :       if (do_substep) then
      67    20062061 :         if (nretries == maxretries) then
      68           0 :           if (do_print) write(LUNOPRT,1) trim(gasname(igas)), iz, xc, yc, gc(iz,igas), gasprod(igas), &
      69           0 :             supsati(iz,igas), supsatl(iz,igas), t(iz)
      70           0 :           if (do_print) write(LUNOPRT,2) gcl(iz,igas), supsatiold(iz,igas), supsatlold(iz,igas), told(iz), d_gc(iz,igas), d_t(iz)
      71             :         end if
      72             :       else
      73           0 :         if (do_print) write(LUNOPRT,1) trim(gasname(igas)), iz, xc, yc, gc(iz,igas), gasprod(igas), &
      74           0 :           supsati(iz, igas), supsatl(iz,igas), t(iz)
      75             :       end if
      76             : 
      77    20062061 :       rc = RC_WARNING_RETRY
      78             :     end if
      79             : 
      80             :     ! If gas changes by too much, then retry the calculation.
      81   426543036 :     threshold = dgc_threshold(igas) / scale_threshold
      82             : 
      83   639814554 :     if (threshold /= 0._f) then
      84           0 :       if ((dtime * gasprod(igas) / gc(iz,igas)) > threshold) then
      85           0 :         if (do_substep) then
      86           0 :           if (nretries == maxretries) then
      87           0 :             if (do_print) write(LUNOPRT,3) trim(gasname(igas)), iz, xc, yc, dtime * gasprod(igas) / gc(iz,igas)
      88           0 :             if (do_print) write(LUNOPRT,2) gcl(iz,igas), supsatiold(iz,igas), supsatlold(iz,igas), told(iz), d_gc(iz,igas), d_t(iz)
      89             :           end if
      90             :         else
      91           0 :           if (do_print) write(LUNOPRT,3) trim(gasname(igas)), iz, xc, yc, dtime * gasprod(igas) / gc(iz,igas)
      92             :         end if
      93             : 
      94           0 :         rc = RC_WARNING_RETRY
      95             :       end if
      96             :     end if
      97             :   end do
      98             : 
      99             :   ! Return to caller with new gas concentrations.
     100   213271518 :   return
     101   213271518 : end