Changeset 5453


Ignore:
Timestamp:
Dec 23, 2024, 8:19:39 PM (24 hours ago)
Author:
aborella
Message:

Added water emissions and IO routines for contrails.
Work to be done: contrails initial cross section and radiative transfer

Location:
LMDZ6/branches/contrails/libf/phylmd
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90

    r5452 r5453  
    77CONTAINS
    88
     9SUBROUTINE aviation_water_emissions( &
     10      klon, klev, dtime, paprs, pplay, temp, qtot, cell_area, &
     11      flight_h2o, d_q_avi &
     12      )
     13
     14USE lmdz_lscp_ini, ONLY: RD, RG
     15
     16IMPLICIT NONE
     17
     18INTEGER,                      INTENT(IN)  :: klon, klev ! number of horizontal grid points and vertical levels
     19REAL,                         INTENT(IN)  :: dtime      ! time step [s]
     20REAL, DIMENSION(klon,klev+1), INTENT(IN)  :: paprs      ! inter-layer pressure [Pa]
     21REAL, DIMENSION(klon,klev),   INTENT(IN)  :: pplay      ! mid-layer pressure [Pa]
     22REAL, DIMENSION(klon,klev),   INTENT(IN)  :: temp       ! temperature (K)
     23REAL, DIMENSION(klon,klev),   INTENT(IN)  :: qtot       ! total specific humidity (in vapor phase) [kg/kg]
     24REAL, DIMENSION(klon),        INTENT(IN)  :: cell_area  ! area of each cell [m2]
     25REAL, DIMENSION(klon,klev),   INTENT(IN)  :: flight_h2o ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
     26REAL, DIMENSION(klon,klev),   INTENT(OUT) :: d_q_avi    ! water vapor tendency from aviation [kg/kg]
     27! Local
     28INTEGER :: i, k
     29REAL :: rho, rhodz, dz, M_cell
     30
     31DO i=1, klon
     32  DO k=1, klev
     33    !--Dry density [kg/m3]
     34    rho = pplay(i,k) / temp(i,k) / RD
     35    !--Dry air mass [kg/m2]
     36    rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
     37    !--Cell thickness [m]
     38    dz = rhodz / rho
     39    !--Cell dry air mass [kg]
     40    M_cell = rhodz * cell_area(i)
     41
     42    !--q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
     43    ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 
     44    !      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
     45    !--The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
     46    !--flight_h2O is in kg H2O / s / mesh
     47   
     48    !d_q_avi(i,k) = ( M_cell * qtot(i,k) + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
     49    !             / ( M_cell             + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
     50    !             - qtot(i,k)
     51    !--Same formula, more computationally effective but less readable
     52    d_q_avi(i,k) = flight_h2o(i,k) * ( 1. - qtot(i,k) ) &
     53                 / ( M_cell / dtime / ( 1. - qtot(i,k) ) + flight_h2o(i,k) )
     54  ENDDO
     55ENDDO
     56
     57END SUBROUTINE aviation_water_emissions
     58
    959
    1060!**********************************************************************************
    1161SUBROUTINE contrails_formation_evolution( &
    1262      dtime, pplay, temp, qsat, qsatl, gamma_cond, rcont_seri, flight_dist, &
    13       cldfra, qvc, V_cell, M_cell, pdf_loc, pdf_scale, pdf_alpha, &
     63      cldfra, qvc, dz, V_cell, pdf_loc, pdf_scale, pdf_alpha, &
    1464      Tcritcont, qcritcont, potcontfraP, potcontfraNP, contfra, &
    1565      dcf_avi, dqvc_avi, dqi_avi &
    … …  
    3585REAL, INTENT(IN)  :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
    3686REAL, INTENT(IN)  :: rcont_seri   ! ratio of contrails fraction to total cloud fraction [-]
    37 REAL, INTENT(IN)  :: flight_dist  !
     87REAL, INTENT(IN)  :: flight_dist  ! aviation distance flown within the mesh [m/s/mesh]
    3888REAL, INTENT(IN)  :: cldfra       ! cloud fraction [-]
    3989REAL, INTENT(IN)  :: qvc          ! gridbox-mean vapor in the cloud [kg/kg]
     90REAL, INTENT(IN)  :: dz           ! cell width [m]
    4091REAL, INTENT(IN)  :: V_cell       ! cell volume [m3]
    41 REAL, INTENT(IN)  :: M_cell       ! cell mass [kg]
    4292REAL, INTENT(IN)  :: pdf_loc      ! location parameter of the clear sky PDF [%]
    4393REAL, INTENT(IN)  :: pdf_scale    ! scale parameter of the clear sky PDF [%]
    … …  
    4696! Output
    4797!
    48 REAL, INTENT(OUT) :: Tcritcont    !
    49 REAL, INTENT(OUT) :: qcritcont    !
    50 REAL, INTENT(OUT) :: potcontfraP  !
    51 REAL, INTENT(OUT) :: potcontfraNP !
    52 REAL, INTENT(OUT) :: contfra      !
    53 REAL, INTENT(OUT) :: dcf_avi      !
    54 REAL, INTENT(OUT) :: dqvc_avi     !
    55 REAL, INTENT(OUT) :: dqi_avi      !
     98REAL, INTENT(OUT) :: Tcritcont    ! critical temperature for contrail formation [K]
     99REAL, INTENT(OUT) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
     100REAL, INTENT(OUT) :: potcontfraP  ! potential persistent contrail fraction [-]
     101REAL, INTENT(OUT) :: potcontfraNP ! potential non-persistent contrail fraction [-]
     102REAL, INTENT(OUT) :: contfra      ! contrail fraction [-]
     103REAL, INTENT(OUT) :: dcf_avi      ! cloud fraction tendency because of aviation [s-1]
     104REAL, INTENT(OUT) :: dqvc_avi     ! specific ice content tendency because of aviation [kg/kg/s]
     105REAL, INTENT(OUT) :: dqi_avi      ! specific cloud water vapor tendency because of aviation [kg/kg/s]
    56106!
    57107! Local
    … …  
    65115REAL :: qpotcontP
    66116!
    67 !
     117! for new contrail formation
    68118REAL :: contrail_cross_section, contfra_new
    69119
    70120qzero(:) = 0.
    71 
    72 !--more local variables for diagnostics
    73 !--values from Schumann, Meteorol Zeitschrift, 1996
    74 !--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen
    75 !--Qheat = 43.  /  50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen
    76 !REAL, PARAMETER :: EiH2O=1.25  !--emission index of water vapour for kerosene (kg kg-1)
    77 !REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1)
    78 !REAL, PARAMETER :: eta=0.3     !--average propulsion efficiency of the aircraft
    79 !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
    80 !--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010.
    81 !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2)
    82121
    83122!---------------------------------
    … …  
    86125!--Revised by Schumann (1996) and Rap et al. (2010)
    87126
     127!--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
     128!--in Pa / K. See Rap et al. (2010) in JGR.
    88129!--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
    89 !--in Pa / K
    90130Gcont = EI_H2O_aviation * RCPD * pplay &
    91131       / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
    92 !--critical T_LM below which no liquid contrail can form in exhaust
     132!--critical temperature below which no liquid contrail can form in exhaust
     133!--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2
    93134!--in Kelvins
    94135Tcritcont = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
    … …  
    105146
    106147
    107 IF ( temp .LT. Tcritcont ) THEN !--contrail formation is possible
     148IF ( ( ( 1. - cldfra ) .GT. eps ) .AND. ( temp .LT. Tcritcont ) ) THEN
    108149
    109150  pdf_x = qcritcont / qsatl * 100.
    … …  
    132173
    133174  potcontfraNP = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
    134   potcontfraP = MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
    135                     pdf_fra_above_qcritcont - pdf_fra_above_qnuc)
    136   qpotcontP = MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
    137                   pdf_q_above_qcritcont - pdf_q_above_qnuc)
     175  potcontfraP = MAX(0., MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
     176                    pdf_fra_above_qcritcont - pdf_fra_above_qnuc))
     177  qpotcontP = MAX(0., MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
     178                  pdf_q_above_qcritcont - pdf_q_above_qnuc))
    138179
    139180ELSE
    … …  
    143184
    144185ENDIF ! temp .LT. Tcritcont
    145 
    146 
    147 ! Ajout des émissions de H2O dues à l'aviation
    148 ! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
    149 ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 
    150 !      = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) )
    151 ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q)
    152 ! flight_h2O is in kg H2O / s / cell
    153 !
    154 !IF (ok_plane_h2o) THEN
    155 !  q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) )
    156 !ENDIF
    157186
    158187!--Convert existing contrail fraction into "natural" cirrus cloud fraction
    … …  
    165194dqvc_avi = 0.
    166195IF ( potcontfraP .GT. eps ) THEN
    167   contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
     196  contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA(dz)
    168197  contfra_new = MIN(1., flight_dist * dtime * contrail_cross_section / V_cell)
    169198  dcf_avi = potcontfraP * contfra_new
    … …  
    565594
    566595!**********************************************************************************
    567 FUNCTION contrail_cross_section_onera()
     596FUNCTION contrail_cross_section_onera(dz)
     597
     598USE lmdz_lscp_init, ONLY: initial_width_contrails
    568599
    569600IMPLICIT NONE
    570601
    571 ! input
    572 ! output
     602!
     603! Input
     604!
     605REAL :: dz ! cell width [m]
     606!
     607! Output
     608!
    573609REAL :: contrail_cross_section_onera ! [m2]
    574 ! local
    575 
    576 contrail_cross_section_onera = 200. * 200.
     610!
     611! Local
     612!
     613
     614contrail_cross_section_onera = initial_width_contrails * dz
    577615
    578616END FUNCTION contrail_cross_section_onera
     617
     618
     619SUBROUTINE read_aviation_emissions( &
     620        klon, klev, latitude_deg, longitude_deg, pplay, &
     621        flight_dist, flight_h2o &
     622        )
     623
     624IMPLICIT NONE
     625!
     626! Input
     627!
     628INTEGER, INTENT(IN)                        :: klon, klev       ! number of horizontal grid points and vertical levels
     629REAL,    INTENT(IN),  DIMENSION(klon)      :: latitude_deg     ! latitude of the grid points [deg]
     630REAL,    INTENT(IN),  DIMENSION(klon)      :: longitude_deg    ! longitude of the grid points [deg]
     631REAL,    INTENT(IN),  DIMENSION(klon,klev) :: pplay            ! layer pressure [Pa]
     632!
     633! Output
     634!
     635REAL,    INTENT(OUT), DIMENSION(klon,klev) :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
     636REAL,    INTENT(OUT), DIMENSION(klon,klev) :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
     637!
     638! Local
     639!
     640INTEGER :: i, k
     641
     642!--Initialisation
     643flight_dist(:,:) = 0.
     644flight_h2o(:,:) = 0.
     645
     646!DO i=1, klon
     647! IF ( ( latitude_deg(i) .GE. 42. ) .AND. ( latitude_deg(i) .LE. 48. ) ) THEN
     648!   flight_dist(i,14) = 50000.  !--5000 m of flight/second in grid cell x 10 scaling
     649!   flight_h2o(i,14) = 100.  !--10 kgH2O/second in grid cell x 10 scaling
     650! ENDIF
     651!ENDDO
     652
     653END SUBROUTINE read_aviation_emissions
    579654
    580655END MODULE lmdz_aviation
    … …  
    754829!!
    755830!END SUBROUTINE airplane
    756 !
    757 !!********************************************************************
    758 !! simple routine to initialise flight_m and test a flight corridor
    759 !!--Olivier Boucher - 2021
    760 !!
    761 !SUBROUTINE flight_init()
    762 !  USE dimphy
    763 !  USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
    764 !  IMPLICIT NONE
    765 !  INTEGER :: i
    766 !
    767 !  ALLOCATE(flight_m(klon,klev))
    768 !  ALLOCATE(flight_h2o(klon,klev))
    769 !  !
    770 !  flight_m(:,:) = 0.0    !--initialisation
    771 !  flight_h2o(:,:) = 0.0  !--initialisation
    772 !  !
    773 !  DO i=1, klon
    774 !   IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN
    775 !     flight_m(i,38) = 50000.0  !--5000 m of flight/second in grid cell x 10 scaling
    776 !   ENDIF
    777 !  ENDDO
    778 ! 
    779 !  RETURN
    780 !END SUBROUTINE flight_init
    781 !
    782 !     !--add a source of cirrus from aviation contrails
    783 !     IF (ok_plane_contrail) THEN
    784 !       drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell    !--tendency rneb due to aviation [s-1]
    785 !       drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime)                     !--majoration
    786 !       dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss)                          !--tendency q aviation [kg kg-1 s-1]
    787 !       rneb = rneb + drneb_avi(i,k)*dtime                                     !--add tendency to rneb
    788 !       qcld = qcld + dqss_avi*dtime                                           !--add tendency to qcld
    789 !       rnebss = rnebss - drneb_avi(i,k)*dtime                                 !--add tendency to rnebss
    790 !       qss = qss - dqss_avi*dtime                                             !--add tendency to qss
    791 !     ELSE
    792 !       drneb_avi(i,k)=0.0
    793 !     ENDIF
    794 !
    795 !  RETURN
    796 !END SUBROUTINE ice_sursat

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90

    r5452 r5453  
    175175  !--------------------------------------------------
    176176  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rcont_seri       ! ratio of contrails fraction to total cloud fraction [-]
    177   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
    178   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
     177  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
     178  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
    179179 
    180180  ! OUTPUT variables
    … …  
    10661066          rvc_seri(i,k) = qvc(i) / zq(i)
    10671067        ELSE
    1068           rvc_seri(i,k) = min_ratio
     1068          rvc_seri(i,k) = 0.
    10691069        ENDIF
    10701070        !--The MIN barrier is NEEDED because of:
    … …  
    10771077            rcont_seri(i,k) = contfra(i) / rneb(i,k)
    10781078          ELSE
    1079             rcont_seri(i,k) = min_ratio
     1079            rcont_seri(i,k) = 0.
    10801080          ENDIF
    10811081          !--This barrier should never be activated

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90

    r5452 r5453  
    157157!
    158158REAL,     INTENT(INOUT), DIMENSION(klon) :: rcont_seri    ! ratio of contrails fraction to total cloud fraction [-]
    159 REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   !
    160 REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    !
     159REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   ! aviation distance flown within the mesh [m/s/mesh]
     160REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
    161161!
    162162!  Output
    … …  
    238238!
    239239! for cell properties
    240 REAL :: rho, rhodz, dz
    241 REAL :: V_cell, M_cell
     240REAL :: rho, rhodz, dz, V_cell
    242241
    243242qzero(:) = 0.
    … …  
    336335      !--Cell volume [m3]
    337336      V_cell = dz * cell_area(i)
    338       !--Cell dry air mass [kg]
    339       M_cell = rhodz * cell_area(i)
    340337
    341338
    … …  
    504501        ENDIF
    505502        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
    506         pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3
    507         pdf_alpha = MIN(10., pdf_alpha)
     503        pdf_alpha = EXP( MIN(1000., rhl_clr) / 100. ) * pdf_e3
    508504       
    509505        IF ( ok_warm_cloud ) THEN
    … …  
    790786            dtime, pplay(i), temp(i), qsat(i), qsatl(i), gamma_cond(i), &
    791787            rcont_seri(i), flight_dist(i), cldfra(i), qvc(i), &
    792             V_cell, M_cell, pdf_loc, pdf_scale, pdf_alpha, &
     788            dz, V_cell, pdf_loc, pdf_scale, pdf_alpha, &
    793789            Tcritcont(i), qcritcont(i), potcontfraP(i), potcontfraNP(i), contfra(i), &
    794790            dcf_avi(i), dqvc_avi(i), dqi_avi(i) &

  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90

    r5452 r5453  
    236236  REAL, SAVE, PROTECTED :: linear_contrails_lifetime=10800.  ! [s] timescale of the lifetime of linear contrails
    237237  !$OMP THREADPRIVATE(linear_contrails_lifetime)
     238
     239  REAL, SAVE, PROTECTED :: initial_width_contrails=200.      ! [m] initial width of the linear contrails formed
     240  !$OMP THREADPRIVATE(initial_width_contrails)
    238241  !--End of the parameters for aviation
    239242
    … …  
    475478    CALL getin_p('prop_efficiency_aviation',prop_efficiency_aviation)
    476479    CALL getin_p('linear_contrails_lifetime',linear_contrails_lifetime)
     480    CALL getin_p('initial_width_contrails',initial_width_contrails)
    477481
    478482
    … …  
    566570    WRITE(lunout,*) 'lscp_ini, prop_efficiency_aviation:', prop_efficiency_aviation
    567571    WRITE(lunout,*) 'lscp_ini, linear_contrails_lifetime:', linear_contrails_lifetime
     572    WRITE(lunout,*) 'lscp_ini, initial_width_contrails:', initial_width_contrails
    568573
    569574

  • LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90

    r5452 r5453  
    665665
    666666!-- LSCP - aviation and contrails variables
     667      REAL, SAVE, ALLOCATABLE :: d_q_avi(:,:)
     668      !$OMP THREADPRIVATE(d_q_avi)
    667669      REAL, SAVE, ALLOCATABLE :: rcont_seri(:,:), d_rcont_dyn(:,:)
    668670      !$OMP THREADPRIVATE(rcont_seri, d_rcont_dyn)
    … …  
    12151217
    12161218!-- LSCP - aviation and contrails variables
     1219      ALLOCATE(d_q_avi(klon,klev))
    12171220      ALLOCATE(rcont_seri(klon,klev), d_rcont_dyn(klon,klev))
    12181221      ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
    … …  
    16211624
    16221625!-- LSCP - aviation and contrails variables
    1623       DEALLOCATE(rcont_seri, d_rcont_dyn, flight_dist, flight_h2o)
     1626      DEALLOCATE(d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o)
    16241627      DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
    16251628      DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)

  • LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90

    r5452 r5453  
    21682168
    21692169!-- LSCP - aviation variables
     2170  TYPE(ctrl_out), SAVE :: o_dqavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2171    'dqavi', 'Water vapor emissions from aviation tendency', 'kg/kg/s', (/ ('',i=1,10) /))
     2172  TYPE(ctrl_out), SAVE :: o_rcontseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2173    'rcontseri', 'Contrails fraction to total cloud fraction ratio', '-', (/ ('',i=1,10) /))
     2174  TYPE(ctrl_out), SAVE :: o_drcontdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2175    'drcontdyn', 'Dynamics contrails fraction to total cloud fraction ratio tendency', 's-1', (/ ('',i=1,10) /))
    21702176  TYPE(ctrl_out), SAVE :: o_Tcritcont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21712177    'Tcritcont', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))

  • LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90

    r5452 r5453  
    226226         o_dqvcadj, o_dqvcsub, o_dqvccon, o_dqvcmix, o_qsatl, o_qsati, &
    227227!-- LSCP - aviation variables
     228         o_rcontseri, o_drcontdyn, o_dqavi, &
    228229         o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
    229230         o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
    … …  
    349350         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
    350351         qsatliq, qsatice, &
     352         rcont_seri, d_rcont_dyn, d_q_avi, &
    351353         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    352354         dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
    … …  
    21232125!-- LSCP - aviation variables
    21242126       IF (ok_plane_contrail) THEN
     2127         CALL histwrite_phy(o_rcontseri, rcont_seri)
     2128         CALL histwrite_phy(o_drcontdyn, d_rcont_dyn)
    21252129         CALL histwrite_phy(o_flight_dist, flight_dist)
    21262130         CALL histwrite_phy(o_Tcritcont, Tcritcont)
    … …  
    21342138       IF (ok_plane_h2o) THEN
    21352139         CALL histwrite_phy(o_flight_h2o, flight_h2o)
     2140         CALL histwrite_phy(o_dqavi, d_q_avi)
    21362141       ENDIF
    21372142       

  • LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

    r5452 r5453  
    7575    USE write_field_phy
    7676    use wxios_mod, ONLY: g_ctx, wxios_set_context
     77    USE lmdz_aviation, ONLY : read_aviation_emissions, aviation_water_emissions
    7778    USE lmdz_lscp, ONLY : lscp
    7879    USE lmdz_call_cloud_optics_prop, ONLY : call_cloud_optics_prop
    … …  
    327328       dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    328329       !-- LSCP - aviation and contrails variables
    329        rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
     330       d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
    330331       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    331332       dcf_avi, dqi_avi, dqvc_avi, &
    … …  
    38633864         ratqs,ratqsc,ratqs_inter_,sigma_qtherm)
    38643865
     3866    !--Read the aviation emissions
     3867    IF ( ok_plane_h2o .OR. ok_plane_contrail ) THEN
     3868      CALL read_aviation_emissions(klon, klev, latitude_deg, longitude_deg, pplay, &
     3869                                   flight_dist, flight_h2o)
     3870    ENDIF
     3871
     3872    !--Add the water emissions from aviation
     3873    IF ( ok_plane_h2o ) THEN
     3874       CALL aviation_water_emissions(klon, klev, phys_tstep, paprs, pplay, &
     3875            t_seri, q_seri, cell_area, flight_h2o, d_q_avi)
     3876       CALL add_phys_tend(du0, dv0, dt0, d_q_avi, dql0, dqi0, dqbs0, paprs, &
     3877            'avi', abortphy, flag_inhib_tend, itap, 0)
     3878       d_q_avi = d_q_avi / phys_tstep
     3879    ENDIF
     3880
    38653881    !
    38663882    ! Appeler le processus de condensation a grande echelle
    … …  
    38753891
    38763892    IF (ok_new_lscp) THEN
    3877 
    38783893 
    38793894    DO k = 1, klev
    … …  
    38833898      ENDDO
    38843899    ENDDO
    3885 
    3886 
    3887     !--mise à jour de flight_m et flight_h2o dans leur module
    3888     !IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
    3889     !  CALL airplane(debut,pphis,pplay,paprs,t_seri)
    3890     !ENDIF
    38913900
    38923901    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
Note: See TracChangeset for help on using the changeset viewer.