[4059] | 1 | MODULE ice_sursat_mod |
---|
| 2 | |
---|
| 3 | IMPLICIT NONE |
---|
| 4 | |
---|
[4225] | 5 | !--flight inventories |
---|
| 6 | ! |
---|
[4059] | 7 | REAL, SAVE, ALLOCATABLE :: flight_m(:,:) !--flown distance m s-1 per cell |
---|
| 8 | !$OMP THREADPRIVATE(flight_m) |
---|
| 9 | REAL, SAVE, ALLOCATABLE :: flight_h2o(:,:) !--emitted kg H2O s-1 per cell |
---|
| 10 | !$OMP THREADPRIVATE(flight_h2o) |
---|
| 11 | ! |
---|
[4225] | 12 | !--Fixed Parameters |
---|
[4059] | 13 | ! |
---|
| 14 | !--safety parameters for ERF function |
---|
| 15 | REAL, PARAMETER :: erf_lim = 5., eps = 1.e-10 |
---|
| 16 | ! |
---|
[4099] | 17 | !--Tuning parameters (and their default values) |
---|
[4059] | 18 | ! |
---|
[4099] | 19 | !--chi gère la répartition statistique de la longueur des frontières |
---|
| 20 | ! entre les zones nuages et ISSR/ciel clair sous-saturé. Gamme de valeur : |
---|
| 21 | ! chi > 1, je n'ai pas regardé de limite max (pour chi = 1, la longueur de |
---|
| 22 | ! la frontière entre ne nuage et l'ISSR est proportionnelle à la |
---|
| 23 | ! répartition ISSR/ciel clair sous-sat dans la maille, i.e. il n'y a pas |
---|
| 24 | ! de favorisation de la localisation de l'ISSR près de nuage. Pour chi = inf, |
---|
| 25 | ! le nuage n'est en contact qu'avec de l'ISSR, quelle que soit la taille |
---|
| 26 | ! de l'ISSR dans la maille.) |
---|
| 27 | ! |
---|
| 28 | !--l_turb est la longueur de mélange pour la turbulence. |
---|
| 29 | ! dans les tests, ça n'a jamais été modifié pour l'instant. |
---|
| 30 | ! |
---|
| 31 | !--tun_N est le paramètre qui contrôle l'importance relative de N_2 par rapport à N_1. |
---|
| 32 | ! La valeur est comprise entre 1 et 2 (tun_N = 1 => N_1 = N_2) |
---|
| 33 | ! |
---|
| 34 | !--tun_ratqs : paramètre qui modifie ratqs en fonction de la valeur de |
---|
| 35 | ! alpha_cld selon la formule ratqs_new = ratqs_old / ( 1 + tun_ratqs * |
---|
| 36 | ! alpha_cld ). Dans le rapport il est appelé beta. Il varie entre 0 et 5 |
---|
| 37 | ! (tun_ratqs = 0 => pas de modification de ratqs). |
---|
| 38 | ! |
---|
[4225] | 39 | !--gamma0 and Tgamma: define RHcrit limit above which heterogeneous freezing occurs as a function of T |
---|
| 40 | !--Karcher and Lohmann (2002) |
---|
| 41 | !--gamma = 2.583 - t / 207.83 |
---|
| 42 | !--Ren and MacKenzie (2005) reused by Kärcher |
---|
| 43 | !--gamma = 2.349 - t / 259.0 |
---|
| 44 | ! |
---|
| 45 | !--N_cld: number of clouds in cell (needs to be parametrized at some point) |
---|
| 46 | ! |
---|
| 47 | !--contrail cross section: typical value found in Freudenthaler et al, GRL, 22, 3501-3504, 1995 |
---|
[4059] | 48 | !--in m2, 1000x200 = 200 000 m2 after 15 min |
---|
[4225] | 49 | ! |
---|
| 50 | REAL, SAVE :: chi=1.1, l_turb=50.0, tun_N=1.3, tun_ratqs=3.0 |
---|
| 51 | REAL, SAVE :: gamma0=2.349, Tgamma=259.0, N_cld=100, contrail_cross_section=200000.0 |
---|
| 52 | !$OMP THREADPRIVATE(chi,l_turb,tun_N,tun_ratqs,gamma0,Tgamma,N_cld,contrail_cross_section) |
---|
[4059] | 53 | |
---|
| 54 | CONTAINS |
---|
| 55 | |
---|
| 56 | !******************************************************************* |
---|
[4099] | 57 | ! |
---|
| 58 | SUBROUTINE ice_sursat_init() |
---|
[4059] | 59 | |
---|
[4099] | 60 | USE print_control_mod, ONLY: lunout |
---|
| 61 | USE ioipsl_getin_p_mod, ONLY : getin_p |
---|
| 62 | |
---|
| 63 | IMPLICIT NONE |
---|
| 64 | |
---|
| 65 | CALL getin_p('flag_chi',chi) |
---|
| 66 | CALL getin_p('flag_l_turb',l_turb) |
---|
| 67 | CALL getin_p('flag_tun_N',tun_N) |
---|
| 68 | CALL getin_p('flag_tun_ratqs',tun_ratqs) |
---|
[4225] | 69 | CALL getin_p('gamma0',gamma0) |
---|
| 70 | CALL getin_p('Tgamma',Tgamma) |
---|
| 71 | CALL getin_p('N_cld',N_cld) |
---|
| 72 | CALL getin_p('contrail_cross_section',contrail_cross_section) |
---|
[4099] | 73 | |
---|
[4225] | 74 | WRITE(lunout,*) 'Parameters for ice_sursat param' |
---|
[4099] | 75 | WRITE(lunout,*) 'flag_chi = ', chi |
---|
| 76 | WRITE(lunout,*) 'flag_l_turb = ', l_turb |
---|
| 77 | WRITE(lunout,*) 'flag_tun_N = ', tun_N |
---|
| 78 | WRITE(lunout,*) 'flag_tun_ratqs = ', tun_ratqs |
---|
[4225] | 79 | WRITE(lunout,*) 'gamma0 = ', gamma0 |
---|
| 80 | WRITE(lunout,*) 'Tgamma = ', Tgamma |
---|
| 81 | WRITE(lunout,*) 'N_cld = ', N_cld |
---|
| 82 | WRITE(lunout,*) 'contrail_cross_section = ', contrail_cross_section |
---|
[4099] | 83 | |
---|
| 84 | END SUBROUTINE ice_sursat_init |
---|
| 85 | |
---|
| 86 | !******************************************************************* |
---|
| 87 | ! |
---|
[4059] | 88 | SUBROUTINE airplane(debut,pphis,pplay,paprs,t_seri) |
---|
| 89 | |
---|
| 90 | USE dimphy |
---|
| 91 | USE mod_grid_phy_lmdz, ONLY: klon_glo |
---|
| 92 | USE geometry_mod, ONLY: cell_area |
---|
| 93 | USE phys_cal_mod, ONLY : mth_cur |
---|
| 94 | USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root |
---|
| 95 | USE mod_phys_lmdz_omp_data, ONLY: is_omp_root |
---|
| 96 | USE mod_phys_lmdz_para, ONLY: scatter, bcast |
---|
| 97 | USE print_control_mod, ONLY: lunout |
---|
| 98 | |
---|
| 99 | IMPLICIT NONE |
---|
| 100 | |
---|
| 101 | INCLUDE "YOMCST.h" |
---|
[5084] | 102 | INCLUDE 'netcdf.inc' |
---|
[4059] | 103 | |
---|
| 104 | !-------------------------------------------------------- |
---|
| 105 | !--input variables |
---|
| 106 | !-------------------------------------------------------- |
---|
| 107 | LOGICAL, INTENT(IN) :: debut |
---|
| 108 | REAL, INTENT(IN) :: pphis(klon), pplay(klon,klev), paprs(klon,klev+1), t_seri(klon,klev) |
---|
| 109 | |
---|
| 110 | !-------------------------------------------------------- |
---|
| 111 | ! ... Local variables |
---|
| 112 | !-------------------------------------------------------- |
---|
| 113 | |
---|
| 114 | CHARACTER (LEN=20) :: modname='airplane_mod' |
---|
| 115 | INTEGER :: i, k, kori, iret, varid, error, ncida, klona |
---|
| 116 | INTEGER,SAVE :: nleva, ntimea |
---|
| 117 | !$OMP THREADPRIVATE(nleva,ntimea) |
---|
| 118 | REAL, ALLOCATABLE :: pkm_airpl_glo(:,:,:) !--km/s |
---|
| 119 | REAL, ALLOCATABLE :: ph2o_airpl_glo(:,:,:) !--molec H2O/cm3/s |
---|
| 120 | REAL, ALLOCATABLE, SAVE :: zmida(:), zinta(:) |
---|
| 121 | REAL, ALLOCATABLE, SAVE :: pkm_airpl(:,:,:) |
---|
| 122 | REAL, ALLOCATABLE, SAVE :: ph2o_airpl(:,:,:) |
---|
| 123 | !$OMP THREADPRIVATE(pkm_airpl,ph2o_airpl,zmida,zinta) |
---|
| 124 | REAL :: zalt(klon,klev+1) |
---|
| 125 | REAL :: zrho, zdz(klon,klev), zfrac |
---|
| 126 | |
---|
| 127 | ! |
---|
| 128 | IF (debut) THEN |
---|
| 129 | !-------------------------------------------------------------------------------- |
---|
| 130 | ! ... Open the file and read airplane emissions |
---|
| 131 | !-------------------------------------------------------------------------------- |
---|
| 132 | ! |
---|
| 133 | IF (is_mpi_root .AND. is_omp_root) THEN |
---|
| 134 | ! |
---|
| 135 | iret = nf_open('aircraft_phy.nc', 0, ncida) |
---|
| 136 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to open aircraft_phy.nc file',1) |
---|
| 137 | ! ... Get lengths |
---|
| 138 | iret = nf_inq_dimid(ncida, 'time', varid) |
---|
| 139 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimid in aircraft_phy.nc file',1) |
---|
| 140 | iret = nf_inq_dimlen(ncida, varid, ntimea) |
---|
| 141 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get time dimlen aircraft_phy.nc file',1) |
---|
| 142 | iret = nf_inq_dimid(ncida, 'vector', varid) |
---|
| 143 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimid aircraft_phy.nc file',1) |
---|
| 144 | iret = nf_inq_dimlen(ncida, varid, klona) |
---|
| 145 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get vector dimlen aircraft_phy.nc file',1) |
---|
| 146 | iret = nf_inq_dimid(ncida, 'lev', varid) |
---|
| 147 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) |
---|
| 148 | iret = nf_inq_dimlen(ncida, varid, nleva) |
---|
| 149 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimlen aircraft_phy.nc file',1) |
---|
| 150 | ! |
---|
| 151 | IF ( klona /= klon_glo ) THEN |
---|
[4099] | 152 | WRITE(lunout,*) 'klona & klon_glo =', klona, klon_glo |
---|
[4059] | 153 | CALL abort_physic(modname,'problem klon in aircraft_phy.nc file',1) |
---|
| 154 | ENDIF |
---|
| 155 | ! |
---|
| 156 | IF ( ntimea /= 12 ) THEN |
---|
[4099] | 157 | WRITE(lunout,*) 'ntimea=', ntimea |
---|
[4059] | 158 | CALL abort_physic(modname,'problem ntime<>12 in aircraft_phy.nc file',1) |
---|
| 159 | ENDIF |
---|
| 160 | ! |
---|
| 161 | ALLOCATE(zmida(nleva), STAT=error) |
---|
| 162 | IF (error /= 0) CALL abort_physic(modname,'problem to allocate zmida',1) |
---|
| 163 | ALLOCATE(pkm_airpl_glo(klona,nleva,ntimea), STAT=error) |
---|
| 164 | IF (error /= 0) CALL abort_physic(modname,'problem to allocate pkm_airpl_glo',1) |
---|
| 165 | ALLOCATE(ph2o_airpl_glo(klona,nleva,ntimea), STAT=error) |
---|
| 166 | IF (error /= 0) CALL abort_physic(modname,'problem to allocate ph2o_airpl_glo',1) |
---|
| 167 | ! |
---|
| 168 | iret = nf_inq_varid(ncida, 'lev', varid) |
---|
| 169 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get lev dimid aircraft_phy.nc file',1) |
---|
[5084] | 170 | iret = nf_get_var_double(ncida, varid, zmida) |
---|
[4059] | 171 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read zmida file',1) |
---|
| 172 | ! |
---|
| 173 | iret = nf_inq_varid(ncida, 'emi_co2_aircraft', varid) !--CO2 as a proxy for m flown - |
---|
| 174 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_distance dimid aircraft_phy.nc file',1) |
---|
[5084] | 175 | iret = nf_get_var_double(ncida, varid, pkm_airpl_glo) |
---|
[4059] | 176 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read pkm_airpl file',1) |
---|
| 177 | ! |
---|
| 178 | iret = nf_inq_varid(ncida, 'emi_h2o_aircraft', varid) |
---|
| 179 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to get emi_h2o_aircraft dimid aircraft_phy.nc file',1) |
---|
[5084] | 180 | iret = nf_get_var_double(ncida, varid, ph2o_airpl_glo) |
---|
[4059] | 181 | IF (iret /= NF_NOERR) CALL abort_physic(modname,'problem to read ph2o_airpl file',1) |
---|
| 182 | ! |
---|
| 183 | ENDIF !--is_mpi_root and is_omp_root |
---|
| 184 | ! |
---|
| 185 | CALL bcast(nleva) |
---|
| 186 | CALL bcast(ntimea) |
---|
| 187 | ! |
---|
| 188 | IF (.NOT.ALLOCATED(zmida)) ALLOCATE(zmida(nleva), STAT=error) |
---|
| 189 | IF (.NOT.ALLOCATED(zinta)) ALLOCATE(zinta(nleva+1), STAT=error) |
---|
| 190 | ! |
---|
| 191 | ALLOCATE(pkm_airpl(klon,nleva,ntimea)) |
---|
| 192 | ALLOCATE(ph2o_airpl(klon,nleva,ntimea)) |
---|
| 193 | ! |
---|
| 194 | ALLOCATE(flight_m(klon,klev)) |
---|
| 195 | ALLOCATE(flight_h2o(klon,klev)) |
---|
| 196 | ! |
---|
| 197 | CALL bcast(zmida) |
---|
| 198 | zinta(1)=0.0 !--surface |
---|
| 199 | DO k=2, nleva |
---|
| 200 | zinta(k) = (zmida(k-1)+zmida(k))/2.0*1000.0 !--conversion from km to m |
---|
| 201 | ENDDO |
---|
| 202 | zinta(nleva+1)=zinta(nleva)+(zmida(nleva)-zmida(nleva-1))*1000.0 !--extrapolation for last interface |
---|
| 203 | !print *,'zinta=', zinta |
---|
| 204 | ! |
---|
| 205 | CALL scatter(pkm_airpl_glo,pkm_airpl) |
---|
| 206 | CALL scatter(ph2o_airpl_glo,ph2o_airpl) |
---|
| 207 | ! |
---|
| 208 | !$OMP MASTER |
---|
| 209 | IF (is_mpi_root .AND. is_omp_root) THEN |
---|
| 210 | DEALLOCATE(pkm_airpl_glo) |
---|
| 211 | DEALLOCATE(ph2o_airpl_glo) |
---|
| 212 | ENDIF !--is_mpi_root |
---|
| 213 | !$OMP END MASTER |
---|
| 214 | |
---|
| 215 | ENDIF !--debut |
---|
| 216 | ! |
---|
| 217 | !--compute altitude of model level interfaces |
---|
| 218 | ! |
---|
| 219 | DO i = 1, klon |
---|
| 220 | zalt(i,1)=pphis(i)/RG !--in m |
---|
| 221 | ENDDO |
---|
| 222 | ! |
---|
| 223 | DO k=1, klev |
---|
| 224 | DO i = 1, klon |
---|
| 225 | zrho=pplay(i,k)/t_seri(i,k)/RD |
---|
| 226 | zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho/RG |
---|
| 227 | zalt(i,k+1)=zalt(i,k)+zdz(i,k) !--in m |
---|
| 228 | ENDDO |
---|
| 229 | ENDDO |
---|
| 230 | ! |
---|
| 231 | !--vertical reprojection |
---|
| 232 | ! |
---|
| 233 | flight_m(:,:)=0.0 |
---|
| 234 | flight_h2o(:,:)=0.0 |
---|
| 235 | ! |
---|
| 236 | DO k=1, klev |
---|
| 237 | DO kori=1, nleva |
---|
| 238 | DO i=1, klon |
---|
| 239 | !--fraction of layer kori included in layer k |
---|
| 240 | zfrac=max(0.0,min(zalt(i,k+1),zinta(kori+1))-max(zalt(i,k),zinta(kori)))/(zinta(kori+1)-zinta(kori)) |
---|
| 241 | !--reproject |
---|
| 242 | flight_m(i,k)=flight_m(i,k) + pkm_airpl(i,kori,mth_cur)*zfrac |
---|
| 243 | !--reproject |
---|
| 244 | flight_h2o(i,k)=flight_h2o(i,k) + ph2o_airpl(i,kori,mth_cur)*zfrac |
---|
| 245 | ENDDO |
---|
| 246 | ENDDO |
---|
| 247 | ENDDO |
---|
| 248 | ! |
---|
| 249 | DO k=1, klev |
---|
| 250 | DO i=1, klon |
---|
| 251 | !--molec.cm-3.s-1 / (molec/mol) * kg CO2/mol * m2 * m * cm3/m3 / (kg CO2/m) => m s-1 per cell |
---|
| 252 | flight_m(i,k)=flight_m(i,k)/RNAVO*44.e-3*cell_area(i)*zdz(i,k)*1.e6/16.37e-3 |
---|
| 253 | flight_m(i,k)=flight_m(i,k)*100.0 !--x100 to augment signal to noise |
---|
| 254 | !--molec.cm-3.s-1 / (molec/mol) * kg H2O/mol * m2 * m * cm3/m3 => kg H2O s-1 per cell |
---|
| 255 | flight_h2o(i,k)=flight_h2o(i,k)/RNAVO*18.e-3*cell_area(i)*zdz(i,k)*1.e6 |
---|
| 256 | ENDDO |
---|
| 257 | ENDDO |
---|
| 258 | ! |
---|
| 259 | END SUBROUTINE airplane |
---|
| 260 | |
---|
| 261 | !******************************************************************** |
---|
[4099] | 262 | ! simple routine to initialise flight_m and test a flight corridor |
---|
| 263 | !--Olivier Boucher - 2021 |
---|
| 264 | ! |
---|
[4059] | 265 | SUBROUTINE flight_init() |
---|
| 266 | USE dimphy |
---|
| 267 | USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg |
---|
| 268 | IMPLICIT NONE |
---|
| 269 | INTEGER :: i |
---|
| 270 | |
---|
| 271 | ALLOCATE(flight_m(klon,klev)) |
---|
| 272 | ALLOCATE(flight_h2o(klon,klev)) |
---|
| 273 | ! |
---|
| 274 | flight_m(:,:) = 0.0 !--initialisation |
---|
| 275 | flight_h2o(:,:) = 0.0 !--initialisation |
---|
| 276 | ! |
---|
| 277 | DO i=1, klon |
---|
[5084] | 278 | IF (latitude_deg(i).GE.42.0.AND.latitude_deg(i).LE.48.0) THEN |
---|
[4059] | 279 | flight_m(i,38) = 50000.0 !--5000 m of flight/second in grid cell x 10 scaling |
---|
| 280 | ENDIF |
---|
| 281 | ENDDO |
---|
| 282 | |
---|
| 283 | RETURN |
---|
| 284 | END SUBROUTINE flight_init |
---|
| 285 | |
---|
| 286 | !******************************************************************* |
---|
| 287 | !--Routine to deal with ice supersaturation |
---|
| 288 | !--Determines the respective fractions of unsaturated clear sky, ice supersaturated clear sky and cloudy sky |
---|
| 289 | !--Diagnoses regions prone for non-persistent and persistent contrail formation |
---|
| 290 | ! |
---|
| 291 | !--Audran Borella - 2021 |
---|
| 292 | ! |
---|
| 293 | SUBROUTINE ice_sursat(pplay, dpaprs, dtime, i, k, t, q, gamma_ss, & |
---|
| 294 | qsat, t_actuel, rneb_seri, ratqs, rneb, qincld, & |
---|
| 295 | rnebss, qss, Tcontr, qcontr, qcontr2, fcontrN, fcontrP) |
---|
| 296 | ! |
---|
| 297 | USE dimphy |
---|
| 298 | USE print_control_mod, ONLY: prt_level, lunout |
---|
| 299 | USE phys_state_var_mod, ONLY: pbl_tke, t_ancien |
---|
| 300 | USE phys_local_var_mod, ONLY: N1_ss, N2_ss |
---|
| 301 | USE phys_local_var_mod, ONLY: drneb_sub, drneb_con, drneb_tur, drneb_avi |
---|
| 302 | !! USE phys_local_var_mod, ONLY: Tcontr, qcontr, fcontrN, fcontrP |
---|
| 303 | USE indice_sol_mod, ONLY: is_ave |
---|
| 304 | USE geometry_mod, ONLY: cell_area |
---|
| 305 | ! |
---|
| 306 | IMPLICIT NONE |
---|
| 307 | INCLUDE "YOMCST.h" |
---|
| 308 | INCLUDE "YOETHF.h" |
---|
| 309 | INCLUDE "FCTTRE.h" |
---|
[4062] | 310 | INCLUDE "clesphys.h" |
---|
| 311 | |
---|
[4059] | 312 | ! |
---|
| 313 | ! Input |
---|
| 314 | ! Beware: this routine works on a gridpoint! |
---|
| 315 | ! |
---|
| 316 | REAL, INTENT(IN) :: pplay ! layer pressure (Pa) |
---|
| 317 | REAL, INTENT(IN) :: dpaprs ! layer delta pressure (Pa) |
---|
| 318 | REAL, INTENT(IN) :: dtime ! intervalle du temps (s) |
---|
| 319 | REAL, INTENT(IN) :: t ! température advectée (K) |
---|
| 320 | REAL, INTENT(IN) :: qsat ! vapeur de saturation |
---|
| 321 | REAL, INTENT(IN) :: t_actuel ! temperature actuelle de la maille (K) |
---|
| 322 | REAL, INTENT(IN) :: rneb_seri ! fraction nuageuse en memoire |
---|
| 323 | INTEGER, INTENT(IN) :: i, k |
---|
| 324 | ! |
---|
| 325 | ! Input/output |
---|
| 326 | ! |
---|
| 327 | REAL, INTENT(INOUT) :: q ! vapeur de la maille (=zq) |
---|
| 328 | REAL, INTENT(INOUT) :: ratqs ! determine la largeur de distribution de vapeur |
---|
| 329 | REAL, INTENT(INOUT) :: Tcontr, qcontr, qcontr2, fcontrN, fcontrP |
---|
| 330 | ! |
---|
| 331 | ! Output |
---|
| 332 | ! |
---|
| 333 | REAL, INTENT(OUT) :: gamma_ss ! |
---|
| 334 | REAL, INTENT(OUT) :: rneb ! cloud fraction |
---|
| 335 | REAL, INTENT(OUT) :: qincld ! in-cloud total water |
---|
| 336 | REAL, INTENT(OUT) :: rnebss ! ISSR fraction |
---|
| 337 | REAL, INTENT(OUT) :: qss ! in-ISSR total water |
---|
| 338 | ! |
---|
| 339 | ! Local |
---|
| 340 | ! |
---|
| 341 | REAL PI |
---|
| 342 | PARAMETER (PI=4.*ATAN(1.)) |
---|
| 343 | REAL rnebclr, gamma_prec |
---|
| 344 | REAL qclr, qvc, qcld, qi |
---|
| 345 | REAL zrho, zdz, zrhodz |
---|
| 346 | REAL pdf_N, pdf_N1, pdf_N2 |
---|
| 347 | REAL pdf_a, pdf_b |
---|
| 348 | REAL pdf_e1, pdf_e2, pdf_k |
---|
| 349 | REAL drnebss, drnebclr, dqss, dqclr, sum_rneb_rnebss, dqss_avi |
---|
| 350 | REAL V_cell !--volume of the cell |
---|
| 351 | REAL M_cell !--dry mass of the cell |
---|
| 352 | REAL tke, sig, L_tur, b_tur, q_eq |
---|
| 353 | REAL V_env, V_cld, V_ss, V_clr |
---|
| 354 | REAL zcor |
---|
| 355 | ! |
---|
| 356 | !--more local variables for diagnostics |
---|
| 357 | !--imported from YOMCST.h |
---|
| 358 | !--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1) |
---|
| 359 | !--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air |
---|
| 360 | !--values from Schumann, Meteorol Zeitschrift, 1996 |
---|
| 361 | !--EiH2O = 1.25 / 2.24 / 8.94 kg H2O / kg fuel for kerosene / methane / dihydrogen |
---|
| 362 | !--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen |
---|
| 363 | REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1) |
---|
| 364 | REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) |
---|
| 365 | REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft |
---|
| 366 | !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute |
---|
| 367 | !--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010. |
---|
| 368 | REAL :: Gcontr |
---|
| 369 | !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) |
---|
| 370 | !--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1) |
---|
| 371 | REAL :: qsatliqcontr |
---|
| 372 | |
---|
| 373 | ! Initialisations |
---|
| 374 | zrho = pplay / t / RD !--dry density kg m-3 |
---|
| 375 | zrhodz = dpaprs / RG !--dry air mass kg m-2 |
---|
| 376 | zdz = zrhodz / zrho !--cell thickness m |
---|
| 377 | V_cell = zdz * cell_area(i) !--cell volume m3 |
---|
| 378 | M_cell = zrhodz * cell_area(i) !--cell dry air mass kg |
---|
| 379 | ! |
---|
| 380 | ! Recuperation de la memoire sur la couverture nuageuse |
---|
| 381 | rneb = rneb_seri |
---|
| 382 | ! |
---|
| 383 | ! Ajout des émissions de H2O dues à l'aviation |
---|
| 384 | ! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q |
---|
| 385 | ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) |
---|
| 386 | ! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) |
---|
| 387 | ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) |
---|
| 388 | ! flight_h2O is in kg H2O / s / cell |
---|
| 389 | ! |
---|
[4062] | 390 | IF (ok_plane_h2o) THEN |
---|
[4059] | 391 | q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) |
---|
| 392 | ENDIF |
---|
| 393 | ! |
---|
| 394 | !--Estimating gamma |
---|
[4225] | 395 | gamma_ss = MAX(1.0, gamma0 - t_actuel/Tgamma) |
---|
| 396 | !gamma_prec = MAX(1.0, gamma0 - t_ancien(i,k)/Tgamma) !--formulation initiale d Audran |
---|
| 397 | gamma_prec = MAX(1.0, gamma0 - t/Tgamma) !--autre formulation possible basée sur le T du pas de temps |
---|
[4059] | 398 | ! |
---|
| 399 | ! Initialisation de qvc : q_sat du pas de temps precedent |
---|
| 400 | !qvc = R2ES*FOEEW(t_ancien(i,k),1.)/pplay !--formulation initiale d Audran |
---|
| 401 | qvc = R2ES*FOEEW(t,1.)/pplay !--autre formulation possible basée sur le T du pas de temps |
---|
| 402 | qvc = min(0.5,qvc) |
---|
| 403 | zcor = 1./(1.-RETV*qvc) |
---|
| 404 | qvc = qvc*zcor |
---|
| 405 | ! |
---|
| 406 | ! Modification de ratqs selon formule proposee : ksi_new = ksi_old/(1+beta*alpha_cld) |
---|
| 407 | ratqs = ratqs / (tun_ratqs*rneb_seri + 1.) |
---|
| 408 | ! |
---|
| 409 | ! Calcul de N |
---|
| 410 | pdf_k = -sqrt(log(1.+ratqs**2.)) |
---|
| 411 | pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) |
---|
| 412 | pdf_b = pdf_k/(2.*sqrt(2.)) |
---|
| 413 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 414 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 415 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 416 | pdf_N = max(0.,sign(rneb,pdf_e1)) |
---|
| 417 | ELSE |
---|
| 418 | pdf_e1 = erf(pdf_e1) |
---|
| 419 | pdf_e1 = 0.5*(1.+pdf_e1) |
---|
| 420 | pdf_N = max(0.,rneb/pdf_e1) |
---|
| 421 | ENDIF |
---|
| 422 | ! |
---|
| 423 | ! On calcule ensuite N1 et N2. Il y a deux cas : N > 1 et N <= 1 |
---|
| 424 | ! Cas 1 : N > 1. N'arrive en theorie jamais, c'est une barriere |
---|
| 425 | ! On perd la memoire sur la temperature (sur qvc) pour garder |
---|
| 426 | ! celle sur alpha_cld |
---|
[5084] | 427 | IF (pdf_N.GT.1.) THEN |
---|
[4059] | 428 | ! On inverse alpha_cld = int_qvc^infty P(q) dq |
---|
| 429 | ! pour determiner qvc = f(alpha_cld) |
---|
| 430 | ! On approxime en serie entiere erf-1(x) |
---|
| 431 | qvc = 2.*rneb-1. |
---|
[4260] | 432 | qvc = qvc + PI/12.*qvc**3 + 7.*PI**2/480.*qvc**5 & |
---|
| 433 | + 127.*PI**3/40320.*qvc**7 + 4369.*PI**4/5806080.*qvc**9 & |
---|
| 434 | + 34807.*PI**5/182476800.*qvc**11 |
---|
[4059] | 435 | qvc = sqrt(PI)/2.*qvc |
---|
| 436 | qvc = (qvc-pdf_b)*pdf_k*sqrt(2.) |
---|
| 437 | qvc = q*exp(qvc) |
---|
| 438 | |
---|
| 439 | ! On met a jour rneb avec la nouvelle valeur de qvc |
---|
| 440 | ! La maj est necessaire a cause de la serie entiere approximative |
---|
| 441 | pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) |
---|
| 442 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 443 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 444 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 445 | ELSE |
---|
| 446 | pdf_e1 = erf(pdf_e1) |
---|
| 447 | ENDIF |
---|
| 448 | pdf_e1 = 0.5*(1.+pdf_e1) |
---|
| 449 | rneb = pdf_e1 |
---|
| 450 | |
---|
| 451 | ! Si N > 1, N1 et N2 = 1 |
---|
| 452 | pdf_N1 = 1. |
---|
| 453 | pdf_N2 = 1. |
---|
| 454 | |
---|
| 455 | ! Cas 2 : N <= 1 |
---|
| 456 | ELSE |
---|
| 457 | ! D'abord on calcule N2 avec le tuning |
---|
| 458 | pdf_N2 = min(1.,pdf_N*tun_N) |
---|
| 459 | |
---|
| 460 | ! Puis N1 pour assurer la conservation de alpha_cld |
---|
| 461 | pdf_a = log(qvc*gamma_prec/q)/(pdf_k*sqrt(2.)) |
---|
| 462 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 463 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 464 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 465 | ELSE |
---|
| 466 | pdf_e2 = erf(pdf_e2) |
---|
| 467 | ENDIF |
---|
| 468 | pdf_e2 = 0.5*(1.+pdf_e2) ! integrale sous P pour q > gamma qsat |
---|
| 469 | |
---|
[5084] | 470 | IF (abs(pdf_e1-pdf_e2).LT.eps) THEN |
---|
[4059] | 471 | pdf_N1 = pdf_N2 |
---|
| 472 | ELSE |
---|
| 473 | pdf_N1 = (rneb-pdf_N2*pdf_e2)/(pdf_e1-pdf_e2) |
---|
| 474 | ENDIF |
---|
| 475 | |
---|
| 476 | ! Barriere qui traite le cas gamma_prec = 1. |
---|
[5084] | 477 | IF (pdf_N1.LE.0.) THEN |
---|
[4059] | 478 | pdf_N1 = 0. |
---|
[5084] | 479 | IF (pdf_e2.GT.eps) THEN |
---|
[4059] | 480 | pdf_N2 = rneb/pdf_e2 |
---|
| 481 | ELSE |
---|
| 482 | pdf_N2 = 0. |
---|
| 483 | ENDIF |
---|
| 484 | ENDIF |
---|
| 485 | ENDIF |
---|
| 486 | |
---|
| 487 | ! Physique 1 |
---|
| 488 | ! Sublimation |
---|
[5084] | 489 | IF (qvc.LT.qsat) THEN |
---|
[4059] | 490 | pdf_a = log(qvc/q)/(pdf_k*sqrt(2.)) |
---|
| 491 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 492 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 493 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 494 | ELSE |
---|
| 495 | pdf_e1 = erf(pdf_e1) |
---|
| 496 | ENDIF |
---|
| 497 | |
---|
| 498 | pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 499 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 500 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 501 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 502 | ELSE |
---|
| 503 | pdf_e2 = erf(pdf_e2) |
---|
| 504 | ENDIF |
---|
| 505 | |
---|
| 506 | pdf_e1 = 0.5*pdf_N1*(pdf_e1-pdf_e2) |
---|
| 507 | |
---|
| 508 | ! Calcul et ajout de la tendance |
---|
| 509 | drneb_sub(i,k) = - pdf_e1/dtime !--unit [s-1] |
---|
| 510 | rneb = rneb + drneb_sub(i,k)*dtime |
---|
| 511 | ELSE |
---|
| 512 | drneb_sub(i,k) = 0. |
---|
| 513 | ENDIF |
---|
| 514 | |
---|
| 515 | ! NOTE : verifier si ca marche bien pour gamma proche de 1. |
---|
| 516 | |
---|
| 517 | ! Condensation |
---|
[5084] | 518 | IF (gamma_ss*qsat.LT.gamma_prec*qvc) THEN |
---|
[4059] | 519 | |
---|
| 520 | pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 521 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 522 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 523 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 524 | ELSE |
---|
| 525 | pdf_e1 = erf(pdf_e1) |
---|
| 526 | ENDIF |
---|
| 527 | |
---|
| 528 | pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.)) |
---|
| 529 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 530 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 531 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 532 | ELSE |
---|
| 533 | pdf_e2 = erf(pdf_e2) |
---|
| 534 | ENDIF |
---|
| 535 | |
---|
| 536 | pdf_e1 = 0.5*(1.-pdf_N1)*(pdf_e1-pdf_e2) |
---|
| 537 | pdf_e2 = 0.5*(1.-pdf_N2)*(1.+pdf_e2) |
---|
| 538 | |
---|
| 539 | ! Calcul et ajout de la tendance |
---|
| 540 | drneb_con(i,k) = (pdf_e1 + pdf_e2)/dtime !--unit [s-1] |
---|
| 541 | rneb = rneb + drneb_con(i,k)*dtime |
---|
| 542 | |
---|
| 543 | ELSE |
---|
| 544 | |
---|
| 545 | pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 546 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 547 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 548 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 549 | ELSE |
---|
| 550 | pdf_e1 = erf(pdf_e1) |
---|
| 551 | ENDIF |
---|
| 552 | pdf_e1 = 0.5*(1.-pdf_N2)*(1.+pdf_e1) |
---|
| 553 | |
---|
| 554 | ! Calcul et ajout de la tendance |
---|
| 555 | drneb_con(i,k) = pdf_e1/dtime !--unit [s-1] |
---|
| 556 | rneb = rneb + drneb_con(i,k)*dtime |
---|
| 557 | |
---|
| 558 | ENDIF |
---|
| 559 | |
---|
| 560 | ! Calcul des grandeurs diagnostiques |
---|
| 561 | ! Determination des grandeurs ciel clair |
---|
| 562 | pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 563 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 564 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 565 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 566 | ELSE |
---|
| 567 | pdf_e1 = erf(pdf_e1) |
---|
| 568 | ENDIF |
---|
| 569 | pdf_e1 = 0.5*(1.-pdf_e1) |
---|
| 570 | |
---|
| 571 | pdf_e2 = pdf_a-pdf_b |
---|
[5084] | 572 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 573 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 574 | ELSE |
---|
| 575 | pdf_e2 = erf(pdf_e2) |
---|
| 576 | ENDIF |
---|
| 577 | pdf_e2 = 0.5*q*(1.-pdf_e2) |
---|
| 578 | |
---|
| 579 | rnebclr = pdf_e1 |
---|
| 580 | qclr = pdf_e2 |
---|
| 581 | |
---|
| 582 | ! Determination de q_cld |
---|
| 583 | ! Partie 1 |
---|
| 584 | pdf_a = log(max(qsat,qvc)/q)/(pdf_k*sqrt(2.)) |
---|
| 585 | pdf_e1 = pdf_a-pdf_b |
---|
[5084] | 586 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 587 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 588 | ELSE |
---|
| 589 | pdf_e1 = erf(pdf_e1) |
---|
| 590 | ENDIF |
---|
| 591 | |
---|
| 592 | pdf_a = log(min(gamma_ss*qsat,gamma_prec*qvc)/q)/(pdf_k*sqrt(2.)) |
---|
| 593 | pdf_e2 = pdf_a-pdf_b |
---|
[5084] | 594 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 595 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 596 | ELSE |
---|
| 597 | pdf_e2 = erf(pdf_e2) |
---|
| 598 | ENDIF |
---|
| 599 | |
---|
| 600 | pdf_e1 = 0.5*q*pdf_N1*(pdf_e1-pdf_e2) |
---|
| 601 | |
---|
| 602 | qcld = pdf_e1 |
---|
| 603 | |
---|
| 604 | ! Partie 2 (sous condition) |
---|
[5084] | 605 | IF (gamma_ss*qsat.GT.gamma_prec*qvc) THEN |
---|
[4059] | 606 | pdf_a = log(gamma_ss*qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 607 | pdf_e1 = pdf_a-pdf_b |
---|
[5084] | 608 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 609 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 610 | ELSE |
---|
| 611 | pdf_e1 = erf(pdf_e1) |
---|
| 612 | ENDIF |
---|
| 613 | |
---|
| 614 | pdf_e2 = 0.5*q*pdf_N2*(pdf_e2-pdf_e1) |
---|
| 615 | |
---|
| 616 | qcld = qcld + pdf_e2 |
---|
| 617 | |
---|
| 618 | pdf_e2 = pdf_e1 |
---|
| 619 | ENDIF |
---|
| 620 | |
---|
| 621 | ! Partie 3 |
---|
| 622 | pdf_e2 = 0.5*q*(1.+pdf_e2) |
---|
| 623 | |
---|
| 624 | qcld = qcld + pdf_e2 |
---|
| 625 | |
---|
| 626 | ! Fin du calcul de q_cld |
---|
| 627 | |
---|
| 628 | ! Determination des grandeurs ISSR via les equations de conservation |
---|
| 629 | rneb=MIN(rneb, 1. - rnebclr - eps) !--ajout OB - barrière |
---|
| 630 | rnebss = MAX(0.0, 1. - rnebclr - rneb) !--ajout OB |
---|
| 631 | qss = MAX(0.0, q - qclr - qcld) !--ajout OB |
---|
| 632 | |
---|
| 633 | ! Physique 2 : Turbulence |
---|
[5084] | 634 | IF (rneb.GT.eps.AND.rneb.LT.1.-eps) THEN ! rneb != 0 and != 1 |
---|
[4059] | 635 | ! |
---|
| 636 | tke = pbl_tke(i,k,is_ave) |
---|
| 637 | ! A MODIFIER formule a revoir |
---|
| 638 | L_tur = min(l_turb, sqrt(tke)*dtime) |
---|
| 639 | |
---|
| 640 | ! On fait pour l'instant l'hypothese a = 3b. V = 4/3 pi a b**2 = alpha_cld |
---|
| 641 | ! donc b = alpha_cld/4pi **1/3. |
---|
| 642 | b_tur = (rneb*V_cell/4./PI/N_cld)**(1./3.) |
---|
| 643 | ! On verifie que la longeur de melange n'est pas trop grande |
---|
[5084] | 644 | IF (L_tur.GT.b_tur) THEN |
---|
[4059] | 645 | L_tur = b_tur |
---|
| 646 | ENDIF |
---|
| 647 | |
---|
| 648 | V_env = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. + 3.*b_tur*(L_tur**2.)) |
---|
| 649 | V_cld = N_cld*4.*PI*(3.*(b_tur**2.)*L_tur + L_tur**3. - 3.*b_tur*(L_tur**2.)) |
---|
| 650 | V_cld = abs(V_cld) |
---|
| 651 | |
---|
| 652 | ! Repartition statistique des zones de contact nuage-ISSR et nuage-ciel clair |
---|
| 653 | sig = rnebss/(chi*rnebclr+rnebss) |
---|
| 654 | V_ss = MIN(sig*V_env,rnebss*V_cell) |
---|
| 655 | V_clr = MIN((1.-sig)*V_env,rnebclr*V_cell) |
---|
| 656 | V_cld = MIN(V_cld,rneb*V_cell) |
---|
| 657 | V_env = V_ss + V_clr |
---|
| 658 | |
---|
| 659 | ! ISSR => rneb |
---|
| 660 | drnebss = -1.*V_ss/V_cell |
---|
| 661 | dqss = drnebss*qss/MAX(eps,rnebss) |
---|
| 662 | |
---|
| 663 | ! Clear sky <=> rneb |
---|
| 664 | q_eq = V_env*qclr/MAX(eps,rnebclr) + V_cld*qcld/MAX(eps,rneb) |
---|
| 665 | q_eq = q_eq/(V_env + V_cld) |
---|
| 666 | |
---|
[5084] | 667 | IF (q_eq.GT.qsat) THEN |
---|
[4059] | 668 | drnebclr = - V_clr/V_cell |
---|
| 669 | dqclr = drnebclr*qclr/MAX(eps,rnebclr) |
---|
| 670 | ELSE |
---|
| 671 | drnebclr = V_cld/V_cell |
---|
| 672 | dqclr = drnebclr*qcld/MAX(eps,rneb) |
---|
| 673 | ENDIF |
---|
| 674 | |
---|
| 675 | ! Maj des variables avec les tendances |
---|
| 676 | rnebclr = MAX(0.0,rnebclr + drnebclr) !--OB ajout d'un max avec eps (il faudrait modified drnebclr pour le diag) |
---|
| 677 | qclr = MAX(0.0, qclr + dqclr) !--OB ajout d'un max avec 0 |
---|
| 678 | |
---|
| 679 | rneb = rneb - drnebclr - drnebss |
---|
| 680 | qcld = qcld - dqclr - dqss |
---|
| 681 | |
---|
| 682 | rnebss = MAX(0.0,rnebss + drnebss) !--OB ajout d'un max avec eps (il faudrait modifier drnebss pour le diag) |
---|
| 683 | qss = MAX(0.0, qss + dqss) !--OB ajout d'un max avec 0 |
---|
| 684 | |
---|
| 685 | ! Tendances pour le diagnostic |
---|
| 686 | drneb_tur(i,k) = (drnebclr + drnebss)/dtime !--unit [s-1] |
---|
| 687 | |
---|
| 688 | ENDIF ! rneb |
---|
| 689 | |
---|
| 690 | !--add a source of cirrus from aviation contrails |
---|
[4062] | 691 | IF (ok_plane_contrail) THEN |
---|
[4059] | 692 | drneb_avi(i,k) = rnebss*flight_m(i,k)*contrail_cross_section/V_cell !--tendency rneb due to aviation [s-1] |
---|
| 693 | drneb_avi(i,k) = MIN(drneb_avi(i,k), rnebss/dtime) !--majoration |
---|
| 694 | dqss_avi = qss*drneb_avi(i,k)/MAX(eps,rnebss) !--tendency q aviation [kg kg-1 s-1] |
---|
| 695 | rneb = rneb + drneb_avi(i,k)*dtime !--add tendency to rneb |
---|
| 696 | qcld = qcld + dqss_avi*dtime !--add tendency to qcld |
---|
| 697 | rnebss = rnebss - drneb_avi(i,k)*dtime !--add tendency to rnebss |
---|
| 698 | qss = qss - dqss_avi*dtime !--add tendency to qss |
---|
| 699 | ELSE |
---|
| 700 | drneb_avi(i,k)=0.0 |
---|
| 701 | ENDIF |
---|
| 702 | |
---|
| 703 | ! Barrieres |
---|
| 704 | ! ISSR trop petite |
---|
[5084] | 705 | IF (rnebss.LT.eps) THEN |
---|
[4059] | 706 | rneb = MIN(rneb + rnebss,1.0-eps) !--ajout OB barriere |
---|
| 707 | qcld = qcld + qss |
---|
| 708 | rnebss = 0. |
---|
| 709 | qss = 0. |
---|
| 710 | ENDIF |
---|
| 711 | |
---|
| 712 | ! le nuage est trop petit |
---|
[5084] | 713 | IF (rneb.LT.eps) THEN |
---|
[4059] | 714 | ! s'il y a une ISSR on met tout dans l'ISSR, sinon dans le |
---|
| 715 | ! clear sky |
---|
[5084] | 716 | IF (rnebss.LT.eps) THEN |
---|
[4059] | 717 | rnebclr = 1. |
---|
| 718 | rnebss = 0. !--ajout OB |
---|
| 719 | qclr = q |
---|
| 720 | ELSE |
---|
| 721 | rnebss = MIN(rnebss + rneb,1.0-eps) !--ajout OB barriere |
---|
| 722 | qss = qss + qcld |
---|
| 723 | ENDIF |
---|
| 724 | rneb = 0. |
---|
| 725 | qcld = 0. |
---|
| 726 | qincld = qsat * gamma_ss |
---|
| 727 | ELSE |
---|
| 728 | qincld = qcld / rneb |
---|
| 729 | ENDIF |
---|
| 730 | |
---|
| 731 | !--OB ajout borne superieure |
---|
| 732 | sum_rneb_rnebss=rneb+rnebss |
---|
| 733 | rneb=rneb*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss) |
---|
| 734 | rnebss=rnebss*MIN(1.-eps,sum_rneb_rnebss)/MAX(eps,sum_rneb_rnebss) |
---|
| 735 | |
---|
| 736 | ! On ecrit dans la memoire |
---|
| 737 | N1_ss(i,k) = pdf_N1 |
---|
| 738 | N2_ss(i,k) = pdf_N2 |
---|
| 739 | |
---|
| 740 | !--Diagnostics only used from last iteration |
---|
| 741 | !--test |
---|
| 742 | !!Tcontr(i,k)=200. |
---|
| 743 | !!fcontrN(i,k)=1.0 |
---|
| 744 | !!fcontrP(i,k)=0.5 |
---|
| 745 | ! |
---|
| 746 | !--slope of dilution line in exhaust |
---|
| 747 | !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1) |
---|
| 748 | Gcontr = EiH2O * RCPD * pplay / (eps_w*Qheat*(1.-eta)) !--Pa K-1 |
---|
| 749 | !--critical T_LM below which no liquid contrail can form in exhaust |
---|
| 750 | !Tcontr(i,k) = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K |
---|
[5084] | 751 | IF (Gcontr .GT. 0.1) THEN |
---|
[4059] | 752 | ! |
---|
[4074] | 753 | Tcontr = 226.69+9.43*log(Gcontr-0.053)+0.72*(log(Gcontr-0.053))**2 !--K |
---|
| 754 | !print *,'Tcontr=',iter,i,k,eps_w,pplay,Gcontr,Tcontr(i,k) |
---|
| 755 | !--Psat with index 0 in FOEEW to get saturation wrt liquid water corresponding to Tcontr |
---|
| 756 | !qsatliqcontr = RESTT*FOEEW(Tcontr(i,k),0.) !--Pa |
---|
| 757 | qsatliqcontr = RESTT*FOEEW(Tcontr,0.) !--Pa |
---|
| 758 | !--Critical water vapour above which there is contrail formation for ambiant temperature |
---|
| 759 | !qcontr(i,k) = Gcontr*(t-Tcontr(i,k)) + qsatliqcontr !--Pa |
---|
| 760 | qcontr = Gcontr*(t-Tcontr) + qsatliqcontr !--Pa |
---|
| 761 | !--Conversion of qcontr in specific humidity - method 1 |
---|
| 762 | !qcontr(i,k) = RD/RV*qcontr(i,k)/pplay !--so as to return to something similar to R2ES*FOEEW/pplay |
---|
| 763 | qcontr2 = RD/RV*qcontr/pplay !--so as to return to something similar to R2ES*FOEEW/pplay |
---|
| 764 | !qcontr(i,k) = min(0.5,qcontr(i,k)) !--and then we apply the same correction term as for qsat |
---|
| 765 | qcontr2 = min(0.5,qcontr2) !--and then we apply the same correction term as for qsat |
---|
| 766 | !zcor = 1./(1.-RETV*qcontr(i,k)) !--for consistency with qsat but is it correct at all? |
---|
| 767 | zcor = 1./(1.-RETV*qcontr2) !--for consistency with qsat but is it correct at all as p is dry? |
---|
| 768 | !zcor = 1./(1.+qcontr2) !--for consistency with qsat |
---|
| 769 | !qcontr(i,k) = qcontr(i,k)*zcor |
---|
| 770 | qcontr2 = qcontr2*zcor |
---|
| 771 | qcontr2=MAX(1.e-10,qcontr2) !--eliminate negative values due to extrapolation on dilution curve |
---|
| 772 | !--Conversion of qcontr in specific humidity - method 2 |
---|
| 773 | !qcontr(i,k) = eps_w*qcontr(i,k) / (pplay+eps_w*qcontr(i,k)) |
---|
| 774 | !qcontr=MAX(1.E-10,qcontr) |
---|
| 775 | !qcontr2 = eps_w*qcontr / (pplay+eps_w*qcontr) |
---|
| 776 | ! |
---|
[5084] | 777 | IF (t .LT. Tcontr) THEN !--contrail formation is possible |
---|
[4074] | 778 | ! |
---|
| 779 | !--compute fractions of persistent (P) and non-persistent(N) contrail potential regions |
---|
| 780 | !!IF (qcontr(i,k).GE.qsat) THEN |
---|
[5084] | 781 | IF (qcontr2.GE.qsat) THEN |
---|
[4059] | 782 | !--none of the unsaturated clear sky is prone for contrail formation |
---|
| 783 | !!fcontrN(i,k) = 0.0 |
---|
| 784 | fcontrN = 0.0 |
---|
| 785 | ! |
---|
| 786 | !--integral of P(q) from qsat to qcontr in ISSR |
---|
| 787 | pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 788 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 789 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 790 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 791 | ELSE |
---|
| 792 | pdf_e1 = erf(pdf_e1) |
---|
| 793 | ENDIF |
---|
| 794 | ! |
---|
| 795 | !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.)) |
---|
| 796 | pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.)) |
---|
| 797 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 798 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 799 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 800 | ELSE |
---|
| 801 | pdf_e2 = erf(pdf_e2) |
---|
| 802 | ENDIF |
---|
| 803 | ! |
---|
| 804 | !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
---|
| 805 | fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
---|
| 806 | ! |
---|
| 807 | pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 808 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 809 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 810 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 811 | ELSE |
---|
| 812 | pdf_e1 = erf(pdf_e1) |
---|
| 813 | ENDIF |
---|
| 814 | ! |
---|
| 815 | !!pdf_a = log(MIN(qcontr(i,k),qvc)/q)/(pdf_k*sqrt(2.)) |
---|
| 816 | pdf_a = log(MIN(qcontr2,qvc)/q)/(pdf_k*sqrt(2.)) |
---|
| 817 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 818 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 819 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 820 | ELSE |
---|
| 821 | pdf_e2 = erf(pdf_e2) |
---|
| 822 | ENDIF |
---|
| 823 | ! |
---|
| 824 | !!fcontrP(i,k) = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
---|
| 825 | fcontrP = MAX(0., 0.5*(pdf_e1-pdf_e2)) |
---|
| 826 | ! |
---|
| 827 | pdf_a = log(MAX(qsat,qvc)/q)/(pdf_k*sqrt(2.)) |
---|
| 828 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 829 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 830 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 831 | ELSE |
---|
| 832 | pdf_e1 = erf(pdf_e1) |
---|
| 833 | ENDIF |
---|
| 834 | ! |
---|
| 835 | !!pdf_a = log(MIN(qcontr(i,k),MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.)) |
---|
| 836 | pdf_a = log(MIN(qcontr2,MIN(gamma_prec*qvc,gamma_ss*qsat))/q)/(pdf_k*sqrt(2.)) |
---|
| 837 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 838 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 839 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 840 | ELSE |
---|
| 841 | pdf_e2 = erf(pdf_e2) |
---|
| 842 | ENDIF |
---|
| 843 | ! |
---|
| 844 | !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2)) |
---|
| 845 | fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N1)*(pdf_e1-pdf_e2)) |
---|
| 846 | ! |
---|
| 847 | pdf_a = log(gamma_prec*qvc/q)/(pdf_k*sqrt(2.)) |
---|
| 848 | pdf_e1 = pdf_a+pdf_b |
---|
[5084] | 849 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 850 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 851 | ELSE |
---|
| 852 | pdf_e1 = erf(pdf_e1) |
---|
| 853 | ENDIF |
---|
| 854 | ! |
---|
| 855 | !!pdf_a = log(MIN(qcontr(i,k),gamma_ss*qsat)/q)/(pdf_k*sqrt(2.)) |
---|
| 856 | pdf_a = log(MIN(qcontr2,gamma_ss*qsat)/q)/(pdf_k*sqrt(2.)) |
---|
| 857 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 858 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 859 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 860 | ELSE |
---|
| 861 | pdf_e2 = erf(pdf_e2) |
---|
| 862 | ENDIF |
---|
| 863 | ! |
---|
| 864 | !!fcontrP(i,k) = fcontrP(i,k) + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2)) |
---|
| 865 | fcontrP = fcontrP + MAX(0., 0.5*(1-pdf_N2)*(pdf_e1-pdf_e2)) |
---|
| 866 | ! |
---|
[4074] | 867 | ELSE !--qcontr LT qsat |
---|
[4059] | 868 | ! |
---|
| 869 | !--all of ISSR is prone for contrail formation |
---|
| 870 | !!fcontrP(i,k) = rnebss |
---|
| 871 | fcontrP = rnebss |
---|
| 872 | ! |
---|
| 873 | !--integral of zq from qcontr to qsat in unsaturated clear-sky region |
---|
| 874 | !!pdf_a = log(qcontr(i,k)/q)/(pdf_k*sqrt(2.)) |
---|
| 875 | pdf_a = log(qcontr2/q)/(pdf_k*sqrt(2.)) |
---|
| 876 | pdf_e1 = pdf_a+pdf_b !--normalement pdf_b est deja defini |
---|
[5084] | 877 | IF (abs(pdf_e1).GE.erf_lim) THEN |
---|
[4059] | 878 | pdf_e1 = sign(1.,pdf_e1) |
---|
| 879 | ELSE |
---|
| 880 | pdf_e1 = erf(pdf_e1) |
---|
| 881 | ENDIF |
---|
| 882 | ! |
---|
| 883 | pdf_a = log(qsat/q)/(pdf_k*sqrt(2.)) |
---|
| 884 | pdf_e2 = pdf_a+pdf_b |
---|
[5084] | 885 | IF (abs(pdf_e2).GE.erf_lim) THEN |
---|
[4059] | 886 | pdf_e2 = sign(1.,pdf_e2) |
---|
| 887 | ELSE |
---|
| 888 | pdf_e2 = erf(pdf_e2) |
---|
| 889 | ENDIF |
---|
| 890 | ! |
---|
| 891 | !!fcontrN(i,k) = 0.5*(pdf_e1-pdf_e2) |
---|
| 892 | fcontrN = 0.5*(pdf_e1-pdf_e2) |
---|
| 893 | !!fcontrN=2.0 |
---|
| 894 | ! |
---|
[4074] | 895 | ENDIF |
---|
| 896 | ! |
---|
| 897 | ENDIF !-- t < Tcontr |
---|
[4059] | 898 | ! |
---|
[4074] | 899 | ENDIF !-- Gcontr > 0.1 |
---|
[4059] | 900 | |
---|
| 901 | RETURN |
---|
| 902 | END SUBROUTINE ice_sursat |
---|
| 903 | ! |
---|
| 904 | !******************************************************************* |
---|
| 905 | ! |
---|
| 906 | END MODULE ice_sursat_mod |
---|