Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 (revision 5487) +++ LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90 (revision 5488) @@ -678,158 +678,4 @@ ! 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 Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90 (revision 5487) +++ LMDZ6/branches/contrails/libf/phylmd/lmdz_call_cloud_optics_prop.f90 (revision 5488) @@ -10,5 +10,5 @@ reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, & reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, & - icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon) + icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail) ! Interface between the LMDZ physics monitor and the cloud properties calculation routines @@ -47,4 +47,6 @@ LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection LOGICAL, INTENT(IN) :: ok_newmicro, ok_aie + + REAL, INTENT(IN) :: rcontrail(klon, klev) ! ratio of contrail to total cloud fraction, used only if ok_plane_contrail=y [-] ! inout: @@ -116,5 +118,5 @@ reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, & reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, & - icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon) + icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail) ELSE CALL nuage (paprs, pplay, & Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90 (revision 5487) +++ LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90 (revision 5488) @@ -9,5 +9,5 @@ reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, & reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, & - icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon) + icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, rcontrail) USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc @@ -28,4 +28,5 @@ USE lmdz_cloud_optics_prop_ini , ONLY : ok_icefra_lscp, rei_max, rei_min USE lmdz_cloud_optics_prop_ini , ONLY : zepsec, novlp, iflag_ice_thermo, ok_new_lscp + USE lmdz_cloud_optics_prop_ini , ONLY : ok_plane_contrail, re_ice_crystals_contrails @@ -69,4 +70,6 @@ LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection + + REAL, INTENT(IN) :: rcontrail(klon, klev) ! ratio of contrails to total cloud fraction, used only if ok_plane_contrail=y [-] ! inout: @@ -333,4 +336,10 @@ IF (zfiwp_var==0. .OR. rei<=0.) rei = 1. + IF ( ok_plane_contrail ) THEN + !--If contrails are activated, rei is a weighted average between the natural + !--rei and the contrails rei, with the weights being the fraction of natural + !--vs contrail cirrus in the gridbox + rei = rei * ( 1. - rcontrail(i,k) ) + re_ice_crystals_contrails * rcontrail(i,k) + ENDIF pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + & zfiwp_var*(3.448E-03+2.431/rei) @@ -442,4 +451,10 @@ IF (zflwp_var==0.) rel = 1. IF (zfiwp_var==0. .OR. rei<=0.) rei = 1. + IF ( ok_plane_contrail ) THEN + !--If contrails are activated, rei is a weighted average between the natural + !--rei and the contrails rei, with the weights being the fraction of natural + !--vs contrail cirrus in the gridbox + rei = rei * ( 1. - rcontrail(i,k) ) + re_ice_crystals_contrails * rcontrail(i,k) + ENDIF pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ & rei) Index: LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop_ini.f90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop_ini.f90 (revision 5487) +++ LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop_ini.f90 (revision 5488) @@ -11,4 +11,5 @@ LOGICAL, PROTECTED :: ok_cdnc LOGICAL, PROTECTED :: ok_icefra_lscp, ok_new_lscp + LOGICAL, PROTECTED :: ok_plane_contrail REAL, PROTECTED :: bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula REAL, ALLOCATABLE :: latitude_deg(:) @@ -22,4 +23,5 @@ REAL, PROTECTED :: rei_max, rei_min REAL, PROTECTED :: zepsec + REAL, PROTECTED :: re_ice_crystals_contrails=7.5E-6 ! [m] effective radius of ice crystals in contrails REAL, PARAMETER :: thres_tau=0.3, thres_neb=0.001 REAL, PARAMETER :: prmhc=440.*100. ! Pressure between medium and high level cloud in Pa @@ -39,4 +41,5 @@ !$OMP THREADPRIVATE(rad_chau1, rad_chau2, rad_froid, rei_max, rei_min) !$OMP THREADPRIVATE(zepsec) +!$OMP THREADPRIVATE(re_ice_crystals_contrails) @@ -46,5 +49,6 @@ & ok_cdnc_in, bl95_b0_in, & & bl95_b1_in, latitude_deg_in, rpi_in, rg_in, rd_in, zepsec_in, novlp_in, & - & iflag_ice_thermo_in, ok_new_lscp_in) + & iflag_ice_thermo_in, ok_new_lscp_in, & + & ok_plane_contrail_in) USE ioipsl_getin_p_mod, ONLY : getin_p @@ -56,4 +60,5 @@ INTEGER, INTENT(IN) :: novlp_in, iflag_ice_thermo_in LOGICAL, INTENT(IN) :: ok_cdnc_in, ok_new_lscp_in + LOGICAL, INTENT(IN) :: ok_plane_contrail_in REAL, INTENT(IN) :: bl95_b0_in, bl95_b1_in REAL, INTENT(IN) :: latitude_deg_in(klon) @@ -77,4 +82,5 @@ iflag_ice_thermo = iflag_ice_thermo_in ok_new_lscp = ok_new_lscp_in + ok_plane_contrail = ok_plane_contrail_in call getin_p('cdnc_min',cdnc_min) @@ -98,4 +104,5 @@ rei_max = 61.29 CALL getin_p('rei_max',rei_max) + CALL getin_p('re_ice_crystals_contrails', re_ice_crystals_contrails) Index: LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90 (revision 5487) +++ LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90 (revision 5488) @@ -1877,5 +1877,6 @@ & ok_cdnc, bl95_b0, & & bl95_b1, latitude_deg, rpi, rg, rd, & - & zepsec, novlp, iflag_ice_thermo, ok_new_lscp) + & zepsec, novlp, iflag_ice_thermo, ok_new_lscp, & + & ok_plane_contrail) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -4465,5 +4466,5 @@ ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, & reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, & - zfice, dNovrN, ptconv, rnebcon, clwcon) + zfice, dNovrN, ptconv, rnebcon, clwcon, rcont_seri) !