Ignore:
Timestamp:
Dec 17, 2024, 9:48:07 AM (43 hours ago)
Author:
evignon
Message:

fin de l'externalisation des precips dans lscp.

  1. Borella, E Vignon
File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp_precip.f90

    r5411 r5413  
    295295    ELSE
    296296      ! if no precip, we reinitialize the cloud fraction used for the precip to 0
    297       ! AB note that this assignment is useless, as znebprecip is not re-used
    298297      znebprecip(i)=0.
    299298
    … …  
    305304
    306305END SUBROUTINE histprecip_precld
     306
     307
     308SUBROUTINE histprecip_postcld( &
     309        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, &
     310        zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
     311        rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
     312        zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
     313        zradocond, zradoice, dqrauto, dqsauto &
     314        )
     315
     316USE lmdz_lscp_ini, ONLY : RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT
     317USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, ffallv_con, &
     318                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, ffallv_lsc, &
     319                          seuil_neb, rain_int_min, cice_velo, dice_velo
     320USE lmdz_lscp_ini, ONLY : iflag_evap_prec, iflag_cloudth_vert, iflag_rain_incloud_vol, &
     321                          iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
     322                          niter_lscp
     323USE lmdz_lscp_tools, ONLY : fallice_velocity
     324
     325IMPLICIT NONE
     326
     327
     328INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
     329REAL,    INTENT(IN)                     :: dtime          !--time step [s]
     330LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
     331
     332REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
     333REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
     334REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
     335REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
     336LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
     337REAL,    INTENT(IN),    DIMENSION(klon) :: zdqsdT_raw     !--derivative of qsat w.r.t. temperature times L/Cp [SI]
     338
     339REAL,    INTENT(INOUT), DIMENSION(klon) :: zt             !--current temperature [K]
     340REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity [kg/kg]
     341REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliq          !--current liquid+ice water specific humidity [kg/kg]
     342REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliql         !--current liquid water specific humidity [kg/kg]
     343REAL,    INTENT(INOUT), DIMENSION(klon) :: zoliqi         !--current ice water specific humidity [kg/kg]
     344REAL,    INTENT(INOUT), DIMENSION(klon) :: zcond          !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg]
     345REAL,    INTENT(IN),    DIMENSION(klon) :: zfice          !--ice fraction AFTER CONDENSATION [-]
     346REAL,    INTENT(IN),    DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
     347
     348REAL,    INTENT(IN),    DIMENSION(klon) :: rneb           !--cloud fraction [-]
     349REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
     350REAL,    INTENT(INOUT), DIMENSION(klon) :: znebprecipcld  !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-]
     351REAL,    INTENT(INOUT), DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
     352REAL,    INTENT(INOUT), DIMENSION(klon) :: tot_zneb       !--total cloud cover above the current mesh [-]
     353REAL,    INTENT(INOUT), DIMENSION(klon) :: zrho_up        !--air density IN THE LAYER ABOVE [kg/m3]
     354REAL,    INTENT(INOUT), DIMENSION(klon) :: zvelo_up       !--ice fallspeed velocity IN THE LAYER ABOVE [m/s]
     355
     356REAL,    INTENT(INOUT), DIMENSION(klon) :: zrfl           !--flux of rain gridbox-mean [kg/s/m2]
     357REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflclr        !--flux of rain gridbox-mean in clear sky [kg/s/m2]
     358REAL,    INTENT(INOUT), DIMENSION(klon) :: zrflcld        !--flux of rain gridbox-mean in cloudy air [kg/s/m2]
     359REAL,    INTENT(INOUT), DIMENSION(klon) :: zifl           !--flux of snow gridbox-mean [kg/s/m2]
     360REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
     361REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
     362
     363REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
     364REAL,    INTENT(OUT),   DIMENSION(klon) :: zradoice       !--condensed ice water used in the radiation scheme [kg/kg]
     365REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
     366REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
     367
     368
     369! Local variables for precip fraction update
     370REAL :: smallestreal
     371REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb
     372REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld
     373REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
     374REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
     375
     376! Local variables for autoconversion
     377REAL :: zct, zcl, zexpo, ffallv
     378REAL :: zchau, zfroi
     379REAL :: zrain, zsnow, zprecip
     380REAL :: effective_zneb
     381REAL, DIMENSION(klon) :: zrho, zvelo
     382REAL, DIMENSION(klon) :: zdz, iwc
     383
     384! Local variables for Bergeron process
     385REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
     386REAL, DIMENSION(klon) :: zqpreci, zqprecl
     387
     388! Local variables for calculation of quantity of condensates seen by the radiative scheme
     389REAL, DIMENSION(klon) :: zradoliq
     390REAL, DIMENSION(klon) :: ziflprev
     391
     392! Misc
     393INTEGER :: i, n
     394
     395! Initialisation
     396smallestreal=1.e-9
     397zqprecl(:) = 0.
     398zqpreci(:) = 0.
     399ziflprev(:) = 0.
     400
     401
     402IF (iflag_evap_prec .GE. 4) THEN
     403
     404  !Partitionning between precipitation coming from clouds and that coming from CS
     405
     406  !0) Calculate tot_zneb, total cloud fraction above the cloud
     407  !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
     408 
     409  DO i=1, klon
     410    tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i),zneb(i))) &
     411            /(1.-min(zneb(i),1.-smallestreal))
     412    d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
     413    tot_zneb(i) = tot_znebn(i)
     414
     415
     416    !1) Cloudy to clear air
     417    d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i),znebprecipcld(i))
     418    IF (znebprecipcld(i) .GT. 0.) THEN
     419      d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
     420      d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
     421    ELSE
     422      d_zrfl_cld_clr(i) = 0.
     423      d_zifl_cld_clr(i) = 0.
     424    ENDIF
     425
     426    !2) Clear to cloudy air
     427    d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i) - d_tot_zneb(i) - zneb(i)))
     428    IF (znebprecipclr(i) .GT. 0) THEN
     429      d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
     430      d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
     431    ELSE
     432      d_zrfl_clr_cld(i) = 0.
     433      d_zifl_clr_cld(i) = 0.
     434    ENDIF
     435
     436    !Update variables
     437    znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
     438    znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
     439    zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
     440    ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
     441    zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
     442    ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
     443
     444    ! c_iso  do the same thing for isotopes precip
     445  ENDDO
     446ENDIF
     447
     448
     449! Autoconversion
     450! -------------------------------------------------------------------------------
     451
     452! Initialisation
     453DO i = 1, klon
     454  zrho(i)  = pplay(i) / zt(i) / RD
     455  zdz(i)   = (paprsdn(i)-paprsup(i)) / (zrho(i)*RG)
     456  iwc(i)   = 0.
     457  zneb(i)  = MAX(rneb(i),seuil_neb)
     458
     459  ! variables for calculation of quantity of condensates seen by the radiative scheme
     460  ! NB. only zradocond and zradoice are outputed, but zradoliq is used if ok_radocond_snow
     461  ! is TRUE
     462  zradocond(i) = zoliq(i)/REAL(niter_lscp+1)
     463  zradoliq(i)  = zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
     464  zradoice(i)  = zoliq(i)*zfice(i)/REAL(niter_lscp+1)
     465ENDDO
     466
     467   
     468DO n = 1, niter_lscp
     469
     470  DO i=1,klon
     471    IF (rneb(i).GT.0.0) THEN
     472      iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
     473    ENDIF
     474  ENDDO
     475
     476  CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo)
     477
     478  DO i = 1, klon
     479
     480    IF (rneb(i).GT.0.0) THEN
     481
     482      ! Initialization of zrain, zsnow and zprecip:
     483      zrain=0.
     484      zsnow=0.
     485      zprecip=0.
     486      ! c_iso same init for isotopes. Externalisation?
     487
     488      IF (zneb(i).EQ.seuil_neb) THEN
     489        zprecip = 0.0
     490        zsnow = 0.0
     491        zrain= 0.0
     492      ELSE
     493
     494        IF (ptconv(i)) THEN ! if convective point
     495          zcl=cld_lc_con
     496          zct=1./cld_tau_con
     497          zexpo=cld_expo_con
     498          ffallv=ffallv_con
     499        ELSE
     500          zcl=cld_lc_lsc
     501          zct=1./cld_tau_lsc
     502          zexpo=cld_expo_lsc
     503          ffallv=ffallv_lsc
     504        ENDIF
     505
     506
     507        ! if vertical heterogeneity is taken into account, we use
     508        ! the "true" volume fraction instead of a modified
     509        ! surface fraction (which is larger and artificially
     510        ! reduces the in-cloud water).
     511        IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
     512          effective_zneb=ctot_vol(i)
     513        ELSE
     514          effective_zneb=zneb(i)
     515        ENDIF
     516
     517
     518        ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
     519        ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
     520        !.........................................................
     521        IF (iflag_autoconversion .EQ. 2) THEN
     522        ! two-steps resolution with niter_lscp=1 sufficient
     523        ! we first treat the second term (with exponential) in an explicit way
     524        ! and then treat the first term (-q/tau) in an exact way
     525          zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
     526            *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
     527        ELSE
     528        ! old explicit resolution with subtimesteps
     529          zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
     530            *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
     531        ENDIF
     532
     533
     534        ! Ice water quantity to remove (Zender & Kiehl, 1997)
     535        ! dqice/dt=1/rho*d(rho*wice*qice)/dz
     536        !....................................
     537        IF (iflag_autoconversion .EQ. 2) THEN
     538        ! exact resolution, niter_lscp=1 is sufficient but works only
     539        ! with iflag_vice=0
     540          IF (zoliqi(i) .GT. 0.) THEN
     541            zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
     542              +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
     543          ELSE
     544            zfroi=0.
     545          ENDIF
     546        ELSE
     547        ! old explicit resolution with subtimesteps
     548          zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*zvelo(i)
     549        ENDIF
     550
     551        zrain   = MIN(MAX(zchau,0.0),zoliql(i))
     552        zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
     553        zprecip = MAX(zrain + zsnow,0.0)
     554
     555      ENDIF
     556
     557
     558      IF (iflag_autoconversion .GE. 1) THEN
     559        ! debugged version with phase conservation through the autoconversion process
     560        zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
     561        zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
     562        zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
     563      ELSE
     564        ! bugged version with phase resetting
     565        zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
     566        zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
     567        zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
     568      ENDIF
     569
     570      ! c_iso: call isotope_conversion (for readibility) or duplicate
     571
     572      ! variables for calculation of quantity of condensates seen by the radiative scheme
     573      zradocond(i) = zradocond(i) + zoliq(i)/REAL(niter_lscp+1)
     574      zradoliq(i)  = zradoliq(i) + zoliql(i)/REAL(niter_lscp+1)
     575      zradoice(i)  = zradoice(i) + zoliqi(i)/REAL(niter_lscp+1)
     576
     577      !--Diagnostics
     578      dqrauto(i) = dqrauto(i) + zrain / dtime
     579      dqsauto(i) = dqsauto(i) + zsnow / dtime
     580
     581    ENDIF ! rneb >0
     582
     583  ENDDO  ! i = 1,klon
     584
     585ENDDO      ! n = 1,niter
     586
     587! Precipitation flux calculation
     588
     589DO i = 1, klon
     590       
     591  IF (iflag_evap_prec.GE.4) THEN
     592    ziflprev(i)=ziflcld(i)
     593  ELSE
     594    ziflprev(i)=zifl(i)*zneb(i)
     595  ENDIF
     596
     597  IF (rneb(i) .GT. 0.0) THEN
     598
     599    ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
     600    ! If T<0C, liquid precip are converted into ice, which leads to an increase in
     601    ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
     602    ! taken into account through a linearization of the equations and by approximating
     603    ! the liquid precipitation process with a "threshold" process. We assume that
     604    ! condensates are not modified during this operation. Liquid precipitation is
     605    ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
     606    ! Water vapor increases as well
     607    ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
     608
     609    IF (ok_bug_phase_lscp) THEN
     610      zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
     611      zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
     612    ELSE
     613      zqpreci(i)=zcond(i)*zfice(i)-zoliqi(i)
     614      zqprecl(i)=zcond(i)*(1.-zfice(i))-zoliql(i)
     615    ENDIF
     616    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     617    coef1 = rneb(i)*RLSTT/zcp*zdqsdT_raw(i)
     618    ! Computation of DT if all the liquid precip freezes
     619    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
     620    ! T should not exceed the freezing point
     621    ! that is Delta > RTT-zt(i)
     622    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
     623    zt(i)  = zt(i) + DeltaT
     624    ! water vaporization due to temp. increase
     625    Deltaq = rneb(i)*zdqsdT_raw(i)*DeltaT
     626    ! we add this vaporized water to the total vapor and we remove it from the precip
     627    zq(i) = zq(i) +  Deltaq
     628    ! The three "max" lines herebelow protect from rounding errors
     629    zcond(i) = max( zcond(i)- Deltaq, 0. )
     630    ! liquid precipitation converted to ice precip
     631    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
     632    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
     633    ! iced water budget
     634    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
     635
     636    ! c_iso : mv here condensation of isotopes + redispatchage en precip
     637
     638    IF (iflag_evap_prec.GE.4) THEN
     639      zrflcld(i) = zrflcld(i)+zqprecl(i) &
     640        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     641      ziflcld(i) = ziflcld(i)+ zqpreci(i) &
     642        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     643      znebprecipcld(i) = rneb(i)
     644      zrfl(i) = zrflcld(i) + zrflclr(i)
     645      zifl(i) = ziflcld(i) + ziflclr(i)
     646    ELSE
     647      zrfl(i) = zrfl(i)+ zqprecl(i)        &
     648        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     649      zifl(i) = zifl(i)+ zqpreci(i)        &
     650        *(paprsdn(i)-paprsup(i))/(RG*dtime)
     651    ENDIF
     652    ! c_iso : same for isotopes
     653
     654  ENDIF ! rneb>0
     655
     656ENDDO
     657
     658! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
     659! if iflag_evap_prec>=4
     660IF (iflag_evap_prec.GE.4) THEN
     661
     662  DO i=1,klon
     663
     664    IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
     665      znebprecipclr(i) = min(znebprecipclr(i),max( &
     666        zrflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), &
     667        ziflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
     668    ELSE
     669      znebprecipclr(i)=0.0
     670    ENDIF
     671
     672    IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
     673      znebprecipcld(i) = min(znebprecipcld(i), max( &
     674        zrflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), &
     675        ziflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
     676    ELSE
     677      znebprecipcld(i)=0.0
     678    ENDIF
     679  ENDDO
     680
     681ENDIF
     682
     683! we recalculate zradoice to account for contributions from the precipitation flux
     684! if ok_radocond_snow is true
     685IF ( ok_radocond_snow ) THEN
     686  IF ( .NOT. iftop ) THEN
     687    ! for the solid phase (crystals + snowflakes)
     688    ! we recalculate zradoice by summing
     689    ! the ice content calculated in the mesh
     690    ! + the contribution of the non-evaporated snowfall
     691    ! from the overlying layer
     692    DO i=1,klon
     693      IF (ziflprev(i) .NE. 0.0) THEN
     694        zradoice(i)=zoliqi(i)+zqpreci(i)+ziflprev(i)/zrho_up(i)/zvelo_up(i)
     695      ELSE
     696        zradoice(i)=zoliqi(i)+zqpreci(i)
     697      ENDIF
     698      zradocond(i)=zradoliq(i)+zradoice(i)
     699    ENDDO
     700  ENDIF
     701  ! keep in memory air density and ice fallspeed velocity
     702  zrho_up(:) = zrho(:)
     703  zvelo_up(:) = zvelo(:)
     704ENDIF
     705
     706END SUBROUTINE histprecip_postcld
    307707
    308708
Note: See TracChangeset for help on using the changeset viewer.