Ignore:
Timestamp:
Jun 27, 2025, 7:30:06 PM (5 weeks ago)
Author:
aborella
Message:

Correction to rare floating point errors + added support for the sedimentation of contrails

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5717 r5728  
    1818SUBROUTINE histprecip_precld( &
    1919           klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
    20            zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
     20           zmqc, zneb, znebprecip, znebprecipclr, flsed, flsed_lincont, flsed_circont, &
    2121           zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
    2222           )
    … …  
    4444REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
    4545REAL,    INTENT(INOUT), DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
    46 REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice flux [kg/s/m2]
     46REAL,    INTENT(IN),    DIMENSION(klon) :: flsed          !--sedimentated ice flux [kg/s/m2]
     47REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_lincont  !--linear contrails sedimentated ice flux [kg/s/m2]
     48REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_circont  !--contrail cirrus sedimentated ice flux [kg/s/m2]
    4749
    4850REAL,    INTENT(IN),    DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
    … …  
    132134    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    133135    !--added to the total water content of the gridbox
    134     IF ( icesed_flux(i) .GT. 0. ) THEN
    135       qice_sedim = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     136    IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
     137      qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
     138          / ( paprsdn(i) - paprsup(i) ) * RG * dtime
    136139
    137140      !--Thermalisation of sedimentated ice
    … …  
    339342SUBROUTINE histprecip_postcld( &
    340343        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, pt_pron_clds, &
    341         zdqsdT_raw, zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
     344        zdqsdT_raw, zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
    342345        rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
    343346        zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
    … …  
    394397REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
    395398REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
    396 REAL,    INTENT(INOUT), DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    397399
    398400REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
    … …  
    438440zqpreci(:) = 0.
    439441ziflprev(:) = 0.
    440 
    441 ! AB ICE SEDIM
    442 !icesed_flux(:) = 0.
    443442
    444443
    … …  
    700699
    701700ENDDO
    702 
    703 !---------------------------------------------------------
    704 !--                  ICE SEDIMENTATION
    705 !---------------------------------------------------------
    706 !IF ( ok_ice_sedim ) THEN
    707 IF ( .FALSE. ) THEN ! AB ICE SEDIM
    708 
    709   DO i = 1, klon
    710     IF ( ( zoliqi(i) .GT. 0. ) .AND. ( rneb(i) .GT. eps ) .AND. pt_pron_clds(i) ) THEN
    711 
    712       iwcg = MAX( zrho(i) * zoliqi(i) / zneb(i) * 1000., eps ) ! iwc in g/m3
    713       !tempc = zt(i) - RTT ! celcius temp
    714       !! so-called 'empirical parameterization' in Stubenrauch et al. 2019
    715       !IF ( tempc .GE. -60. ) THEN
    716       !  icevelo = EXP( 1.9714 + 0.00216078 * tempc &
    717       !          - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
    718       !ELSE
    719       !  icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
    720       !ENDIF
    721       !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
    722 
    723       icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16
    724 
    725       qice_sedim = zoliqi(i) * MIN(1., icevelo * dtime / zdz(i))
    726       !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
    727       !--times in the same timestep)
    728       qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
    729       !--Add tendencies
    730       zoliqi(i) = zoliqi(i) - qice_sedim
    731       zoliq(i)  = zoliq(i)  - qice_sedim
    732       !--Convert to flux
    733       icesed_flux(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
    734 
    735       !--Diagnostic tendencies
    736       dqised(i) = dqised(i) - qice_sedim / dtime
    737     ENDIF
    738   ENDDO
    739 ENDIF
    740701
    741702! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
    … …  
    796757SUBROUTINE poprecip_precld( &
    797758           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
    798            qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, icesed_flux, &
     759           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
     760           flsed, flsed_lincont, flsed_circont, &
    799761           cldfra, qvc, qliq, qice, &
    800762           rain, rainclr, raincld, snow, snowclr, snowcld, &
    … …  
    831793REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
    832794REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
    833 REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice water flux [kg/s/m2]
     795REAL,    INTENT(IN),    DIMENSION(klon) :: flsed          !--sedimentated ice water flux [kg/s/m2]
     796REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_lincont  !--sedimentated ice water flux [kg/s/m2]
     797REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_circont  !--sedimentated ice water flux [kg/s/m2]
    834798
    835799REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
    … …  
    933897  cpw = RCPD * RVTMP2
    934898  DO i = 1, klon
    935     IF ( icesed_flux(i) .GT. 0. ) THEN
    936       qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
     899    IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
     900      qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) / dhum_to_dflux(i)
    937901      !--No condensed water so cp=cp(vapor+dry air)
    938902      !-- RVTMP2=rcpv/rcpd-1
    … …  
    10421006    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    10431007    !--added to the total water content of the gridbox
    1044     IF ( icesed_flux(i) .GT. 0. ) THEN
    1045       qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
     1008    IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
     1009      qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) / dhum_to_dflux(i)
    10461010      !--Vapor is updated after evaporation/sublimation (it is increased)
    10471011      qvap(i) = qvap(i) + qice_sedim
    … …  
    14081372SUBROUTINE poprecip_postcld( &
    14091373           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
    1410            temp, qvap, qliq, qice, icefrac, cldfra, icesed_flux, &
     1374           temp, qvap, qliq, qice, icefrac, cldfra, &
    14111375           precipfracclr, precipfraccld, &
    14121376           rain, rainclr, raincld, snow, snowclr, snowcld, &
    … …  
    14641428REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
    14651429
    1466 REAL,    INTENT(INOUT), DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    14671430REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
    14681431REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
    … …  
    15281491
    15291492qzero(:)    = 0.
    1530 ! AB ICE SEDIM
    1531 !icesed_flux(:) = 0.
    15321493
    15331494dqrcol(:)   = 0.
    … …  
    20001961
    20011962
    2002   !---------------------------------------------------------
    2003   !--                  ICE SEDIMENTATION
    2004   !---------------------------------------------------------
    2005   !IF ( ok_ice_sedim ) THEN
    2006   IF ( .FALSE. ) THEN ! AB ICE SEDIM
    2007 
    2008     IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. eps ) ) THEN
    2009 
    2010       rho = pplay(i) / temp(i) / RD
    2011       iwcg = MAX( rho * qice(i) / cldfra(i) * 1000., eps ) ! iwc in g/m3
    2012       !tempc = temp(i) - RTT ! celcius temp
    2013       !! so-called 'empirical parameterization' in Stubenrauch et al. 2019
    2014       !IF ( tempc .GE. -60. ) THEN
    2015       !  icevelo = EXP( 1.9714 + 0.00216078 * tempc &
    2016       !          - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
    2017       !ELSE
    2018       !  icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
    2019       !ENDIF
    2020       !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
    2021 
    2022       icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16
    2023 
    2024       dz = ( paprsdn(i) - paprsup(i) ) / RG / rho
    2025       qice_sedim = qice(i) * MIN(1., icevelo * dtime / dz)
    2026       !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
    2027       !--times in the same timestep)
    2028       qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
    2029       !--Add tendencies
    2030       qice(i) = qice(i) - qice_sedim
    2031       !--Convert to flux
    2032       icesed_flux(i) = qice_sedim * dhum_to_dflux(i)
    2033 
    2034       !--Diagnostic tendencies
    2035       dqised(i) = dqised(i) - qice_sedim / dtime
    2036     ENDIF
    2037   ENDIF
    2038 
    2039 
    20401963  !--If the local flux of rain+snow in clear air is lower than min_precip,
    20411964  !--we reduce the precipiration fraction in the clear air so that the new
Note: See TracChangeset for help on using the changeset viewer.