Changeset 4226 for LMDZ6/trunk


Ignore:
Timestamp:
Jul 23, 2022, 6:04:45 PM (2 years ago)
Author:
oboucher
Message:

Cleaning up the routine

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/phylmd/lscp_mod.F90

    r4114 r4226  
    2323!--------------------------------------------------------------------------------
    2424! Aim: Large Scale Clouds and Precipitation (LSCP)
    25 !     
     25!
    2626! This code is a new version of the fisrtilp.F90 routine, which is itself a
    2727! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
     
    4646! Code structure:
    4747!
    48 ! P0>     Thermalization of the  precipitation coming from the overlying layer
     48! P0>     Thermalization of the precipitation coming from the overlying layer
    4949! P1>     Evaporation of the precipitation (falling from the k+1 level)
    5050! P2>     Cloud formation (at the k level)
     
    8686! (which is not trivial)
    8787!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    88  
     88!
    8989USE dimphy
    9090USE print_control_mod, ONLY: prt_level, lunout
     
    112112include "nuage.h"
    113113
    114  
    115114! INPUT VARIABLES:
    116115!-----------------
     
    125124  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
    126125  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
    127                                                                    ! CR: if iflag_ice_thermo=2, only convection is active   
     126                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
    128127  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated
    129128
    130129  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
    131130
    132  
     131
    133132  !Inputs associated with thermal plumes
    134  
     133
    135134  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv             ! virtual potential temperature [K]
    136135  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta            ! specific humidity within thermals [kg/kg]
     
    144143
    145144  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs            ! function of pressure that sets the large-scale
    146                                                                     ! cloud PDF (sigma=ratqs*qt)
    147  
    148145
    149146  ! Input sursaturation en glace
     
    166163  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
    167164
    168  
    169165  ! cofficients of scavenging fraction (for off-line)
    170166 
     
    178174  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl        ! scavenging fraction due tu nucleation [-]
    179175 
    180 
    181 
    182 
    183176  ! PROGRAM PARAMETERS:
    184177  !--------------------
    185178 
    186  
    187   REAL, SAVE :: seuil_neb=0.001                                     ! cloud fraction threshold: a cloud really exists when exceeded
     179  REAL, SAVE :: seuil_neb=0.001                 ! cloud fraction threshold: a cloud really exists when exceeded
    188180  !$OMP THREADPRIVATE(seuil_neb)
    189181
    190   INTEGER, SAVE :: ninter=5                                         ! number of iterations to calculate autoconversion to precipitation
     182  INTEGER, SAVE :: ninter=5                     ! number of iterations to calculate autoconversion to precipitation
    191183  !$OMP THREADPRIVATE(ninter)
    192184
    193   INTEGER,SAVE :: iflag_evap_prec=1                                 ! precipitation evap. flag. 0: nothing, 1: "old way", 2: Max cloud fraction above to calculate the max of reevaporation
    194                                                                     ! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
     185  INTEGER,SAVE :: iflag_evap_prec=1             ! precipitation evaporation flag. 0: nothing, 1: "old way",
     186                                                ! 2: Max cloud fraction above to calculate the max of reevaporation
     187                                                ! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
    195188  !$OMP THREADPRIVATE(iflag_evap_prec)
    196                                                                    
    197 
    198   REAL t_coup                                                       ! temperature threshold which determines the phase
    199   PARAMETER (t_coup=234.0)                                          ! for which the saturation vapor pressure is calculated
    200 
    201   REAL DDT0                                                         ! iteration parameter
     189
     190  REAL t_coup                                   ! temperature threshold which determines the phase
     191  PARAMETER (t_coup=234.0)                      ! for which the saturation vapor pressure is calculated
     192
     193  REAL DDT0                                     ! iteration parameter
    202194  PARAMETER (DDT0=.01)
    203195
    204   REAL ztfondue                                                     ! parameter to calculate melting fraction of precipitation
     196  REAL ztfondue                                 ! parameter to calculate melting fraction of precipitation
    205197  PARAMETER (ztfondue=278.15)
    206198
    207   REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction
     199  REAL, SAVE    :: rain_int_min=0.001           ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction
    208200  !$OMP THREADPRIVATE(rain_int_min)
    209 
    210201
    211202
     
    218209  REAL ctot(klon,klev)
    219210  REAL ctot_vol(klon,klev)
    220   REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
     211  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
    221212  REAL zdqsdT_raw(klon)
    222   REAL gammasat(klon),dgammasatdt(klon)                                 ! coefficient to make cold condensation at the correct RH and derivative wrt T
     213  REAL gammasat(klon),dgammasatdt(klon)                ! coefficient to make cold condensation at the correct RH and derivative wrt T
    223214  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
    224215  REAL cste
    225216  REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    226217  REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
    227   REAL erf   
     218  REAL erf
    228219  REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon)
    229220  REAL zrfl(klon), zrfln(klon), zqev, zqevt
     
    248239  ! temporary: alpha parameter for scavenging
    249240  ! 4 possible scavenging processes
    250   REAL a_tr_sca(4)
    251   save a_tr_sca
    252   !$OMP THREADPRIVATE(a_tr_sca) 
     241  REAL,SAVE :: a_tr_sca(4)
     242  !$OMP THREADPRIVATE(a_tr_sca)
    253243  REAL zalpha_tr
    254244  REAL zfrac_lessi
     
    257247  REAL zmair(klon), zcpair, zcpeau
    258248  REAL zlh_solid(klon), zm_solid         ! for liquid -> solid conversion
    259   REAL zrflclr(klon), zrflcld(klon)                                                                   
    260   REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)                                                   
     249  REAL zrflclr(klon), zrflcld(klon)
     250  REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
    261251  REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon)
    262   REAL ziflclr(klon), ziflcld(klon)                                                                   
    263   REAL znebprecipclr(klon), znebprecipcld(klon)                                                       
    264   REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)                                             
     252  REAL ziflclr(klon), ziflcld(klon)
     253  REAL znebprecipclr(klon), znebprecipcld(klon)
     254  REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
    265255  REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
    266256  REAL velo(klon,klev), vr
     
    270260  INTEGER i, k, n, kk
    271261  INTEGER n_i(klon), iter
    272   INTEGER ncoreczq 
     262  INTEGER ncoreczq
    273263  INTEGER mpc_bl_points(klon,klev)
    274264  INTEGER,SAVE :: itap=0
     
    277267  LOGICAL lognormale(klon)
    278268  LOGICAL convergence(klon)
    279   LOGICAL appel1er
    280   SAVE appel1er
     269  LOGICAL,SAVE :: appel1er
    281270  !$OMP THREADPRIVATE(appel1er)
    282271
     
    290279!===============================================================================
    291280
    292 ! Few initial checksS
     281! Few initial checks
    293282
    294283IF (iflag_t_glace.EQ.0) THEN
     
    303292
    304293IF (.NOT.thermcep) THEN
    305      abort_message = 'lscp cannot be used when thermcep=false'
    306      CALL abort_physic(modname,abort_message,1)
     294    abort_message = 'lscp cannot be used when thermcep=false'
     295    CALL abort_physic(modname,abort_message,1)
    307296ENDIF
    308297
    309298IF (iflag_fisrtilp_qsat .LT. 0) THEN
    310         abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
    311      CALL abort_physic(modname,abort_message,1)
     299    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
     300    CALL abort_physic(modname,abort_message,1)
    312301ENDIF
    313302
    314 
    315303! Few initialisations
    316304
    317305itap=itap+1
    318 znebprecip(:)=0.
     306znebprecip(:)=0.0
    319307ctot_vol(1:klon,1:klev)=0.0
    320308rneblsvol(1:klon,1:klev)=0.0
    321 smallestreal=1.e-9                                                                                 
    322 znebprecipclr(:)=0.   
    323 znebprecipcld(:)=0.                                                                           
     309smallestreal=1.e-9
     310znebprecipclr(:)=0.0
     311znebprecipcld(:)=0.0
    324312mpc_bl_points(:,:)=0
    325313
     
    341329        WRITE(lunout,*) 'I would prefer a 6 min sub-timestep'
    342330    ENDIF
    343      
     331
    344332    !AA Temporary initialisation
    345333    a_tr_sca(1) = -0.5
     
    351339    DO k = 1, klev
    352340        DO i = 1, klon
    353             pfrac_nucl(i,k)=1.
    354             pfrac_1nucl(i,k)=1.
    355             pfrac_impa(i,k)=1.
    356             beta(i,k)=0.       
     341            pfrac_nucl(i,k)=1.0
     342            pfrac_1nucl(i,k)=1.0
     343            pfrac_impa(i,k)=1.0
     344            beta(i,k)=0.0
    357345        ENDDO
    358346    ENDDO
     
    362350    appel1er = .FALSE.
    363351
    364 ENDIF     !(appel1er) 
     352ENDIF     !(appel1er)
    365353
    366354! AA for 'safety' reasons
     
    379367radliq(:,:) = 0.0
    380368radicefrac(:,:) = 0.0
    381 frac_nucl(:,:) = 1.
    382 frac_impa(:,:) = 1.
     369frac_nucl(:,:) = 1.0
     370frac_impa(:,:) = 1.0
    383371
    384372rain(:) = 0.0
     
    393381ziflprev(:)=0.0
    394382zneb(:) = seuil_neb
    395 zrflclr(:) = 0.0                                                                           
    396 ziflclr(:) = 0.0                                                                           
    397 zrflcld(:) = 0.0                                                                               
    398 ziflcld(:) = 0.0                                                                                 
    399 tot_zneb(:) = 0.0                                                                                 
    400 tot_znebn(:) = 0.0                                                                               
    401 d_tot_zneb(:) = 0.0   
    402 qzero(:)=0.0
     383zrflclr(:) = 0.0
     384ziflclr(:) = 0.0
     385zrflcld(:) = 0.0
     386ziflcld(:) = 0.0
     387tot_zneb(:) = 0.0
     388tot_znebn(:) = 0.0
     389d_tot_zneb(:) = 0.0
     390qzero(:) = 0.0
    403391
    404392!--ice sursaturation
    405 gamma_ss(:,:) = 1.
    406 qss(:,:) = 0.
    407 rnebss(:,:) = 0.
     393gamma_ss(:,:) = 1.0
     394qss(:,:) = 0.0
     395rnebss(:,:) = 0.0
    408396Tcontr(:,:) = missing_val
    409397qcontr(:,:) = missing_val
     
    419407
    420408DO k = klev, 1, -1
    421 
    422409
    423410    ! Initialisation temperature and specific humidity
     
    447434            ! RVTMP2=rcpv/rcpd-1
    448435            zcpair=RCPD*(1.0+RVTMP2*zq(i))
    449             zcpeau=RCPD*RVTMP2             
    450                                          
     436            zcpeau=RCPD*RVTMP2
     437
    451438            ! zmqc: precipitation mass that has to be thermalized with
    452439            ! layer's air so that precipitation at the ground has the
     
    456443            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
    457444             / (zcpair + zmqc(i)*zcpeau)
    458            
     445
    459446        ENDDO
    460447 
    461     ELSE 
     448    ELSE
    462449 
    463450        DO i = 1, klon
     
    474461    ! --------------------------------------------------------------------
    475462
    476 
    477463    IF (iflag_evap_prec.GE.1) THEN
    478464
     
    486472
    487473        DO i = 1, klon
    488      
     474
    489475            ! if precipitation
    490476            IF (zrfl(i)+zifl(i).GT.0.) THEN
    491                    
     477
    492478            ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
    493                 IF (iflag_evap_prec.EQ.4) THEN                                                         
    494                     zrfl(i) = zrflclr(i)                                                         
    495                     zifl(i) = ziflclr(i)                                                           
    496                 ENDIF                     
     479                IF (iflag_evap_prec.EQ.4) THEN
     480                    zrfl(i) = zrflclr(i)
     481                    zifl(i) = ziflclr(i)
     482                ENDIF
    497483
    498484                IF (iflag_evap_prec.EQ.1) THEN
     
    501487                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
    502488                ENDIF
    503      
    504                                                                                                                                                                                                                  
    505                 IF (iflag_evap_prec.EQ.4) THEN                                                                 
    506                 ! Max evaporation not to saturate the whole mesh                                             
    507                     zqev0 = MAX(0.0, zqs(i)-zq(i))                                                             
    508                 ELSE       
    509                 ! Max evap not to saturate the fraction below the cloud                                                   
    510                     zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))                                                                                 
    511                 ENDIF   
    512 
     489
     490                IF (iflag_evap_prec.EQ.4) THEN
     491                ! Max evaporation not to saturate the whole mesh
     492                    zqev0 = MAX(0.0, zqs(i)-zq(i))
     493                ELSE
     494                ! Max evap not to saturate the fraction below the cloud
     495                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
     496                ENDIF
    513497
    514498                ! Evaporation of liquid precipitation coming from above
     
    530514                ENDIF
    531515
    532 
    533516                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
    534                 *RG*dtime/(paprs(i,k)-paprs(i,k+1))
    535          
     517                    *RG*dtime/(paprs(i,k)-paprs(i,k+1))
     518
    536519                ! sublimation of the solid precipitation coming from above
    537520                IF (iflag_evap_prec.EQ.3) THEN
     
    547530                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    548531                ENDIF
    549            
     532
    550533                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
    551                 *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
     534                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
    552535
    553536                ! A. JAM
     
    556539                ! does not reach saturation. In this case, we
    557540                ! redistribute zqev0 conserving the ratio liquid/ice
    558    
     541
    559542                IF (zqevt+zqevti.GT.zqev0) THEN
    560543                    zqev=zqev0*zqevt/(zqevt+zqevti)
     
    589572                zifl(i) = zifln(i)
    590573
    591                 IF (iflag_evap_prec.EQ.4) THEN                                                                 
    592                     zrflclr(i) = zrfl(i)                                                                 
    593                     ziflclr(i) = zifl(i)   
    594                     IF(zrflclr(i) + ziflclr(i).LE.0) THEN                                               
    595                         znebprecipclr(i) = 0.                                                         
    596                     ENDIF                                                                               
    597                     zrfl(i) = zrflclr(i) + zrflcld(i)                                                     
    598                     zifl(i) = ziflclr(i) + ziflcld(i)                                                     
    599                 ENDIF     
     574                IF (iflag_evap_prec.EQ.4) THEN
     575                    zrflclr(i) = zrfl(i)
     576                    ziflclr(i) = zifl(i)
     577                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
     578                        znebprecipclr(i) = 0.0
     579                    ENDIF
     580                    zrfl(i) = zrflclr(i) + zrflcld(i)
     581                    zifl(i) = ziflclr(i) + ziflcld(i)
     582                ENDIF
    600583
    601584
     
    605588                zmelt = MIN(MAX(zmelt,0.),1.)
    606589                ! Ice melting
    607                 IF (iflag_evap_prec.EQ.4) THEN                                                             
    608                     zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)                                           
    609                     zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)                                             
    610                     zrfl(i)=zrflclr(i)+zrflcld(i)                                                       
    611                 ELSE                                                                                       
    612                     zrfl(i)=zrfl(i)+zmelt*zifl(i)                                                     
    613                 ENDIF                 
     590                IF (iflag_evap_prec.EQ.4) THEN
     591                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
     592                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
     593                    zrfl(i)=zrflclr(i)+zrflcld(i)
     594                ELSE
     595                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
     596                ENDIF
    614597
    615598                ! Latent heat of melting with precipitation thermalization
     
    617600                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
    618601               
    619                 IF (iflag_evap_prec.EQ.4) THEN                                                             
    620                     ziflclr(i)=ziflclr(i)*(1.-zmelt)                                                   
    621                     ziflcld(i)=ziflcld(i)*(1.-zmelt)                                                   
    622                     zifl(i)=ziflclr(i)+ziflcld(i)                                                       
    623                 ELSE                                                                                     
    624                     zifl(i)=zifl(i)*(1.-zmelt)                                                         
     602                IF (iflag_evap_prec.EQ.4) THEN
     603                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
     604                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
     605                    zifl(i)=ziflclr(i)+ziflcld(i)
     606                ELSE
     607                    zifl(i)=zifl(i)*(1.-zmelt)
    625608                ENDIF
    626609
     
    632615
    633616        ENDDO
    634          
     617
    635618    ENDIF ! (iflag_evap_prec>=1)
    636619
     
    638621    ! End precip evaporation
    639622    ! --------------------------------------------------------------------
    640 
    641 
    642 
    643623   
    644624    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
    645625    !-------------------------------------------------------
    646      
     626
    647627     qtot(:)=zq(:)+zmqc(:)
    648628     CALL CALC_QSAT_ECMWF(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
     
    657637     ENDDO
    658638   
    659 
    660639    ! --------------------------------------------------------------------
    661640    ! P2> Cloud formation
     
    670649    ! -------------------------------------------------------------------
    671650
    672            
    673651        ! P2.1> With the PDFs (log-normal, bigaussian)
    674652        ! cloud propertues calculation with the initial values of t and q
     
    697675                         ratqs,zqs,t)
    698676
    699 
    700677                !Jean Jouhaud's version, Decembre 2018
    701                 !-------------------------------------             
    702              
     678                !-------------------------------------
     679
    703680                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
    704681
     
    740717
    741718        ELSE
    742                 ! When iflag_cld_th=5, we always assume 
     719                ! When iflag_cld_th=5, we always assume
    743720                ! bi-gaussian distribution
    744721            lognormale = .FALSE.
     
    746723        ENDIF
    747724
    748 
    749725        DT(:) = 0.
    750726        n_i(:)=0
     
    755731        ! condensation with qsat(T) variation (adaptation)
    756732        ! Iterative Loop to converge towards qsat
    757        
    758        
     733
    759734        DO iter=1,iflag_fisrtilp_qsat+1
    760                
     735
    761736                DO i=1,klon
    762737
    763738                    ! convergence = .true. until when convergence is not satisfied
    764739                    convergence(i)=ABS(DT(i)).GT.DDT0
    765                    
     740
    766741                    IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
    767                          
     742
    768743                        ! if not convergence:
    769744
     
    781756                        ! Rneb, qzn and zcond for lognormal PDFs
    782757                        qtot(i)=zq(i)+zmqc(i)
    783                    
     758
    784759                      ENDIF
    785760
     
    787762
    788763                  ! Calculation of saturation specific humidity and ce fraction
    789                   CALL CALC_QSAT_ECMWF(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
    790                   CALL CALC_GAMMASAT(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
     764                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
     765                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
    791766                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
    792767                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
    793                   CALL ICEFRAC_LSCP(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
    794 
    795 
     768                  CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
    796769
    797770                  DO i=1,klon
     
    810783                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
    811784                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
    812              
     785
    813786                        !--ice sursaturation by Audran
    814787                        IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN
     
    826799                          fcontrP(i,k)=0.0  !--idem
    827800                          qss(i,k)=0.0      !--idem
    828                        
     801
    829802                        ELSE
    830                
     803
    831804                        !------------------------------------
    832805                        ! ICE SUPERSATURATION
     
    834807
    835808                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
    836                              gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),               
    837                              rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                 &
     809                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
     810                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   &
    838811                             Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
    839812
     
    862835                        qlbef(i)=max(0.,zqn(i)-zqs(i))
    863836                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
    864                         &   +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
     837                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
    865838                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
    866839                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
    867                         &     *qlbef(i)*dzfice(i)
     840                              *qlbef(i)*dzfice(i)
    868841                        DT(i)=num/denom
    869842                        n_i(i)=n_i(i)+1
     
    874847
    875848        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
    876 
    877 
    878849
    879850        ! P2.3> Final quantities calculation
     
    886857        ! ----------------------------------------------------------------
    887858
    888 
    889 
    890859        ! Water vapor update, Phase determination and subsequent latent heat exchange
    891860
    892861        ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
    893         CALL ICEFRAC_LSCP(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
     862        CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
    894863
    895864        DO i=1, klon
     
    923892                ENDIF
    924893
    925 
    926894                ! water vapor update
    927                 zq(i) = zq(i) - zcond(i)       
    928            
     895                zq(i) = zq(i) - zcond(i)
     896
    929897            ENDIF
    930898
     
    935903        ENDDO
    936904
    937 
    938905        ! If vertical heterogeneity, change volume fraction
    939906        IF (iflag_cloudth_vert .GE. 3) THEN
     
    945912    ! P3> Precipitation formation
    946913    ! ----------------------------------------------------------------
    947    
     914
    948915    ! LTP:
    949916    IF (iflag_evap_prec==4) THEN
     
    972939
    973940                !2) Clear to cloudy air
    974                 d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) &
    975                         - d_tot_zneb(i) - zneb(i)))
     941                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
    976942                IF (znebprecipclr(i) .GT. 0) THEN
    977943                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
     
    983949
    984950                !Update variables
    985                 znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) 
     951                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
    986952                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
    987953                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
     
    1021987        ENDDO
    1022988
    1023         CALL FALLICE_VELOCITY(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
     989        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
    1024990
    1025991        DO i = 1, klon
     
    10681034
    10691035                ENDIF
    1070            
     1036
    10711037                zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
    10721038                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
     
    10761042                radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1)
    10771043                radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
    1078            
     1044
    10791045            ENDIF ! rneb >0
    10801046
     
    10831049    ENDDO      ! n = 1,niter
    10841050
    1085    
    1086 
    10871051    ! Precipitation flux calculation
    10881052
    10891053    DO i = 1, klon
    10901054
    1091 
    10921055            ziflprev(i)=zifl(i)
    10931056
    10941057            IF (rneb(i,k) .GT. 0.0) THEN
    1095                
    10961058
    10971059               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
     
    11261088                    ! iced water budget
    11271089                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
    1128            
    1129                 IF (iflag_evap_prec.EQ.4) THEN                                                           
    1130                     zrflcld(i) = zrflcld(i)+zqprecl(i) &                                                   
    1131                     *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                               
    1132                     ziflcld(i) = ziflcld(i)+ zqpreci(i) &                                                 
    1133                         *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)                                           
    1134                     znebprecipcld(i) = rneb(i,k)                                                         
    1135                     zrfl(i) = zrflcld(i) + zrflclr(i)                                                     
    1136                     zifl(i) = ziflcld(i) + ziflclr(i)                                                   
    1137                 ELSE                                 
     1090
     1091                IF (iflag_evap_prec.EQ.4) THEN
     1092                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
     1093                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
     1094                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
     1095                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
     1096                    znebprecipcld(i) = rneb(i,k)
     1097                    zrfl(i) = zrflcld(i) + zrflclr(i)
     1098                    zifl(i) = ziflcld(i) + ziflclr(i)
     1099                ELSE
    11381100                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
    11391101                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    11401102                    zifl(i) = zifl(i)+ zqpreci(i)        &
    1141                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 
     1103                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    11421104                ENDIF
    11431105 
    11441106           ENDIF ! rneb>0
    11451107
    1146 
    11471108   ENDDO
    11481109
    1149 
    1150 
    1151 
    1152     ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min                                                                             
     1110    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
    11531111    ! if iflag_evap_pre=4
    1154     IF (iflag_evap_prec.EQ.4) THEN 
    1155 
    1156         DO i=1,klon                                       
    1157                
    1158             IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN     
     1112    IF (iflag_evap_prec.EQ.4) THEN
     1113
     1114        DO i=1,klon
     1115
     1116            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
    11591117                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
    11601118                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
    11611119            ELSE
    1162                 znebprecipclr(i)=0.                                                                   
    1163             ENDIF 
    1164                                                                                                          
    1165             IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN                                                 
     1120                znebprecipclr(i)=0.0
     1121            ENDIF
     1122
     1123            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
    11661124                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
    11671125                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
    1168             ELSE                                                                                     
    1169                 znebprecipcld(i)=0.                                                                   
    1170             ENDIF     
    1171 
    1172         ENDDO       
    1173 
    1174     ENDIF 
     1126            ELSE
     1127                znebprecipcld(i)=0.0
     1128            ENDIF
     1129
     1130        ENDDO
     1131
     1132    ENDIF
    11751133
    11761134    ! End of precipitation formation
     
    11941152    ! we recaulate radliqi to account for contributions from the precipitation flux
    11951153
    1196 
    11971154    IF ((ok_radliq_snow) .AND. (k .LT. klev)) THEN
    11981155        ! for the solid phase (crystals + snowflakes)
     
    12051162              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
    12061163           ELSE
    1207               radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i) 
     1164              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
    12081165           ENDIF
    12091166           radliq(i,k)=radliql(i,k)+radliqi(i,k)
     
    12181175    ENDDO
    12191176
    1220 
    1221 
    12221177    ! ----------------------------------------------------------------
    12231178    ! P4> Wet scavenging
    12241179    ! ----------------------------------------------------------------
    1225 
    12261180
    12271181    !Scavenging through nucleation in the layer
     
    12481202            pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
    12491203            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
    1250    
     1204
    12511205            ! Nucleation with a factor of -1 instead of -0.5
    12521206            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
     
    12541208
    12551209        ENDIF
    1256        
    1257     ENDDO     
    1258 
    1259      
     1210
     1211    ENDDO
     1212
    12601213    ! Scavenging through impaction in the underlying layer
    12611214
     
    12811234
    12821235    ENDDO
    1283      
     1236
    12841237    !--save some variables for ice supersaturation
    12851238    !
     
    12941247        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--added by OB because of pathologic cases with the lognormal
    12951248        qcld(i,k) = qvc(i,k) + zcond(i)
    1296    END DO
     1249   ENDDO
    12971250   !q_sat
    1298    CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:))
    1299    CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:))     
    1300 
    1301 END DO
     1251   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:))
     1252   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:))
     1253
     1254ENDDO
    13021255
    13031256!======================================================================
Note: See TracChangeset for help on using the changeset viewer.