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 17:57:11 |  | 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