Changeset 5589 for LMDZ6


Ignore:
Timestamp:
Mar 26, 2025, 6:05:40 PM (2 months ago)
Author:
aborella
Message:

Multiple changes:

  • added new radiative diagnostics for contrails
  • added ok_rad_contrail option to allow for a double call of RRTM (w/ and w/o contrails)
  • transformed resuspension of snow into ice sedimentation (poprecip)
  • some modifications in poprecip in line with the ones from EV
  • cleaned sublimation of ice clouds in lmdz_lscp_condensation, option ok_ice_supersat
  • aviation emissions can now be read with IOIPSL (in lon/lat mode)
Location:
LMDZ6/branches/contrails
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml

    r5575 r5589  
    912912        <field id="qsati"      long_name="Saturation with respect to ice"    unit="kg/kg" />
    913913
    914         <field id="rcontseri"  long_name="Linear contrail fraction to cloud fraction ratio" unit="-" />
    915         <field id="drcontdyn"  long_name="Dynamics linear contrail fraction to cloud fraction ratio tendency" unit="s-1" />
    916         <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
    917         <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
     914            <field id="rcontseri"  long_name="Linear contrail fraction to cloud fraction ratio" unit="-" />
     915            <field id="drcontdyn"  long_name="Dynamics linear contrail fraction to cloud fraction ratio tendency" unit="s-1" />
     916            <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
     917            <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
    918918        <field id="potcontfraP"  long_name="Fraction with pontential persistent contrail"    unit="-" />
    919919        <field id="potcontfraNP" long_name="Fraction with potential non-persistent contrail"    unit="-" />
    … …  
    925925        <field id="flightdist" long_name="Aviation flown distance concentration"    unit="m/s/m3" />
    926926        <field id="flighth2o"  long_name="Aviation emitted H2O concentration"    unit="kg H2O/s/m3" />
    927         <field id="dqavi"      long_name="Water vapor emissions from aviation tendency"    unit="kg/kg/s" />
     927            <field id="dqavi"      long_name="Water vapor emissions from aviation tendency"    unit="kg/kg/s" />
     928            <field id="cldfra_nocont" long_name="Cloud fraction w/o contrails"    unit="-" />
     929            <field id="cldtau_nocont" long_name="Cloud optical thickness w/o contrails"    unit="1" />
     930            <field id="cldemi_nocont" long_name="Cloud optical emissivity w/o contrails"    unit="1" />
     931            <field id="cldh_nocont"   long_name="High-level cloudiness w/o contrails"    unit="-" />
     932            <field id="contcov"       long_name="Total contrails cover"    unit="-" />
     933        <field id="iwp_nocont"    long_name="Cloud ice water path w/o contrails"    unit="kg/m2" />
     934        <field id="iwc_nocont"    long_name="Cloud ice water content seen by radiation w/o contrails"    unit="kg/m3" />
     935            <field id="ref_ice_nocont" long_name="Effective radius of ice crystals w/o contrails"    unit="microns" />
     936        <field id="tops_nocont"   long_name="Solar rad. at TOA w/o contrails"    unit="W/m2" />
     937        <field id="topl_nocont"   long_name="IR rad. at TOA w/o contrails"    unit="W/m2" />
     938        <field id="sols_nocont"   long_name="Solar rad. at surf. w/o contrails"    unit="W/m2" />
     939        <field id="soll_nocont"   long_name="IR rad. at surface w/o contrails"    unit="W/m2" />
     940        <field id="nettop_nocont" long_name="Net dn radiatif flux at TOA w/o contrails"    unit="W/m2" />
    928941
    929942        <field id="fluxt"     long_name="flux h"     unit="W/m2" />

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

    r5575 r5589  
    4545          , iflag_rrtm, ok_strato, ok_hines, ok_qch4                    &
    4646          , iflag_ice_thermo, ok_ice_supersat, ok_no_issr_strato       &
    47           , ok_plane_h2o, ok_plane_contrail                            &
     47          , ok_plane_h2o, ok_plane_contrail, ok_rad_contrail           &
    4848          , ok_gwd_rando, NSW, iflag_albedo                            &
    4949          , ok_chlorophyll, ok_conserv_q, adjust_tropopause             &
    … …  
    142142  INTEGER :: ip_ebil_phy, iflag_rrtm, iflag_ice_thermo, NSW, iflag_albedo
    143143  LOGICAL :: ok_ice_supersat, ok_no_issr_strato
    144   LOGICAL :: ok_plane_h2o, ok_plane_contrail
     144  LOGICAL :: ok_plane_h2o, ok_plane_contrail, ok_rad_contrail
    145145  LOGICAL :: ok_chlorophyll
    146146  LOGICAL :: ok_strato
    … …  
    203203  !$OMP      , iflag_rrtm, ok_strato, ok_hines, ok_qch4                    &
    204204  !$OMP      , iflag_ice_thermo, ok_ice_supersat, ok_no_issr_strato       &
    205   !$OMP      , ok_plane_h2o, ok_plane_contrail                            &
     205  !$OMP      , ok_plane_h2o, ok_plane_contrail, ok_rad_contrail           &
    206206  !$OMP      , ok_gwd_rando, NSW, iflag_albedo                            &
    207207  !$OMP      , ok_chlorophyll, ok_conserv_q, adjust_tropopause             &

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

    r5575 r5589  
    173173    INTEGER,SAVE :: iflag_ice_thermo_omp
    174174    LOGICAL,SAVE :: ok_ice_supersat_omp, ok_no_issr_strato_omp
    175     LOGICAL,SAVE :: ok_plane_h2o_omp, ok_plane_contrail_omp
     175    LOGICAL,SAVE :: ok_plane_h2o_omp, ok_plane_contrail_omp, ok_rad_contrail_omp
    176176    REAL,SAVE :: rad_froid_omp, rad_chau1_omp, rad_chau2_omp
    177177    INTEGER,SAVE :: iflag_sic_omp, iflag_inertie_omp
    … …  
    13781378    ok_plane_contrail_omp = .FALSE.
    13791379    CALL getin('ok_plane_contrail',ok_plane_contrail_omp)
     1380
     1381    !Config Key  = ok_rad_contrail
     1382    !Config Desc = double call to radiative routine w/ and w/o contrails
     1383    !Config Def  = .FALSE.
     1384    !Config Help =
     1385    !
     1386    ok_rad_contrail_omp = .FALSE.
     1387    CALL getin('ok_rad_contrail',ok_rad_contrail_omp)
    13801388
    13811389    !
    … …  
    23572365    ok_plane_h2o = ok_plane_h2o_omp
    23582366    ok_plane_contrail = ok_plane_contrail_omp
     2367    ok_rad_contrail = ok_rad_contrail_omp
    23592368    top_height = top_height_omp
    23602369    overlap = overlap_omp
    … …  
    27832792    WRITE(lunout,*) ' ok_plane_h2o = ',ok_plane_h2o
    27842793    WRITE(lunout,*) ' ok_plane_contrail = ',ok_plane_contrail
     2794    WRITE(lunout,*) ' ok_rad_contrail = ',ok_rad_contrail
    27852795    WRITE(lunout,*) ' overlap = ',overlap
    27862796    WRITE(lunout,*) ' cdmmax = ',cdmmax

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

    r5579 r5589  
    635635    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, nbp_lev_=>nbp_lev, &
    636636                                 klon_glo, grid2Dto1D_glo, grid_type, unstructured
     637    USE iophy, ONLY : io_lon, io_lat
    637638    USE lmdz_xios
    638639    USE print_control_mod, ONLY: lunout
     640    USE readaerosol_mod, ONLY: check_err
    639641    USE lmdz_lscp_ini, ONLY: EI_H2O_aviation
    640642    IMPLICIT NONE
    … …  
    648650    INTEGER :: ierr
    649651
     652    ! FOR REGULAR LON LAT
     653    CHARACTER(len=30)     :: fname
     654    INTEGER               :: nid, dimid, varid
     655    INTEGER               :: imth, i, j, k
     656    REAL                  :: npole, spole
     657    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D
     658    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_h2o_glo2D
     659    REAL, ALLOCATABLE, DIMENSION(:)       :: vartmp_lev
     660    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: vartmp
     661    REAL, ALLOCATABLE, DIMENSION(:)       :: lon_src              ! longitudes in file
     662    REAL, ALLOCATABLE, DIMENSION(:)       :: lat_src, lat_src_inv ! latitudes in file
     663    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
     664    INTEGER                               :: nbp_lon, nbp_lat
     665
     666
     667IF (grid_type==unstructured) THEN
     668
    650669    ! Get number of vertical levels and level values
    651670    IF (is_omp_master) CALL xios_get_axis_attr( "aviation_lev", n_glo=nleva )
    … …  
    654673    ! Allocation of arrays
    655674    ALLOCATE(flight_dist_read(klon, nleva,1), STAT=ierr)
    656     IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist',1)
     675    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1)
    657676    ALLOCATE(flight_h2o_read(klon, nleva,1), STAT=ierr)
    658     IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_h2o',1)
     677    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_h2o',1)
    659678    ALLOCATE(aviation_lev(nleva), STAT=ierr)
    660     IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate aviation_lev',1)
     679    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1)
    661680
    662681    ! Read the data from the file
    … …  
    675694
    676695    ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev)
    677     ! (klon_mpi,klon) = (200,50) avec 80 MPI, 4 OMP, nbp40
    678696    CALL scatter_omp(flight_dist_mpi, flight_dist_read)
    679697    CALL scatter_omp(flight_h2o_mpi, flight_h2o_read)
    680698    CALL bcast_omp(aviation_lev)
     699
     700ELSE
     701
     702    nbp_lon=nbp_lon_
     703    nbp_lat=nbp_lat_
     704   
     705    IF (is_mpi_root) THEN
     706      ALLOCATE(lon_src(nbp_lon))
     707      ALLOCATE(lat_src(nbp_lat))
     708      ALLOCATE(lat_src_inv(nbp_lat))
     709    ELSE
     710      ALLOCATE(flight_dist_glo2D(0,0,0,0))
     711      ALLOCATE(flight_h2o_glo2D(0,0,0,0))
     712    ENDIF
     713
     714    IF (is_mpi_root .AND. is_omp_root) THEN
     715
     716! 1) Open file
     717!****************************************************************************************
     718      fname = 'aviation.nc'
     719      CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, nid), "pb open "//TRIM(fname) )
     720
     721
     722! Test for equal longitudes and latitudes in file and model
     723!****************************************************************************************
     724      ! Read and test longitudes
     725      CALL check_err( nf90_inq_varid(nid, 'longitude', varid), "pb inq lon" )
     726      CALL check_err( nf90_get_var(nid, varid, lon_src(:)), "pb get lon" )
     727     
     728      IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
     729         WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
     730         WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
     731         WRITE(lunout,*) 'longitudes in model :', io_lon
     732       
     733         CALL abort_physic('lmdz_aviation', 'longitudes are not the same in file and model',1)
     734      END IF
     735
     736      ! Read and test latitudes
     737      CALL check_err( nf90_inq_varid(nid, 'latitude', varid),"pb inq lat" )
     738      CALL check_err( nf90_get_var(nid, varid, lat_src(:)),"pb get lat" )
     739
     740      ! Invert source latitudes
     741      DO j = 1, nbp_lat
     742         lat_src_inv(j) = lat_src(nbp_lat+1-j)
     743      END DO
     744
     745      IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
     746         ! Latitudes are the same
     747         invert_lat=.FALSE.
     748      ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
     749         ! Inverted source latitudes correspond to model latitudes
     750         WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
     751         invert_lat=.TRUE.
     752      ELSE
     753         WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
     754         WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
     755         WRITE(lunout,*) 'latitudes in model :', io_lat
     756         CALL abort_physic('lmdz_aviation', 'latitudes do not correspond between file and model',1)
     757      END IF
     758
     759       
     760! 2) Find vertical dimension nleva
     761!****************************************************************************************
     762       ierr = nf90_inq_dimid(nid, 'pressure_Pa', dimid)
     763       CALL check_err( nf90_inquire_dimension(nid, dimid, len = nleva), "pb inq dim for PRESNIVS or lev" )
     764       
     765     ! Allocate variables depending on the number of vertical levels
     766       ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_h2o_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
     767       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1)
     768
     769       ALLOCATE(aviation_lev(nleva), stat=ierr)
     770       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 2',1)
     771
     772! 3) Read all variables from file
     773!****************************************************************************************
     774
     775       ! Get variable id
     776       CALL check_err( nf90_inq_varid(nid, 'seg_length_m', varid),"pb inq var seg_length_m" )
     777       ! Get the variable
     778       CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m"  )
     779
     780!       ! Get variable id
     781!       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
     782!       ! Get the variable
     783!       CALL check_err( nf90_get_var(nid, varid, flight_h2o_glo2D),"pb get var fuel_burn"  )
     784       flight_h2o_glo2D(:,:,:,:) = 0.
     785
     786       ! Get variable id
     787       CALL check_err( nf90_inq_varid(nid, "pressure_Pa", varid),"pb inq var pressure_Pa" )
     788       ! Get the variable
     789       CALL check_err( nf90_get_var(nid, varid, aviation_lev),"pb get var pressure_Pa" )
     790         
     791
     792! 4) Close file 
     793!****************************************************************************************
     794       CALL check_err( nf90_close(nid), "pb in close" )
     795     
     796
     797! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
     798!****************************************************************************************
     799! Test if vertical levels have to be inversed
     800
     801       IF (aviation_lev(1) < aviation_lev(nleva)) THEN
     802         
     803          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva), vartmp_lev(nleva))
     804
     805          ! Inverse vertical levels for flight_dist_glo2D
     806          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1)
     807          DO k=1, nleva
     808             DO j=1, nbp_lat
     809                DO i=1, nbp_lon
     810                   flight_dist_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
     811                END DO
     812             END DO
     813          END DO
     814
     815          ! Inverse vertical levels for flight_h2o_glo2D
     816          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1)
     817          DO k=1, nleva
     818             DO j=1, nbp_lat
     819                DO i=1, nbp_lon
     820                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
     821                END DO
     822             END DO
     823          END DO
     824           
     825          ALLOCATE(vartmp_lev(nleva))
     826          ! Inverte vertical axes for aviation_lev
     827          vartmp_lev(:) = aviation_lev(:)
     828          DO k=1, nleva
     829             aviation_lev(k) = vartmp_lev(nleva+1-k)
     830          END DO
     831
     832          DEALLOCATE(vartmp, vartmp_lev)
     833
     834       ELSE
     835          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
     836       END IF
     837
     838
     839!     - Invert latitudes if necessary
     840       IF (invert_lat) THEN
     841
     842          ALLOCATE(vartmp(nbp_lon, nbp_lat, nleva))
     843
     844          ! Invert latitudes for the variable
     845          vartmp(:,:,:) = flight_dist_glo2D(:,:,:,1) ! use varmth temporarly
     846          DO k=1,nleva
     847             DO j=1,nbp_lat
     848                DO i=1,nbp_lon
     849                   flight_dist_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
     850                END DO
     851             END DO
     852          END DO
     853
     854          ! Invert latitudes for the variable
     855          vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1) ! use varmth temporarly
     856          DO k=1,nleva
     857             DO j=1,nbp_lat
     858                DO i=1,nbp_lon
     859                   flight_h2o_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
     860                END DO
     861             END DO
     862          END DO
     863
     864          DEALLOCATE(vartmp)
     865       END IF ! invert_lat
     866       
     867       ! Do zonal mead at poles and distribut at whole first and last latitude
     868       DO k=1, nleva
     869          npole=0.  ! North pole, j=1
     870          spole=0.  ! South pole, j=nbp_lat       
     871          DO i=1,nbp_lon
     872             npole = npole + flight_dist_glo2D(i,1,k,1)
     873             spole = spole + flight_dist_glo2D(i,nbp_lat,k,1)
     874          END DO
     875          npole = npole/REAL(nbp_lon)
     876          spole = spole/REAL(nbp_lon)
     877          flight_dist_glo2D(:,1,    k,1) = npole
     878          flight_dist_glo2D(:,nbp_lat,k,1) = spole
     879       END DO
     880
     881       ! Do zonal mead at poles and distribut at whole first and last latitude
     882       DO k=1, nleva
     883          npole=0.  ! North pole, j=1
     884          spole=0.  ! South pole, j=nbp_lat       
     885          DO i=1,nbp_lon
     886             npole = npole + flight_h2o_glo2D(i,1,k,1)
     887             spole = spole + flight_h2o_glo2D(i,nbp_lat,k,1)
     888          END DO
     889          npole = npole/REAL(nbp_lon)
     890          spole = spole/REAL(nbp_lon)
     891          flight_h2o_glo2D(:,1,    k,1) = npole
     892          flight_h2o_glo2D(:,nbp_lat,k,1) = spole
     893       END DO
     894     
     895       ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_h2o_mpi(klon_glo, nleva, 1), stat=ierr)
     896       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1)
     897     
     898       ! Transform from 2D to 1D field
     899       CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi)
     900       CALL grid2Dto1D_glo(flight_h2o_glo2D, flight_h2o_mpi)
     901   
     902    ELSE
     903        ALLOCATE(flight_dist_mpi(0,0,0), flight_h2o_mpi(0,0,0))
     904    END IF ! is_mpi_root .AND. is_omp_root
     905
     906!$OMP BARRIER
     907 
     908! 6) Distribute to all processes
     909!    Scatter global field(klon_glo) to local process domain(klon)
     910!    and distribute nleva to all processes
     911!****************************************************************************************
     912
     913    ! Distribute nleva
     914    CALL bcast(nleva)
     915
     916    ! Allocate and distribute pt_ap and pt_b
     917    IF (.NOT. ALLOCATED(aviation_lev)) THEN  ! if pt_ap is allocated also pt_b is allocated
     918       ALLOCATE(aviation_lev(nleva), stat=ierr)
     919       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 4',1)
     920    END IF
     921    CALL bcast(aviation_lev)
     922
     923    ! Allocate space for output pointer variable at local process
     924    ALLOCATE(flight_dist_read(klon, nleva, 1), flight_h2o_read(klon, nleva, 1), stat=ierr)
     925    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1)
     926
     927    ! Scatter global field to local domain at local process
     928    CALL scatter(flight_dist_mpi, flight_dist_read)
     929    CALL scatter(flight_h2o_mpi, flight_h2o_read)
     930
     931ENDIF
    681932
    682933END SUBROUTINE read_aviation_emissions

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

    r5488 r5589  
    1010    reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
    1111    reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
    12     icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail)
     12    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
     13    !--AB contrails
     14    rcontrail, pclc_nocont, pcltau_nocont, pclemi_nocont, pch_nocont, &
     15    pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
    1316
    1417  ! Interface between the LMDZ physics monitor and the cloud properties calculation routines
    … …  
    4750  LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection
    4851  LOGICAL, INTENT(IN) :: ok_newmicro, ok_aie
    49 
    50   REAL, INTENT(IN) :: rcontrail(klon, klev) ! ratio of contrail to total cloud fraction, used only if ok_plane_contrail=y [-]
    5152
    5253  ! inout:
    … …  
    9293  REAL, INTENT(INOUT) :: icefrac_optics(klon, klev)! ice fraction in clouds seen by radiation [-]
    9394
     95  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
     96  REAL, INTENT(IN) :: rcontrail(klon, klev)      ! ratio of contrails to total cloud fraction [-]
     97  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
     98  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
     99  REAL, INTENT(OUT) :: xfiwp_nocont(klon)        ! ice water path (seen by radiation) without contrails [kg/m2]
     100  REAL, INTENT(OUT) :: xfiwc_nocont(klon, klev)  ! ice water content seen by radiation without contrails [kg/kg]
     101  REAL, INTENT(OUT) :: pclc_nocont(klon, klev)   ! cloud fraction for radiation without contrails [-]
     102  REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [m]
     103  REAL, INTENT(OUT) :: pclemi_nocont(klon, klev) ! cloud emissivity without contrails [-]
     104  REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [m]
     105  !--AB
     106
    94107  ! Local variables
    95108  !----------------
    … …  
    118131    reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
    119132    reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
    120     icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, rcontrail)
     133    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
     134    !--AB for contrails
     135    rcontrail, pclc_nocont, pcltau_nocont, pclemi_nocont, pch_nocont, &
     136    pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
    121137  ELSE
    122138    CALL nuage (paprs, pplay, &

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

    r5576 r5589  
    99    reliq_pi, reice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
    1010    reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
    11     icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, rcontrail)
     11    icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, &
     12    !--AB contrails
     13    rcontrail, pclc_nocont, pcltau_nocont, pclemi_nocont, pch_nocont, &
     14    pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
    1215
    1316  USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc
    … …  
    7275  LOGICAL, INTENT(IN) :: ptconv(klon, klev) ! flag for grid points affected by deep convection
    7376
    74   REAL, INTENT(IN) :: rcontrail(klon, klev) ! ratio of contrails to total cloud fraction, used only if ok_plane_contrail=y [-]
    75 
    7677  ! inout:
    7778  REAL, INTENT(INOUT) :: pclc(klon, klev) ! cloud fraction for radiation [-]
    … …  
    116117  REAL, INTENT(INOUT) :: icefrac_optics(klon, klev)! ice fraction in clouds seen by radiation [-]
    117118
     119  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
     120  REAL, INTENT(IN) :: rcontrail(klon, klev)      ! ratio of contrails to total cloud fraction [-]
     121  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
     122  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
     123  REAL, INTENT(OUT) :: xfiwp_nocont(klon)        ! ice water path (seen by radiation) without contrails [kg/m2]
     124  REAL, INTENT(OUT) :: xfiwc_nocont(klon, klev)  ! ice water content seen by radiation without contrails [kg/kg]
     125  REAL, INTENT(OUT) :: pclc_nocont(klon, klev)   ! cloud fraction for radiation without contrails [-]
     126  REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [m]
     127  REAL, INTENT(OUT) :: pclemi_nocont(klon, klev) ! cloud emissivity without contrails [-]
     128  REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [m]
     129  !--AB
     130
    118131  ! Local variables
    119132  !----------------
    … …  
    173186  xflwc = 0.D0
    174187  xfiwc = 0.D0
     188  !--AB
     189  IF (ok_plane_contrail) THEN
     190    xfiwp_nocont = 0.D0
     191    xfiwc_nocont = 0.D0
     192    reice_nocont = 0.
     193  ENDIF
    175194
    176195  reliq = 0.
    … …  
    222241      ENDDO
    223242    ENDDO
     243
     244    !--AB
     245    IF (ok_plane_contrail) THEN
     246      !--If contrails are activated, we diagnose iwc without contrails
     247      DO k = 1, klev
     248        DO i = 1, klon
     249          pclc_nocont(i,k) = pclc(i,k) * ( 1. - rcontrail(i,k) )
     250          xfiwc_nocont(i, k) = icefrac_optics(i, k)*radocond(i, k)*( 1. - rcontrail(i,k) )
     251        ENDDO
     252      ENDDO
     253    ENDIF
    224254  ENDIF
    225255
    … …  
    363393            !--vs contrail cirrus in the gridbox
    364394            !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
    365             rei = rei * ( 1. - rcontrail(i,k) ) + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
     395            rei = rei * ( 1. - rcontrail(i,k) ) &
     396                + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
    366397          ENDIF
    367398          pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
    … …  
    415446        pcltau(i, k) = 0.0
    416447        pclemi(i, k) = 0.0
     448
     449        !--AB
     450        IF ( ok_plane_contrail ) THEN
     451          reice_nocont(i,k) = 0.
     452          pclc_nocont(i,k) = 0.
     453          pcltau_nocont(i,k) = 0.
     454          pclemi_nocont(i,k) = 0.
     455        ENDIF
    417456
    418457      ELSE
    … …  
    495534        IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
    496535        IF ( ok_plane_contrail ) THEN
     536          !--Diagnostics of clouds emissivity, optical depth and ice crystal radius
     537          !--without contrails
     538          IF ( pclc_nocont(i,k) .LE. seuil_neb ) THEN
     539            reice_nocont(i,k) = 0.
     540            pclc_nocont(i,k) = 0.
     541            pcltau_nocont(i,k) = 0.
     542            pclemi_nocont(i, k) = 0.
     543          ELSE
     544            reice_nocont(i,k) = rei
     545            pcltau_nocont(i,k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/rei)
     546            k_ice = k_ice0 + 1.0/rei
     547            pclemi_nocont(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
     548          ENDIF
     549
    497550          !--If contrails are activated, rei is a weighted average between the natural
    498551          !--rei and the contrails rei, with the weights being the fraction of natural
    499552          !--vs contrail cirrus in the gridbox
    500553          !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
    501           rei = rei * ( 1. - rcontrail(i,k) ) + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
    502         ENDIF
     554          rei = rei * ( 1. - rcontrail(i,k) ) &
     555              + re_ice_crystals_contrails * 1.E6 * rcontrail(i,k)
     556        ENDIF
     557
    503558        pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
    504559          rei)
    … …  
    519574      xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
    520575      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
     576      !--AB
     577      IF ( ok_plane_contrail ) THEN
     578          xfiwp_nocont(i) = xfiwp_nocont(i) + xfiwc_nocont(i, k)*rhodz(i, k)
     579      ENDIF
    521580
    522581    ENDDO
    … …  
    630689    pcl(i) = 1. - pcl(i)
    631690  ENDDO
     691
     692  !--AB we recompute high cloud cover with contrails only
     693  !--and without contrails
     694  IF ( ok_plane_contrail ) THEN
     695    DO i = 1, klon
     696      zclear(i) = 1.
     697      zcloud(i) = 0.
     698      zcloudh(i) = 0.
     699      pch_nocont(i) = 1.
     700      pct_cont(i) = 1.
     701    ENDDO
     702
     703    DO k = klev, 1, -1
     704      DO i = 1, klon
     705        zclear(i) = zclear(i)*(1.-max(pclc(i,k)*rcontrail(i,k),zcloud(i)))/(1.-min(real( &
     706          zcloud(i),kind=8),1.-zepsec))
     707        pct_cont(i) = 1. - zclear(i)
     708        zcloud(i) = pclc(i, k)*rcontrail(i,k)
     709        IF (paprs(i,k)<prmhc) THEN
     710          pch_nocont(i) = pch_nocont(i)*(1.-max(pclc_nocont(i,k),zcloudh(i)))/(1.-min(real(zcloudh &
     711            (i),kind=8),1.-zepsec))
     712          zcloudh(i) = pclc_nocont(i, k)
     713        ENDIF
     714      ENDDO
     715    ENDDO
     716
     717    DO i = 1, klon
     718      pch_nocont(i) = 1. - pch_nocont(i)
     719    ENDDO
     720  ENDIF ! ok_plane_contrail
    632721
    633722  ! ========================================================

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

    r5580 r5589  
    3131     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
    3232     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
    33      dqsmelt, dqsfreez, dcfres, dqsres, dqvcres)
     33     dqsmelt, dqsfreez, dqised, dcfsed, dqvcsed)
    3434
    3535!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    … …  
    264264  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
    265265  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
    266   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcfres         !--cloud fraction tendency due to resuspension of ice crystals [kg/kg/s]
    267   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsres         !--snow tendency due to resuspension of ice crystals [kg/kg/s]
    268   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvcres        !--cloud water vapor tendency due to resuspension of ice crystals [kg/kg/s]
     266  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqised         !--ice water content tendency due to sedmentation of ice crystals [kg/kg/s]
     267  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [kg/kg/s]
     268  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
    269269
    270270  ! for thermals
    … …  
    306306  REAL, DIMENSION(klon) :: ztupnew
    307307  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
     308  REAL, DIMENSION(klon) :: zqice_sedim ! for sedimentation of ice crystals
    308309  REAL, DIMENSION(klon) :: zrflclr, zrflcld
    309310  REAL, DIMENSION(klon) :: ziflclr, ziflcld
    … …  
    474475zqupnew(:)    = 0.
    475476zqvapclr(:)   = 0.
     477zqice_sedim(:)= 0.
    476478
    477479
    … …  
    529531      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    530532                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
    531                         zqvapclr, zqupnew, &
     533                        zqvapclr, zqupnew, zqice_sedim, &
    532534                        cldfra_in, qvc_in, qliq_in, qice_in, &
    533535                        zrfl, zrflclr, zrflcld, &
    534536                        zifl, ziflclr, ziflcld, &
    535537                        dqreva(:,k), dqssub(:,k), &
    536                         dqsres(:,k), dcfres(:,k), dqvcres(:,k) &
     538                        dcfsed(:,k), dqvcsed(:,k) &
    537539                        )
    538540
    … …  
    974976                            ctot_vol(:,k), ptconv(:,k), &
    975977                            zt, zq, zoliql, zoliqi, zfice, &
    976                             rneb(:,k), znebprecipclr, znebprecipcld, &
     978                            rneb(:,k), zqice_sedim, znebprecipclr, znebprecipcld, &
    977979                            zrfl, zrflclr, zrflcld, &
    978980                            zifl, ziflclr, ziflcld, &
    … …  
    980982                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
    981983                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
    982                             dqsmelt(:,k), dqsfreez(:,k) &
     984                            dqsmelt(:,k), dqsfreez(:,k), dqised(:,k) &
    983985                            )
    984986      DO i = 1, klon

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

    r5580 r5589  
    119119USE lmdz_lscp_ini,        ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
    120120USE lmdz_lscp_ini,        ONLY: eps, temp_nowater, ok_weibull_warm_clouds
    121 USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds, iflag_cloud_sublim_pdf
     121USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds
    122122USE lmdz_lscp_ini,        ONLY: ok_plane_contrail, ok_no_issr_strato
    123123
    124124USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, &
    125                                 mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, &
     125                                beta_pdf_lscp, temp_thresh_pdf_lscp, &
    126126                                std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, &
    127127                                coef_mixing_lscp, coef_shear_lscp, chi_mixing_lscp
    … …  
    361361          qvc(i) = 0.
    362362
     363          !--If all the ice has been sublimated, we sublimate
     364          !--completely the cloud and do not activate the sublimation
     365          !--process
     366          qvapincld_new = 0.
     367
    363368        !--Else, the cloud is adjusted and sublimated
    364369        ELSE
    … …  
    406411          !------------------------------------
    407412
    408           !--If the vapor in cloud is below vapor needed for the cloud to survive
    409           IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN
    410             !--Sublimation of the subsaturated cloud
    411             !--iflag_cloud_sublim_pdf selects the PDF of the ice water content
    412             !--to use.
    413             !--iflag = 1 --> uniform distribution
    414             !--iflag = 2 --> exponential distribution
    415             !--iflag = 3 --> gamma distribution (Karcher et al 2018)
    416 
    417             IF ( iflag_cloud_sublim_pdf .EQ. 1 ) THEN
    418               !--Uniform distribution starting at qvapincld
    419               pdf_e1 = 1. / ( 2. * qiceincld )
    420 
    421               dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1
    422               dqt_sub    = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) &
    423                                        * pdf_e1 / 2.
    424 
    425             ELSEIF ( iflag_cloud_sublim_pdf .EQ. 2 ) THEN
    426               !--Exponential distribution starting at qvapincld
    427               pdf_alpha = 1. / qiceincld
    428               pdf_e1 = EXP( - pdf_alpha * ( qvapincld_new - qvapincld ) )
    429 
    430               dcf_sub(i) = - cldfra(i) * ( 1. - pdf_e1 )
    431               dqt_sub    = - cldfra(i) * ( ( 1. - pdf_e1 ) / pdf_alpha &
    432                                        + qvapincld - qvapincld_new * pdf_e1 )
    433 
    434             ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN
    435               !--Gamma distribution starting at qvapincld
    436               pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld
    437               pdf_y = pdf_alpha * ( qvapincld_new - qvapincld )
    438               pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
    439               pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
    440 
    441               dcf_sub(i) = - cldfra(i) * pdf_e1
    442               dqt_sub    = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 )
    443 
    444             ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN
    445               !--Normal distribution around qvapincld + qiceincld with width
    446               !--constant in the RHi space
    447               pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
    448                     / ( std_subl_pdf_lscp / 100. * qsat(i))
    449               pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
    450               pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
    451 
    452               dcf_sub(i) = - cldfra(i) * pdf_e1
    453               dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
    454                                          - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
    455                                                                                        
    456             ENDIF
     413          !--If the dissipation process must be activated
     414          IF ( qvapincld_new .GT. 0. ) THEN
     415            !--Normal distribution around qvapincld + qiceincld with width
     416            !--constant in the RHi space
     417            pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
     418                  / ( std_subl_pdf_lscp / 100. * qsat(i))
     419            pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
     420            pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
     421
     422            dcf_sub(i) = - cldfra(i) * pdf_e1
     423            dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
     424                                       - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
    457425
    458426            !--Tendencies and diagnostics
    … …  
    464432            qvc(i)    = MAX(0., qvc(i) + dqvc_sub(i))
    465433
    466           ENDIF ! qvapincld .LT. qvapincld_new
     434          ENDIF ! qvapincld_new .GT. 0.
    467435
    468436        ENDIF   ! qiceincld .LT. eps

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

    r5580 r5589  
    162162  !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
    163163
    164   INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf=4       ! iflag for the distribution of water inside ice clouds
    165   !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)
    166 
    167164  REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
    168165  !$OMP THREADPRIVATE(depo_coef_cirrus)
    … …  
    173170  REAL, SAVE, PROTECTED :: std_subl_pdf_lscp=2.              ! [%] standard deviation of the gaussian distribution of water inside ice clouds
    174171  !$OMP THREADPRIVATE(std_subl_pdf_lscp)
    175 
    176   REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] shape factor of the gamma distribution of water inside ice clouds
    177   !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
    178172 
    179173  REAL, SAVE, PROTECTED :: beta_pdf_lscp=1.E-3               ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
    … …  
    260254  !$OMP THREADPRIVATE(ok_corr_vap_evasub)
    261255
     256  LOGICAL, SAVE, PROTECTED :: ok_growth_precip_deposition=.FALSE.
     257  !$OMP THREADPRIVATE(ok_growth_precip_deposition)
     258
    262259  REAL, SAVE, PROTECTED :: cld_lc_lsc_snow=2.e-5            ! snow autoconversion coefficient, stratiform. default from  Chaboureau and PInty 2006
    263260  !$OMP THREADPRIVATE(cld_lc_lsc_snow)
    … …  
    311308  !$OMP THREADPRIVATE(r_snow)
    312309
     310  REAL, SAVE, PROTECTED :: expo_tau_auto_snow=0.1
     311  !$OMP THREADPRIVATE(expo_tau_auto_snow)
     312
    313313  REAL, SAVE, PROTECTED :: tau_auto_snow_min=100.           ! Snow autoconversion minimal timescale (when liquid) [s]
    314314  !$OMP THREADPRIVATE(tau_auto_snow_min)
    … …  
    344344  !$OMP THREADPRIVATE(snow_fallspeed_cld)
    345345
    346   LOGICAL, SAVE, PROTECTED :: ok_precip_resuspension=.FALSE.! Flag to activate the resuspension of ice crystals
    347   !$OMP THREADPRIVATE(ok_precip_resuspension)
    348 
    349   REAL, SAVE, PROTECTED :: snow_thresh_resuspension=1.e-5   ! [kg/m2/s] Threshold flux below which part of the snow flux will be resuspended
    350   !$OMP THREADPRIVATE(snow_thresh_resuspension)
     346  LOGICAL, SAVE, PROTECTED :: ok_ice_sedim=.FALSE.          ! Flag to activate the sedimentation of ice crystals
     347  !$OMP THREADPRIVATE(ok_ice_sedim)
     348
     349  REAL, SAVE, PROTECTED :: ice_fallspeed_sedim=0.1          ! Ice fallspeed velocity for sedimentation [m/s]
     350  !$OMP THREADPRIVATE(ice_fallspeed_sedim)
    351351  !--End of the parameters for poprecip
    352352
    … …  
    448448    CALL getin_p('ok_poprecip',ok_poprecip)
    449449    CALL getin_p('ok_corr_vap_evasub',ok_corr_vap_evasub)
     450    CALL getin_p('ok_growth_precip_deposition',ok_growth_precip_deposition)
    450451    CALL getin_p('rain_int_min',rain_int_min)
    451452    CALL getin_p('gamma_agg',gamma_agg)
    … …  
    467468    CALL getin_p('snow_fallspeed_clr',snow_fallspeed_clr)
    468469    CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
    469     CALL getin_p('ok_precip_resuspension',ok_precip_resuspension)
    470     CALL getin_p('snow_thresh_resuspension',snow_thresh_resuspension)
     470    CALL getin_p('ok_ice_sedim',ok_ice_sedim)
     471    CALL getin_p('ice_fallspeed_sedim',ice_fallspeed_sedim)
    471472    ! for condensation and ice supersaturation
    472473    CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
    473474    CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds)
    474     CALL getin_p('iflag_cloud_sublim_pdf',iflag_cloud_sublim_pdf)
    475475    CALL getin_p('depo_coef_cirrus',depo_coef_cirrus)
    476476    CALL getin_p('capa_cond_cirrus',capa_cond_cirrus)
    477477    CALL getin_p('std_subl_pdf_lscp',std_subl_pdf_lscp)
    478     CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
    479478    CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
    480479    CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
    … …  
    553552    WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
    554553    WRITE(lunout,*) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
     554    WRITE(lunout,*) 'lscp_ini, ok_growth_precip_deposition', ok_growth_precip_deposition
    555555    WRITE(lunout,*) 'lscp_ini, rain_int_min:', rain_int_min
    556556    WRITE(lunout,*) 'lscp_ini, gamma_agg:', gamma_agg
    … …  
    566566    WRITE(lunout,*) 'lscp_ini, snow_fallspeed_clr:', snow_fallspeed_clr
    567567    WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
    568     WRITE(lunout,*) 'lscp_ini, ok_precip_resuspension:', ok_precip_resuspension
    569     WRITE(lunout,*) 'lscp_ini, snow_thresh_resuspension:', snow_thresh_resuspension
     568    WRITE(lunout,*) 'lscp_ini, ok_ice_sedim:', ok_ice_sedim
     569    WRITE(lunout,*) 'lscp_ini, ice_fallspeed_sedim:', ice_fallspeed_sedim
    570570    ! for condensation and ice supersaturation
    571571    WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
    572572    WRITE(lunout,*) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds
    573573    WRITE(lunout,*) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
    574     WRITE(lunout,*) 'lscp_ini, iflag_cloud_sublim_pdf:', iflag_cloud_sublim_pdf
    575574    WRITE(lunout,*) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
    576575    WRITE(lunout,*) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus
    577576    WRITE(lunout,*) 'lscp_ini, std_subl_pdf_lscp:', std_subl_pdf_lscp
    578     WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
    579577    WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
    580578    WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp

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

    r5581 r5589  
    540540          IF (zoliqi(i) .GT. 0.) THEN
    541541            zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
    542               +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
     542              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)/zneb(i)*zrho(i)*ffallv)**(-1./dice_velo))
    543543          ELSE
    544544            zfroi=0.
    … …  
    713713SUBROUTINE poprecip_precld( &
    714714           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
    715            qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
     715           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, qice_sedim, &
    716716           cldfra, qvc, qliq, qice, &
    717717           rain, rainclr, raincld, snow, snowclr, snowcld, &
    718            dqreva, dqssub, dqsres, dcfres, dqvcres &
     718           dqreva, dqssub, dcfsed, dqvcsed &
    719719           )
    720720
    … …  
    722722USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
    723723USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
     724USE lmdz_lscp_ini, ONLY : ok_growth_precip_deposition, ok_ice_sedim
    724725USE lmdz_lscp_ini, ONLY : eps, temp_nowater
    725 USE lmdz_lscp_ini, ONLY : ok_precip_resuspension, snow_thresh_resuspension
    726726USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
    727727
    … …  
    749749REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
    750750REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
     751REAL,    INTENT(IN),    DIMENSION(klon) :: qice_sedim     !--ice water content that sedimentated from the layer above [kg/kg]
    751752
    752753REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
    … …  
    765766REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
    766767REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
    767 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsres         !--snow tendency due to resuspension of ice crystals [kg/kg/s]
    768 REAL,    INTENT(OUT),   DIMENSION(klon) :: dcfres         !--cloud fraction tendency due to resuspension of ice crystals [kg/kg/s]
    769 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqvcres        !--cloud water vapor tendency due to resuspension of ice crystals [kg/kg/s]
     768REAL,    INTENT(OUT),   DIMENSION(klon) :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [s-1]
     769REAL,    INTENT(OUT),   DIMENSION(klon) :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
    770770
    771771
    … …  
    791791!--Specific heat constant
    792792REAL :: cpair, cpw
    793 !--For resuspension of ice crystals
    794 REAL :: dsnowres, qiceinprecip
    795793
    796794!--Initialisation
    … …  
    798796dqreva(:) = 0.
    799797dqssub(:) = 0.
    800 dqrevap   = 0.
    801 dqssubl   = 0.
    802 IF ( ok_precip_resuspension ) THEN
    803   dqsres(:) = 0.
    804   dcfres(:) = 0.
    805   dqvcres(:) = 0.
     798IF ( ok_ice_sedim .AND. ok_ice_supersat ) THEN
     799  dcfsed(:) = 0.
     800  dqvcsed(:) = 0.
    806801ENDIF
    807802
    … …  
    891886DO i = 1, klon
    892887
     888  dqrevap = 0.
     889  dqssubl = 0.
    893890  !--If there is precipitation from the layer above
    894891  IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
    … …  
    924921      IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN
    925922        qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
     923      ELSE
     924        qvapclr = 0.
     925      ENDIF
     926      IF (  precipfraccld(i)  .GT. eps ) THEN
     927        qvapcld = MAX(qtotupnew(i)-qvapclrup(i) , 0.) / qtotupnew(i) * qvap(i) /  precipfraccld(i)
     928      ELSE
     929        qvapcld = 0.
    926930      ENDIF
    927931    ELSE
    … …  
    929933      !--for the evap / subl process
    930934      qvapclr = qvap(i)
    931     ENDIF
     935      qvapcld = qvap(i)
     936    ENDIF
     937
     938
     939    !--------------------
     940    !-- ICE SEDIMENTATION
     941    !--------------------
     942
     943    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
     944    !--added to the total water content of the gridbox
     945    IF ( ok_ice_sedim .AND. ( qice_sedim(i) .GT. 0. ) ) THEN
     946      !--Vapor is updated after evaporation/sublimation (it is increased)
     947      qvap(i) = qvap(i) + qice_sedim(i)
     948      !--Air and precip temperature (i.e., gridbox temperature)
     949      !--is updated due to latent heat cooling
     950      temp(i) = temp(i) &
     951              - qice_sedim(i) * RLSTT / RCPD &
     952              / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
     953   
     954      IF ( ok_ice_supersat ) THEN
     955        !--If ice supersaturation is activated, the cloud properties are prognostic.
     956        !--The falling ice is then considered a new cloud in the gridbox.
     957        !--BEWARE with this parameterisation, we can create a new cloud with a much
     958        !--different ice water content and water vapor content than the existing cloud
     959        !--(if it exists). This implies than unphysical fluxes of ice and vapor
     960        !--occur between the existing cloud and the newly formed cloud (from sedimentation).
     961        !--Note also that currently, the precipitation scheme assume a maximum
     962        !--random overlap, meaning that the new formed clouds will not be affected
     963        !--by vertical wind shear.
     964        ! Maybe we should reduce dcfsed?
     965        dcfsed(i) = precipfracclr_tmp(i)
     966        dqvcsed(i) = qvapclr * dcfsed(i)
     967        !--Add tendencies
     968        cldfra(i) = cldfra(i) + dcfsed(i)
     969        qvc(i) = qvc(i) + dqvcsed(i)
     970        qice(i) = qice(i) + qice_sedim(i)
     971        !--Normalisation
     972        dcfsed(i) = dcfsed(i) / dtime
     973        dqvcsed(i) = dqvcsed(i) / dtime
     974      ENDIF
     975    ENDIF ! ok_ice_sedim
     976
    932977
    933978    !---------------------------
    … …  
    9971042      ENDIF
    9981043
     1044    ELSE
     1045      !--All the precipitation is sublimated if the fraction is zero
     1046      drainclreva = - rainclr_tmp(i)
     1047      dsnowclrsub = - snowclr_tmp(i)
     1048
    9991049    ENDIF ! precipfracclr_tmp .GT. eps
    10001050
    … …  
    10041054    !---------------------------
    10051055
    1006     IF ( ok_unadjusted_clouds .AND. ( temp(i) .LE. temp_nowater ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
     1056    IF ( ok_growth_precip_deposition .AND. ( temp(i) .LE. RTT ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN
    10071057      !--Evaporation of liquid precipitation coming from above
    10081058      !--in the cloud only
    … …  
    10111061      !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly)
    10121062      !--which does not need a barrier on raincld, because included in the formula
    1013       !draincldeva = precipfraccld_tmp(i) * MAX(0., &
    1014       !            - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
    1015       !            + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
    1016       !            )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
     1063      draincldeva = precipfraccld_tmp(i) * MAX(0., &
     1064                  - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) &
     1065                  + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) &
     1066                  )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i)
    10171067               
    10181068      !--Evaporation is limited by 0
    10191069      !--NB. with ok_ice_supersat activated, this barrier should be useless
    1020       !draincldeva = MIN(0., draincldeva)
    1021       draincldeva = 0.
     1070      draincldeva = MIN(0., draincldeva)
    10221071
    10231072
    … …  
    10811130    raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
    10821131    snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)
    1083    
    1084 
    1085     !--If ok_precip_to_clouds is set to .TRUE., some of the precipitations can go back to clouds
    1086     !--Only the snow flux can go back to cloud, and only if there is no rain flux
    1087     !--(mixed-phase clouds)
    1088     IF ( ok_precip_resuspension ) THEN
    1089 
    1090       IF ( (snowclr_tmp(i) .GT. 0.) .AND. (precipfracclr_tmp(i) .GT. eps) &
    1091               .AND. (rainclr_tmp(i) .LT. eps) ) THEN
    1092         dsnowres = snowclr_tmp(i) * EXP( - snowclr_tmp(i) / precipfracclr_tmp(i) &
    1093                                            / snow_thresh_resuspension )
    1094         IF ( dsnowres .GT. 0. ) THEN
    1095           dqsres(i) = dqsres(i) - dsnowres / dhum_to_dflux(i)
    1096           !--The following line determines the in-cloud ice water content of the newly formed cloud
    1097           !--It is a linear combination between the in-precip ice water content (left term) and the
    1098           !--in-existing-cloud ice water content (right term). This is done so that the existing
    1099           !--cloud is not too modified, i.e. unphysical ice fluxes going from the existing cloud
    1100           !--to the new cloud. If the existing cloud is already large, the newly formed cloud
    1101           !--fraction is reduced to match the in-cloud ice water content of the existing cloud.
    1102           !--NB. this is done on qi, not on qvc. Not sure what impact that implies
    1103           qiceinprecip = - dqsres(i) / precipfracclr_tmp(i) * ( 1. - cldfra(i) ) + qice(i)
    1104           dcfres(i) = - dqsres(i) / qiceinprecip
    1105           dqvcres(i) = qvapclr * dcfres(i)
    1106           !--Add tendencies
    1107           snowcld_tmp(i) = snowcld_tmp(i) + dsnowres
    1108         ENDIF
    1109       ENDIF
    1110 
    1111       IF ( (snowcld_tmp(i) .GT. 0.) .AND. (precipfraccld_tmp(i) .GT. 0.) &
    1112               .AND. (raincld_tmp(i) .LT. eps) ) THEN
    1113         dsnowres = snowcld_tmp(i) * EXP( - snowcld_tmp(i) / precipfraccld_tmp(i) &
    1114                                            / snow_thresh_resuspension )
    1115         IF ( dsnowres .GT. 0. ) THEN
    1116           dqsres(i) = dqsres(i) - dsnowres / dhum_to_dflux(i)
    1117           !--Add tendencies
    1118           snowcld_tmp(i) = snowcld_tmp(i) + dsnowres
    1119         ENDIF
    1120       ENDIF
    1121 
    1122       !--Add tendencies
    1123       cldfra(i) = cldfra(i) + dcfres(i)
    1124       qvc(i) = qvc(i) + dqvcres(i)
    1125       qice(i) = qice(i) - dqsres(i)
    1126       !--Recall that qvap is here the total water in the gridbox
    1127       qvap(i) = qvap(i) - dqsres(i)
    1128     ENDIF
    11291132
    11301133
    … …  
    13231326SUBROUTINE poprecip_postcld( &
    13241327           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
    1325            temp, qvap, qliq, qice, icefrac, cldfra, &
     1328           temp, qvap, qliq, qice, icefrac, cldfra, qice_sedim, &
    13261329           precipfracclr, precipfraccld, &
    13271330           rain, rainclr, raincld, snow, snowclr, snowcld, &
    13281331           qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, &
    1329            dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
     1332           dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez, dqised)
    13301333
    13311334USE lmdz_lscp_ini, ONLY : prt_level, lunout
    … …  
    13381341                          rho_rain, r_rain, r_snow, rho_ice,                   &
    13391342                          tau_auto_snow_min, tau_auto_snow_max,                &
    1340                           thresh_precip_frac, eps,                             &
     1343                          expo_tau_auto_snow, thresh_precip_frac, eps,         &
    13411344                          gamma_melt, alpha_freez, beta_freez, temp_nowater,   &
    13421345                          iflag_cloudth_vert, iflag_rain_incloud_vol,          &
    13431346                          cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez,       &
    13441347                          rain_fallspeed_clr, rain_fallspeed_cld,              &
    1345                           snow_fallspeed_clr, snow_fallspeed_cld
     1348                          snow_fallspeed_clr, snow_fallspeed_cld,              &
     1349                          ok_ice_sedim, ice_fallspeed_sedim
    13461350
    13471351
    … …  
    13641368REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
    13651369REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
     1370REAL,    INTENT(INOUT), DIMENSION(klon) :: qice_sedim     !--quantity of sedimentated ice [kg/kg]
    13661371
    13671372REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
    … …  
    13881393REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
    13891394REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
     1395REAL,    INTENT(OUT),   DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
    13901396
    13911397
    … …  
    14281434REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
    14291435REAL :: Eff_rain_ice
     1436
     1437!--Ice sedimentation
     1438REAL :: dz, qice_sedim_tmp
    14301439
    14311440
    … …  
    14441453dqrfreez(:) = 0.
    14451454dqsfreez(:) = 0.
     1455IF ( ok_ice_sedim ) dqised(:) = 0.
    14461456
    14471457
    … …  
    15591569        !--tau for snow depends on the ice fraction in mixed-phase clouds     
    15601570        tau_auto_snow = tau_auto_snow_max &
    1561                       + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
     1571                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
    15621572
    15631573        expo_auto_rain = cld_expo_con
    … …  
    15701580        !--tau for snow depends on the ice fraction in mixed-phase clouds     
    15711581        tau_auto_snow = tau_auto_snow_max &
    1572                       + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
     1582                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow)
    15731583
    15741584        expo_auto_rain = cld_expo_lsc
    … …  
    19031913
    19041914
     1915  !---------------------------------------------------------
     1916  !--                  ICE SEDIMENTATION
     1917  !---------------------------------------------------------
     1918  IF ( ok_ice_sedim ) THEN
     1919    qice_sedim_tmp = 0.
     1920
     1921    IF ( temp(i) .LE. temp_nowater ) THEN
     1922      rho = pplay(i) / temp(i) / RD
     1923      dz = ( paprsdn(i) - paprsup(i) ) / RG / rho
     1924      qice_sedim_tmp = qice(i) * MIN(1., ice_fallspeed_sedim * dtime / dz)
     1925      !--Add tendencies
     1926      qice(i) = qice(i) - qice_sedim_tmp
     1927    ENDIF
     1928
     1929    !--Diagnostic tendencies
     1930    dqised(i) = ( qice_sedim(i) - qice_sedim_tmp ) / dtime
     1931    !--Save quantity that will be sedimentated in the layer below
     1932    qice_sedim(i) = qice_sedim_tmp
     1933  ENDIF
     1934
    19051935
    19061936  !--If the local flux of rain+snow in clear/cloudy air is lower than rain_int_min,
    … …  
    19151945  !--Note that this is physically different than what is proposed in LTP thesis.
    19161946  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
    1917   precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
    19181947
    19191948  !--Calculate outputs

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

    r5579 r5589  
    683683      REAL, SAVE, ALLOCATABLE :: dcf_avi(:,:), dqi_avi(:,:), dqvc_avi(:,:)
    684684      !$OMP THREADPRIVATE(dcf_avi, dqi_avi, dqvc_avi)
     685      REAL, SAVE, ALLOCATABLE :: cldfra_nocont(:,:), cldtau_nocont(:,:), cldemi_nocont(:,:)
     686      !$OMP THREADPRIVATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
     687      REAL, SAVE, ALLOCATABLE :: cldh_nocont(:), contcov(:)
     688      !$OMP THREADPRIVATE(cldh_nocont, contcov)
     689      REAL, SAVE, ALLOCATABLE :: fiwp_nocont(:), fiwc_nocont(:,:), ref_ice_nocont(:,:)
     690      !$OMP THREADPRIVATE(fiwp_nocont, fiwc_nocont, ref_ice_nocont)
     691      REAL, SAVE, ALLOCATABLE :: topsw_nocont(:), solsw_nocont(:)
     692      !$OMP THREADPRIVATE(topsw_nocont, solsw_nocont)
     693      REAL, SAVE, ALLOCATABLE :: toplw_nocont(:), sollw_nocont(:)
     694      !$OMP THREADPRIVATE(toplw_nocont, sollw_nocont)
    685695
    686696!-- LSCP - mixed phase clouds variables
    … …  
    717727      REAL, SAVE, ALLOCATABLE :: dqsfreez(:,:)
    718728      !$OMP THREADPRIVATE(dqsfreez)
    719       REAL, SAVE, ALLOCATABLE :: dcfres(:,:)
    720       !$OMP THREADPRIVATE(dcfres)
    721       REAL, SAVE, ALLOCATABLE :: dqsres(:,:)
    722       !$OMP THREADPRIVATE(dqsres)
    723       REAL, SAVE, ALLOCATABLE :: dqvcres(:,:)
    724       !$OMP THREADPRIVATE(dqvcres)
     729      REAL, SAVE, ALLOCATABLE :: dqised(:,:)
     730      !$OMP THREADPRIVATE(dqised)
     731      REAL, SAVE, ALLOCATABLE :: dcfsed(:,:)
     732      !$OMP THREADPRIVATE(dcfsed)
     733      REAL, SAVE, ALLOCATABLE :: dqvcsed(:,:)
     734      !$OMP THREADPRIVATE(dqvcsed)
    725735
    726736! variables for stratospheric aerosol
    … …  
    12381248      ALLOCATE(contfra(klon,klev), dcontfra_cir(klon,klev))
    12391249      ALLOCATE(dcf_avi(klon,klev), dqi_avi(klon,klev), dqvc_avi(klon,klev))
     1250      ALLOCATE(cldfra_nocont(klon,klev), cldtau_nocont(klon,klev), cldemi_nocont(klon,klev))
     1251      ALLOCATE(cldh_nocont(klon), contcov(klon))
     1252      ALLOCATE(fiwp_nocont(klon), fiwc_nocont(klon,klev), ref_ice_nocont(klon,klev))
     1253      ALLOCATE(topsw_nocont(klon), solsw_nocont(klon))
     1254      ALLOCATE(toplw_nocont(klon), sollw_nocont(klon))
    12401255
    12411256!-- LSCP - POPRECIP variables
    … …  
    12441259      ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
    12451260      ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
    1246       ALLOCATE(dcfres(klon,klev), dqsres(klon,klev), dqvcres(klon,klev))
     1261      ALLOCATE(dqised(klon,klev), dcfsed(klon,klev), dqvcsed(klon,klev))
    12471262
    12481263IF (CPPKEY_STRATAER) THEN
    … …  
    16461661      DEALLOCATE(contfra, dcontfra_cir)
    16471662      DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
     1663      DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
     1664      DEALLOCATE(cldh_nocont, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)
     1665      DEALLOCATE(topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
    16481666
    16491667!-- LSCP - POPRECIP variables
    … …  
    16521670      DEALLOCATE(dqrauto, dqrcol, dqrmelt, dqrfreez)
    16531671      DEALLOCATE(dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
    1654       DEALLOCATE(dcfres, dqsres, dqvcres)
     1672      DEALLOCATE(dqised, dcfsed, dqvcsed)
    16551673
    16561674IF (CPPKEY_STRATAER) THEN

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

    r5579 r5589  
    16111611  TYPE(ctrl_out), SAVE :: o_dqsfreez = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    16121612    'dqsfreez', 'LS snow tendency due to freezing', 'kg/kg/s', (/ ('', i=1, 10) /))
    1613   TYPE(ctrl_out), SAVE :: o_dcfres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    1614     'dcfres', 'LS cloud fraction tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
    1615   TYPE(ctrl_out), SAVE :: o_dqsres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    1616     'dqsres', 'LS snow tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
    1617   TYPE(ctrl_out), SAVE :: o_dqvcres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    1618     'dqvcres', 'LS cloud water vapor tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
     1613  TYPE(ctrl_out), SAVE :: o_dqised = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     1614    'dqised', 'Ice water content tendency due to sedimentation', 'kg/kg/s', (/ ('', i=1, 10) /))
     1615  TYPE(ctrl_out), SAVE :: o_dcfsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     1616    'dcfsed', 'Cloud fraction tendency due to sedimentation', 's-1', (/ ('', i=1, 10) /))
     1617  TYPE(ctrl_out), SAVE :: o_dqvcsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     1618    'dqvcsed', 'Cloud water vapor tendency due to sedimentation', 'kg/kg/s', (/ ('', i=1, 10) /))
    16191619  TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
    16201620    'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))
    … …  
    22142214  TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    22152215    'flighth2o', 'Aviation H2O flight emission', 'kg H2O/s/m^3', (/ ('', i=1, 10) /))
     2216  TYPE(ctrl_out), SAVE :: o_cldfra_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2217    'cldfra_nocont', 'Cloud fraction w/o contrails', '-', (/ ('', i=1, 10) /))
     2218  TYPE(ctrl_out), SAVE :: o_cldtau_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2219    'cldtau_nocont', 'Cloud optical thickness w/o contrails', '1', (/ ('', i=1, 10) /))
     2220  TYPE(ctrl_out), SAVE :: o_cldemi_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2221    'cldemi_nocont', 'Cloud optical emissivity w/o contrails', '1', (/ ('', i=1, 10) /))
     2222  TYPE(ctrl_out), SAVE :: o_cldh_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2223    'cldh_nocont', 'High-level cloudiness w/o contrails', '-', (/ ('', i=1, 10) /))
     2224  TYPE(ctrl_out), SAVE :: o_contcov = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2225    'contcov', 'Total contrails cover', '-', (/ ('', i=1, 10) /))
     2226  TYPE(ctrl_out), SAVE :: o_iwp_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2227    'iwp_nocont', 'Cloud ice water path w/o contrails', 'kg/m2', (/ ('', i=1, 10) /))
     2228  TYPE(ctrl_out), SAVE :: o_iwc_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2229    'iwc_nocont', 'Cloud ice water content seen by radiation w/o contrails', 'kg/kg', (/ ('', i=1, 10) /))
     2230  TYPE(ctrl_out), SAVE :: o_ref_ice_nocont = ctrl_out((/ 4, 10, 10, 10, 10, 10, 11, 11, 11, 11/), &
     2231    'ref_ice_nocont', 'Effective radius of ice crystals w/o contrails', 'microns', (/ ('', i=1, 10) /))
     2232  TYPE(ctrl_out), SAVE :: o_tops_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2233    'tops_nocont', 'Solar rad. at TOA w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
     2234  TYPE(ctrl_out), SAVE :: o_topl_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2235    'topl_nocont', 'IR rad. at TOA w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
     2236  TYPE(ctrl_out), SAVE :: o_sols_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2237    'sols_nocont', 'Solar rad. at surf. w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
     2238  TYPE(ctrl_out), SAVE :: o_soll_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2239    'soll_nocont', 'IR rad. at surface w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
     2240  TYPE(ctrl_out), SAVE :: o_nettop_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2241    'nettop_nocont', 'Net dn radiatif flux at TOA w/o contrails', 'W/m2', (/ ('', i=1, 10) /))
    22162242
    22172243!!!!!!!!!!!!! Sorties niveaux standards de pression NMC

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

    r5579 r5589  
    145145         o_qrainlsc, o_qsnowlsc, o_dqreva, o_dqrauto, o_dqrcol, o_dqrmelt, o_dqrfreez, &
    146146         o_dqssub, o_dqsauto, o_dqsagg, o_dqsrim, o_dqsmelt, o_dqsfreez, &
    147          o_dcfres, o_dqsres, o_dqvcres, &
     147         o_dqised, o_dcfsed, o_dqvcsed, &
    148148         o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, &
    149149         o_dqsphy, o_dqsphy2d, o_dqbsphy, o_dqbsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
    … …  
    232232         o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
    233233         o_dcontfracir, o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
     234         o_cldfra_nocont, o_cldtau_nocont, o_cldemi_nocont, o_cldh_nocont, &
     235         o_contcov, o_iwp_nocont, o_iwc_nocont, o_ref_ice_nocont, &
     236         o_tops_nocont, o_topl_nocont, o_sols_nocont, o_soll_nocont, &
    234237!--interactive CO2
    235238         o_flx_co2_ocean, o_flx_co2_ocean_cor, &
    … …  
    269272         o_SAD_sulfate, o_reff_sulfate, o_sulfmmr, o_nd_mode, o_sulfmmr_mode
    270273
    271     USE lmdz_lscp_ini, ONLY: ok_poprecip, ok_precip_resuspension
     274    USE lmdz_lscp_ini, ONLY: ok_poprecip, ok_ice_sedim
    272275
    273276    USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
    … …  
    358361         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    359362         dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
     363         cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
     364         contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
     365         topsw_nocont, toplw_nocont, solsw_nocont, sollw_nocont, &
    360366         alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
    361367         alp_bl_stat, alp_bl_fluct_tke, slab_wfbils, &
    … …  
    384390         dqrauto,dqrcol,dqrmelt,dqrfreez, &
    385391         dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
    386          dcfres, dqsres, dqvcres, &
     392         dqised, dcfsed, dqvcsed, &
    387393         d_t_dyn,  &
    388394         d_q_dyn,  d_ql_dyn, d_qs_dyn, d_qbs_dyn,  &
    … …  
    20952101           CALL histwrite_phy(o_dqsfreez, dqsfreez)
    20962102           CALL histwrite_phy(o_dqsrim, dqsrim)
    2097            IF ( ok_precip_resuspension ) THEN
    2098             CALL histwrite_phy(o_dcfres, dcfres)
    2099             CALL histwrite_phy(o_dqsres, dqsres)
    2100             CALL histwrite_phy(o_dqvcres, dqvcres)
     2103           IF ( ok_ice_sedim ) THEN
     2104            CALL histwrite_phy(o_dqised, dqised)
     2105            CALL histwrite_phy(o_dcfsed, dcfsed)
     2106            CALL histwrite_phy(o_dqvcsed, dqvcsed)
    21012107           ENDIF
    21022108           ELSE
    … …  
    21542160         CALL histwrite_phy(o_dqiavi, dqi_avi)
    21552161         CALL histwrite_phy(o_dqvcavi, dqvc_avi)
     2162         CALL histwrite_phy(o_cldfra_nocont, cldfra_nocont)
     2163         CALL histwrite_phy(o_cldtau_nocont, cldtau_nocont)
     2164         CALL histwrite_phy(o_cldemi_nocont, cldemi_nocont)
     2165         CALL histwrite_phy(o_cldh_nocont, cldh_nocont)
     2166         CALL histwrite_phy(o_contcov, contcov)
     2167         CALL histwrite_phy(o_iwp_nocont, fiwp_nocont)
     2168         CALL histwrite_phy(o_iwc_nocont, fiwc_nocont)
     2169         CALL histwrite_phy(o_ref_ice_nocont, ref_ice_nocont)
     2170         IF (ok_rad_contrail) THEN
     2171           IF (vars_defined) zx_tmp_fi2d = topsw_nocont * swradcorr
     2172           CALL histwrite_phy(o_tops_nocont, topsw_nocont)
     2173           CALL histwrite_phy(o_topl_nocont, toplw_nocont)
     2174           IF (vars_defined) zx_tmp_fi2d = topsw_nocont * swradcorr - toplw_nocont
     2175           CALL histwrite_phy(o_nettop_nocont, zx_tmp_fi2d)
     2176           IF (vars_defined) zx_tmp_fi2d = solsw_nocont * swradcorr
     2177           CALL histwrite_phy(o_sols_nocont, solsw_nocont)
     2178           CALL histwrite_phy(o_soll_nocont, sollw_nocont)
     2179         ENDIF
    21562180       ENDIF
    21572181       IF (ok_plane_h2o) THEN

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

    r5579 r5589  
    187187       ! proprecip
    188188       qraindiag, qsnowdiag, &
    189        dqreva, dqssub, dcfres, dqsres, dqvcres, &
     189       dqreva, dqssub, dqised, dcfsed, dqvcsed, &
    190190       dqrauto,dqrcol,dqrmelt,dqrfreez, &
    191191       dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
    … …  
    331331       contfra, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    332332       dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, &
     333       cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
     334       contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
     335       topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont, &
    333336       !
    334337       stratomask, &
    … …  
    12261229    ! "CRF off"
    12271230    REAL, dimension(klon, klev) :: cldfrarad   ! fraction nuageuse
     1231    !--AB contrails
     1232    REAL, dimension(klon, klev) :: cldfrarad_nocont   ! fraction nuageuse without contrails
    12281233
    12291234    REAL :: calday, zxsnow_dummy(klon)
    … …  
    14301435          WRITE (lunout, *) ' ok_plane_contrail=y requires 6 H2O tracers ', &
    14311436              '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c, H2O_a) but nqo=', nqo, '. Might as well stop here.'
     1437          abort_message='see above'
     1438          CALL abort_physic(modname,abort_message,1)
     1439       ENDIF
     1440
     1441       IF (ok_rad_contrail.AND..NOT.ok_plane_contrail) THEN
     1442          WRITE (lunout, *) ' ok_rad_contrail=y requires ok_plane_contrail=y '
    14321443          abort_message='see above'
    14331444          CALL abort_physic(modname,abort_message,1)
    … …  
    21112122       !=============================================================
    21122123
     2124       !--Read the aviation emissions
     2125       IF ( ok_plane_h2o .OR. ok_plane_contrail ) CALL read_aviation_emissions(klon, klev)
     2126
    21132127       IF (using_xios) THEN
    21142128         ! Get "missing_val" value from XML files (from temperature variable)
    … …  
    21232137         IF (is_omp_master) CALL xios_get_field_attr("temp",default_value=missing_val)
    21242138         CALL bcast_omp(missing_val)
    2125 
    2126          !--Read the aviation emissions
    2127          IF ( ok_plane_h2o .OR. ok_plane_contrail ) CALL read_aviation_emissions(klon, klev)
    21282139       !
    21292140       ! Now we activate some double radiation call flags only if some
    … …  
    39243935         qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
    39253936         dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez, &
    3926          dcfres, dqsres, dqvcres)
     3937         dqised, dcfsed, dqvcsed)
    39273938
    39283939
    … …  
    44704481               ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
    44714482               reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  &
    4472                zfice, dNovrN, ptconv, rnebcon, clwcon, rcont_seri)
     4483               zfice, dNovrN, ptconv, rnebcon, clwcon, &
     4484               !--AB contrails
     4485               rcont_seri, cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
     4486               contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)
    44734487
    44744488       !
    … …  
    44794493       cldemirad   = cldemi
    44804494       cldfrarad   = cldfra
     4495       !--AB contrails
     4496       IF (ok_rad_contrail) cldfrarad_nocont = cldfra_nocont
    44814497
    44824498       !
    … …  
    45054521             ENDDO
    45064522          ENDDO
     4523
     4524          !--AB contrails
     4525          IF ( ok_rad_contrail ) THEN
     4526            DO k=1, klev
     4527              DO i=1, klon
     4528                cldfrarad_nocont(i,k) = cldfra_nocont(i,k) * beta(i,k)
     4529              ENDDO
     4530            ENDDO
     4531          ENDIF
    45074532          !
    45084533       ELSE
    … …  
    45334558             ENDDO
    45344559          ENDDO
     4560
     4561          !--AB contrails
     4562          IF ( ok_rad_contrail ) THEN
     4563            DO k=1, klev
     4564              DO i=1, klon
     4565                cldfrarad_nocont(i,k) = cldfra_nocont(i,k) * beta(i,k)
     4566              ENDDO
     4567            ENDDO
     4568          ENDIF
    45354569       !
    45364570       ENDIF
    … …  
    46494683               ZLWFT0_i, ZFLDN0, ZFLUP0, &
    46504684               ZSWFT0_i, ZFSDN0, ZFSUP0, &
    4651                cloud_cover_sw)
     4685               cloud_cover_sw, &
     4686               !--AB contrails radiative effects
     4687               cldfrarad_nocont, fiwc_nocont, ref_ice_nocont, &
     4688               topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
    46524689
    46534690          !lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
    … …  
    46604697                        sollwdown(:))
    46614698          cool = cool + betalwoff * (cool0 - cool)
     4699          !--AB contrails
     4700          IF ( ok_rad_contrail ) THEN
     4701            toplw_nocont = toplw_nocont + betalwoff * (toplw0 - toplw_nocont)
     4702            sollw_nocont = sollw_nocont + betalwoff * (sollw0 - sollw_nocont)
     4703          ENDIF
    46624704
    46634705          IF (.NOT. using_xios) THEN
    … …  
    47274769                     ZLWFT0_i, ZFLDN0, ZFLUP0, &
    47284770                     ZSWFT0_i, ZFSDN0, ZFSUP0, &
    4729                      cloud_cover_sw)
     4771                     cloud_cover_sw, &
     4772                     !--AB contrails radiative effects
     4773                     cldfrarad_nocont, fiwc_nocont, ref_ice_nocont, &
     4774                     topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
    47304775          ENDIF !ok_4xCO2atm
    47314776

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

    r5325 r5589  
    4545       ZLWFT0_i, ZFLDN0, ZFLUP0, &
    4646       ZSWFT0_i, ZFSDN0, ZFSUP0, &
    47        cloud_cover_sw)
     47       cloud_cover_sw, &
     48       !--AB contrails radiative effects
     49       cldfra_nocont, fiwc_nocont, ref_ice_nocont, &
     50       topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)
    4851
    4952    ! Modules necessaires
    … …  
    279282    REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZSWFT0_i
    280283    REAL, DIMENSION(kdlon,kflev+1), INTENT(out) :: ZLWFT0_i
     284    !--AB contrails radiative effects
     285    REAL, DIMENSION(klon,klev), INTENT(in)  :: cldfra_nocont
     286    REAL, DIMENSION(klon,klev), INTENT(in)  :: fiwc_nocont
     287    REAL, DIMENSION(klon,klev), INTENT(in)  :: ref_ice_nocont
     288    REAL, DIMENSION(klon),      INTENT(out) :: topsw_nocont
     289    REAL, DIMENSION(klon),      INTENT(out) :: solsw_nocont
     290    REAL, DIMENSION(klon),      INTENT(out) :: toplw_nocont
     291    REAL, DIMENSION(klon),      INTENT(out) :: sollw_nocont
    281292
    282293    ! Local variables
    … …  
    467478    REAL(KIND=8) ZFLCCDWN_i (klon,klev+1)
    468479    REAL(KIND=8) ZFLCCUP_i (klon,klev+1)
     480    !--AB contrails radiative effects
     481    REAL(KIND=8) cldfra_nocont_i(klon,klev)
     482    REAL(KIND=8) fiwc_nocont_i(klon,klev)
     483    REAL(KIND=8) ref_ice_nocont_i(klon,klev)
     484    REAL(KIND=8) ZTOPSWNOCONT(klon)
     485    REAL(KIND=8) ZSOLSWNOCONT(klon)
     486    REAL(KIND=8) ZTOPLWNOCONT(klon)
     487    REAL(KIND=8) ZSOLLWNOCONT(klon)
    469488    ! 3 lignes suivantes a activer pour CCMVAL (MPL 20100412)
    470489    !      REAL(KIND=8) RSUN(3,2)
    … …  
    900919             ref_liq_pi_i(1:klon,k) =ref_liq_pi(1:klon,klev+1-k)
    901920             ref_ice_pi_i(1:klon,k) =ref_ice_pi(1:klon,klev+1-k)
     921             IF (ok_rad_contrail) THEN
     922               !--AB contrails radiative effects
     923               cldfra_nocont_i(1:klon,k)  = cldfra_nocont(1:klon,klev+1-k)
     924               fiwc_nocont_i(1:klon,k)    = fiwc_nocont(1:klon,klev+1-k)
     925               ref_ice_nocont_i(1:klon,k) = ref_ice_nocont(1:klon,klev+1-k)
     926             ENDIF
    902927          ENDDO
    903928          DO k=1,kflev
    … …  
    9801005               ZLWADAERO, & !--NL
    9811006               volmip_solsw, flag_volc_surfstrat, & !--VOLMIP
    982                ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback) ! flags aerosols
     1007               ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat, flag_aer_feedback, & ! flags aerosols
     1008               !--AB contrails radiative effect
     1009               ok_rad_contrail, cldfra_nocont_i, fiwc_nocont_i, ref_ice_nocont_i, &
     1010               ZTOPSWNOCONT, ZSOLSWNOCONT, ZTOPLWNOCONT, ZSOLLWNOCONT)
    9831011
    9841012          !--OB diagnostics
    … …  
    10801108          ZSOLSWCF_AERO(:,2)=ZSOLSWCF_AERO(:,2)*fract(:)
    10811109          ZSOLSWCF_AERO(:,3)=ZSOLSWCF_AERO(:,3)*fract(:)
     1110
     1111          IF (ok_rad_contrail) THEN
     1112            !--AB contrails radiative effect
     1113            ZTOPSWNOCONT(:) = ZTOPSWNOCONT(:) * fract(:)
     1114            ZSOLSWNOCONT(:) = ZSOLSWNOCONT(:) * fract(:)
     1115          ENDIF
    10821116
    10831117          ! ---------
    … …  
    16551689          ENDDO
    16561690       ENDIF
     1691       !--AB radiative effect from contrails
     1692       IF (ok_rad_contrail) THEN
     1693         DO i = 1, kdlon
     1694           topsw_nocont(iof+i) = ZTOPSWNOCONT(i)
     1695           solsw_nocont(iof+i) = ZSOLSWNOCONT(i)
     1696           toplw_nocont(iof+i) = ZTOPLWNOCONT(i)
     1697           sollw_nocont(iof+i) = ZSOLLWNOCONT(i)
     1698         ENDDO
     1699       ENDIF
    16571700       DO k = 1, kflev
    16581701          DO i = 1, kdlon

  • LMDZ6/branches/contrails/libf/phylmd/rrtm/recmwf_aero.F90

    r5294 r5589  
    4040     !..end
    4141     & ok_ade, ok_aie, ok_volcan, flag_aerosol,flag_aerosol_strat,&
    42      & flag_aer_feedback)
     42     & flag_aer_feedback, &
     43     !--AB contrails radiative effect
     44     & ok_rad_contrail, PCLFR_NOCONT, PQIWP_NOCONT, PREF_ICE_NOCONT, &
     45     & PTOPSWNOCONT, PSOLSWNOCONT, PTOPLWNOCONT, PSOLLWNOCONT)
    4346  !--fin
    4447
    … …  
    266269  REAL(KIND=JPRB)   ,INTENT(OUT)   :: volmip_solsw(KPROMA) ! SW clear sky in the case of VOLMIP
    267270  INTEGER, INTENT(IN)              :: flag_volc_surfstrat !--VOlMIP Modif
     271  !--AB contrails radiative effect
     272  LOGICAL           ,INTENT(IN)    :: ok_rad_contrail
     273  REAL(KIND=JPRB)   ,INTENT(IN)    :: PCLFR_NOCONT(KPROMA,KLEV)
     274  REAL(KIND=JPRB)   ,INTENT(IN)    :: PQIWP_NOCONT(KPROMA,KLEV)
     275  REAL(KIND=JPRB)   ,INTENT(IN)    :: PREF_ICE_NOCONT(KPROMA,KLEV)
     276  REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPSWNOCONT(KPROMA), PSOLSWNOCONT(KPROMA) ! No contrails experiment forcing at TOA and surface (SW)
     277  REAL(KIND=JPRB)   ,INTENT(OUT)   :: PTOPLWNOCONT(KPROMA), PSOLLWNOCONT(KPROMA) ! No contrails experiment forcing at TOA and surface (LW)
    268278
    269279  !     ==== COMPUTED IN RADITE ===
    … …  
    346356  REAL(KIND=JPRB) ::  LWUP0_AERO(KPROMA,KLEV+1,5)
    347357  REAL(KIND=JPRB) ::  LWDN0_AERO(KPROMA,KLEV+1,5)
     358  !--AB contrails radiative effect
     359  REAL(KIND=JPRB) :: ZRCLC_NOCONT(KPROMA,KLEV), ZQIWP_NOCONT(KPROMA,KLEV)
     360  REAL(KIND=JPRB) :: PREF_LIQ_NOCONT(KPROMA,KLEV)
     361  REAL(KIND=JPRB) :: PPIZA_NOCONT(KPROMA,KLEV,NSW)
     362  REAL(KIND=JPRB) :: PCGA_NOCONT(KPROMA,KLEV,NSW)
     363  REAL(KIND=JPRB) :: PTAU_NOCONT(KPROMA,KLEV,NSW)
     364  REAL(KIND=JPRB) :: PTAU_LW_NOCONT(KPROMA,KLEV,NLW)
    348365
    349366#include "radlsw.intfb.h"
    … …  
    673690  ENDIF ! .not. AEROSOLFEEDBACK_ACTIVE
    674691
     692
     693  !--Double call to radiative routine for contrails
     694  !--The calculation are done again WITHOUT contrails
     695  IF (ok_rad_contrail) THEN
     696
     697     !--The same base case is used
     698     !--NB. here we could use pointers...
     699     IF ( flag_aerosol .EQ. 0 ) THEN
     700        PREF_LIQ_NOCONT(:,:) = PREF_LIQ_PI(:,:)
     701        PPIZA_NOCONT(:,:,:) = PPIZA_ZERO(:,:,:)
     702        PCGA_NOCONT(:,:,:) = PCGA_ZERO(:,:,:)
     703        PTAU_NOCONT(:,:,:) = PTAU_ZERO(:,:,:)
     704        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_ZERO(:,:,:)
     705     ELSEIF ( .not. ok_ade .AND. .not. ok_aie ) THEN
     706        PREF_LIQ_NOCONT(:,:) = PREF_LIQ_PI(:,:)
     707        PPIZA_NOCONT(:,:,:) = PPIZA_NAT(:,:,:)
     708        PCGA_NOCONT(:,:,:) = PCGA_NAT(:,:,:)
     709        PTAU_NOCONT(:,:,:) = PTAU_NAT(:,:,:)
     710        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_NAT(:,:,:)
     711     ELSEIF ( .not. ok_ade .AND. ok_aie ) THEN
     712        PREF_LIQ_NOCONT(:,:) = PREF_LIQ(:,:)
     713        PPIZA_NOCONT(:,:,:) = PPIZA_NAT(:,:,:)
     714        PCGA_NOCONT(:,:,:) = PCGA_NAT(:,:,:)
     715        PTAU_NOCONT(:,:,:) = PTAU_NAT(:,:,:)
     716        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_NAT(:,:,:)
     717     ELSEIF ( ok_ade .AND. .not. ok_aie ) THEN
     718        PREF_LIQ_NOCONT(:,:) = PREF_LIQ_PI(:,:)
     719        PPIZA_NOCONT(:,:,:) = PPIZA_TOT(:,:,:)
     720        PCGA_NOCONT(:,:,:) = PCGA_TOT(:,:,:)
     721        PTAU_NOCONT(:,:,:) = PTAU_TOT(:,:,:)
     722        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_TOT(:,:,:)
     723     ELSEIF ( ok_ade .AND. ok_aie ) THEN
     724        PREF_LIQ_NOCONT(:,:) = PREF_LIQ(:,:)
     725        PPIZA_NOCONT(:,:,:) = PPIZA_TOT(:,:,:)
     726        PCGA_NOCONT(:,:,:) = PCGA_TOT(:,:,:)
     727        PTAU_NOCONT(:,:,:) = PTAU_TOT(:,:,:)
     728        PTAU_LW_NOCONT(:,:,:) = PTAU_LW_TOT(:,:,:)
     729     ENDIF
     730
     731     DO JK=1,KLEV
     732         DO JL=IBEG,IEND
     733             ZRCLC_NOCONT(JL,JK)=MAX( 0.0_JPRB ,MIN( 1.0_JPRB ,PCLFR_NOCONT(JL,JK)))
     734             IF (ZRCLC_NOCONT(JL,JK) > REPCLC) THEN
     735                 ZQIWP_NOCONT(JL,JK)=PQIWP_NOCONT(JL,JK)
     736             ELSE
     737                 ZQIWP_NOCONT(JL,JK)=REPH2O*ZRCLC_NOCONT(JL,JK)
     738             ENDIF
     739         ENDDO
     740     ENDDO
     741
     742     CALL RADLSW (&
     743          & IBEG  , IEND   , KPROMA  , KLEV  , KMODE , NAER,&
     744          & ZRII0 ,&
     745          & ZRAER , PALBD  , PALBP   , PAPRS , ZRPR  ,&
     746          & ZCCNL , ZCCNO  ,&
     747          & PCCO2 , ZRCLC_NOCONT  , PDP     , PEMIS , ZEMIW ,PSLM    , ZRMU0 , ZPQO3,&
     748          & ZQ    , ZQIWP_NOCONT  , ZQLWP   , ZQS   , ZQRAIN,ZQRAINT ,&
     749          & PTH   , ZRTI   , PTS     , ZNBAS , ZNTOP ,&
     750          & PREF_LIQ_NOCONT, PREF_ICE_NOCONT,&
     751          & ZEMIT , ZFCT   , ZFLT    , ZFCS    , ZFLS  ,&
     752          & ZFRSOD, ZSUDU  , ZUVDF   , ZPARF   , ZPARCF, ZTINCF, PSFSWDIR,&
     753          & PSFSWDIF,PFSDNN, PFSDNV  ,& 
     754          & LRDUST,PPIZA_NOCONT,PCGA_NOCONT,PTAU_NOCONT,PTAU_LW_NOCONT,PFLUX,PFLUC,&
     755          & PFSDN , PFSUP  , PFSCDN  , PFSCUP )
     756
     757     ! save budgets LW and SW at TOA and surface
     758     PTOPSWNOCONT(:) = PFSDN(:,KLEV+1) - PFSUP(:,KLEV+1)
     759     PSOLSWNOCONT(:) = PFSDN(:,1) - PFSUP(:,1)
     760     PTOPLWNOCONT(:) = PFLUX(:,1,KLEV+1) + PFLUX(:,2,KLEV+1)
     761     PSOLLWNOCONT(:) = - PFLUX(:,1,1) - PFLUX(:,2,1)
     762
     763  ENDIF
     764
    675765  !*         4.2     TRANSFORM FLUXES TO MODEL HISTORICAL VARIABLES
    676766
Note: See TracChangeset for help on using the changeset viewer.