LCOV - code coverage report
Current view: top level - physics/cam - ssatcontrail.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 19 87 21.8 %
Date: 2024-12-17 22:39:59 Functions: 1 1 100.0 %

          Line data    Source code
       1             : module ssatcontrail
       2             : ! contrail parameterization
       3             : ! see Chen et al., 2012: Global contrail coverage simulated
       4             : !     by CAM5 with the inventory of 2006 global aircraft emissions, JAMES
       5             : !     https://doi.org/10.1029/2011MS000105
       6             :     use shr_kind_mod,   only: r8 => shr_kind_r8
       7             :     use ppgrid,         only: pcols, pver
       8             :     use physics_types,  only: physics_state, physics_ptend, physics_ptend_init
       9             :     use physconst,      only: cpair,mwdry,mwh2o, gravit, zvir, rair, pi, rearth, tmelt
      10             :     use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc,  pbuf_old_tim_idx
      11             :     use constituents,   only: cnst_get_ind, pcnst
      12             :     use phys_grid,      only: get_wght_all_p
      13             :     use wv_saturation,  only: qsat_water, qsat_ice
      14             :     use aircraft_emit,  only: get_aircraft
      15             : 
      16             :     implicit none
      17             :     private
      18             :     save
      19             : 
      20             :     public ssatcontrail_d0
      21             : 
      22             :     ! Private data
      23             :     real(r8), parameter :: rhoi = 500.0_r8             ! density of ice (500 kg/m3)
      24             :     real(r8), parameter :: radius = 3.75e-6_r8         ! diameter of ice particle = 7.5 microns
      25             : 
      26             :  
      27             : contains
      28             : 
      29     2978352 :     subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc)
      30             :     implicit none
      31             : 
      32             :     type(physics_state), intent(in) :: state1
      33             :     type(physics_ptend), intent(inout) :: ptend_loc
      34             :     type(physics_buffer_desc), pointer :: pbuf(:)
      35             :     real(r8),            intent(in)    :: dtime      ! time step
      36             : !------------------------Local storage------------------------------------------------------
      37             :     real(r8) :: Ma, Mh2o, epsi, Q, eta, p, G, T_contr, eslTc, eslT, RH_contr, qslTc, qslT
      38             :     real(r8) :: w, esiT, qsiT, ws, RH, ei
      39             :     integer  :: i,k
      40             :     integer  :: lchnk,ncol
      41             :     real(r8) :: contrail(pcols,pver), pcontrail(pcols,pver)
      42     1489176 :     real(r8), pointer, dimension(:,:) ::  cld    ! cloud fraction
      43     1489176 :     real(r8), pointer, dimension(:,:) :: ac_H2O
      44     1489176 :     real(r8), pointer, dimension(:,:) :: ac_SLANT_DIST
      45             :     integer  :: itim, ifld
      46             :     integer :: ixcldice, ixcldliq                  ! indices for CLDICE and CLDLIQ
      47             :     integer :: ixnumice, ixnumliq
      48             :     real(r8):: zi, zm, rog
      49             :     logical :: has_aircraft_H2O
      50             :     logical :: has_aircraft_distance
      51             :     real(r8) :: hkl, hkk, tv
      52             :     real(r8) :: particle_mass
      53             :     real(r8) :: ICIWC0(pcols,pver), ICIWC, rho
      54             :     real(r8) :: qs
      55             :     real(r8) :: wght(pcols)
      56             :     real(r8) :: dz, ratio(pcols,pver)
      57             :     real(r8) :: dcld(pcols,pver)
      58             :     real(r8) :: ac_Q, ac_Q1, ac_Q2
      59             :     real(r8) :: RHcts(pcols,pver)
      60             :     logical :: lq(pcnst)
      61             : 
      62             :     integer :: yr, mon, day, ncsec
      63             :     real(r8) :: curr_factor
      64             :     integer :: aircraft_cnt
      65             :     character(len=16) :: spc_name_list(30)
      66             :  
      67     1489176 :     has_aircraft_H2O = .false.
      68     1489176 :     has_aircraft_distance = .false.
      69             : 
      70             :     ! Update constituents, all schemes use time split q: no tendency kept
      71     1489176 :     call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
      72     1489176 :     call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
      73             :     ! Check for number concentration of cloud liquid and cloud ice (if not present)
      74             :     ! the indices will be set to -1)
      75     1489176 :     call cnst_get_ind('NUMICE', ixnumice, abort=.false.)
      76     1489176 :     call cnst_get_ind('NUMLIQ', ixnumliq, abort=.false.)
      77             : 
      78     1489176 :     call get_aircraft(aircraft_cnt, spc_name_list)
      79             : !-----------------------------------------------------------------------------------------
      80             : ! check if ac_H2O in namelist
      81             : ! if not, then bypass this subroutine
      82             : !-----------------------------------------------------------------------------------------
      83     1489176 :     if(aircraft_cnt>0) then
      84           0 :      do i=1,aircraft_cnt
      85           0 :       if(trim(spc_name_list(i)) == 'ac_H2O') then
      86           0 :         has_aircraft_H2O = .true.
      87             :       endif
      88           0 :       if(trim(spc_name_list(i)) == 'ac_SLANT_DIST' .or. trim(spc_name_list(i)) == 'ac_TRACK_DIST') then
      89           0 :         has_aircraft_distance = .true.
      90             :       endif
      91             :      enddo
      92             :     endif
      93             : !------------------------------------------------------------------------------------------
      94     1489176 :     lq(:) = .FALSE.
      95     1489176 :     lq(1) = .TRUE.
      96     1489176 :     lq(ixcldice) = .TRUE.
      97     1489176 :     lq(ixnumice) = .TRUE.
      98             : 
      99     1489176 :     call physics_ptend_init(ptend_loc, state1%psetcols,'ssatcontrail',ls=.true.,lq = lq)
     100             : !-----------------------------------------------------------------------------------------
     101     1489176 :     if(.not. has_aircraft_H2O)  return
     102           0 :     if(.not. has_aircraft_distance) return
     103             : 
     104           0 :     particle_mass = 4._r8/3._r8*pi*rhoi*radius**3   ! mass of ice particle
     105             :  
     106           0 :     rog = rair/gravit                               ! Rd/Cp
     107             : 
     108             : 
     109           0 :     contrail(:,:) = 0.0_r8
     110           0 :     pcontrail(:,:) = 0.0_r8
     111           0 :     ICIWC0(:,:) = 0.0_r8
     112           0 :     RHcts(:,:) = 0.0_r8
     113             : 
     114           0 :     lchnk = state1%lchnk
     115           0 :     ncol  = state1%ncol
     116             :  
     117           0 :     call get_wght_all_p(lchnk, ncol, wght)
     118             : 
     119           0 :     itim = pbuf_old_tim_idx()
     120           0 :     ifld = pbuf_get_index('CLD')
     121           0 :     call pbuf_get_field(pbuf, ifld, cld, (/1,1,itim/),(/pcols,pver,1/))
     122             : 
     123           0 :     ifld = pbuf_get_index('ac_H2O')
     124           0 :     call pbuf_get_field(pbuf,ifld,ac_H2O)
     125           0 :     ifld = pbuf_get_index('ac_SLANT_DIST')
     126           0 :     call pbuf_get_field(pbuf,ifld,ac_SLANT_DIST)
     127             :  
     128             : ! adjust h2o to volume mixing ratio (mass adjustment and conversion from g/kg to kg/kg)
     129           0 :     Ma = mwdry
     130           0 :     Mh2o = mwh2o
     131             : 
     132             : ! contrail paramter (G = CF*p/epi)
     133             : ! and Schumann 1996 DOI: 10.1127/metz/5/1996/4, reprinted by Ponater 2002, JGR (eq 6-8) DOI: 10.1029/2011MS000105
     134             : 
     135           0 :     epsi = Mh2o/Ma
     136           0 :     ei = 1.21_r8      ! water vapor emision index (g) h2o per kg fuel (Schumann 96)
     137           0 :     Q = 43.e6_r8      ! specific combustion heat Schummann 1996, Q = 43 MJ/kg
     138           0 :     eta = 0.3_r8      ! propulsion effieciency (Ponater 2002)
     139             :  
     140           0 :     ratio = 0._r8
     141           0 :     dcld = 0._r8
     142             : 
     143           0 :     do i=1,ncol
     144           0 :       do k=1,pver
     145           0 :         p = (state1%pint(i,k)+state1%pint(i,k+1))/2.0_r8
     146             :  
     147           0 :         G = (ei*cpair*p)/(epsi*Q*(1.0_r8-eta))   ! eq 7, Ponater JGR 2002
     148             : 
     149           0 :         if( G > 0.053_r8 ) then
     150           0 :             T_contr = -46.46_r8+9.43_r8*log(G-0.053_r8)+0.72_r8*log(G-0.053_r8)*log(G-0.053_r8) ! eq 6, Ponater JGR 2002
     151           0 :             T_contr = T_contr + tmelt  ! convert to Kelvin
     152             :  
     153             :             ! compute saturation pressure
     154           0 :             call qsat_water(T_contr, p, eslTc, qslTc)
     155           0 :             call qsat_water(state1%t(i,k), p, eslT, qslT)
     156             : 
     157           0 :             RH_contr = (G*(state1%t(i,k)-T_contr)+eslTc)/eslT
     158             :             ! RH_contr ranges between 0 and 1
     159           0 :             if(RH_contr>1.0_r8) RH_contr = 1.0_r8
     160           0 :             if(RH_contr<0.0_r8) RH_contr = 0.0_r8
     161             :  
     162           0 :             w = state1%q(i,k,1)/(1.0_r8-state1%q(i,k,1))  ! mixing ratio from specific humidity
     163           0 :             call qsat_ice(state1%t(i,k), p, esiT, qsiT)
     164           0 :             ws = epsi*esiT/(p-esiT)  ! saturation mixing ration with respect to ice
     165           0 :             qs = ws/(1.0_r8+ws)
     166             : 
     167           0 :             RH = w/ws  ! relative humidity with respect to ice
     168           0 :             if( RH>=1.0_r8 ) RHcts(i,k) = 1.0_r8
     169             : 
     170             : ! Schumann, U. “Contrail Cirrus.” In Cirrus, edited by D. K. Lynch and others, 231–55. Oxford University Press, 2002
     171             : !                 IWC(g/m3) = exp(6.97+0.103*T(C))*1e-3
     172             : !                 IWC(kg/m3) = exp(6.97+0.103*T(C))*1e-6
     173             : 
     174           0 :             ICIWC0(i,k) = exp(6.97_r8+0.103_r8*(state1%t(i,k)-tmelt)) ! in mg/m3
     175           0 :             rho = p/(rair*state1%t(i,k))
     176           0 :             ICIWC = ICIWC0(i,k)/rho*1.0e-6_r8
     177             : 
     178             : 
     179             : ! persistent contrail condition
     180           0 :             if( (state1%t(i,k)<T_contr).and.(RH>RH_contr).and.(RH>1.0_r8).and.(ac_H2O(i,k)>0.0_r8) ) then
     181             : 
     182             : ! if persistent contrail, H2O emitted from aircraft turns into cloud ice
     183           0 :                  dz = state1%zi(i,k)-state1%zi(i,k+1)
     184           0 :                  ratio(i,k) = (ac_SLANT_DIST(i,k)*dtime*1.e4_r8)/(dz*rearth*rearth*wght(i))
     185             : 
     186           0 :                  ac_Q = min(ac_H2O(i,k)*dtime + (state1%q(i,k,1)-qs)*ratio(i,k),ratio(i,k)*ICIWC)
     187           0 :                  ptend_loc%q(i,k,ixcldice) = ac_Q/dtime
     188             : 
     189             : ! take out water vapor from q
     190           0 :                  ptend_loc%q(i,k,1) = -(ac_Q-ac_H2O(i,k)*dtime)/dtime
     191             : 
     192             : ! modify cloud fraction
     193             : ! by a prescribed ICIWC, we may deduce the new cloud fraction
     194             :  
     195           0 :                  cld(i,k) = min(1._r8, cld(i,k)+ac_Q/ICIWC)
     196             :  
     197             : ! modify cloud ice number concentration,
     198             : ! by assuming the particle size, the number of ice particles may be obtained
     199           0 :                  ptend_loc%q(i,k,ixnumice) = ac_Q/particle_mass/dtime
     200             : 
     201             :            else
     202             : ! if not persistent contrail, just add ac_H2O to state1%q(1) (vapor phase)
     203             : 
     204           0 :                  ptend_loc%q(i,k,1) = ac_H2O(i,k)
     205             :  
     206             :            endif
     207             :        
     208             :        else
     209           0 :            ptend_loc%q(i,k,1) = ac_H2O(i,k)
     210             :        end if
     211             : 
     212             :       enddo
     213             : 
     214             : ! modify dry static energy if water field is added to any grid cell
     215             : ! this bypasses geopotential_t which assumes dry static energy conservation
     216             : ! water vapor added to the system is assumed to increase dry static energy
     217             : ! conservation of dry static energy by geopotential_t will lower temperature to compensate
     218             :  
     219             :       zi = 0.0_r8
     220           0 :       do k=pver,1,-1
     221           0 :          hkl = state1%lnpint(i,k+1)-state1%lnpint(i,k)
     222           0 :          hkk = 1._r8 - state1%pint(i,k) * hkl * state1%rpdel(i,k)
     223             : 
     224           0 :          tv = state1%t(i,k) * (1._r8 + zvir*(state1%q(i,k,1)+ptend_loc%q(i,k,1)*dtime))
     225             : 
     226           0 :          zm = zi + rog * tv * hkk
     227           0 :          zi = zi + rog * tv * hkl
     228             : 
     229           0 :          ptend_loc%s(i,k) = (cpair*state1%t(i,k)+gravit*zm + state1%phis(i) - state1%s(i,k) )/dtime
     230             :       enddo
     231             : 
     232             :      enddo
     233             : 
     234     1489176 :     end subroutine ssatcontrail_d0
     235             : 
     236             : end module ssatcontrail

Generated by: LCOV version 1.14