Changeset 2086


Ignore:
Timestamp:
Jul 9, 2014, 11:19:07 PM (10 years ago)
Author:
fhourdin
Message:

nclusion de la thermodynamique de la glace
Ice thermodynamics
(Catherine Rio)

Location:
LMDZ5/trunk/libf
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/trunk/libf/dyn3d_common/infotrac.F90

    r1945 r2086  
    55! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
    66  INTEGER, SAVE :: nqtot
     7!CR: on ajoute le nombre de traceurs de l eau
     8  INTEGER, SAVE :: nqo
    79
    810! nbtr : number of tracers not including higher order of moment or water vapor or liquid
     
    228230         endif ! of if (planet_type=="earth")
    229231       END IF
     232
     233!CR: nombre de traceurs de l eau
     234       if (tnom_0(3) == 'H2Oi') then
     235          nqo=3
     236       else
     237          nqo=2
     238       endif
    230239       
    231240       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
  • LMDZ5/trunk/libf/phylmd/add_pbl_tend.F90

    r1998 r2086  
    1 SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, paprs, text)
     1SUBROUTINE add_pbl_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text)
    22  ! ======================================================================
    33  ! Ajoute les tendances de couche limite, soit determinees par la
     
    2727  ! ------------
    2828  REAL zdu(klon, klev), zdv(klon, klev)
    29   REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev)
     29  REAL zdt(klon, klev), zdq(klon, klev), zdql(klon, klev), zdqi(klon, klev)
    3030  CHARACTER *(*) text
    3131  REAL paprs(klon,klev+1)
     
    4646    PRINT *, ' add_pbl_tend, zzdt ', zzdt
    4747    PRINT *, ' add_pbl_tend, zzdq ', zzdq
    48     CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, paprs, text)
     48    CALL add_phys_tend(zdu, zdv, zzdt, zzdq, zdql, zdqi, paprs, text)
    4949  ELSE
    50     CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, paprs, text)
     50    CALL add_phys_tend(zdu, zdv, zdt, zdq, zdql, zdqi, paprs, text)
    5151  END IF
    5252
  • LMDZ5/trunk/libf/phylmd/add_phys_tend.F90

    r2007 r2086  
    22! $Id$
    33!
    4 SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,paprs,text)
     4SUBROUTINE add_phys_tend (zdu,zdv,zdt,zdq,zdql,zdqi,paprs,text)
    55!======================================================================
    66! Ajoute les tendances des variables physiques aux variables
     
    2525!------------
    2626REAL zdu(klon,klev),zdv(klon,klev)
    27 REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev)
     27REAL zdt(klon,klev),zdq(klon,klev),zdql(klon,klev),zdqi(klon,klev)
    2828REAL paprs(klon,klev+1)
    2929CHARACTER*(*) text
     
    6666     v_seri(:,:)=v_seri(:,:)+zdv(:,:)
    6767     ql_seri(:,:)=ql_seri(:,:)+zdql(:,:)
     68     qs_seri(:,:)=qs_seri(:,:)+zdqi(:,:)
    6869
    6970!======================================================================
  • LMDZ5/trunk/libf/phylmd/fisrtilp.F90

    r2077 r2086  
    44!
    55SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
    6      d_t, d_q, d_ql, rneb, radliq, rain, snow,          &
     6     d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow,          &
    77     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
    88     frac_impa, frac_nucl, beta,                        &
     
    4141  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
    4242  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
     43  REAL d_qi(klon,klev) ! incrementation de l'eau glace
    4344  REAL rneb(klon,klev) ! fraction nuageuse
    4445  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
     
    107108  REAL DDT0
    108109  PARAMETER (DDT0=.01)
    109   INTEGER n_i, iter
     110  INTEGER n_i(klon), iter
     111  REAL cste
    110112 
    111113  REAL zrfl(klon), zrfln(klon), zqev, zqevt
     
    123125  REAL zmelt, zpluie, zice, zcondold
    124126  PARAMETER (ztfondue=278.15)
     127  REAL dzfice(klon)
    125128  !
    126129  LOGICAL appel1er
     
    165168  DATA appel1er /.TRUE./
    166169  !ym
    167 
    168   ice_thermo = iflag_ice_thermo .GE. 1
     170!CR: pour iflag_ice_thermo=2, on active que la convection
     171!  ice_thermo = iflag_ice_thermo .GE. 1
     172   ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3)
    169173  zdelq=0.0
    170174
     
    209213  !  nexpo regle la raideur de la transition eau liquide / eau glace.
    210214  !
     215!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut
    211216  IF (iflag_t_glace.EQ.0) THEN
    212217!   ztglace = RTT - 15.0
     
    220225      exposant_glace_old = 6
    221226    ENDIF
     227   
    222228  ENDIF
    223229 
     
    243249        d_q(i,k) = 0.0
    244250        d_ql(i,k) = 0.0
     251        d_qi(i,k) = 0.0
    245252        rneb(i,k) = 0.0
    246253        radliq(i,k) = 0.0
     
    451458         zifl(i) = zifln(i)
    452459
     460!CR ATTENTION: deplacement de la fonte de la glace
     461           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
     462           zmelt = MIN(MAX(zmelt,0.),1.)
     463           zrfl(i)=zrfl(i)+zmelt*zifl(i)
     464           zifl(i)=zifl(i)*(1.-zmelt)
     465!           print*,zt(i),'octavio1'
     466           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
     467                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
     468!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
     469!fin CR
     470
     471
     472
    453473           ENDIF ! (zrfl(i)+zifl(i).GT.0.)
    454474        ENDDO
     
    487507     !
    488508!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
    489        if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
    490          write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
    491         " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
    492          stop
    493        endif
     509!       if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
     510!         write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
     511!        " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
     512!         stop
     513!       endif
    494514
    495515     IF (cpartiel) THEN
     
    550570           end if
    551571
     572!CR: variation de qsat avec T en présence de glace ou non
     573!initialisations
    552574           do i=1,klon
     575              DT(i) = 0.
     576              n_i(i)=0
    553577              Tbef(i)=zt(i)
    554               qlbef(i)=0.
    555               if (lognormale(i)) then
     578              qlbef(i)=0.
     579           enddo
     580
     581
     582!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici
     583           if (iflag_fisrtilp_qsat.ge.0) then
     584             do iter=1,iflag_fisrtilp_qsat+1
     585               
     586               do i=1,klon
     587!                 do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0))
     588                 convergence(i)=abs(DT(i)).gt.DDT0
     589                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
     590                 Tbef(i)=Tbef(i)+DT(i)
     591                 if (.not.ice_thermo) then   
     592                    zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i)))
     593                 else
     594                    if (iflag_t_glace.eq.0) then
     595                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
     596                    else if (iflag_t_glace.eq.1) then
     597                    zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
     598                    endif
     599                 endif
     600                 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
     601                 zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
     602                 zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k)
     603                 zqs(i) = MIN(0.5,zqs(i))
     604                 zcor = 1./(1.-RETV*zqs(i))
     605                 zqs(i) = zqs(i)*zcor
     606                 zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
    556607                 zpdf_sig(i)=ratqs(i,k)*zq(i)
    557608                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
     
    574625                 endif
    575626
     627                 endif !convergence
     628               enddo ! boucle en i
     629
     630                 if (.not. ice_thermo) then
     631
     632                 do i=1,klon
     633                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
     634
    576635                 qlbef(i)=max(0.,zqn(i)-zqs(i))
    577636                 num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
    578637                 denom=1.+rneb(i,k)*zdqs(i)
    579638                 DT(i)=num/denom
    580               endif
    581            enddo
     639                 n_i(i)=n_i(i)+1
     640                 endif
     641                 enddo
     642
     643                 else
     644
     645!calcul de la fraction de glace
     646!CR: on utilise la nouvelle fonction de JBM pour l ancien calcul
     647!                 zfice(i) = icefrac_lsc(Tbef(i), t_glace_min, &
     648!                                     t_glace_max, exposant_glace)
     649!                 zfice(i) = 1.0 - (Tbef(i)-ztglace) / (RTT-ztglace)
     650!                 zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
     651!                 zfice(i) = zfice(i)**nexpo
     652                 if (iflag_t_glace.eq.1) then
     653                 CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
     654                      t_glace_min,t_glace_max,exposant_glace,zfice(:))
     655                 endif
     656
     657                 do i=1,klon
     658                 if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then
     659                 
     660                 if (iflag_t_glace.eq.0) then
     661                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
     662                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
     663                    zfice(i) = zfice(i)**exposant_glace_old
     664                    dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) / (t_glace_min_old - RTT)
     665                 endif
     666                 
     667                 if (iflag_t_glace.eq.1) then
     668                 dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) / (t_glace_min - t_glace_max)
     669                 endif
     670               
     671                 if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then
     672                    dzfice(i)=0.
     673                 endif
     674
     675                 if (zfice(i).lt.1) then
     676                    cste=RLVTT
     677                 else
     678                    cste=RLSTT
     679                 endif
     680
     681                 qlbef(i)=max(0.,zqn(i)-zqs(i))
     682                 num=-Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
     683                 denom=1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
     684                         -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i)
     685                 DT(i)=num/denom
     686                 n_i(i)=n_i(i)+1
     687
     688                 endif ! fin convergence
     689                 enddo ! fin boucle i
     690
     691                 endif !ice_thermo
     692
     693!                 endif
     694!               enddo
    582695 
    583            n_i=1
    584 !               do while ((abs(DT(i)).gt.DDT0).and.((zqn(i)-zqs(i)).gt.0.))
    585 ! iflag_fisrtilp_qsat=0: qsat ne varie pas avec T
    586 ! iflag_fisrtilp_qsat > 1 : nombre d iterations max pour converger sur le calcul de qsat(T)
    587 !               do while (n_i.le.iflag_fisrtilp_qsat)
    588            if (iflag_fisrtilp_qsat.ge.1) then
    589            do iter=1,iflag_fisrtilp_qsat
    590               do i=1,klon
    591                  convergence(i)=abs(DT(i)).gt.DDT0
    592                  if (convergence(i).and.lognormale(i)) then
    593                  Tbef(i)=Tbef(i)+DT(i)                 
    594                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef(i)))
    595                  zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
    596                  zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
    597                  zqs(i)= R2ES * FOEEW(Tbef(i),zdelta)/pplay(i,k)
    598                  zqs(i)=MIN(0.5,zqs(i))
    599                  zcor=1./(1.-retv*zqs(i))
    600                  zqs(i)=zqs(i)*zcor
    601                  zdqs(i)=FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
    602 
    603                  zpdf_sig(i)=ratqs(i,k)*zq(i)
    604                  zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
    605                  zpdf_delta(i)=log(zq(i)/zqs(i))
    606                  zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
    607                  zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
    608                  zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
    609                  zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
    610                  zpdf_e1(i)=1.-erf(zpdf_e1(i))
    611                  zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
    612                  zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
    613                  zpdf_e2(i)=1.-erf(zpdf_e2(i))
    614              
    615                  if (zpdf_e1(i).lt.1.e-10) then
    616                     rneb(i,k)=0.
    617                     zqn(i)=zqs(i)
    618                  else
    619                     rneb(i,k)=0.5*zpdf_e1(i)
    620                     zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
    621                  endif
    622 
    623                  qlbef(i)=max(0.,zqn(i)-zqs(i))   
    624                  num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i)
    625                  denom=1.+rneb(i,k)*zdqs(i)
    626                  DT(i)=num/denom
    627                  n_i=n_i+1
    628                  endif
    629               enddo
    630               enddo
    631 
     696         
     697             enddo
    632698           endif
    633699
     
    635701        endif ! iflag_pdf
    636702
    637         if (iflag_fisrtilp_qsat.eq.-1) then
    638 !CR: ATTENTION option fausse mais a existe
    639         DO i=1,klon
     703
     704!        if (iflag_fisrtilp_qsat.eq.-1) then
     705!CR: ATTENTION option fausse mais a existe: pour la re-activer, prendre iflag_fisrtilp_qsat=0 et activer les lignes suivantes:
     706       IF (1.eq.0) THEN
     707       DO i=1,klon
    640708           IF (rneb(i,k) .LE. 0.0) THEN
    641709              zqn(i) = 0.0
     
    653721           ENDIF
    654722        ENDDO
    655 
    656         ELSE
     723        ENDIF
     724
     725!        ELSE
    657726
    658727        DO i=1,klon
     
    674743
    675744
    676         ENDIF
     745!        ENDIF
    677746
    678747        !         do i=1,klon
     
    723792        endif
    724793     ELSE
    725        IF (iflag_t_glace.EQ.0) THEN
     794         if (iflag_t_glace.eq.1) then
     795            CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
     796                 t_glace_min,t_glace_max,exposant_glace,zfice(:))
     797         endif
    726798         if (iflag_fisrtilp_qsat.lt.1) then
    727799           DO i = 1, klon
    728               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.15-t_glace_min_old)
    729               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    730               zfice(i) = zfice(i)**exposant_glace_old
    731 !             zfice(i) = zfice(i)**nexpo
     800! JBM: icefrac_lsc is now a function contained in microphys_mod
     801!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
     802!                                     t_glace_max, exposant_glace)
     803              if (iflag_t_glace.eq.0) then
     804                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
     805                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
     806                    zfice(i) = zfice(i)**exposant_glace_old
     807              endif
    732808              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
    733809                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
     
    735811         else
    736812           DO i=1, klon
    737               zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.15-t_glace_min_old)
    738               zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    739               zfice(i) = zfice(i)**exposant_glace_old
    740 !             zfice(i) = zfice(i)**nexpo
    741 !CR: ATTENTION zt different de Tbef: à corriger
     813              if (lognormale(i)) then
     814                 zt(i)=Tbef(i)
     815              else
     816! JBM: icefrac_lsc is now a function contained in microphys_mod
     817!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
     818!                                     t_glace_max, exposant_glace)
     819              if (iflag_t_glace.eq.0) then
     820                    zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old)
     821                    zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
     822                    zfice(i) = zfice(i)**exposant_glace_old
     823              endif
    742824              zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
    743825                       +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
    744            ENDDO
    745          endif
    746 !         print*,zt(i),zrfl(i),zifl(i),'temp1'
    747        ELSE ! of IF (iflag_t_glace.EQ.0)
    748          CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
    749 &               t_glace_min,t_glace_max,exposant_glace,zfice(:))
    750          if (iflag_fisrtilp_qsat.lt.1) then
    751            DO i = 1, klon
    752 ! JBM: icefrac_lsc is now a function contained in microphys_mod
    753 !             zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
    754 !                                    t_glace_max, exposant_glace)
    755               zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
    756                        +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i))
    757            ENDDO
    758          else
    759            DO i=1, klon
    760 ! JBM: icefrac_lsc is now a function contained in microphys_mod
    761 !             zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
    762 !                                    t_glace_max, exposant_glace)
    763 !CR: ATTENTION zt different de Tbef: à corriger
    764               zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
    765                        +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
     826              endif
    766827           ENDDO
    767828         endif
    768829!         print*,zt(i),zrfl(i),zifl(i),'temp1'
    769830       ENDIF
    770      ENDIF
    771831!>AJ
    772832     !
     
    794854       ELSE ! of IF (iflag_t_glace.EQ.0)
    795855         CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1), &
    796 &               t_glace_min,t_glace_max,exposant_glace,zfice(:))
    797 !        DO i = 1, klon
    798 !           IF (rneb(i,k).GT.0.0) THEN
     856               t_glace_min,t_glace_max,exposant_glace,zfice(:))
     857!         DO i = 1, klon
     858!            IF (rneb(i,k).GT.0.0) THEN
    799859! JBM: icefrac_lsc is now a function contained in microphys_mod
    800 !             zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
    801 !                                    t_glace_max, exposant_glace)
    802 !           ENDIF
    803 !        ENDDO
     860!              zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
     861!                                     t_glace_max, exposant_glace)
     862!            ENDIF
     863!         ENDDO
    804864       ENDIF
    805865     ENDIF
     
    874934       DO i = 1, klon
    875935         IF (rneb(i,k).GT.0.0) THEN
     936!CR on prend en compte la phase glace
     937           if (.not.ice_thermo) then
    876938           d_ql(i,k) = zoliq(i)
     939           d_qi(i,k) = 0.
     940           else
     941           d_ql(i,k) = (1-zfice(i))*zoliq(i)
     942           d_qi(i,k) = zfice(i)*zoliq(i)
     943           endif
    877944!AJ<
    878945           zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
     
    889956     ENDIF
    890957
    891      IF (ice_thermo) THEN
    892        DO i = 1, klon
    893            zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
    894            zmelt = MIN(MAX(zmelt,0.),1.)
    895            zrfl(i)=zrfl(i)+zmelt*zifl(i)
    896            zifl(i)=zifl(i)*(1.-zmelt)
     958!CR: la fonte est faite au debut
     959!      IF (ice_thermo) THEN
     960!       DO i = 1, klon
     961!           zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
     962!           zmelt = MIN(MAX(zmelt,0.),1.)
     963!           zrfl(i)=zrfl(i)+zmelt*zifl(i)
     964!           zifl(i)=zifl(i)*(1.-zmelt)
    897965!           print*,zt(i),'octavio1'
    898            zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
    899                       *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
     966!           zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
     967!                      *RLMLT/RCPD/(1.0+RVTMP2*zq(i))
    900968!           print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
    901        ENDDO
    902      ENDIF
     969!       ENDDO
     970!     ENDIF
    903971
    904972       
     
    10031071  ! Pluie ou neige au sol selon la temperature de la 1ere couche
    10041072  !
     1073!CR: si la thermo de la glace est active, on calcule zifl directement
     1074  IF (.NOT.ice_thermo) THEN
    10051075  DO i = 1, klon
    10061076     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
    10071077!AJ<
    1008 !!        snow(i) = zrfl(i)
     1078!        snow(i) = zrfl(i)
    10091079        snow(i) = zrfl(i)+zifl(i)
    10101080!>AJ
     
    10151085     ENDIF
    10161086  ENDDO
     1087
     1088  ELSE
     1089     DO i = 1, klon
     1090        snow(i) = zifl(i)
     1091        rain(i) = zrfl(i)
     1092     ENDDO
     1093   
     1094   ENDIF
    10171095  !
    10181096  ! For energy conservation : when snow is present, the solification
    10191097  ! latent heat is considered.
     1098!CR: si thermo de la glace, neige deja prise en compte
     1099  IF (.not.ice_thermo) THEN
    10201100  DO k = 1, klev
    10211101     DO i = 1, klon
     
    10261106     END DO
    10271107  END DO
     1108  ENDIF
    10281109  !
    10291110
  • LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90

    r2003 r2086  
    3333      REAL, SAVE, ALLOCATABLE :: d_t_wake(:,:),d_q_wake(:,:)
    3434      !$OMP THREADPRIVATE( d_t_wake,d_q_wake)
    35       REAL, SAVE, ALLOCATABLE :: d_t_lsc(:,:),d_q_lsc(:,:),d_ql_lsc(:,:)
     35      REAL, SAVE, ALLOCATABLE :: d_t_lsc(:,:),d_q_lsc(:,:),d_ql_lsc(:,:),d_qi_lsc(:,:)
    3636      !$OMP THREADPRIVATE(d_t_lsc,d_q_lsc,d_ql_lsc)
    3737      REAL, SAVE, ALLOCATABLE :: d_t_ajsb(:,:), d_q_ajsb(:,:)
     
    307307      allocate(d_t_wake(klon,klev),d_q_wake(klon,klev))
    308308      allocate(d_t_lsc(klon,klev),d_q_lsc(klon,klev))
    309       allocate(d_ql_lsc(klon,klev))
     309      allocate(d_ql_lsc(klon,klev),d_qi_lsc(klon,klev))
    310310      allocate(d_t_ajsb(klon,klev),d_q_ajsb(klon,klev))
    311311      allocate(d_t_ajs(klon,klev),d_q_ajs(klon,klev))
     
    467467      deallocate(d_t_wake,d_q_wake)
    468468      deallocate(d_t_lsc,d_q_lsc)
    469       deallocate(d_ql_lsc)
     469      deallocate(d_ql_lsc,d_qi_lsc)
    470470      deallocate(d_t_ajsb,d_q_ajsb)
    471471      deallocate(d_t_ajs,d_q_ajs)
  • LMDZ5/trunk/libf/phylmd/phys_output_mod.F90

    r2037 r2086  
    399399#endif
    400400
    401       IF (nqtot>=3) THEN
    402             DO iq=3,nqtot 
     401!CR: ajout d'une variable eau
     402!      IF (nqtot>=3) THEN
     403       IF (nqtot>=nqo+1) THEN
     404!            DO iq=3,nqtot 
     405            DO iq=nqo+1,nqtot
    403406            iiq=niadv(iq)
    404             o_trac(iq-2) = ctrl_out((/ 4, 5, 1, 1, 1, 10, 11, 11, 11 /), &
     407            o_trac(iq-nqo) = ctrl_out((/ 4, 5, 1, 1, 1, 10, 11, 11, 11 /), &
    405408                           tname(iiq),'Tracer '//ttext(iiq), "-",  &
    406409                           (/ '', '', '', '', '', '', '', '', '' /))
    407410
    408             o_dtr_vdf(iq-2) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     411            o_dtr_vdf(iq-nqo) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    409412                              'd'//trim(tname(iq))//'_vdf',  &
    410413                              'Tendance tracer '//ttext(iiq), "-" , &
    411414                              (/ '', '', '', '', '', '', '', '', '' /))
    412415
    413             o_dtr_the(iq-2) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     416            o_dtr_the(iq-nqo) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    414417                              'd'//trim(tname(iq))//'_the', &
    415418                              'Tendance tracer '//ttext(iiq), "-", &
    416419                              (/ '', '', '', '', '', '', '', '', '' /))
    417420
    418             o_dtr_con(iq-2) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     421            o_dtr_con(iq-nqo) = ctrl_out((/ 5, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    419422                              'd'//trim(tname(iq))//'_con', &
    420423                              'Tendance tracer '//ttext(iiq), "-", &
    421424                              (/ '', '', '', '', '', '', '', '', '' /))
    422425
    423             o_dtr_lessi_impa(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     426            o_dtr_lessi_impa(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    424427                                     'd'//trim(tname(iq))//'_lessi_impa', &
    425428                                     'Tendance tracer '//ttext(iiq), "-", &
    426429                                     (/ '', '', '', '', '', '', '', '', '' /))
    427430
    428             o_dtr_lessi_nucl(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     431            o_dtr_lessi_nucl(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    429432                                     'd'//trim(tname(iq))//'_lessi_nucl', &
    430433                                     'Tendance tracer '//ttext(iiq), "-", &
    431434                                     (/ '', '', '', '', '', '', '', '', '' /))
    432435
    433             o_dtr_insc(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     436            o_dtr_insc(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    434437                               'd'//trim(tname(iq))//'_insc', &
    435438                               'Tendance tracer '//ttext(iiq), "-", &
    436439                               (/ '', '', '', '', '', '', '', '', '' /))
    437440
    438             o_dtr_bcscav(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     441            o_dtr_bcscav(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    439442                                 'd'//trim(tname(iq))//'_bcscav', &
    440443                                 'Tendance tracer '//ttext(iiq), "-", &
    441444                                 (/ '', '', '', '', '', '', '', '', '' /))
    442445
    443             o_dtr_evapls(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     446            o_dtr_evapls(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    444447                                 'd'//trim(tname(iq))//'_evapls', &
    445448                                 'Tendance tracer '//ttext(iiq), "-", &
    446449                                 (/ '', '', '', '', '', '', '', '', '' /))
    447450
    448             o_dtr_ls(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     451            o_dtr_ls(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    449452                             'd'//trim(tname(iq))//'_ls', &
    450453                             'Tendance tracer '//ttext(iiq), "-", &
    451454                             (/ '', '', '', '', '', '', '', '', '' /))
    452455
    453             o_dtr_trsp(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     456            o_dtr_trsp(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    454457                               'd'//trim(tname(iq))//'_trsp', &
    455458                               'Tendance tracer '//ttext(iiq), "-", &
    456459                               (/ '', '', '', '', '', '', '', '', '' /))
    457460
    458             o_dtr_sscav(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     461            o_dtr_sscav(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    459462                                'd'//trim(tname(iq))//'_sscav', &
    460463                                'Tendance tracer '//ttext(iiq), "-", &
    461464                                (/ '', '', '', '', '', '', '', '', '' /))
    462465
    463             o_dtr_sat(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     466            o_dtr_sat(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    464467                               'd'//trim(tname(iq))//'_sat', &
    465468                               'Tendance tracer '//ttext(iiq), "-", &
    466469                               (/ '', '', '', '', '', '', '', '', '' /))
    467470
    468             o_dtr_uscav(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     471            o_dtr_uscav(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    469472                                'd'//trim(tname(iq))//'_uscav', &
    470473                                'Tendance tracer '//ttext(iiq), "-", &
    471474                                 (/ '', '', '', '', '', '', '', '', '' /))
    472475
    473             o_dtr_dry(iq-2) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
     476            o_dtr_dry(iq-nqo) = ctrl_out((/ 7, 7, 7, 7, 10, 10, 11, 11, 11 /), &
    474477                              'cum'//'d'//trim(tname(iq))//'_dry', &
    475478                              'tracer tendency dry deposition'//ttext(iiq), "-", &
    476479                              (/ '', '', '', '', '', '', '', '', '' /))
    477480
    478             o_trac_cum(iq-2) = ctrl_out((/ 3, 4, 10, 10, 10, 10, 11, 11, 11 /), &
     481            o_trac_cum(iq-nqo) = ctrl_out((/ 3, 4, 10, 10, 10, 10, 11, 11, 11 /), &
    479482                               'cum'//tname(iiq),&
    480483                               'Cumulated tracer '//ttext(iiq), "-", &
  • LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90

    r2057 r2086  
    222222    USE ocean_slab_mod, only: tslab, slab_bils
    223223    USE indice_sol_mod, only: nbsrf
    224     USE infotrac, only: nqtot
     224    USE infotrac, only: nqtot, nqo
    225225    USE comgeomphy, only: airephy
    226226    USE surface_data, only: type_ocean, ok_veget, ok_snow
     
    12881288       ENDDO !nfiles
    12891289!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1290        IF (nqtot.GE.3) THEN
    1291           DO iq=3,nqtot
    1292              CALL histwrite_phy(o_trac(iq-2), qx(:,:,iq))
    1293              CALL histwrite_phy(o_dtr_vdf(iq-2),d_tr_cl(:,:,iq-2))
    1294              CALL histwrite_phy(o_dtr_the(iq-2),d_tr_th(:,:,iq-2))
    1295              CALL histwrite_phy(o_dtr_con(iq-2),d_tr_cv(:,:,iq-2))
    1296              CALL histwrite_phy(o_dtr_lessi_impa(iq-2),d_tr_lessi_impa(:,:,iq-2))
    1297              CALL histwrite_phy(o_dtr_lessi_nucl(iq-2),d_tr_lessi_nucl(:,:,iq-2))
    1298              CALL histwrite_phy(o_dtr_insc(iq-2),d_tr_insc(:,:,iq-2))
    1299              CALL histwrite_phy(o_dtr_bcscav(iq-2),d_tr_bcscav(:,:,iq-2))
    1300              CALL histwrite_phy(o_dtr_evapls(iq-2),d_tr_evapls(:,:,iq-2))
    1301              CALL histwrite_phy(o_dtr_ls(iq-2),d_tr_ls(:,:,iq-2))
    1302              CALL histwrite_phy(o_dtr_trsp(iq-2),d_tr_trsp(:,:,iq-2))
    1303              CALL histwrite_phy(o_dtr_sscav(iq-2),d_tr_sscav(:,:,iq-2))
    1304              CALL histwrite_phy(o_dtr_sat(iq-2),d_tr_sat(:,:,iq-2))
    1305              CALL histwrite_phy(o_dtr_uscav(iq-2),d_tr_uscav(:,:,iq-2))
     1290!       IF (nqtot.GE.3) THEN
     1291!          DO iq=3,nqtot
     1292          IF (nqtot.GE.nqo+1) THEN
     1293            DO iq=nqo+1,nqtot
     1294             CALL histwrite_phy(o_trac(iq-nqo), qx(:,:,iq))
     1295             CALL histwrite_phy(o_dtr_vdf(iq-nqo),d_tr_cl(:,:,iq-nqo))
     1296             CALL histwrite_phy(o_dtr_the(iq-nqo),d_tr_th(:,:,iq-nqo))
     1297             CALL histwrite_phy(o_dtr_con(iq-nqo),d_tr_cv(:,:,iq-nqo))
     1298             CALL histwrite_phy(o_dtr_lessi_impa(iq-nqo),d_tr_lessi_impa(:,:,iq-nqo))
     1299             CALL histwrite_phy(o_dtr_lessi_nucl(iq-nqo),d_tr_lessi_nucl(:,:,iq-nqo))
     1300             CALL histwrite_phy(o_dtr_insc(iq-nqo),d_tr_insc(:,:,iq-nqo))
     1301             CALL histwrite_phy(o_dtr_bcscav(iq-nqo),d_tr_bcscav(:,:,iq-nqo))
     1302             CALL histwrite_phy(o_dtr_evapls(iq-nqo),d_tr_evapls(:,:,iq-nqo))
     1303             CALL histwrite_phy(o_dtr_ls(iq-nqo),d_tr_ls(:,:,iq-nqo))
     1304             CALL histwrite_phy(o_dtr_trsp(iq-nqo),d_tr_trsp(:,:,iq-nqo))
     1305             CALL histwrite_phy(o_dtr_sscav(iq-nqo),d_tr_sscav(:,:,iq-nqo))
     1306             CALL histwrite_phy(o_dtr_sat(iq-nqo),d_tr_sat(:,:,iq-nqo))
     1307             CALL histwrite_phy(o_dtr_uscav(iq-nqo),d_tr_uscav(:,:,iq-nqo))
    13061308             zx_tmp_fi2d=0.
    13071309             IF(vars_defined) THEN
     
    13101312                ENDDO
    13111313             ENDIF
    1312              CALL histwrite_phy(o_trac_cum(iq-2), zx_tmp_fi2d)
     1314             CALL histwrite_phy(o_trac_cum(iq-nqo), zx_tmp_fi2d)
    13131315          ENDDO
    13141316       ENDIF
  • LMDZ5/trunk/libf/phylmd/physiq.F90

    r2076 r2086  
    201201  INTEGER iliq          ! indice de traceurs pour eau liquide
    202202  PARAMETER (iliq=2)
    203 
     203!CR: on ajoute la phase glace
     204  INTEGER isol          ! indice de traceurs pour eau glace
     205  PARAMETER (isol=3)
    204206  !
    205207  !
     
    596598
    597599  ! tendance nulles
    598   REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0
     600  REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
    599601
    600602  !
     
    13771379  dq0(:,:)=0.
    13781380  dql0(:,:)=0.
     1381  dqi0(:,:)=0.
    13791382  !
    13801383  ! Mettre a zero des variables de sortie (pour securite)
     
    14221425        q_seri(i,k)  = qx(i,k,ivap)
    14231426        ql_seri(i,k) = qx(i,k,iliq)
    1424         qs_seri(i,k) = 0.
     1427!CR: ATTENTION, on rajoute la variable glace
     1428        if (nqo.eq.2) then
     1429           qs_seri(i,k) = 0.
     1430        else if (nqo.eq.3) then
     1431           qs_seri(i,k) = qx(i,k,isol)
     1432        endif
    14251433     ENDDO
    14261434  ENDDO
    14271435  tke0(:,:)=pbl_tke(:,:,is_ave)
    1428   IF (nqtot.GE.3) THEN
    1429      DO iq = 3, nqtot
     1436!CR:Nombre de traceurs de l'eau: nqo
     1437!  IF (nqtot.GE.3) THEN
     1438   IF (nqtot.GE.(nqo+1)) THEN
     1439!     DO iq = 3, nqtot       
     1440     DO iq = nqo+1, nqtot 
    14301441        DO  k = 1, klev
    14311442           DO  i = 1, klon
    1432               tr_seri(i,k,iq-2) = qx(i,k,iq)
     1443!              tr_seri(i,k,iq-2) = qx(i,k,iq)
     1444              tr_seri(i,k,iq-nqo) = qx(i,k,iq)
    14331445           ENDDO
    14341446        ENDDO
     
    16171629        ENDIF
    16181630        !>jyg
     1631     
     1632        if (iflag_ice_thermo.eq.0) then   
     1633!pas necessaire a priori
    16191634
    16201635        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
     
    16271642        d_t_eva(i,k) = za
    16281643        d_q_eva(i,k) = zb
     1644
     1645        else
     1646
     1647!CR: on ré-évapore eau liquide et glace
     1648
     1649!        zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
     1650!        zb = MAX(0.0,ql_seri(i,k))
     1651!        za = - MAX(0.0,ql_seri(i,k)) &
     1652!             * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
     1653        zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
     1654        za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
     1655             - MAX(0.0,qs_seri(i,k))*zlsdcp
     1656        t_seri(i,k) = t_seri(i,k) + za
     1657        q_seri(i,k) = q_seri(i,k) + zb
     1658        ql_seri(i,k) = 0.0
     1659!on évapore la glace
     1660        qs_seri(i,k) = 0.0
     1661        d_t_eva(i,k) = za
     1662        d_q_eva(i,k) = zb
     1663        endif
     1664
    16291665     ENDDO
    16301666  ENDDO
     
    17641800     IF (klon_glo==1) THEN
    17651801        CALL add_pbl_tend &
    1766              (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
     1802             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
    17671803     ELSE
    17681804        CALL add_phys_tend &
    1769              (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,paprs,'vdf')
     1805             (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
    17701806     ENDIF
    17711807     !--------------------------------------------------------------------
     
    21352171  !-----------------------------------------------------------------------------------------
    21362172  ! ajout des tendances de la diffusion turbulente
    2137   CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,paprs,'con')
     2173  CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,dqi0,paprs,'con')
    21382174  !-----------------------------------------------------------------------------------------
    21392175
     
    22522288     d_t_wake(:,:)=dt_wake(:,:)*dtime
    22532289     d_q_wake(:,:)=dq_wake(:,:)*dtime
    2254      CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,paprs,'wake')
     2290     CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake')
    22552291     !-----------------------------------------------------------------------------------------
    22562292
     
    25122548        !-----------------------------------------------------------------------------------------
    25132549        ! ajout des tendances de l'ajustement sec ou des thermiques
    2514         CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,paprs,'ajsb')
     2550        CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs,'ajsb')
    25152551        d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
    25162552        d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
     
    25582594  CALL fisrtilp(dtime,paprs,pplay, &
    25592595       t_seri, q_seri,ptconv,ratqs, &
    2560        d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
     2596       d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
    25612597       rain_lsc, snow_lsc, &
    25622598       pfrac_impa, pfrac_nucl, pfrac_1nucl, &
     
    25702606  !-----------------------------------------------------------------------------------------
    25712607  ! ajout des tendances de la diffusion turbulente
    2572   CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,paprs,'lsc')
     2608  CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs,'lsc')
    25732609  !-----------------------------------------------------------------------------------------
    25742610  DO k = 1, klev
    25752611     DO i = 1, klon
    25762612        cldfra(i,k) = rneb(i,k)
     2613!CR: a quoi ca sert? Faut-il ajouter qs_seri?
    25772614        IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
    25782615     ENDDO
     
    28452882        zx_t = t_seri(i,k)
    28462883        IF (thermcep) THEN
     2884           if (iflag_ice_thermo.eq.0) then
    28472885           zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
     2886           else
     2887           zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
     2888           endif
    28482889           zx_qs  = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
    28492890           zx_qs  = MIN(0.5,zx_qs)
     
    32883329     !-----------------------------------------------------------------------------------------
    32893330     ! ajout des tendances de la trainee de l'orographie
    3290      CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,paprs,'oro')
     3331     CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro')
    32913332     !-----------------------------------------------------------------------------------------
    32923333     !
     
    33343375     !-----------------------------------------------------------------------------------------
    33353376     ! ajout des tendances de la portance de l'orographie
    3336      CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,paprs,'lif')
     3377     CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif')
    33373378     !-----------------------------------------------------------------------------------------
    33383379     !
     
    33483389     !
    33493390     !  ajout des tendances
    3350      CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,paprs,'hin')
     3391     CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin')
    33513392
    33523393  ENDIF
     
    33563397          rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
    33573398          du_gwd_rando, dv_gwd_rando)
    3358      CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,paprs, &
     3399     CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
    33593400          'flott_gwd_rando')
    33603401  end if
     
    36233664        d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
    36243665        d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
     3666!CR: on ajoute le contenu en glace
     3667        if (nqo.eq.3) then
     3668        d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
     3669        endif
    36253670     ENDDO
    36263671  ENDDO
    36273672  !
    3628   IF (nqtot.GE.3) THEN
    3629      DO iq = 3, nqtot
     3673!CR: nb de traceurs eau: nqo
     3674!  IF (nqtot.GE.3) THEN
     3675   IF (nqtot.GE.(nqo+1)) THEN
     3676!     DO iq = 3, nqtot
     3677     DO iq = nqo+1, nqtot
    36303678        DO  k = 1, klev
    36313679           DO  i = 1, klon
    3632               d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
     3680!              d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
     3681               d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
    36333682           ENDDO
    36343683        ENDDO
     
    36533702
    36543703!!! RomP >>>
    3655   IF (nqtot.GE.3) THEN
    3656      DO iq = 3, nqtot
     3704!CR: nb de traceurs eau: nqo
     3705!  IF (nqtot.GE.3) THEN
     3706   IF (nqtot.GE.(nqo+1)) THEN
     3707!     DO iq = 3, nqtot
     3708     DO iq = nqo+1, nqtot
    36573709        DO k = 1, klev
    36583710           DO i = 1, klon
    3659               tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
     3711!              tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
     3712              tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
    36603713           ENDDO
    36613714        ENDDO
Note: See TracChangeset for help on using the changeset viewer.