MODULE ice_sursat_mod IMPLICIT NONE REAL, SAVE, ALLOCATABLE :: flight_m(:,:) !--flown distance m s-1 per cell !$OMP THREADPRIVATE(flight_m) REAL, SAVE, ALLOCATABLE :: flight_h2o(:,:) !--emitted kg H2O s-1 per cell !$OMP THREADPRIVATE(flight_h2o) !--parameters for gamma function !--Karcher and Lohmann (2002) !--gamma = 2.583 - t / 207.83 ! REAL, SAVE, PARAMETER :: gamma0=2.583, Tgamma=207.83 !--Ren and MacKenzie (2005) reused by Kärcher !--gamma = 2.349 - t / 259.0 REAL, PARAMETER :: gamma0=2.349, Tgamma=259.0 ! !--number of clouds in cell (needs to be parametrized at some point) REAL, PARAMETER :: N_cld = 100. ! !--safety parameters for ERF function REAL, PARAMETER :: erf_lim = 5., eps = 1.e-10 ! !--Tuning parameters (and their default values) ! !--chi gère la répartition statistique de la longueur des frontières ! entre les zones nuages et ISSR/ciel clair sous-saturé. Gamme de valeur : ! chi > 1, je n'ai pas regardé de limite max (pour chi = 1, la longueur de ! la frontière entre ne nuage et l'ISSR est proportionnelle à la ! répartition ISSR/ciel clair sous-sat dans la maille, i.e. il n'y a pas ! de favorisation de la localisation de l'ISSR près de nuage. Pour chi = inf, ! le nuage n'est en contact qu'avec de l'ISSR, quelle que soit la taille ! de l'ISSR dans la maille.) ! !--l_turb est la longueur de mélange pour la turbulence. ! dans les tests, ça n'a jamais été modifié pour l'instant. ! !--tun_N est le paramètre qui contrôle l'importance relative de N_2 par rapport à N_1. ! La valeur est comprise entre 1 et 2 (tun_N = 1 => N_1 = N_2) ! !--tun_ratqs : paramètre qui modifie ratqs en fonction de la valeur de ! alpha_cld selon la formule ratqs_new = ratqs_old / ( 1 + tun_ratqs * ! alpha_cld ). Dans le rapport il est appelé beta. Il varie entre 0 et 5 ! (tun_ratqs = 0 => pas de modification de ratqs). REAL, SAVE :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0 !$OMP THREADPRIVATE(chi,l_turb,tun_N,tun_ratqs) ! !--contrail cross section, typical value found in Freudenthaler et al, GRL, 22, 3501-3504, 1995 !--in m2, 1000x200 = 200 000 m2 after 15 min REAL, PARAMETER :: contrail_cross_section=200000.0 CONTAINS !******************************************************************* ! SUBROUTINE ice_sursat_init() USE print_control_mod, ONLY: lunout USE ioipsl_getin_p_mod, ONLY : getin_p IMPLICIT NONE CALL getin_p('flag_chi',chi) CALL getin_p('flag_l_turb',l_turb) CALL getin_p('flag_tun_N',tun_N) CALL getin_p('flag_tun_ratqs',tun_ratqs) WRITE(lunout,*) 'flag_chi = ', chi WRITE(lunout,*) 'flag_l_turb = ', l_turb WRITE(lunout,*) 'flag_tun_N = ', tun_N WRITE(lunout,*) 'flag_tun_ratqs = ', tun_ratqs END SUBROUTINE ice_sursat_init !******************************************************************* ! SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri) USE dimphy USE mod_grid_phy_lmdz, ONLY: klon_glo USE geometry_mod, ONLY: cell_area USE phys_cal_mod, ONLY : mth_cur USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root USE mod_phys_lmdz_omp_data, ONLY: is_omp_root USE mod_phys_lmdz_para, ONLY: scatter, bcast USE print_control_mod, ONLY: lunout IMPLICIT NONE INCLUDE "YOMCST.h" INCLUDE 'netcdf.inc' !-------------------------------------------------------- !--input variables !-------------------------------------------------------- LOGICAL, INTENT(IN) :: debut REAL, INTENT(IN) :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev) !-------------------------------------------------------- ! ... Local variables !-------------------------------------------------------- CHARACTER (LEN=20) :: modname='airplane_mod' INTEGER :: i, k, kori, iret, varid, error, ncida, klona INTEGER,SAVE :: nleva, ntimea !$OMP THREADPRIVATE(nleva,ntimea) REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:) !--km/s REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:) !--molec H2O/cm3/s REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:) REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:) REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:) !$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta) REAL :: zalt(klon,klev+1) REAL :: zrho, zdz(klon,klev), zfrac ! IF (debut) THEN !-------------------------------------------------------------------------------- ! ... Open the file and read airplane emissions !-------------------------------------------------------------------------------- ! IF (is_mpi_root .AND. is_omp_root) THEN ! iret = nf_open('aircraft_phy.nc', 0, ncida) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1) ! ... Get lengths iret = nf_inq_dimid(ncida, 'time', varid) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1) iret = nf_inq_dimlen(ncida, varid, ntimea) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1) iret = nf_inq_dimid(ncida, 'vector', varid) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1) iret = nf_inq_dimlen(ncida, varid, klona) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1) iret = nf_inq_dimid(ncida, 'lev', varid) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) iret = nf_inq_dimlen(ncida, varid, nleva) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1) ! IF ( klona /= klon_glo ) THEN WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1) ENDIF ! IF ( ntimea /= 12 ) THEN WRITE(lunout,*) 'ntimea=', ntimea CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1) ENDIF ! ALLOCATE(zmida(nleva), STAT=error) IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1) ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error) IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1) ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error) IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1) ! iret = nf_inq_varid(ncida, 'lev', varid) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) iret = nf_get_var_double(ncida, varid, zmida) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1) ! iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid) !--CO2 as a proxy for m flown - IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1) iret = nf_get_var_double(ncida, varid, pkm_airpl_glo) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1) ! iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1) iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo) IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1) ! ENDIF !--is_mpi_root and is_omp_root ! CALL bcast(nleva) CALL bcast(ntimea) ! IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error) IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error) ! ALLOCATE(pkm_airpl(klon,nleva,ntimea)) ALLOCATE(ph2o_airpl(klon,nleva,ntimea)) ! ALLOCATE(flight_m(klon,klev)) ALLOCATE(flight_h2o(klon,klev)) ! CALL bcast(zmida) zinta(1)=0.0 !--surface DO k=2, nleva zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0 !--conversion from km to m ENDDO zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface !print *,'zinta=', zinta ! CALL scatter(pkm_airpl_glo,pkm_airpl) CALL scatter(ph2o_airpl_glo,ph2o_airpl) ! !$OMP MASTER IF (is_mpi_root .AND. is_omp_root) THEN DEALLOCATE(pkm_airpl_glo) DEALLOCATE(ph2o_airpl_glo) ENDIF !--is_mpi_root !$OMP END MASTER ENDIF !--debut ! !--compute altitude of model level interfaces ! DO i = 1, klon zalt(i,1)=pphis(i)/RG !--in m ENDDO ! DO k=1, klev DO i = 1, klon zrho=pplay(i,k)/t_seri(i,k)/RD zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG zalt(i,k+1)=zalt(i,k)+zdz(i,k) !--in m ENDDO ENDDO ! !--vertical reprojection ! flight_m(:,:)=0.0 flight_h2o(:,:)=0.0 ! DO k=1, klev DO kori=1, nleva DO i=1, klon !--fraction of layer kori included in layer k zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori)) !--reproject flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac !--reproject flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac ENDDO ENDDO ENDDO ! DO k=1, klev DO i=1, klon !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3 flight_m(i,k)=flight_m(i,k)*100.0 !--x100 to augment signal to noise !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6 ENDDO ENDDO ! END SUBROUTINE airplane !******************************************************************** ! simple routine to initialise flight_m and test a flight corridor !--Olivier Boucher - 2021 ! SUBROUTINE flight_init() USE dimphy USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg IMPLICIT NONE INTEGER :: i ALLOCATE(flight_m(klon,klev)) ALLOCATE(flight_h2o(klon,klev)) ! flight_m(:,:) = 0.0 !--initialisation flight_h2o(:,:) = 0.0 !--initialisation ! DO i=1, klon IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN flight_m(i,38) = 50000.0 !--5000 m of flight/second in grid cell x 10 scaling ENDIF ENDDO RETURN END SUBROUTINE flight_init !******************************************************************* !--Routine to deal with ice supersaturation !--Determines the respective fractions of unsaturated clear sky, ice supersaturated clear sky and cloudy sky !--Diagnoses regions prone for non-persistent and persistent contrail formation ! !--Audran Borella - 2021 ! SUBROUTINE ice_sursat(pplay, dpaprs, dtime, i, k, t, q, gamma_ss, & qsat, t_actuel, rneb_seri, ratqs, rneb, qincld, & rnebss, qss, Tcontr, qcontr, qcontr2, fcontrN, fcontrP) ! USE dimphy USE print_control_mod, ONLY: prt_level, lunout USE phys_state_var_mod, ONLY: pbl_tke, t_ancien USE phys_local_var_mod, ONLY: N1_ss, N2_ss USE phys_local_var_mod, ONLY: drneb_sub, drneb_con, drneb_tur, drneb_avi !! USE phys_local_var_mod, ONLY: Tcontr, qcontr, fcontrN, fcontrP USE indice_sol_mod, ONLY: is_ave USE geometry_mod, ONLY: cell_area ! IMPLICIT NONE INCLUDE "YOMCST.h" INCLUDE "YOETHF.h" INCLUDE "FCTTRE.h" INCLUDE "fisrtilp.h" INCLUDE "clesphys.h" ! ! Input ! Beware: this routine works on a gridpoint! ! REAL, INTENT(IN) :: pplay ! layer pressure (Pa) REAL, INTENT(IN) :: dpaprs ! layer delta pressure (Pa) REAL, INTENT(IN) :: dtime ! intervalle du temps (s) REAL, INTENT(IN) :: t ! température advectée (K) REAL, INTENT(IN) :: qsat ! vapeur de saturation REAL, INTENT(IN) :: t_actuel ! temperature actuelle de la maille (K) REAL, INTENT(IN) :: rneb_seri ! fraction nuageuse en memoire INTEGER, INTENT(IN) :: i, k ! ! Input/output ! REAL, INTENT(INOUT) :: q ! vapeur de la maille (=zq) REAL, INTENT(INOUT) :: ratqs ! determine la largeur de distribution de vapeur REAL, INTENT(INOUT) :: Tcontr, qcontr, qcontr2, fcontrN, fcontrP ! ! Output ! REAL, INTENT(OUT) :: gamma_ss ! REAL, INTENT(OUT) :: rneb ! cloud fraction REAL, INTENT(OUT) :: qincld ! in-cloud total water REAL, INTENT(OUT) :: rnebss ! ISSR fraction REAL, INTENT(OUT) :: qss ! in-ISSR total water ! ! Local ! REAL PI PARAMETER (PI=4.*ATAN(1.)) REAL rnebclr, gamma_prec REAL qclr, qvc, qcld, qi REAL zrho, zdz, zrhodz REAL pdf_N, pdf_N1, pdf_N2 REAL pdf_a, pdf_b REAL pdf_e1, pdf_e2, pdf_k REAL drnebss, drnebclr, dqss, dqclr, sum_rneb_rnebss, dqss_avi REAL V_cell !--volume of the cell REAL M_cell !--dry mass of the cell REAL tke, sig, L_tur, b_tur, q_eq REAL V_env, V_cld, V_ss, V_clr REAL zcor ! !--more local variables for diagnostics !--imported from YOMCST.h !--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1) !--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air !--values from Schumann, Meteorol Zeitschrift, 1996 !--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen !--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1) REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute !--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010. REAL :: Gcontr !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) !--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1) REAL :: qsatliqcontr ! Initialisations zrho = pplay / t / RD !--dry density kg m-3 zrhodz = dpaprs / RG !--dry air mass kg m-2 zdz = zrhodz / zrho !--cell thickness m V_cell = zdz * cell_area(i) !--cell volume m3 M_cell = zrhodz * cell_area(i) !--cell dry air mass kg ! ! Recuperation de la memoire sur la couverture nuageuse rneb = rneb_seri ! ! Ajout des émissions de H2O dues à l'aviation ! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) ! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) ! flight_h2O is in kg H2O / s / cell ! IF (ok_plane_h2o) THEN q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) ENDIF ! !--Estimating gamma gamma_ss = gamma0 - t_actuel/Tgamma !gamma_prec = gamma0 - t_ancien(i,k)/Tgamma !--formulation initiale d Audran gamma_prec = gamma0 - t/Tgamma !--autre formulation possible basée sur le T du pas de temps ! ! Initialisation de qvc : q_sat du pas de temps precedent !qvc = R2ES*FOEEW(t_ancien(i,k),1.)/pplay !--formulation initiale d Audran qvc = R2ES*FOEEW(t,1.)/pplay !--autre formulation possible basée sur le T du pas de temps qvc = min(0.5,qvc) zcor = 1./(1.-RETV*qvc) qvc = qvc*zcor ! ! Modification de ratqs selon formule proposee : ksi_new = ksi_old/(1+beta*alpha_cld) ratqs = ratqs / (tun_ratqs*rneb_seri + 1.) ! ! Calcul de N pdf_k = -sqrt(log(1.+ratqs**2.)) pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) pdf_b = pdf_k/(2.*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) pdf_N = max(0.,sign(rneb,pdf_e1)) ELSE pdf_e1 = erf(pdf_e1) pdf_e1 = 0.5*(1.+pdf_e1) pdf_N = max(0.,rneb/pdf_e1) ENDIF ! ! On calcule ensuite N1 et N2. Il y a deux cas : N > 1 et N <= 1 ! Cas 1 : N > 1. N'arrive en theorie jamais, c'est une barriere ! On perd la memoire sur la temperature (sur qvc) pour garder ! celle sur alpha_cld IF (pdf_N.GT.1.) THEN ! On inverse alpha_cld = int_qvc^infty P(q) dq ! pour determiner qvc = f(alpha_cld) ! On approxime en serie entiere erf-1(x) qvc = 2.*rneb-1. qvc = qvc + PI/12.*qvc**3 + 7.*PI**2/480.*qvc**5 + 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 + 34807.*PI**5/182476800.*qvc**11 qvc = sqrt(PI)/2.*qvc qvc = (qvc-pdf_b)*pdf_k*sqrt(2.) qvc = q*exp(qvc) ! On met a jour rneb avec la nouvelle valeur de qvc ! La maj est necessaire a cause de la serie entiere approximative pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF pdf_e1 = 0.5*(1.+pdf_e1) rneb = pdf_e1 ! Si N > 1, N1 et N2 = 1 pdf_N1 = 1. pdf_N2 = 1. ! Cas 2 : N <= 1 ELSE ! D'abord on calcule N2 avec le tuning pdf_N2 = min(1.,pdf_N*tun_N) ! Puis N1 pour assurer la conservation de alpha_cld pdf_a = log(qvc*gamma_prec/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF pdf_e2 = 0.5*(1.+pdf_e2) ! integrale sous P pour q > gamma qsat IF (abs(pdf_e1-pdf_e2).LT.eps) THEN pdf_N1 = pdf_N2 ELSE pdf_N1 = (rneb-pdf_N2*pdf_e2)/(pdf_e1-pdf_e2) ENDIF ! Barriere qui traite le cas gamma_prec = 1. IF (pdf_N1.LE.0.) THEN pdf_N1 = 0. IF (pdf_e2.GT.eps) THEN pdf_N2 = rneb/pdf_e2 ELSE pdf_N2 = 0. ENDIF ENDIF ENDIF ! Physique 1 ! Sublimation IF (qvc.LT.qsat) THEN pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF pdf_e1 = 0.5*pdf_N1*(pdf_e1-pdf_e2) ! Calcul et ajout de la tendance drneb_sub(i,k) = - pdf_e1/dtime !--unit [s-1] rneb = rneb + drneb_sub(i,k)*dtime ELSE drneb_sub(i,k) = 0. ENDIF ! NOTE : verifier si ca marche bien pour gamma proche de 1. ! Condensation IF (gamma_ss*qsat.LT.gamma_prec*qvc) THEN pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF pdf_e1 = 0.5*(1.-pdf_N1)*(pdf_e1-pdf_e2) pdf_e2 = 0.5*(1.-pdf_N2)*(1.+pdf_e2) ! Calcul et ajout de la tendance drneb_con(i,k) = (pdf_e1 + pdf_e2)/dtime !--unit [s-1] rneb = rneb + drneb_con(i,k)*dtime ELSE pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF pdf_e1 = 0.5*(1.-pdf_N2)*(1.+pdf_e1) ! Calcul et ajout de la tendance drneb_con(i,k) = pdf_e1/dtime !--unit [s-1] rneb = rneb + drneb_con(i,k)*dtime ENDIF ! Calcul des grandeurs diagnostiques ! Determination des grandeurs ciel clair pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF pdf_e1 = 0.5*(1.-pdf_e1) pdf_e2 = pdf_a-pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF pdf_e2 = 0.5*q*(1.-pdf_e2) rnebclr = pdf_e1 qclr = pdf_e2 ! Determination de q_cld ! Partie 1 pdf_a = log(max(qsat,qvc)/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a-pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF pdf_a = log(min(gamma_ss*qsat,gamma_prec*qvc)/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a-pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF pdf_e1 = 0.5*q*pdf_N1*(pdf_e1-pdf_e2) qcld = pdf_e1 ! Partie 2 (sous condition) IF (gamma_ss*qsat.GT.gamma_prec*qvc) THEN pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a-pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF pdf_e2 = 0.5*q*pdf_N2*(pdf_e2-pdf_e1) qcld = qcld + pdf_e2 pdf_e2 = pdf_e1 ENDIF ! Partie 3 pdf_e2 = 0.5*q*(1.+pdf_e2) qcld = qcld + pdf_e2 ! Fin du calcul de q_cld ! Determination des grandeurs ISSR via les equations de conservation rneb=MIN(rneb, 1. - rnebclr - eps) !--ajout OB - barrière rnebss = MAX(0.0, 1. - rnebclr - rneb) !--ajout OB qss = MAX(0.0, q - qclr - qcld) !--ajout OB ! Physique 2 : Turbulence IF (rneb.GT.eps.AND.rneb.LT.1.-eps) THEN ! rneb != 0 and != 1 ! tke = pbl_tke(i,k,is_ave) ! A MODIFIER formule a revoir L_tur = min(l_turb, sqrt(tke)*dtime) ! On fait pour l'instant l'hypothese a = 3b. V = 4/3 pi a b**2 = alpha_cld ! donc b = alpha_cld/4pi **1/3. b_tur = (rneb*V_cell/4./PI/N_cld)**(1./3.) ! On verifie que la longeur de melange n'est pas trop grande IF (L_tur.GT.b_tur) THEN L_tur = b_tur ENDIF V_env = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. + 3.*b_tur*(L_tur**2.)) V_cld = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. - 3.*b_tur*(L_tur**2.)) V_cld = abs(V_cld) ! Repartition statistique des zones de contact nuage-ISSR et nuage-ciel clair sig = rnebss/(chi*rnebclr+rnebss) V_ss = MIN(sig*V_env,rnebss*V_cell) V_clr = MIN((1.-sig)*V_env,rnebclr*V_cell) V_cld = MIN(V_cld,rneb*V_cell) V_env = V_ss + V_clr ! ISSR => rneb drnebss = -1.*V_ss/V_cell dqss = drnebss*qss/MAX(eps,rnebss) ! Clear sky <=> rneb q_eq = V_env*qclr/MAX(eps,rnebclr) + V_cld*qcld/MAX(eps,rneb) q_eq = q_eq/(V_env + V_cld) IF (q_eq.GT.qsat) THEN drnebclr = - V_clr/V_cell dqclr = drnebclr*qclr/MAX(eps,rnebclr) ELSE drnebclr = V_cld/V_cell dqclr = drnebclr*qcld/MAX(eps,rneb) ENDIF ! Maj des variables avec les tendances rnebclr = MAX(0.0,rnebclr + drnebclr) !--OB ajout d'un max avec eps (il faudrait modified drnebclr pour le diag) qclr = MAX(0.0, qclr + dqclr) !--OB ajout d'un max avec 0 rneb = rneb - drnebclr - drnebss qcld = qcld - dqclr - dqss rnebss = MAX(0.0,rnebss + drnebss) !--OB ajout d'un max avec eps (il faudrait modifier drnebss pour le diag) qss = MAX(0.0, qss + dqss) !--OB ajout d'un max avec 0 ! Tendances pour le diagnostic drneb_tur(i,k) = (drnebclr + drnebss)/dtime !--unit [s-1] ENDIF ! rneb !--add a source of cirrus from aviation contrails IF (ok_plane_contrail) THEN drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell !--tendency rneb due to aviation [s-1] drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime) !--majoration dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss) !--tendency q aviation [kg kg-1 s-1] rneb = rneb + drneb_avi(i,k)*dtime !--add tendency to rneb qcld = qcld + dqss_avi*dtime !--add tendency to qcld rnebss = rnebss - drneb_avi(i,k)*dtime !--add tendency to rnebss qss = qss - dqss_avi*dtime !--add tendency to qss ELSE drneb_avi(i,k)=0.0 ENDIF ! Barrieres ! ISSR trop petite IF (rnebss.LT.eps) THEN rneb = MIN(rneb + rnebss,1.0-eps) !--ajout OB barriere qcld = qcld + qss rnebss = 0. qss = 0. ENDIF ! le nuage est trop petit IF (rneb.LT.eps) THEN ! s'il y a une ISSR on met tout dans l'ISSR, sinon dans le ! clear sky IF (rnebss.LT.eps) THEN rnebclr = 1. rnebss = 0. !--ajout OB qclr = q ELSE rnebss = MIN(rnebss + rneb,1.0-eps) !--ajout OB barriere qss = qss + qcld ENDIF rneb = 0. qcld = 0. qincld = qsat * gamma_ss ELSE qincld = qcld / rneb ENDIF !--OB ajout borne superieure sum_rneb_rnebss=rneb+rnebss rneb=rneb*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss) rnebss=rnebss*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss) ! On ecrit dans la memoire N1_ss(i,k) = pdf_N1 N2_ss(i,k) = pdf_N2 !--Diagnostics only used from last iteration !--test !!Tcontr(i,k)=200. !!fcontrN(i,k)=1.0 !!fcontrP(i,k)=0.5 ! !--slope of dilution line in exhaust !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1) Gcontr = EiH2O * RCPD * pplay / (eps_w*Qheat*(1.-eta)) !--Pa K-1 !--critical T_LM below which no liquid contrail can form in exhaust !Tcontr(i,k) = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K IF (Gcontr .GT. 0.1) THEN ! Tcontr = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K !print *,'Tcontr=',iter,i,k,eps_w,pplay,Gcontr,Tcontr(i,k) !--Psat with index 0 in FOEEW to get saturation wrt liquid water corresponding to Tcontr !qsatliqcontr = RESTT*FOEEW(Tcontr(i,k),0.) !--Pa qsatliqcontr = RESTT*FOEEW(Tcontr,0.) !--Pa !--Critical water vapour above which there is contrail formation for ambiant temperature !qcontr(i,k) = Gcontr*(t-Tcontr(i,k)) + qsatliqcontr !--Pa qcontr = Gcontr*(t-Tcontr) + qsatliqcontr !--Pa !--Conversion of qcontr in specific humidity - method 1 !qcontr(i,k) = RD/RV*qcontr(i,k)/pplay !--so as to return to something similar to R2ES*FOEEW/pplay qcontr2 = RD/RV*qcontr/pplay !--so as to return to something similar to R2ES*FOEEW/pplay !qcontr(i,k) = min(0.5,qcontr(i,k)) !--and then we apply the same correction term as for qsat qcontr2 = min(0.5,qcontr2) !--and then we apply the same correction term as for qsat !zcor = 1./(1.-RETV*qcontr(i,k)) !--for consistency with qsat but is it correct at all? zcor = 1./(1.-RETV*qcontr2) !--for consistency with qsat but is it correct at all as p is dry? !zcor = 1./(1.+qcontr2) !--for consistency with qsat !qcontr(i,k) = qcontr(i,k)*zcor qcontr2 = qcontr2*zcor qcontr2=MAX(1.e-10,qcontr2) !--eliminate negative values due to extrapolation on dilution curve !--Conversion of qcontr in specific humidity - method 2 !qcontr(i,k) = eps_w*qcontr(i,k) / (pplay+eps_w*qcontr(i,k)) !qcontr=MAX(1.E-10,qcontr) !qcontr2 = eps_w*qcontr / (pplay+eps_w*qcontr) ! IF (t .LT. Tcontr) THEN !--contrail formation is possible ! !--compute fractions of persistent (P) and non-persistent(N) contrail potential regions !!IF (qcontr(i,k).GE.qsat) THEN IF (qcontr2.GE.qsat) THEN !--none of the unsaturated clear sky is prone for contrail formation !!fcontrN(i,k) = 0.0 fcontrN = 0.0 ! !--integral of P(q) from qsat to qcontr in ISSR pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF ! !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.)) pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF ! !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2)) fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2)) ! pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF ! !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.)) pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF ! !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2)) fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2)) ! pdf_a = log(MAX(qsat,qvc)/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF ! !!pdf_a = log(MIN(qcontr(i,k),MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.)) pdf_a = log(MIN(qcontr2,MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF ! !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2)) fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2)) ! pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF ! !!pdf_a = log(MIN(qcontr(i,k),gamma_ss*qsat)/q)/(pdf_k*sqrt(2.)) pdf_a = log(MIN(qcontr2,gamma_ss*qsat)/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF ! !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2)) fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2)) ! ELSE !--qcontr LT qsat ! !--all of ISSR is prone for contrail formation !!fcontrP(i,k) = rnebss fcontrP = rnebss ! !--integral of zq from qcontr to qsat in unsaturated clear-sky region !!pdf_a = log(qcontr(i,k)/q)/(pdf_k*sqrt(2.)) pdf_a = log(qcontr2/q)/(pdf_k*sqrt(2.)) pdf_e1 = pdf_a+pdf_b !--normalement pdf_b est deja defini IF (abs(pdf_e1).GE.erf_lim) THEN pdf_e1 = sign(1.,pdf_e1) ELSE pdf_e1 = erf(pdf_e1) ENDIF ! pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) pdf_e2 = pdf_a+pdf_b IF (abs(pdf_e2).GE.erf_lim) THEN pdf_e2 = sign(1.,pdf_e2) ELSE pdf_e2 = erf(pdf_e2) ENDIF ! !!fcontrN(i,k) = 0.5*(pdf_e1-pdf_e2) fcontrN = 0.5*(pdf_e1-pdf_e2) !!fcontrN=2.0 ! ENDIF ! ENDIF !-- t < Tcontr ! ENDIF !-- Gcontr > 0.1 RETURN END SUBROUTINE ice_sursat ! !******************************************************************* ! END MODULE ice_sursat_mod