Changeset 2500


Ignore:
Timestamp:
Apr 28, 2016, 4:58:11 PM (8 years ago)
Author:
jbmadeleine
Message:

Ajouts de commentaires dans fisrtilp pour clarifier la lecture.
Ajout de flags P1> P2> P3> pour se situer dans le code plus facilement.
Aucune incidence sur les resultats du modele.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2466 r2500  
    2121  ! Objet: condensation et precipitation stratiforme.
    2222  !        schema de nuage
     23  ! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval)
     24  !             et ilp (il pleut, L. Li)
     25  ! Principales parties:
     26  ! P1> Evaporation de la precipitation (qui vient du niveau k+1)
     27  ! P2> Formation du nuage (en k)
     28  ! P3> Formation de la precipitation (en k)
    2329  !======================================================================
    2430  !======================================================================
    … …  
    2834
    2935  !
    30   ! Arguments:
     36  ! Principaux inputs:
    3137  !
    3238  REAL dtime ! intervalle du temps (s)
    … …  
    3541  REAL t(klon,klev) ! temperature (K)
    3642  REAL q(klon,klev) ! humidite specifique (kg/kg)
     43  !
     44  ! Principaux outputs:
     45  !
    3746  REAL d_t(klon,klev) ! incrementation de la temperature (K)
    3847  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
    … …  
    4655  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    4756  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     57  !
     58  ! Autres arguments
     59  !
    4860  REAL ztv(klon,klev)
    4961  REAL zqta(klon,klev),fraca(klon,klev)
    … …  
    155167  !     Pour la conversion eau-neige
    156168  REAL zlh_solid(klon), zm_solid
    157   !IM
    158   !ym      INTEGER klevm1
    159169  !---------------------------------------------------------------
    160170  !
    161171  ! Fonctions en ligne:
    162172  !
    163   REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
     173  REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace
     174                     ! (Heymsfield & Donner, 1990)
    164175  REAL zzz
    165176  include "YOETHF.h"
    … …  
    278289  zfrac_lessi = 0.
    279290
    280   !AA----------------------------------------------------------
     291  !AA==================================================================
    281292  !
    282293  ncoreczq=0
    283   ! Boucle verticale (du haut vers le bas)
    284   !
    285   !IM : klevm1
    286   !ym      klevm1=klev-1
     294  ! BOUCLE VERTICALE (DU HAUT VERS LE BAS)
     295  !
    287296  DO k = klev, 1, -1
    288297     !
    289      !AA----------------------------------------------------------
    290      !
     298     !AA===============================================================
     299     !
     300     ! Initialisation temperature et vapeur
    291301     DO i = 1, klon
    292302        zt(i)=t(i,k)
    … …  
    312322        ENDDO
    313323     ENDIF
    314      !
    315      !
    316      ! Calculer l'evaporation de la precipitation
    317      !
    318 
    319 
    320      ! Calculer l'evaporation de la precipitation
    321      !
    322 
    323 
     324     ! ----------------------------------------------------------------
     325     ! P1> Debut evaporation de la precipitation
     326     ! ----------------------------------------------------------------
    324327     IF (evap_prec) THEN
    325328        DO i = 1, klon
    … …  
    328331           IF (zrfl(i)+zifl(i).GT.0.) THEN
    329332!>AJ
     333              ! Calcul du qsat
    330334              IF (thermcep) THEN
    331335                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
    … …  
    350354           IF (zrfl(i)+zifl(i).GT.0.) THEN
    351355!>AJ
     356                ! Evap max pour ne pas saturer la fraction sous le nuage
    352357                zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
     358                ! Calcul de l'evaporation du flux de precip herite
     359                !   d'au-dessus
    353360                zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
    354361                     * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    355362                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
    356363                     * RG*dtime/(paprs(i,k)-paprs(i,k+1))
     364                ! Seuil pour ne pas saturer la fraction sous le nuage
    357365                zqev = MIN (zqev, zqevt)
     366                ! Nouveau flux de precip
    358367                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
    359368                     /RG/dtime
    360        
    361                 ! pour la glace, on ré-évapore toute la précip dans la
    362                 ! couche du dessous
    363                 ! la glace venant de la couche du dessus est simplement
    364                 ! dans la couche du dessous.
    365        
     369                ! Aucun flux liquide pour T < t_coup
    366370                IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
    367        
     371                ! Nouvelle vapeur
    368372                zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
    369373                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     374                ! Nouvelle temperature (chaleur latente)
    370375                zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
    371376                     * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
    … …  
    384389!>AJ
    385390     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    386      ! Modification de l'évaporation avec la glace
    387      ! Différentiation entre précipitation liquide et solide
    388      ! On suppose que coef_evai=2*coef_eva
     391     ! Modification de l'evaporation avec la glace
     392     ! Differentiation entre precipitation liquide et solide
    389393     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    390394     
     395     ! Evap max pour ne pas saturer la fraction sous le nuage
    391396         zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
    392397     !    zqev0 = MAX (0.0, zqs(i)-zq(i) )
    393398
    394399     !JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    395      ! On différencie qsat pour l'eau et la glace
     400     ! On differencie qsat pour l'eau et la glace
    396401     ! Si zdelta=1. --> glace
    397402     ! Si zdelta=0. --> eau liquide
    398403     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    399          
     404       
     405         ! Calcul du qsat par rapport a l'eau liquide
    400406         qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
    401407         qsl= MIN(0.5,qsl)
    … …  
    403409         qsl= qsl*zcor
    404410         
     411         ! Calcul de l'evaporation du flux de precip herite
     412         !   d'au-dessus
     413         ! Formulation en racine du flux de precip
     414         ! (Klemp & Wilhelmson, 1978; Sundqvist, 1988)
    405415         zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
    406416              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    … …  
    408418              *RG*dtime/(paprs(i,k)-paprs(i,k+1))
    409419
     420         
     421         ! Calcul du qsat par rapport a la glace
    410422         qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
    411423         qsi= MIN(0.5,qsi)
    … …  
    413425         qsi= qsi*zcor
    414426
     427         ! Calcul de la sublimation du flux de precip solide herite
     428         !   d'au-dessus
    415429         zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
    416430              *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    … …  
    420434             
    421435     !JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    422      ! Vérification sur l'évaporation
     436     ! Verification sur l'evaporation
     437     ! On s'assure qu'on ne sature pas
     438     ! la fraction sous le nuage sinon on
     439     ! repartit zqev0 en gardant la proportion
     440     ! liquide / glace
    423441     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    424442     
    … …  
    436454             ENDIF
    437455         ENDIF
    438      
     456         ! Nouveaux flux de precip liquide et solide
    439457         zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
    440458                                 /RG/dtime)
    441459         zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
    442460                                 /RG/dtime)
    443          
    444      ! Pour la glace, on révapore toute la précip dans la couche du dessous
    445      ! la glace venant de la couche du dessus est simplement dans la couche
    446      ! du dessous.
    447      
    448      !    IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
    449 !         print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
     461
     462         ! Mise a jour de la vapeur, temperature et flux de precip
    450463         zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
    451464                  * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    … …  
    466479           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))             ! jyg
    467480           zmelt = MIN(MAX(zmelt,0.),1.)
     481           ! Fusion de la glace
    468482           zrfl(i)=zrfl(i)+zmelt*zifl(i)
    469483           zifl(i)=zifl(i)*(1.-zmelt)
    470484!           print*,zt(i),'octavio1'
     485           ! Chaleur latente de fusion
    471486           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
    472487                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
    … …  
    481496      ENDIF ! (.NOT. ice_thermo)
    482497     
     498     ! ----------------------------------------------------------------
     499     ! Fin evaporation de la precipitation
     500     ! ----------------------------------------------------------------
    483501     ENDIF ! (evap_prec)
    484502     !
    … …  
    518536!       endif
    519537
     538     ! ----------------------------------------------------------------
     539     ! P2> Formation du nuage
     540     ! ----------------------------------------------------------------
    520541     IF (cpartiel) THEN
    521542
    … …  
    534555        if (iflag_pdf.eq.0) then
    535556
     557           ! version creneau de (Li, 1998)
    536558           do i=1,klon
    537559              zdelq = min(ratqs(i,k),0.99) * zq(i)
    … …  
    575597           end if
    576598
    577 !CR: variation de qsat avec T en présence de glace ou non
     599!CR: variation de qsat avec T en presence de glace ou non
    578600!initialisations
    579601           do i=1,klon
    … …  
    587609!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
    588610           if (iflag_fisrtilp_qsat.ge.0) then
     611             ! Iteration pour condensation avec variation de qsat(T)
     612             ! -----------------------------------------------------
    589613             do iter=1,iflag_fisrtilp_qsat+1
    590614               
    … …  
    603627                    endif
    604628                 endif
     629                 ! Calcul des PDF lognormales
    605630                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
    606631                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
    … …  
    647672
    648673                 else
    649 
    650 !calcul de la fraction de glace
    651 !CR: on utilise la nouvelle fonction de JBM pour l ancien calcul
    652 !                 zfice(i) = icefrac_lsc(Tbef(i), t_glace_min, &
    653 !                                     t_glace_max, exposant_glace)
    654 !                 zfice(i) = 1.0 - (Tbef(i)-ztglace) / (RTT-ztglace)
    655 !                 zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    656 !                 zfice(i) = zfice(i)**nexpo
     674                 ! Iteration pour convergence avec qsat(T)
    657675                 if (iflag_t_glace.eq.1) then
    658676                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:))
    … …  
    699717 
    700718         
    701              enddo
     719             enddo ! iter=1,iflag_fisrtilp_qsat+1
     720             ! Fin d'iteration pour condensation avec variation de qsat(T)
     721             ! -----------------------------------------------------------
    702722           endif
    703723
    … …  
    707727
    708728!        if (iflag_fisrtilp_qsat.eq.-1) then
    709 !CR: ATTENTION option fausse mais a existe: pour la re-activer, prendre iflag_fisrtilp_qsat=0 et activer les lignes suivantes:
     729!------------------------------------------
     730!CR: ATTENTION option fausse mais a existe:
     731! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et
     732! activer les lignes suivantes:
    710733       IF (1.eq.0) THEN
    711734       DO i=1,klon
    … …  
    724747              rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
    725748           ENDIF
    726         ENDDO
    727         ENDIF
     749       ENDDO
     750       ENDIF
     751!------------------------------------------
    728752
    729753!        ELSE
    730754
     755        ! Calcul de l'eau in-cloud (zqn),
     756        ! moyenne dans la maille (zcond),
     757        ! fraction nuageuse (rneb) et
     758        ! humidite relative ciel-clair (rhcl)
    731759        DO i=1,klon
    732760           IF (rneb(i,k) .LE. 0.0) THEN
    … …  
    749777!        ENDIF
    750778
    751         !         do i=1,klon
    752         !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
    753         !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
    754         !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
    755         !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
    756         !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
    757         !c  la convection.
    758         !c  ATTENTION !!! Il va falloir verifier tout ca.
    759         !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
    760         !c           print*,'ZDQS ',zdqs(i)
    761         !c--Olivier
    762         !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
    763         !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
    764         !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
    765         !c--fin
    766         !           ENDDO
    767      ELSE
     779     ELSE ! de IF (cpartiel)
     780        ! Cas "tout ou rien"
    768781        DO i = 1, klon
    769782           IF (zq(i).GT.zqs(i)) THEN
    … …  
    775788        ENDDO
    776789     ENDIF
    777      !
     790     ! ----------------------------------------------------------------
     791     ! Fin de formation du nuage
     792     ! ----------------------------------------------------------------
     793     !
     794     ! Mise a jour vapeur d'eau
    778795     DO i = 1, klon
    779796        zq(i) = zq(i) - zcond(i)
    … …  
    781798     ENDDO
    782799!AJ<
     800     ! Chaleur latente apres formation nuage
     801     ! -------------------------------------
    783802     IF (.NOT. ice_thermo) THEN
    784803        if (iflag_fisrtilp_qsat.lt.1) then
    … …  
    825844       ENDIF
    826845!>AJ
     846     ! ----------------------------------------------------------------
     847     ! P3> Formation des precipitations
     848     ! ----------------------------------------------------------------
    827849     !
    828850     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
    829851     !
     852
     853     ! Initialisation de zoliq (eau condensee moyenne dans la maille)
    830854     DO i = 1, klon
    831855        IF (rneb(i,k).GT.0.0) THEN
    … …  
    858882       ENDIF
    859883     ENDIF
     884
     885     ! Calcul de radliq (eau condensee pour le rayonnement)
     886     ! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip
     887     ! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une
     888     ! eau moyenne restante dans le nuage sur la duree du pas de temps qui est
     889     ! transmise au rayonnement;
     890     ! ----------------------------------------------------------------
    860891     DO i = 1, klon
    861892        IF (rneb(i,k).GT.0.0) THEN
    … …  
    878909                 ztot = 0.0
    879910              ELSE
    880                  !  quantite d'eau a eliminer: zchau
    881                  !  meme chose pour la glace: zfroi
     911                 !  quantite d'eau a eliminer: zchau (Sundqvist, 1978)
     912                 !  meme chose pour la glace: zfroi (Zender & Kiehl, 1997)
    882913                 if (ptconv(i,k)) then
    883914                    zcl   =cld_lc_con
    … …  
    916947        ENDDO  ! i = 1,klon
    917948     ENDDO     ! n = 1,ninter
     949     ! ----------------------------------------------------------------
    918950     !
    919951     IF (.NOT. ice_thermo) THEN
    … …  
    10281060     ELSE
    10291061     ! JAM*************************************************
    1030      ! Revoir partie ci-dessous: à quoi servent psfl et prfl?
     1062     ! Revoir partie ci-dessous: a quoi servent psfl et prfl?
    10311063     ! *****************************************************
    10321064
    … …  
    10401072       ENDDO
    10411073     ENDIF
     1074     ! ----------------------------------------------------------------
     1075     ! Fin de formation des precipitations
     1076     ! ----------------------------------------------------------------
    10421077     !
    10431078     !
    … …  
    11101145     ENDDO
    11111146     !
    1112      !AA----------------------------------------------------------
    1113      !                     FIN DE BOUCLE SUR K  
     1147     !AA===============================================================
     1148     !                     FIN DE LA BOUCLE VERTICALE 
    11141149  end DO
    11151150  !
    1116   !AA-----------------------------------------------------------
     1151  !AA==================================================================
    11171152  !
    11181153  ! Pluie ou neige au sol selon la temperature de la 1ere couche
Note: See TracChangeset for help on using the changeset viewer.