Changeset 2807 for LMDZ5


Ignore:
Timestamp:
Mar 6, 2017, 12:10:54 AM (8 years ago)
Author:
jyg
Message:

Various corrections in fisrtilp.F90 yielding an
energy conserving scheme. Energy is conserved
provided: fl_cor_ebil>0, ok_bad_ecmwf_thermo=n
and iflag_bergeron=2

Location:
LMDZ5/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/fisrtilp.F90

    r2703 r2807  
    1717  USE cloudth_mod
    1818  USE ioipsl_getin_p_mod, ONLY : getin_p
     19  USE phys_local_var_mod, ONLY: ql_seri,qs_seri
     20  ! flag to include modifications to ensure energy conservation (if flag >0)
     21  USE add_phys_tend_mod, only : fl_cor_ebil
    1922  IMPLICIT none
    2023  !======================================================================
    … …  
    2629  !             et ilp (il pleut, L. Li)
    2730  ! Principales parties:
     31  ! P0> Thermalisation des precipitations venant de la couche du dessus
    2832  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
    2933  ! P2> Formation du nuage (en k)
     34  ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
     35  ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
     36  ! les valeurs de T et Q initiales
     37  ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
     38  ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
     39  ! P2.B> Nuage "tout ou rien"
     40  ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
    3041  ! P3> Formation de la precipitation (en k)
    3142  !======================================================================
     43  ! JLD:
     44  ! * Routine probablement fausse (au moins incoherente) si thermcep = .false.
     45  ! * fl_cor_ebil doit etre > 0 ;
     46  !   fl_cor_ebil= 0 pour reproduire anciens bugs
    3247  !======================================================================
    3348  include "YOMCST.h"
    … …  
    115130  INTEGER i, k, n, kk
    116131  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 
     132  REAL zdqsdT_raw(klon)
    117133  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
    118134  LOGICAL convergence(klon)
    … …  
    140156!!!!
    141157!  Variables pour Bergeron
    142   REAL zcp, coef1, DeltaT
     158  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
    143159  REAL zqpreci(klon), zqprecl(klon)
     160! Variable pour conservation enegie des precipitations
     161  REAL zmqc(klon)
    144162  !
    145163  LOGICAL appel1er
    … …  
    170188  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
    171189! RomP <<<
    172   REAL zmair, zcpair, zcpeau
     190  REAL zmair(klon), zcpair, zcpeau
    173191  !     Pour la conversion eau-neige
    174192  REAL zlh_solid(klon), zm_solid
    … …  
    180198                     ! (Heymsfield & Donner, 1990)
    181199  REAL zzz
     200
    182201  include "YOETHF.h"
    183202  include "FCTTRE.h"
    … …  
    312331     ENDDO
    313332     !
     333     ! ----------------------------------------------------------------
     334     ! P0> Thermalisation des precipitations venant de la couche du dessus
     335     ! ----------------------------------------------------------------
    314336     ! Calculer la varition de temp. de l'air du a la chaleur sensible
    315      ! transporter par la pluie.
    316      ! Il resterait a rajouter cet effet de la chaleur sensible sur les
    317      ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
    318      ! surface.
     337     ! transporter par la pluie. On thermalise la pluie avec l'air de la couche.
     338     ! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors
     339     ! des differentes transformations thermodynamiques. Cette masse d'eau doit
     340     ! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation
     341     ! de l'enthalpie  de la couche avec la temperature
     342     ! Variables calculees ou modifiees:
     343     !   -  zt: temperature de la cocuhe
     344     !   - zmqc: masse de precip qui doit etre thermalisee
    319345     !
    320346     IF(k.LE.klevm1) THEN         
    321347        DO i = 1, klon
    322348           !IM
    323            zmair=(paprs(i,k)-paprs(i,k+1))/RG
     349           zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
     350           ! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit
    324351           zcpair=RCPD*(1.0+RVTMP2*zq(i))
    325352           zcpeau=RCPD*RVTMP2
     353         if (fl_cor_ebil .GT. 0) then
     354           ! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm
     355           ! pour s'assurer que la precip arrivant au sol aura bien la temperature de la
     356           ! derniere couche
     357           zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
     358           ! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus
     359           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
     360                 / (zcpair + zmqc(i)*zcpeau)
     361         else ! si on maintient les anciennes erreurs
    326362           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
    327                 + zmair*zcpair*zt(i) ) &
    328                 / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
    329            !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
     363                + zmair(i)*zcpair*zt(i) ) &
     364                / (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau)
     365         end if
    330366        ENDDO
    331      ENDIF
    332      ! ----------------------------------------------------------------
    333      ! P1> Debut evaporation de la precipitation
    334      ! ----------------------------------------------------------------
     367     ENDIF ! end IF(k.LE.klevm1)
     368     !
     369     ! ----------------------------------------------------------------
     370     ! P1> Calcul de l'evaporation de la precipitation
     371     ! ----------------------------------------------------------------
     372     ! On evapore une partie des precipitations venant de la maille du dessus.
     373     ! On calcule l'evaporation et la sublimation des precipitations, jusqu'a
     374     ! ce que la fraction de cette couche qui est sous le nuage soit saturee.
     375     ! Variables calculees ou modifiees:
     376     !   - zrfl et zifl: flux de precip liquide et glace
     377     !   - zq, zt: humidite et temperature de la cocuhe
     378     !   - zmqc: masse de precip qui doit etre thermalisee
     379     !
    335380     IF (evap_prec) THEN
    336381        DO i = 1, klon
    337 !AJ<
    338 !!           IF (zrfl(i) .GT.0.) THEN
     382!          S'il y a des precipitations
    339383           IF (zrfl(i)+zifl(i).GT.0.) THEN
    340 !>AJ
    341384              ! Calcul du qsat
    342385              IF (thermcep) THEN
    … …  
    356399        ENDDO
    357400!AJ<
     401
    358402       IF (.NOT. ice_thermo) THEN
    359403        DO i = 1, klon
    360 !AJ<
    361 !!           IF (zrfl(i) .GT.0.) THEN
     404!          S'il y a des precipitations
    362405           IF (zrfl(i)+zifl(i).GT.0.) THEN
    363 !>AJ
    364406                ! Evap max pour ne pas saturer la fraction sous le nuage
     407                ! Evap max jusqu'à atteindre la saturation dans la partie
     408                ! de la maille qui est sous le nuage de la couche du dessus
     409                !!! On ne tient compte de cette fraction que sous une seule
     410                !!! couche sous le nuage
    365411                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
     412             ! Ajout de la prise en compte des precip a thermiser
     413             ! avec petite reecriture
     414             if  (fl_cor_ebil .GT. 0) then ! nouveau
     415                ! Calcul de l'evaporation du flux de precip herite
     416                !   d'au-dessus
     417                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
     418                     * zmair(i)/pplay(i,k)*zt(i)*RD
     419                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i)
     420
     421                ! Seuil pour ne pas saturer la fraction sous le nuage
     422                zqev = MIN (zqev, zqevt)
     423                ! Nouveau flux de precip
     424                zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime
     425                ! Aucun flux liquide pour T < t_coup, on reevapore tout.
     426                IF (zt(i) .LT. t_coup.and.reevap_ice) THEN
     427                  zrfln(i)=0.
     428                  zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime
     429                END IF
     430                ! Nouvelle vapeur
     431                zq(i) = zq(i) + zqev
     432                zmqc(i) = zmqc(i)-zqev
     433                ! Nouvelle temperature (chaleur latente)
     434                zt(i) = zt(i) - zqev &
     435                     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     436!!JLD debut de partie a supprimer a therme
     437            else ! if  (fl_cor_ebil .GT. 0)
    366438                ! Calcul de l'evaporation du flux de precip herite
    367439                !   d'au-dessus
    … …  
    384456                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    385457                     * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     458              end if ! if  (fl_cor_ebil .GT. 0)
     459!!JLD fin de partie a supprimer a therme
    386460                zrfl(i) = zrfln(i)
    387461                zifl(i) = 0.
    … …  
    390464!
    391465       ELSE ! (.NOT. ice_thermo)
    392 !
     466!      ================================
     467!      Avec thermodynamique de la glace
     468!      ================================
    393469        DO i = 1, klon
    394470!AJ<
    395 !!           IF (zrfl(i) .GT.0.) THEN
    396            IF (zrfl(i)+zifl(i).GT.0.) THEN
    397 !>AJ
    398      !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    399      ! Modification de l'evaporation avec la glace
    400      ! Differentiation entre precipitation liquide et solide
    401      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     471!        S'il y a des precipitations
     472         IF (zrfl(i)+zifl(i).GT.0.) THEN
    402473     
    403      ! Evap max pour ne pas saturer la fraction sous le nuage
     474        ! Evap max pour ne pas saturer la fraction sous le nuage
    404475         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
    405      !    zqev0 = MAX (0.0, zqs(i)-zq(i) )
    406 
    407      !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    408      ! On differencie qsat pour l'eau et la glace
    409      ! Si zdelta=1. --> glace
    410      ! Si zdelta=0. --> eau liquide
    411      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     476
     477         !JAM
     478         ! On differencie qsat pour l'eau et la glace
     479         ! Si zdelta=1. --> glace
     480         ! Si zdelta=0. --> eau liquide
    412481       
    413482         ! Calcul du qsat par rapport a l'eau liquide
    … …  
    417486         qsl= qsl*zcor
    418487         
    419          ! Calcul de l'evaporation du flux de precip herite
    420          !   d'au-dessus
     488         ! Calcul de l'evaporation du flux de precip venant du dessus
    421489         ! Formulation en racine du flux de precip
    422490         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
    … …  
    425493         zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
    426494              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
    427 
    428495         
    429496         ! Calcul du qsat par rapport a la glace
    … …  
    440507              *RG*dtime/(paprs(i,k)-paprs(i,k+1))   
    441508
    442              
    443      !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    444      ! Verification sur l'evaporation
    445      ! On s'assure qu'on ne sature pas
    446      ! la fraction sous le nuage sinon on
    447      ! repartit zqev0 en gardant la proportion
    448      ! liquide / glace
    449      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     509        !JAM
     510        ! Limitation de l'evaporation. On s'assure qu'on ne sature pas
     511        ! la fraction de la couche sous le nuage sinon on repartit zqev0
     512        ! en conservant la proportion liquide / glace
    450513     
    451514         IF (zqevt+zqevti.GT.zqev0) THEN
    452                   zqev=zqev0*zqevt/(zqevt+zqevti)
    453                   zqevi=zqev0*zqevti/(zqevt+zqevti)
    454              
     515            zqev=zqev0*zqevt/(zqevt+zqevti)
     516            zqevi=zqev0*zqevti/(zqevt+zqevti)
    455517         ELSE
     518!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips
     519!       liquides et solides meme si on ne sature pas la couche.
     520!       A mon avis, le test est inutile, et il faudrait tout remplacer par:
     521!            zqev=zqevt
     522!            zqevi=zqevti
    456523             IF (zqevt+zqevti.GT.0.) THEN
    457                   zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
    458                   zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
     524                zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
     525                zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
    459526             ELSE
    460527             zqev=0.
    … …  
    471538         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
    472539                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     540       if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
     541         zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
     542                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     543         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
     544                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     545                  * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
     546                  + (zifln(i)-zifl(i)) &
     547                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     548                  * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     549       else ! sans correction thermalisation des precips
    473550         zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
    474551                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    … …  
    477554                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    478555                  * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
    479      
     556       end if
     557        ! Nouvelles vaeleurs des precips liquides et solides
    480558         zrfl(i) = zrfln(i)
    481559         zifl(i) = zifln(i)
    … …  
    486564!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg
    487565           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
     566           ! fraction de la precip solide qui est fondue
    488567           zmelt = MIN(MAX(zmelt,0.),1.)
    489568           ! Fusion de la glace
    490569           zrfl(i)=zrfl(i)+zmelt*zifl(i)
    491            zifl(i)=zifl(i)*(1.-zmelt)
    492 !           print*,zt(i),'octavio1'
     570           if (fl_cor_ebil .LE. 0) then
     571             ! the following line should not be here. Indeed, if zifl is modified
     572             ! now, zifl(i)*zmelt is no more the amount of ice that has melt
     573             ! and therefore the change in temperature computed below is wrong
     574             zifl(i)=zifl(i)*(1.-zmelt)
     575           end if
    493576           ! Chaleur latente de fusion
     577        if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips
     578           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
     579                      *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     580        else ! sans correction thermalisation des precips
    494581           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
    495582                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
    496 !           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
    497 !fin CR
    498 
    499 
     583        end if
     584           if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente
     585             zifl(i)=zifl(i)*(1.-zmelt)
     586           end if
    500587
    501588           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
    … …  
    515602           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
    516603           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
     604       if  (fl_cor_ebil .GT. 0) then ! nouveau
     605           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     606       else   
    517607           zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
     608       endif
    518609           zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
    519610           zqs(i) = MIN(0.5,zqs(i))
    … …  
    521612           zqs(i) = zqs(i)*zcor
    522613           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
     614           zdqsdT_raw(i) = zdqs(i)*  &
     615         &         RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
    523616        ENDDO
    524617     ELSE
    … …  
    547640     ! P2> Formation du nuage
    548641     ! ----------------------------------------------------------------
     642     ! Variables calculees:
     643     !   rneb  : fraction nuageuse
     644     !   zcond : eau condensee moyenne dans la maille.
     645     !   rhcl: humidite relative ciel-clair
     646     !   zt : temperature de la maille
     647     ! ----------------------------------------------------------------
     648     !
    549649     IF (cpartiel) THEN
    550 
    551         !        print*,'Dans partiel k=',k
     650        ! -------------------------
     651        ! P2.A> Nuage fractionnaire
     652        ! -------------------------
    552653        !
    553654        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
    … …  
    561662        !   Version avec les raqts
    562663
     664        ! ----------------------------------------------------------------
     665        ! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau
     666        ! ----------------------------------------------------------------
    563667        if (iflag_pdf.eq.0) then
    564668
    … …  
    570674           enddo
    571675
    572         else
    573            !
    574            !   Version avec les nouvelles PDFs.
     676        else !  if (iflag_pdf.eq.0)
     677           ! ----------------------------------------------------------------
     678           ! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour
     679           ! les valeurs de T et Q initiales
     680           ! ----------------------------------------------------------------
    575681           do i=1,klon
    576682              if(zq(i).lt.1.e-15) then
    … …  
    620726           enddo
    621727
     728        ! ----------------------------------------------------------------
     729        ! P2.A.2> Prise en compte du couplage entre eau condensee et T.
     730        ! Calcul des grandeurs nuageuses en tenant compte de l'effet de
     731        ! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses
     732        ! qui en dependent. Ce changement de temperature est provisoire, et
     733        ! la valeur definitive sera calcule plus tard.
     734        ! Variables calculees:
     735        !   rneb : nebulosite
     736        !   zcond: eau condensee en moyenne dans la maille
     737        ! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble
     738        ! pas clair, il n'y a probablement pas de prise en compte de l'effet de
     739        ! T sur qsat
     740        ! ----------------------------------------------------------------
    622741
    623742!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
    … …  
    629748               do i=1,klon
    630749!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
     750                 ! !! convergence = .true. tant que l'on n'a pas converge !!
     751                 ! ------------------------------
    631752                 convergence(i)=abs(DT(i)).gt.DDT0
    632753                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
    633                  Tbef(i)=Tbef(i)+DT(i)
     754                 ! si on n'a pas converge
     755                 !
     756                 ! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee
     757                 ! ---------------------------------------------------------------
     758                 ! Variables calculees:
     759                 ! rneb : nebulosite
     760                 ! zqn : eau condensee, dans le nuage (in cloud water content)
     761                 ! zcond: eau condensee en moyenne dans la maille
     762                 ! rhcl: humidite relative ciel-clair
     763                 !
     764                 Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature
    634765                 if (.not.ice_thermo) then   
    635766                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
    … …  
    641772                          zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
    642773                       else
    643 !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
     774                          !avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
    644775                          zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
    645776                       endif
    646777                    endif
    647778                 endif
    648                  ! Calcul des PDF lognormales
     779                 ! Calcul de rneb, qzn et zcond pour les PDF lognormales
    649780                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
     781               if (fl_cor_ebil .GT. 0) then
     782                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     783               else
    650784                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
     785               end if
    651786                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
    652787                 zqs(i) = MIN(0.5,zqs(i))
    … …  
    677812               enddo ! boucle en i
    678813
     814                 ! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT
     815                 !         due a la condensation.
     816                 ! ---------------------------------------------------------------
     817                 ! Variables calculees:
     818                 ! DT : variation de temperature due a la condensation
     819
    679820                 if (.not. ice_thermo) then
    680 
     821                 ! --------------------------
    681822                 do i=1,klon
    682823                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
    683824
    684825                 qlbef(i)=max(0.,zqn(i)-zqs(i))
     826               if (fl_cor_ebil .GT. 0) then
     827                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
     828               else
    685829                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
     830               end if
    686831                 denom=1.+rneb(i,k)*zdqs(i)
    687832                 DT(i)=num/denom
    … …  
    690835                 enddo
    691836
    692                  else
    693                  ! Iteration pour convergence avec qsat(T)
     837                 else ! if (.not. ice_thermo)
     838                 ! --------------------------
    694839                 if (iflag_t_glace.ge.1) then
    695840                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
    … …  
    703848                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    704849                    zfice(i) = zfice(i)**exposant_glace_old
    705                     dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
     850                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) &
     851          &                     / (t_glace_min_old - RTT)
    706852                 endif
    707853                 
    708854                 if (iflag_t_glace.ge.1) then
    709                  dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
     855                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) &
     856          &                    / (t_glace_min - t_glace_max)
    710857                 endif
    711858               
    … …  
    721868
    722869                 qlbef(i)=max(0.,zqn(i)-zqs(i))
    723                  num=-Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
    724                  denom=1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
     870               if (fl_cor_ebil .GT. 0) then
     871                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
     872           &          +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
     873                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
     874                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)    &
     875           &               *qlbef(i)*dzfice(i)
     876               else
     877                 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
     878           &         +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
     879                 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
    725880                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
     881               end if
    726882                 DT(i)=num/denom
    727883                 n_i(i)=n_i(i)+1
    … …  
    732888                 endif !ice_thermo
    733889
    734 !                 endif
    735 !               enddo
    736  
    737          
    738890             enddo ! iter=1,iflag_fisrtilp_qsat+1
    739891             ! Fin d'iteration pour condensation avec variation de qsat(T)
    740892             ! -----------------------------------------------------------
    741            endif
    742 
     893           endif !  if (iflag_fisrtilp_qsat.ge.0)
     894     ! ----------------------------------------------------------------
     895     ! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T
     896     ! ----------------------------------------------------------------
    743897
    744898        endif ! iflag_pdf
    745 
    746899
    747900!        if (iflag_fisrtilp_qsat.eq.-1) then
    … …  
    771924
    772925!        ELSE
    773 
    774         ! Calcul de l'eau in-cloud (zqn),
    775         ! moyenne dans la maille (zcond),
    776         ! fraction nuageuse (rneb) et
    777         ! humidite relative ciel-clair (rhcl)
     926        ! ----------------------------------------------------------------
     927        ! P2.A.3> Calcul des valeures finales associees a la formation des nuages
     928        ! Variables calculees:
     929        !   rneb : nebulosite
     930        !   zcond: eau condensee en moyenne dans la maille
     931        !   zq : eau vapeur dans la maille
     932        !   zt : temperature de la maille
     933        !   rhcl: humidite relative ciel-clair
     934        ! ----------------------------------------------------------------
     935        !
     936        ! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb)
     937        ! Calcule de l'eau condensee moyenne dans la maille (zcond),
     938        ! et de l'humidite relative ciel-clair (rhcl)
    778939        DO i=1,klon
    779940           IF (rneb(i,k) .LE. 0.0) THEN
    … …  
    793954        ENDDO
    794955
    795 
    796956!        ENDIF
    797957
    798958     ELSE ! de IF (cpartiel)
    799         ! Cas "tout ou rien"
     959        ! -------------------------
     960        ! P2.B> Nuage "tout ou rien"
     961        ! -------------------------
     962        ! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences?
    800963        DO i = 1, klon
    801964           IF (zq(i).GT.zqs(i)) THEN
    … …  
    806969           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
    807970        ENDDO
    808      ENDIF
    809      ! ----------------------------------------------------------------
    810      ! Fin de formation du nuage
    811      ! ----------------------------------------------------------------
     971     ENDIF ! de IF (cpartiel)
    812972     !
    813973     ! Mise a jour vapeur d'eau
     974     ! -------------------------
    814975     DO i = 1, klon
    815976        zq(i) = zq(i) - zcond(i)
    … …  
    817978     ENDDO
    818979!AJ<
    819      ! Chaleur latente apres formation nuage
     980     ! ------------------------------------
     981     ! P2.C> Prise en compte de la Chaleur latente apres formation nuage
    820982     ! -------------------------------------
     983     ! Variable calcule:
     984     !   zt : temperature de la maille
     985     !
    821986     IF (.NOT. ice_thermo) THEN
    822987        if (iflag_fisrtilp_qsat.lt.1) then
    … …  
    826991        else if (iflag_fisrtilp_qsat.gt.0) then
    827992           DO i= 1, klon
     993    if (fl_cor_ebil .GT. 0) then
     994             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     995    else
    828996             zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
     997    end if
    829998           ENDDO
    830999        endif
    … …  
    8561025                    zfice(i) = zfice(i)**exposant_glace_old
    8571026              endif
     1027        if (fl_cor_ebil .GT. 0) then
     1028              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
     1029           &             * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
     1030                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     1031        else
    8581032              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
    859                        +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
     1033                      +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
     1034        end if
    8601035           ENDDO
    8611036         endif
    … …  
    8631038       ENDIF
    8641039!>AJ
     1040
    8651041     ! ----------------------------------------------------------------
    8661042     ! P3> Formation des precipitations
    … …  
    9581134     !         zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie   , 0.0)
    9591135     !         zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice   , 0.0)
     1136!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip
     1137!      temporaires et ne doivent pas etre calcule (alors qu'elles le sont
     1138!      si iflag_bergeron <> 2
     1139!      A SUPPRIMER A TERME
    9601140              zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
    9611141              zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
    … …  
    9661146        ENDDO  ! i = 1,klon
    9671147     ENDDO     ! n = 1,ninter
     1148
    9681149     ! ----------------------------------------------------------------
    9691150     !
    … …  
    9941175             zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
    9951176             zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
     1177           if (fl_cor_ebil .GT. 0) then
     1178             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
     1179             coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
     1180!            Calcul de DT si toute les precips liquides congelent
     1181             DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
     1182!            On ne veut pas que T devienne superieur a la temp. de congelation.
     1183!            donc que Delta > RTT-zt(i
     1184             DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
     1185             zt(i)      = zt(i)      + DeltaT
     1186!            Eau vaporisee du fait de l'augmentation de T
     1187             Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
     1188!            on reajoute cette eau vaporise a la vapeur et on l'enleve des precips
     1189             zq(i) = zq(i) +  Deltaq
     1190!            Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies
     1191             zcond(i)   = max( zcond(i)- Deltaq, 0. )
     1192!            precip liquide qui congele ou qui s'evapore
     1193             Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
     1194             zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
     1195!            bilan eau glacee
     1196             zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
     1197           else ! if (fl_cor_ebil .GT. 0)
     1198!            ancien calcul
    9961199             zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i)))
    9971200             coef1 = RLMLT*zdqs(i)/RLVTT
    … …  
    10021205             zq(i)      = zq(i)      + zcp/RLVTT*zdqs(i)*DeltaT
    10031206             zt(i)      = zt(i)      + DeltaT
     1207           end if ! if (fl_cor_ebil .GT. 0)
    10041208           ENDIF  ! rneb(i,k) .GT. 0.0
    10051209         ENDDO
    … …  
    10211225         IF (rneb(i,k).GT.0.0) THEN
    10221226!CR on prend en compte la phase glace
    1023            if (.not.ice_thermo) then
    1024            d_ql(i,k) = zoliq(i)
    1025            d_qi(i,k) = 0.
    1026            else
     1227!JLD inutile car on ne passe jamais ici si .not.ice_thermo
     1228!           if (.not.ice_thermo) then
     1229!           d_ql(i,k) = zoliq(i)
     1230!           d_qi(i,k) = 0.
     1231!           else
    10271232           d_ql(i,k) = (1-zfice(i))*zoliq(i)
    10281233           d_qi(i,k) = zfice(i)*zoliq(i)
    1029            endif
     1234!           endif
    10301235!AJ<
    10311236           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
    … …  
    10431248              zifl(i) = zifl(i)+zrfl(i)
    10441249              zrfl(i) = 0.
     1250           if (fl_cor_ebil .GT. 0) then
     1251              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
     1252                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
     1253           else
    10451254              zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
    10461255                      *(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))
     1256           end if
    10471257           ENDIF  ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)
    10481258!RC   
    … …  
    10941304     ! Fin de formation des precipitations
    10951305     ! ----------------------------------------------------------------
    1096      !
    10971306     !
    10981307     ! Calculer les tendances de q et de t:
    … …  
    12021411     DO i = 1, klon
    12031412        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
    1204         zmair=(paprs(i,k)-paprs(i,k+1))/RG
     1413        zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
    12051414        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
    1206         d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
     1415        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i))
    12071416     END DO
    12081417  END DO

  • LMDZ5/trunk/libf/phylmd/reevap.F90

    r2799 r2807  
    2424    ! Re-evaporer l'eau liquide nuageuse
    2525    !
    26 
     26print *,'rrevap ; fl_cor_ebil:',fl_cor_ebil,' iflag_ice_thermo:',iflag_ice_thermo,' RVTMP2',RVTMP2
    2727    DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
    2828       DO i = 1, klon
    … …  
    4646
    4747             zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
     48  zdelta = 0.
    4849             zb = MAX(0.0,ql_seri(i,k))
    4950             za = - MAX(0.0,ql_seri(i,k)) &
Note: See TracChangeset for help on using the changeset viewer.