Index: LMDZ6/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- LMDZ6/trunk/DefLists/field_def_lmdz.xml (revision 4058)
+++ LMDZ6/trunk/DefLists/field_def_lmdz.xml (revision 4059)
@@ -758,4 +758,27 @@
pres*ovap*461.5 / (287.04*(1.+ (10.9491/18.0153)*ovap))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: LMDZ6/trunk/DefLists/file_def_histday_lmdz.xml
===================================================================
--- LMDZ6/trunk/DefLists/file_def_histday_lmdz.xml (revision 4058)
+++ LMDZ6/trunk/DefLists/file_def_histday_lmdz.xml (revision 4059)
@@ -638,4 +638,27 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: LMDZ6/trunk/DefLists/file_def_histhf_lmdz.xml
===================================================================
--- LMDZ6/trunk/DefLists/file_def_histhf_lmdz.xml (revision 4058)
+++ LMDZ6/trunk/DefLists/file_def_histhf_lmdz.xml (revision 4059)
@@ -485,5 +485,4 @@
-
@@ -500,5 +499,5 @@
-
+
@@ -654,4 +653,25 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: LMDZ6/trunk/libf/phylmd/clesphys.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/clesphys.h (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/clesphys.h (revision 4059)
@@ -85,5 +85,5 @@
LOGICAL :: ok_cosp,ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP
LOGICAL :: ok_airs
- INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
+ INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, iflag_ice_sursat, NSW, iflag_albedo
LOGICAL :: ok_chlorophyll
LOGICAL :: ok_strato
@@ -123,5 +123,5 @@
! THEN INTEGER AND LOGICALS
& , top_height &
- & , iflag_cycle_diurne, soil_model, new_oliq &
+ & , iflag_cycle_diurne, soil_model, new_oliq &
& , ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad &
& , iflag_con, nbapp_cv, nbapp_wk &
@@ -139,5 +139,6 @@
& , ok_lic_melt, ok_lic_cond, aer_type &
& , iflag_rrtm, ok_strato,ok_hines, ok_qch4 &
- & , iflag_ice_thermo, ok_gwd_rando, NSW, iflag_albedo &
+ & , iflag_ice_thermo, iflag_ice_sursat &
+ & , ok_gwd_rando, NSW, iflag_albedo &
& , ok_chlorophyll,ok_conserv_q, adjust_tropopause &
& , ok_daily_climoz, ok_all_xml, ok_lwoff &
Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90 (revision 4059)
@@ -179,4 +179,5 @@
INTEGER,SAVE :: iflag_pdf_omp
INTEGER,SAVE :: iflag_ice_thermo_omp
+ INTEGER,SAVE :: iflag_ice_sursat_omp
INTEGER,SAVE :: iflag_t_glace_omp
INTEGER,SAVE :: iflag_cloudth_vert_omp
@@ -1493,6 +1494,4 @@
CALL getin('iflag_vice',iflag_vice_omp)
-
-
!
!Config Key = iflag_ice_thermo
@@ -1504,4 +1503,14 @@
CALL getin('iflag_ice_thermo',iflag_ice_thermo_omp)
+ !
+ !Config Key = iflag_ice_sursat
+ !Config Desc =
+ !Config Def = 0
+ !Config Help =
+ !
+ iflag_ice_sursat_omp = 0
+ CALL getin('iflag_ice_sursat',iflag_ice_sursat_omp)
+
+ !
!Config Key = rei_min
!Config Desc =
@@ -2491,4 +2500,5 @@
iflag_vice=iflag_vice_omp
iflag_ice_thermo = iflag_ice_thermo_omp
+ iflag_ice_sursat = iflag_ice_sursat_omp
rei_min = rei_min_omp
rei_max = rei_max_omp
@@ -2918,4 +2928,5 @@
WRITE(lunout,*) ' iflag_vice = ',iflag_vice
WRITE(lunout,*) ' iflag_ice_thermo = ',iflag_ice_thermo
+ WRITE(lunout,*) ' iflag_ice_sursat = ',iflag_ice_sursat
WRITE(lunout,*) ' rei_min = ',rei_min
WRITE(lunout,*) ' rei_max = ',rei_max
Index: LMDZ6/trunk/libf/phylmd/ice_sursat_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/ice_sursat_mod.F90 (revision 4059)
+++ LMDZ6/trunk/libf/phylmd/ice_sursat_mod.F90 (revision 4059)
@@ -0,0 +1,847 @@
+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
+!
+!--Parametres de tuning
+REAL, PARAMETER :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0
+!
+! Contrail cirrus
+LOGICAL, PARAMETER :: contrail_cirrus=.FALSE., h2o_emissions=.FALSE.
+!
+!--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 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
+
+!********************************************************************
+! temporary routine to initialise flight_m and test a flight corridor
+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"
+ !
+ ! 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 (h2o_emissions) 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 (contrail_cirrus) 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
+ 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
+
+ RETURN
+END SUBROUTINE ice_sursat
+!
+!*******************************************************************
+!
+END MODULE ice_sursat_mod
Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (revision 4059)
@@ -6,12 +6,13 @@
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE LSCP(dtime,paprs,pplay,t,q,ptconv,ratqs, &
- d_t, d_q, d_ql, d_qi, rneb, radliq, radicefrac, &
- rain, snow, &
+SUBROUTINE LSCP(dtime,missing_val, &
+ paprs,pplay,t,q,ptconv,ratqs, &
+ d_t, d_q, d_ql, d_qi, rneb, rneb_seri, &
+ radliq, radicefrac, rain, snow, &
pfrac_impa, pfrac_nucl, pfrac_1nucl, &
frac_impa, frac_nucl, beta, &
prfl, psfl, rhcl, zqta, fraca, &
ztv, zpspsk, ztla, zthl, iflag_cld_th, &
- iflag_ice_thermo)
+ iflag_ice_thermo, iflag_ice_sursat)
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -93,9 +94,11 @@
USE phys_local_var_mod, ONLY: rneblsvol
USE lscp_tools_mod, ONLY : CALC_QSAT_ECMWF, ICEFRAC_LSCP, CALC_GAMMASAT, FALLICE_VELOCITY
-
+USE ice_sursat_mod
+!--ice supersaturation
+USE phys_local_var_mod, ONLY: zqsats, zqsatl
+USE phys_local_var_mod, ONLY: qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss
+USE phys_local_var_mod, ONLY: Tcontr, qcontr, qcontr2, fcontrN, fcontrP
IMPLICIT NONE
-
-
!===============================================================================
@@ -114,4 +117,6 @@
REAL, INTENT(IN) :: dtime ! time step [s]
+ REAL, INTENT(IN) :: missing_val ! missing value for output
+
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
@@ -121,4 +126,6 @@
INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
! CR: if iflag_ice_thermo=2, only convection is active
+ INTEGER, INTENT(IN) :: iflag_ice_sursat ! 0 = sursat desativee, 1 = sursat activee
+
LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
@@ -138,4 +145,7 @@
REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale
! cloud PDF (sigma=ratqs*qt)
+
+ ! Input sursaturation en glace
+ REAL, DIMENSION(klon,klev), INTENT(INOUT):: rneb_seri ! fraction nuageuse en memoire
! OUTPUT variables
@@ -388,4 +398,13 @@
d_tot_zneb(:) = 0.0
+!--ice sursaturation
+gamma_ss(:,:) = 1.
+qss(:,:) = 0.
+rnebss(:,:) = 0.
+Tcontr(:,:) = missing_val
+qcontr(:,:) = missing_val
+qcontr2(:,:) = missing_val
+fcontrN(:,:) = 0.0
+fcontrP(:,:) = 0.0
!===============================================================================
@@ -645,6 +664,5 @@
qtot=zq(i)+zmqc(i)
CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
- zdqsdT_raw(i) = zdqs(i)* &
- & RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
+ zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
IF (zq(i) .LT. 1.e-15) THEN
@@ -677,6 +695,4 @@
qincloud_mpc(:)=0.
-
-
IF (iflag_cld_th.GE.5) THEN
@@ -778,5 +794,4 @@
! new temperature:
Tbef(i)=Tbef(i)+DT(i)
-
! Rneb, qzn and zcond for lognormal PDFs
@@ -800,11 +815,31 @@
zpdf_e2(i)=1.-erf(zpdf_e2(i))
- IF (zpdf_e1(i).LT.1.e-10) THEN
- rneb(i,k)=0.
- zqn(i)=gammasat(i)*zqs(i)
+ !--ice sursaturation by Audran
+ IF ((iflag_ice_sursat.EQ.0).OR.(Tbef(i).GT.t_glace_min)) THEN
+
+ IF (zpdf_e1(i).LT.1.e-10) THEN
+ rneb(i,k)=0.
+ zqn(i)=gammasat(i)*zqs(i)
+ ELSE
+ rneb(i,k)=0.5*zpdf_e1(i)
+ zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
+ ENDIF
+
+ rnebss(i,k)=0.0 !--ajout OB (necessaire car boucle de convergence sur le temps)
+ fcontrN(i,k)=0.0 !--idem
+ fcontrP(i,k)=0.0 !--idem
+ qss(i,k)=0.0 !--idem
+
ELSE
- rneb(i,k)=0.5*zpdf_e1(i)
- zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
- ENDIF
+ !------------------------------------
+ ! SURSATURATION EN GLACE
+ !------------------------------------
+
+ CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
+ gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), &
+ rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), &
+ Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
+
+ ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
! If vertical heterogeneity, change fraction by volume as well
@@ -823,5 +858,4 @@
! EV: calculation of icefrac in one sole function
CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
-
IF (zfice(i).LT.1) THEN
@@ -985,5 +1019,4 @@
! remaining water in the cloud during the time step that is seen by the radiation
! -------------------------------------------------------------------------------
-
DO n = 1, ninter
@@ -1128,7 +1161,4 @@
ENDDO
-
-
-
! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
! if iflag_evap_pre=4
@@ -1137,5 +1167,4 @@
DO i=1, klon
-
IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
@@ -1144,5 +1173,4 @@
znebprecipclr(i)=0.
ENDIF
-
IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
@@ -1152,5 +1180,4 @@
znebprecipcld(i)=0.
ENDIF
-
ENDDO
@@ -1190,10 +1217,7 @@
ENDIF
- zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
- * (paprs(i,k)-paprs(i,k+1))/RG
-
+ zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
-
IF (t(i,k) .GE. t_glace_min) THEN
@@ -1202,5 +1226,4 @@
zalpha_tr = a_tr_sca(4)
ENDIF
-
zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
@@ -1241,11 +1264,28 @@
ENDDO
-
-END DO
+ !--save some variables for ice sursaturation
+ !
+ DO i = 1, klon
+ ! pour la mémoire
+ rneb_seri(i,k) = rneb(i,k)
+
+ ! pour les diagnostics
+ rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
+
+ qvc(i,k) = zqs(i) * rneb(i,k)
+ qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--ajout OB a cause de cas pathologiques avec lognormale=F
+ qcld(i,k) = qvc(i,k) + zcond(i)
+
+ !q_sat
+ CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,1,.false.,zqsatl(i,k),zdqs(i))
+ CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,2,.false.,zqsats(i,k),zdqs(i))
+
+ ENDDO
+
+ENDDO
!======================================================================
! END OF VERTICAL LOOP
!======================================================================
-
! Rain or snow at the surface (depending on the first layer temperature)
@@ -1254,6 +1294,4 @@
rain(i) = zrfl(i)
ENDDO
-
-
IF (ncoreczq>0) THEN
@@ -1261,7 +1299,4 @@
ENDIF
-
-
-
END SUBROUTINE LSCP
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Index: LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 (revision 4059)
@@ -123,5 +123,5 @@
USE print_control_mod, ONLY: lunout, prt_level
- IMPLICIT none
+ IMPLICIT NONE
@@ -202,5 +202,4 @@
ENDDO
-
RETURN
@@ -215,6 +214,5 @@
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
- IMPLICIT none
+ IMPLICIT NONE
include "YOMCST.h"
@@ -236,7 +234,5 @@
REAL, INTENT(OUT) :: dqs ! derivation of saturation specific humidity wrt T
-
REAL delta, cor, cvm5
-
IF (phase .EQ. 1) THEN
@@ -261,6 +257,4 @@
dqs= FOEDE(temp,delta,cvm5,qs,cor)
-
-
END SUBROUTINE CALC_QSAT_ECMWF
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -278,5 +272,5 @@
- IMPLICIT none
+ IMPLICIT NONE
include "YOMCST.h"
@@ -341,7 +335,4 @@
ENDIF
-
-
-
END SUBROUTINE CALC_GAMMASAT
Index: LMDZ6/trunk/libf/phylmd/phyetat0.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyetat0.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/phyetat0.F90 (revision 4059)
@@ -13,5 +13,5 @@
du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, &
- ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, radpas, radsol, rain_fall, ratqs, &
+ ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, rneb_ancien, radpas, radsol, rain_fall, ratqs, &
rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
@@ -361,4 +361,5 @@
ancien_ok=ancien_ok.AND.phyetat0_get(klev,ql_ancien,"QLANCIEN","QLANCIEN",0.)
ancien_ok=ancien_ok.AND.phyetat0_get(klev,qs_ancien,"QSANCIEN","QSANCIEN",0.)
+ ancien_ok=ancien_ok.AND.phyetat0_get(klev,rneb_ancien,"RNEBANCIEN","RNEBANCIEN",0.)
ancien_ok=ancien_ok.AND.phyetat0_get(klev,u_ancien,"UANCIEN","UANCIEN",0.)
ancien_ok=ancien_ok.AND.phyetat0_get(klev,v_ancien,"VANCIEN","VANCIEN",0.)
@@ -373,4 +374,5 @@
(maxval(ql_ancien).EQ.minval(ql_ancien)) .OR. &
(maxval(qs_ancien).EQ.minval(qs_ancien)) .OR. &
+ (maxval(rneb_ancien).EQ.minval(rneb_ancien)) .OR. &
(maxval(prw_ancien).EQ.minval(prw_ancien)) .OR. &
(maxval(prlw_ancien).EQ.minval(prlw_ancien)) .OR. &
Index: LMDZ6/trunk/libf/phylmd/phyredem.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyredem.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/phyredem.F90 (revision 4059)
@@ -19,5 +19,5 @@
zval, rugoro, t_ancien, q_ancien, &
prw_ancien, prlw_ancien, prsw_ancien, &
- ql_ancien, qs_ancien, u_ancien, &
+ ql_ancien, qs_ancien, rneb_ancien, u_ancien, &
v_ancien, clwcon, rnebcon, ratqs, pbl_tke, &
wake_delta_pbl_tke, zmax0, f0, sig1, w01, &
@@ -239,4 +239,6 @@
CALL put_field(pass,"QSANCIEN", "QSANCIEN", qs_ancien)
+
+ CALL put_field(pass,"RNEBANCIEN", "RNEBANCIEN", rneb_ancien)
CALL put_field(pass,"PRWANCIEN", "PRWANCIEN", prw_ancien)
Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 4059)
@@ -16,4 +16,8 @@
REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
!$OMP THREADPRIVATE(u_seri, v_seri)
+ REAL, SAVE, ALLOCATABLE :: rneb_seri(:,:)
+ !$OMP THREADPRIVATE(rneb_seri)
+ REAL, SAVE, ALLOCATABLE :: d_rneb_dyn(:,:)
+ !$OMP THREADPRIVATE(d_rneb_dyn)
REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),tke_dissip(:,:,:),wprime(:,:,:)
!$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip,wprime)
@@ -473,4 +477,45 @@
!$OMP THREADPRIVATE(zn2mout)
+ REAL, SAVE, ALLOCATABLE :: qclr(:,:)
+ !$OMP THREADPRIVATE(qclr)
+ REAL, SAVE, ALLOCATABLE :: qcld(:,:)
+ !$OMP THREADPRIVATE(qcld)
+ REAL, SAVE, ALLOCATABLE :: qss(:,:)
+ !$OMP THREADPRIVATE(qss)
+ REAL, SAVE, ALLOCATABLE :: qvc(:,:)
+ !$OMP THREADPRIVATE(qvc)
+ REAL, SAVE, ALLOCATABLE :: rnebclr(:,:)
+ !$OMP THREADPRIVATE(rnebclr)
+ REAL, SAVE, ALLOCATABLE :: rnebss(:,:)
+ !$OMP THREADPRIVATE(rnebss)
+ REAL, SAVE, ALLOCATABLE :: gamma_ss(:,:)
+ !$OMP THREADPRIVATE(gamma_ss)
+ REAL, SAVE, ALLOCATABLE :: N1_ss(:,:)
+ !$OMP THREADPRIVATE(N1_ss)
+ REAL, SAVE, ALLOCATABLE :: N2_ss(:,:)
+ !$OMP THREADPRIVATE(N2_ss)
+ REAL, SAVE, ALLOCATABLE :: drneb_sub(:,:)
+ !$OMP THREADPRIVATE(drneb_sub)
+ REAL, SAVE, ALLOCATABLE :: drneb_con(:,:)
+ !$OMP THREADPRIVATE(drneb_con)
+ REAL, SAVE, ALLOCATABLE :: drneb_tur(:,:)
+ !$OMP THREADPRIVATE(drneb_tur)
+ REAL, SAVE, ALLOCATABLE :: drneb_avi(:,:)
+ !$OMP THREADPRIVATE(drneb_avi)
+ REAL, SAVE, ALLOCATABLE :: zqsatl(:,:)
+ !$OMP THREADPRIVATE(zqsatl)
+ REAL, SAVE, ALLOCATABLE :: zqsats(:,:)
+ !$OMP THREADPRIVATE(zqsats)
+ REAL, SAVE, ALLOCATABLE :: Tcontr(:,:)
+ !$OMP THREADPRIVATE(Tcontr)
+ REAL, SAVE, ALLOCATABLE :: qcontr(:,:)
+ !$OMP THREADPRIVATE(qcontr)
+ REAL, SAVE, ALLOCATABLE :: qcontr2(:,:)
+ !$OMP THREADPRIVATE(qcontr2)
+ REAL, SAVE, ALLOCATABLE :: fcontrN(:,:)
+ !$OMP THREADPRIVATE(fcontrN)
+ REAL, SAVE, ALLOCATABLE :: fcontrP(:,:)
+ !$OMP THREADPRIVATE(fcontrP)
+
#ifdef CPP_StratAer
!
@@ -828,4 +873,14 @@
ALLOCATE(zn2mout(klon,6))
+
+! Supersaturation
+ ALLOCATE(rneb_seri(klon,klev))
+ ALLOCATE(d_rneb_dyn(klon,klev))
+ ALLOCATE(qclr(klon,klev), qcld(klon,klev), qss(klon,klev), qvc(klon,klev))
+ ALLOCATE(rnebclr(klon,klev), rnebss(klon,klev), gamma_ss(klon,klev))
+ ALLOCATE(N1_ss(klon,klev), N2_ss(klon,klev))
+ ALLOCATE(drneb_sub(klon,klev), drneb_con(klon,klev), drneb_tur(klon,klev), drneb_avi(klon,klev))
+ ALLOCATE(zqsatl(klon,klev), zqsats(klon,klev))
+ ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev), fcontrN(klon,klev), fcontrP(klon,klev))
#ifdef CPP_StratAer
@@ -1117,4 +1172,14 @@
DEALLOCATE(zn2mout)
+! Supersaturation
+ DEALLOCATE(rneb_seri)
+ DEALLOCATE(d_rneb_dyn)
+ DEALLOCATE(qclr, qcld, qss, qvc)
+ DEALLOCATE(rnebclr, rnebss, gamma_ss)
+ DEALLOCATE(N1_ss, N2_ss)
+ DEALLOCATE(drneb_sub, drneb_con, drneb_tur, drneb_avi)
+ DEALLOCATE(zqsatl, zqsats)
+ DEALLOCATE(Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
+
#ifdef CPP_StratAer
! variables for strat. aerosol CK
Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 4059)
@@ -1917,4 +1917,52 @@
'albslw3', 'Surface albedo LW3', '-', (/ ('', i=1, 10) /))
+!--aviation & supersaturation
+ TYPE(ctrl_out), SAVE :: o_oclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'oclr', 'Clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_ocld = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'ocld', 'Cloudy sky total water', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_oss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'oss', 'ISSR total water', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_ovc = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'ovc', 'In-cloup vapor', 'kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_rnebclr = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'rnebclr', 'Clear sky fraction', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_rnebss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'rnebss', 'ISSR fraction', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_rnebseri = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'rnebseri', 'Cloud fraction', '-', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_gammass = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'gammass', 'Gamma supersaturation', '', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_N1_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'N1ss', 'N1', '', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_N2_ss = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'N2ss', 'N2', '', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_drnebsub = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'drnebsub', 'Cloud fraction change because of sublimation', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_drnebcon = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'drnebcon', 'Cloud fraction change because of condensation', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_drnebtur = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'drnebtur', 'Cloud fraction change because of turbulence', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_drnebavi = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'drnebavi', 'Cloud fraction change because of aviation', 's-1', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qsatl = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'qsatl', 'Saturation with respect to liquid water', '', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qsats = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'qsats', 'Saturation with respect to solid water', '', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_flight_m = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'flightm', 'Flown meters', 'm/s/mesh', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 1, 1, 1, 1, 10, 10, 11, 11, 11, 11/), &
+ 'flighth2o', 'H2O flight emission', 'kg H2O/s/mesh', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_Tcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
+ 'Tcontr', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
+ TYPE(ctrl_out), SAVE :: o_qcontr = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
+ 'qcontr', 'Specific humidity threshold for contrail formation','Pa', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_qcontr2 = ctrl_out((/ 1, 1, 1, 1, 11, 11, 11, 11, 11, 11/),&
+ 'qcontr2', 'Specific humidity threshold for contrail formation','kg/kg', (/ ('', i=1, 10) /))
+ TYPE(ctrl_out), SAVE :: o_fcontrN = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
+ 'fcontrN', 'Fraction with non-persistent contrail in clear-sky', '-', (/ ('', i=1,10)/))
+ TYPE(ctrl_out), SAVE :: o_fcontrP = ctrl_out((/ 2, 2, 2, 2, 2, 2, 11, 11, 11, 11/),&
+ 'fcontrP', 'Fraction with persistent contrail in ISSR', '-', (/ ('', i=1,10)/))
+
!!!!!!!!!!!!! Sorties niveaux standards de pression NMC
TYPE(ctrl_out), SAVE :: o_tnondef = ctrl_out((/ 11, 11, 11, 11, 11, 11, 5, 5, 5, 11/), &
Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 4059)
@@ -209,4 +209,9 @@
o_p_tropopause, o_z_tropopause, o_t_tropopause, &
o_col_O3_strato, o_col_O3_tropo, &
+!--aviation & supersaturation
+ o_oclr, o_ocld, o_oss, o_ovc, o_rnebss, o_rnebclr, o_rnebseri, o_gammass, &
+ o_N1_ss, o_N2_ss, o_qsatl, o_qsats, &
+ o_drnebsub, o_drnebcon, o_drnebtur, o_drnebavi, o_flight_m, o_flight_h2o, &
+ o_Tcontr, o_qcontr, o_qcontr2, o_fcontrN, o_fcontrP, &
!--interactive CO2
o_flx_co2_ocean, o_flx_co2_ocean_cor, &
@@ -230,4 +235,6 @@
o_vsed_aer, o_tau_strat_1020, o_ext_strat_1020, o_f_r_wet
#endif
+
+ USE ice_sursat_mod, ONLY: flight_m, flight_h2o
USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
@@ -280,4 +287,5 @@
cldq, flwp, fiwp, ue, ve, uq, vq, &
uwat, vwat, &
+ rneb_seri, d_rneb_dyn, &
plcl, plfc, wbeff, convoccur, upwd, dnwd, dnwd0, prw, prlw, prsw, &
s_pblh, s_pblt, s_lcl, s_therm, uwriteSTD, &
@@ -293,4 +301,8 @@
wdtrainA, wdtrainS, wdtrainM, n2, s2, proba_notrig, &
random_notrig, &
+ qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
+ N1_ss, N2_ss, zqsatl, zqsats, &
+ Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
+ drneb_sub, drneb_con, drneb_tur, drneb_avi, &
alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
alp_bl_stat, alp_bl_fluct_tke, slab_wfbils, &
@@ -1813,4 +1825,28 @@
ENDIF
+!--aviation & supersaturation
+ CALL histwrite_phy(o_oclr, qclr)
+ CALL histwrite_phy(o_ocld, qcld)
+ CALL histwrite_phy(o_oss, qss)
+ CALL histwrite_phy(o_ovc, qvc)
+ CALL histwrite_phy(o_rnebclr, rnebclr)
+ CALL histwrite_phy(o_rnebss, rnebss)
+ CALL histwrite_phy(o_rnebseri, rneb_seri)
+ CALL histwrite_phy(o_gammass, gamma_ss)
+ CALL histwrite_phy(o_N1_ss, N1_ss)
+ CALL histwrite_phy(o_N2_ss, N2_ss)
+ CALL histwrite_phy(o_drnebsub, drneb_sub)
+ CALL histwrite_phy(o_drnebcon, drneb_con)
+ CALL histwrite_phy(o_drnebtur, drneb_tur)
+ CALL histwrite_phy(o_drnebavi, drneb_avi)
+ CALL histwrite_phy(o_qsatl, zqsatl)
+ CALL histwrite_phy(o_qsats, zqsats)
+ CALL histwrite_phy(o_flight_m, flight_m)
+ CALL histwrite_phy(o_flight_h2o, flight_h2o)
+ CALL histwrite_phy(o_Tcontr, Tcontr)
+ CALL histwrite_phy(o_qcontr, qcontr)
+ CALL histwrite_phy(o_qcontr2, qcontr2)
+ CALL histwrite_phy(o_fcontrN, fcontrN)
+ CALL histwrite_phy(o_fcontrP, fcontrP)
IF (vars_defined) zx_tmp_fi3d = wo(:, :, 1) * dobson_u * 1e3 / zmasse / rmo3 * rmd
Index: LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90 (revision 4059)
@@ -90,4 +90,6 @@
REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:)
!$OMP THREADPRIVATE(clwcon,rnebcon)
+ REAL, ALLOCATABLE, SAVE :: rneb_ancien(:,:)
+!$OMP THREADPRIVATE(rneb_ancien)
REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:)
!$OMP THREADPRIVATE(qtc_cv,sigt_cv)
@@ -514,4 +516,5 @@
!!! Rom P <<<
ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
+ ALLOCATE(rneb_ancien(klon,klev))
ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev))
ALLOCATE(ratqs(klon,klev))
@@ -689,5 +692,5 @@
DEALLOCATE(zthe, zpic, zval)
DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
- DEALLOCATE(qs_ancien, ql_ancien)
+ DEALLOCATE(qs_ancien, ql_ancien, rneb_ancien)
DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien)
DEALLOCATE(qtc_cv,sigt_cv)
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 4058)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 4059)
@@ -73,4 +73,5 @@
USE tracinca_mod, ONLY: config_inca
USE tropopause_m, ONLY: dyn_tropopause
+ USE ice_sursat_mod, ONLY: flight_init, airplane
USE vampir
USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
@@ -127,7 +128,7 @@
! [Variables internes non sauvegardees de la physique]
! Variables locales pour effectuer les appels en serie
- t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
+ t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri,rneb_seri, &
! Dynamic tendencies (diagnostics)
- d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
+ d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn,d_rneb_dyn, &
d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
! Physic tendencies
@@ -444,7 +445,8 @@
INTEGER iliq ! indice de traceurs pour eau liquide
PARAMETER (iliq=2)
- !CR: on ajoute la phase glace
INTEGER isol ! indice de traceurs pour eau glace
PARAMETER (isol=3)
+ INTEGER irneb ! indice de traceurs pour fraction nuageuse LS (optional)
+ PARAMETER (irneb=4)
!
!
@@ -1283,4 +1285,6 @@
#endif
+ !!CALL flight_init
+
print*, '================================================='
!
@@ -1289,4 +1293,17 @@
WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
'(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
+ abort_message='see above'
+ CALL abort_physic(modname,abort_message,1)
+ ENDIF
+
+ IF ((iflag_ice_sursat.GT.0).AND.(iflag_ice_thermo.EQ.0)) THEN
+ WRITE (lunout, *) ' iflag_ice_sursat=1 requires iflag_ice_thermo=1 as well'
+ abort_message='see above'
+ CALL abort_physic(modname,abort_message,1)
+ ENDIF
+
+ IF ((iflag_ice_sursat.GT.0).AND.(nqo.NE.4)) THEN
+ WRITE (lunout, *) ' iflag_ice_sursat=1 requires 4 H2O tracers ', &
+ '(H2Ov, H2Ol, H2Oi, rnebi) but nqo=', nqo, '. Might as well stop here.'
abort_message='see above'
CALL abort_physic(modname,abort_message,1)
@@ -2217,4 +2234,7 @@
ELSE IF (nqo.eq.3) THEN
qs_seri(i,k) = qx(i,k,isol)
+ ELSE IF (nqo.eq.4) THEN
+ qs_seri(i,k) = qx(i,k,isol)
+ rneb_seri(i,k) = qx(i,k,irneb)
ENDIF
ENDDO
@@ -2292,4 +2312,6 @@
IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
! !! RomP <<<
+ !!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
+ d_rneb_dyn(:,:)=0.0
ELSE
d_u_dyn(:,:) = 0.0
@@ -2305,4 +2327,5 @@
IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
! !! RomP <<<
+ d_rneb_dyn(:,:)=0.0
ancien_ok = .TRUE.
ENDIF
@@ -3536,15 +3559,19 @@
IF (ok_new_lscp) THEN
- CALL lscp(phys_tstep,paprs,pplay, &
+ !--mise à jour de flight_m dans son module
+ CALL airplane(debut,pphis,pplay,paprs,t_seri)
+
+ CALL lscp(phys_tstep,missing_val,paprs,pplay, &
t_seri, q_seri,ptconv,ratqs, &
- d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, picefra, &
- rain_lsc, snow_lsc, &
+ d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneb_seri, &
+ cldliq, picefra, rain_lsc, snow_lsc, &
pfrac_impa, pfrac_nucl, pfrac_1nucl, &
frac_impa, frac_nucl, beta_prec_fisrt, &
prfl, psfl, rhcl, &
zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
- iflag_ice_thermo)
+ iflag_ice_thermo, iflag_ice_sursat)
ELSE
+
CALL fisrtilp(phys_tstep,paprs,pplay, &
t_seri, q_seri,ptconv,ratqs, &
@@ -3556,4 +3583,5 @@
zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
iflag_ice_thermo)
+
ENDIF
!
@@ -5060,6 +5088,10 @@
d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
!CR: on ajoute le contenu en glace
- IF (nqo.eq.3) THEN
+ IF (nqo.gt.3) THEN
d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
+ ENDIF
+ !--ice_sursat: nqo=4, on ajoute rneb
+ IF (nqo.eq.4) THEN
+ d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
ENDIF
ENDDO
@@ -5110,4 +5142,5 @@
ql_ancien(:,:) = ql_seri(:,:)
qs_ancien(:,:) = qs_seri(:,:)
+ rneb_ancien(:,:) = rneb_seri(:,:)
CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)