Ignore:
Timestamp:
Jan 15, 2016, 8:27:16 AM (9 years ago)
Author:
emillour
Message:

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

  • in dynlonlat_phylonlat:
  • add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)
  • in misc:
  • follow up updates on wxios.F (add missing_val module variable)
  • in dyn3d_common:
  • pression.F => pression.F90
  • misc_mod.F90: moved from misc to dyn3d_common
  • added new iso_verif_dyn.F
  • covcont.F => covcont.F90
  • infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
  • dynetat0.F => dynetat0.F90 with some code factorization
  • dynredem.F => dynredem.F90 with some code factorization
  • added dynredem_mod.F90: routines used by dynredem
  • iniacademic.F90 : added isotopes-related initialization for Earth case
  • in dyn3d:
  • added check_isotopes.F
  • modified (isotopes) advtrac.F90, caladvtrac.F
  • guide_mod.F90: ported updates
  • leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
  • qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
  • vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
  • vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)
  • in dyn3dpar:
  • leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
  • dynredem_p.F => dynredem_p.F90 and some code factorization
  • and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

EM

Location:
trunk/LMDZ.COMMON/libf/dyn3d
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3d/advtrac.F90

    r1441 r1508  
    99  !            M.A Filiberti (04/2002)
    1010  !
    11   USE infotrac, ONLY: nqtot, iadv
     11  USE infotrac, ONLY: nqtot, iadv, nqperes, ok_iso_verif
    1212  USE control_mod, ONLY: iapp_tracvl, day_step
    1313  USE comconst_mod, ONLY: dtvr
    … …  
    218218     !     Appel des sous programmes d'advection
    219219     !-----------------------------------------------------------
    220      do iq=1,nqtot
     220     if (ok_iso_verif) then
     221           write(*,*) 'advtrac 227'
     222           call check_isotopes_seq(q,ip1jmp1,'advtrac 162')
     223     endif !if (ok_iso_verif) then
     224
     225     do iq=1,nqperes
    221226        !        call clock(t_initial)
    222227        if(iadv(iq) == 0) cycle
    … …  
    225230        !   ----------------------------------------------------------------
    226231        if(iadv(iq).eq.10) THEN
    227            call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
     232           ! CRisi: on fait passer tout q pour avoir acces aux fils
     233           
     234           !write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
     235           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
    228236
    229237           !   ----------------------------------------------------------------
    … …  
    233241        else if(iadv(iq).eq.14) then
    234242           !
    235            CALL vlspltqs( q(1,1,1), 2., massem, wg , &
    236                 pbarug,pbarvg,dtvr,p,pk,teta )
     243           !write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
     244           CALL vlspltqs( q, 2., massem, wg , &
     245                pbarug,pbarvg,dtvr,p,pk,teta,iq)
     246
    237247           !   ----------------------------------------------------------------
    238248           !   Schema de Frederic Hourdin
    … …  
    383393     end DO
    384394
     395     if (ok_iso_verif) then
     396           write(*,*) 'advtrac 402'
     397           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
     398     endif !if (ok_iso_verif) then
    385399
    386400     !------------------------------------------------------------------

  • trunk/LMDZ.COMMON/libf/dyn3d/caladvtrac.F

    r1422 r1508  
    5353      if (planet_type.eq."earth") then
    5454C initialisation
    55         dq(:,:,1:2)=q(:,:,1:2)
     55! CRisi: il faut gérer tous les traceurs si on veut pouvoir faire des
     56! isotopes
     57!        dq(:,:,1:2)=q(:,:,1:2)
     58        dq(:,:,1:nqtot)=q(:,:,1:nqtot)
    5659       
    5760c  test des valeurs minmax
    … …  
    8285          ENDDO
    8386         
    84           CALL qminimum( q, 2, finmasse )
     87          !write(*,*) 'caladvtrac 87'
     88          CALL qminimum( q, nqtot, finmasse )
     89          !write(*,*) 'caladvtrac 89'
    8590
    8691          CALL SCOPY   ( ip1jmp1*llm, masse, 1, finmasse,       1 )
    … …  
    9297          dtvrtrac = iapp_tracvl * dtvr
    9398c
    94            DO iq = 1 , 2
     99           DO iq = 1 , nqtot
    95100            DO l = 1 , llm
    96101             DO ij = 1,ip1jmp1
    … …  
    105110        if (planet_type.eq."earth") then
    106111! Earth-specific treatment for the first 2 tracers (water)
    107           dq(:,:,1:2)=0.
     112          dq(:,:,1:nqtot)=0.
    108113        endif ! of if (planet_type.eq."earth")
    109114      ENDIF ! of IF( iapptrac.EQ.iapp_tracvl )

  • trunk/LMDZ.COMMON/libf/dyn3d/guide_mod.F90

    r1441 r1508  
    7979! Lecture des parametres: 
    8080! ---------------------------------------------
     81    call ini_getparam("nudging_parameters_out.txt")
    8182! Variables guidees
    8283    CALL getpar('guide_u',.true.,guide_u,'guidage de u')
    … …  
    144145    CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
    145146    CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
     147
     148    call fin_getparam
    146149
    147150! ---------------------------------------------

  • trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F

    r1422 r1508  
    1212      use IOIPSL
    1313#endif
    14       USE infotrac, ONLY: nqtot
     14      USE infotrac, ONLY: nqtot,ok_iso_verif
    1515      USE guide_mod, ONLY : guide_main
    1616      USE write_field, ONLY: writefield
    … …  
    305305      jH_cur = jH_cur - int(jH_cur)
    306306
     307        if (ok_iso_verif) then
     308           call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
     309        endif !if (ok_iso_verif) then
    307310
    308311#ifdef CPP_IOIPSL
    … …  
    337340!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
    338341!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
     342
     343        if (ok_iso_verif) then
     344           call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
     345        endif !if (ok_iso_verif) then
    339346
    340347   2  CONTINUE ! Matsuno backward or leapfrog step begins here
    … …  
    381388#endif
    382389
     390        if (ok_iso_verif) then
     391           call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
     392        endif !if (ok_iso_verif) then
     393
    383394c-----------------------------------------------------------------------
    384395c   calcul des tendances dynamiques:
    … …  
    419430c   -------------------------------------------------------------
    420431
     432        if (ok_iso_verif) then
     433           call check_isotopes_seq(q,ip1jmp1,
     434     &           'leapfrog 686: avant caladvtrac')
     435        endif !if (ok_iso_verif) then
     436
    421437      IF( forward. OR . leapf )  THEN
    422438! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
    … …  
    443459c   ----------------------------------
    444460
    445 
    446        CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
     461        if (ok_iso_verif) then
     462           write(*,*) 'leapfrog 720'
     463           call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
     464        endif !if (ok_iso_verif) then
     465       
     466       CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
    447467     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
    448468!     $              finvmaold                                    )
     469
     470       if (ok_iso_verif) then
     471          write(*,*) 'leapfrog 724'
     472           call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
     473        endif !if (ok_iso_verif) then
    449474
    450475       IF ((planet_type.eq."titan").and.(tidal)) then
    … …  
    519544         IF (ip_ebil_dyn.ge.1 ) THEN
    520545          ztit='bil dyn'
    521 ! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
     546! Ehouarn: be careful, diagedyn is Earth-specific!
    522547           IF (planet_type.eq."earth") THEN
    523548            CALL diagedyn(ztit,2,1,1,dtphys
    … …  
    619644        CALL massdair(p,masse)
    620645
     646        if (ok_iso_verif) then
     647           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
     648        endif !if (ok_iso_verif) then
     649
    621650c-----------------------------------------------------------------------
    622651c   dissipation horizontale et verticale  des petites echelles:
    … …  
    717746c   preparation du pas d'integration suivant  ......
    718747
     748        if (ok_iso_verif) then
     749           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
     750        endif !if (ok_iso_verif) then
     751
    719752      IF ( .NOT.purmats ) THEN
    720753c       ........................................................
    … …  
    776809
    777810            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
     811
     812        if (ok_iso_verif) then
     813           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
     814        endif !if (ok_iso_verif) then
    778815
    779816c-----------------------------------------------------------------------
    … …  
    858895      ELSE ! of IF (.not.purmats)
    859896
     897        if (ok_iso_verif) then
     898           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
     899        endif !if (ok_iso_verif) then
     900
    860901c       ........................................................
    861902c       ..............       schema  matsuno        ...............
    … …  
    880921
    881922            ELSE ! of IF(forward) i.e. backward step
     923
     924        if (ok_iso_verif) then
     925           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
     926        endif !if (ok_iso_verif) then 
    882927
    883928              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN

  • trunk/LMDZ.COMMON/libf/dyn3d/qminimum.F

    r1422 r1508  
    22! $Header$
    33!
    4       SUBROUTINE qminimum( q,nq,deltap )
     4      SUBROUTINE qminimum( q,nqtot,deltap )
    55
     6      USE infotrac, ONLY: ok_isotopes,ntraciso,iqiso,ok_iso_verif
    67      IMPLICIT none
    78c
    … …  
    1213#include "paramet.h"
    1314c
    14       INTEGER nq
    15       REAL q(ip1jmp1,llm,nq), deltap(ip1jmp1,llm)
     15      INTEGER nqtot
     16      REAL q(ip1jmp1,llm,nqtot), deltap(ip1jmp1,llm)
    1617c
    1718      INTEGER iq_vap, iq_liq
    … …  
    2930      INTEGER i, k, iq
    3031      REAL zx_defau, zx_abc, zx_pump(ip1jmp1), pompe
     32
     33      real zx_defau_diag(ip1jmp1,llm,2)
     34      real q_follow(ip1jmp1,llm,2)
    3135c
    3236      REAL SSUM
    … …  
    3539      SAVE imprim
    3640      DATA imprim /0/
     41      !INTEGER ijb,ije
     42      !INTEGER Index_pump(ij_end-ij_begin+1)
     43      !INTEGER nb_pump
     44      INTEGER ixt
    3745c
    3846c Quand l'eau liquide est trop petite (ou negative), on prend
    … …  
    4048c (sans changer la temperature !)
    4149c
     50
     51        if (ok_iso_verif) then
     52           call check_isotopes_seq(q,ip1jmp1,'qminimum 52')   
     53        endif !if (ok_iso_verif) then     
     54
     55      zx_defau_diag(:,:,:)=0.0
     56      q_follow(:,:,1:2)=q(:,:,1:2) 
    4257      DO 1000 k = 1, llm
    4358        DO 1040 i = 1, ip1jmp1
    4459          if (seuil_liq - q(i,k,iq_liq) .gt. 0.d0 ) then
     60
     61              if (ok_isotopes) then
     62                 zx_defau_diag(i,k,iq_liq)=AMAX1
     63     :               ( seuil_liq - q(i,k,iq_liq), 0.0 )
     64              endif !if (ok_isotopes) then
     65
    4566             q(i,k,iq_vap) = q(i,k,iq_vap) + q(i,k,iq_liq) - seuil_liq
    4667             q(i,k,iq_liq) = seuil_liq
    … …  
    5879        DO i = 1, ip1jmp1
    5980          if ( seuil_vap - q(i,k,iq) .gt. 0.d0 ) then
     81
     82            if (ok_isotopes) then
     83              zx_defau_diag(i,k,iq)=AMAX1( seuil_vap - q(i,k,iq), 0.0 )
     84            endif !if (ok_isotopes) then
     85
    6086            q(i,k-1,iq) =  q(i,k-1,iq) - ( seuil_vap - q(i,k,iq) ) *
    6187     &                     deltap(i,k) / deltap(i,k-1)
    … …  
    82108         ENDDO
    83109      ENDIF
     110
     111      !write(*,*) 'qminimum 128'
     112      if (ok_isotopes) then
     113      ! CRisi: traiter de même les traceurs d'eau
     114      ! Mais il faut les prendre à l'envers pour essayer de conserver la
     115      ! masse.
     116      ! 1) pompage dans le sol 
     117      ! On suppose que ce pompage se fait sans isotopes -> on ne modifie
     118      ! rien ici et on croise les doigts pour que ça ne soit pas trop
     119      ! génant
     120      DO i = 1,ip1jmp1
     121        if (zx_pump(i).gt.0.0) then
     122          q_follow(i,1,iq_vap)=q_follow(i,1,iq_vap)+zx_pump(i)
     123        endif !if (zx_pump(i).gt.0.0) then
     124      enddo !DO i = 1,ip1jmp1
     125
     126      ! 2) transfert de vap vers les couches plus hautes
     127      !write(*,*) 'qminimum 139'
     128      do k=2,llm
     129        DO i = 1,ip1jmp1
     130          if (zx_defau_diag(i,k,iq_vap).gt.0.0) then             
     131              ! on ajoute la vapeur en k             
     132              do ixt=1,ntraciso
     133               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
     134     :              +zx_defau_diag(i,k,iq_vap)
     135     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
     136               
     137              ! et on la retranche en k-1
     138               q(i,k-1,iqiso(ixt,iq_vap))=q(i,k-1,iqiso(ixt,iq_vap))
     139     :              -zx_defau_diag(i,k,iq_vap)
     140     :              *deltap(i,k)/deltap(i,k-1)
     141     :              *q(i,k-1,iqiso(ixt,iq_vap))/q_follow(i,k-1,iq_vap)
     142
     143              enddo !do ixt=1,niso
     144              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
     145     :               +zx_defau_diag(i,k,iq_vap)
     146              q_follow(i,k-1,iq_vap)=   q_follow(i,k-1,iq_vap)
     147     :               -zx_defau_diag(i,k,iq_vap)
     148     :              *deltap(i,k)/deltap(i,k-1)
     149          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
     150        enddo !DO i = 1, ip1jmp1       
     151       enddo !do k=2,llm
     152
     153        if (ok_iso_verif) then     
     154           call check_isotopes_seq(q,ip1jmp1,'qminimum 168')
     155        endif !if (ok_iso_verif) then
     156       
     157     
     158        ! 3) transfert d'eau de la vapeur au liquide
     159        !write(*,*) 'qminimum 164'
     160        do k=1,llm
     161        DO i = 1,ip1jmp1
     162          if (zx_defau_diag(i,k,iq_liq).gt.0.0) then
     163
     164              ! on ajoute eau liquide en k en k             
     165              do ixt=1,ntraciso
     166               q(i,k,iqiso(ixt,iq_liq))=q(i,k,iqiso(ixt,iq_liq))
     167     :              +zx_defau_diag(i,k,iq_liq)
     168     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)
     169              ! et on la retranche à la vapeur en k
     170               q(i,k,iqiso(ixt,iq_vap))=q(i,k,iqiso(ixt,iq_vap))
     171     :              -zx_defau_diag(i,k,iq_liq)
     172     :              *q(i,k,iqiso(ixt,iq_vap))/q_follow(i,k,iq_vap)   
     173              enddo !do ixt=1,niso
     174              q_follow(i,k,iq_liq)=   q_follow(i,k,iq_liq)
     175     :               +zx_defau_diag(i,k,iq_liq)
     176              q_follow(i,k,iq_vap)=   q_follow(i,k,iq_vap)
     177     :               -zx_defau_diag(i,k,iq_liq)
     178          endif !if (zx_defau_diag(i,k,iq_vap).gt.0.0) then
     179        enddo !DO i = 1, ip1jmp1
     180       enddo !do k=2,llm 
     181
     182        if (ok_iso_verif) then
     183           call check_isotopes_seq(q,ip1jmp1,'qminimum 197')
     184        endif !if (ok_iso_verif) then
     185
     186      endif !if (ok_isotopes) then
     187      !write(*,*) 'qminimum 188'
     188     
    84189c
    85190      RETURN

  • trunk/LMDZ.COMMON/libf/dyn3d/vlsplt.F

    r1422 r1508  
    11c
    2 c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z lguez $
    3 c
    4 
    5       SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt)
     2c $Id: vlsplt.F 2286 2015-05-20 13:27:07Z emillour $
     3c
     4
     5      SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
     6      USE infotrac, ONLY: nqtot,nqdesc,iqfils
    67c
    78c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    2930c      REAL masse(iip1,jjp1,llm),pente_max
    3031      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    31       REAL q(ip1jmp1,llm)
     32      REAL q(ip1jmp1,llm,nqtot)
    3233c      REAL q(iip1,jjp1,llm)
    3334      REAL w(ip1jmp1,llm),pdt
     35      INTEGER iq ! CRisi
    3436c
    3537c      Local
    … …  
    3941      INTEGER ijlqmin,iqmin,jqmin,lqmin
    4042c
    41       REAL zm(ip1jmp1,llm),newmasse
     43      REAL zm(ip1jmp1,llm,nqtot),newmasse
    4244      REAL mu(ip1jmp1,llm)
    4345      REAL mv(ip1jm,llm)
    4446      REAL mw(ip1jmp1,llm+1)
    45       REAL zq(ip1jmp1,llm),zz
     47      REAL zq(ip1jmp1,llm,nqtot),zz
    4648      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
    4749      REAL second,temps0,temps1,temps2,temps3
    … …  
    5254      SAVE temps1,temps2,temps3
    5355      INTEGER iminn,imaxx
     56      INTEGER ifils,iq2 ! CRisi
    5457
    5558      REAL qmin,qmax
    … …  
    7679         mw(ij,llm+1)=0.
    7780      ENDDO
    78      
    79       CALL SCOPY(ijp1llm,q,1,zq,1)
    80       CALL SCOPY(ijp1llm,masse,1,zm,1)
     81           
     82      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
     83      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
     84       
     85      if (nqdesc(iq).gt.0) then 
     86        do ifils=1,nqdesc(iq)
     87          iq2=iqfils(ifils,iq)
     88          CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
     89        enddo 
     90      endif !if (nqfils(iq).gt.0) then
    8191
    8292cprint*,'Entree vlx1'
    8393c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
    84       call vlx(zq,pente_max,zm,mu)
     94      call vlx(zq,pente_max,zm,mu,iq)
    8595cprint*,'Sortie vlx1'
    8696c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
    8797
    8898c print*,'Entree vly1'
    89       call vly(zq,pente_max,zm,mv)
     99
     100      call vly(zq,pente_max,zm,mv,iq)
    90101c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
    91102cprint*,'Sortie vly1'
    92       call vlz(zq,pente_max,zm,mw)
     103      call vlz(zq,pente_max,zm,mw,iq)
    93104c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
    94105
    95106
    96       call vly(zq,pente_max,zm,mv)
     107      call vly(zq,pente_max,zm,mv,iq)
    97108c       call minmaxq(zq,qmin,qmax,'apres vly     ')
    98109
    99110
    100       call vlx(zq,pente_max,zm,mu)
     111      call vlx(zq,pente_max,zm,mu,iq)
    101112c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
    102113       
    … …  
    104115      DO l=1,llm
    105116         DO ij=1,ip1jmp1
    106            q(ij,l)=zq(ij,l)
     117           q(ij,l,iq)=zq(ij,l,iq)
    107118         ENDDO
    108119         DO ij=1,ip1jm+1,iip1
    109             q(ij+iim,l)=q(ij,l)
    110          ENDDO
    111       ENDDO
     120            q(ij+iim,l,iq)=q(ij,l,iq)
     121         ENDDO
     122      ENDDO
     123      ! CRisi: aussi pour les fils
     124      if (nqdesc(iq).gt.0) then
     125      do ifils=1,nqdesc(iq)
     126        iq2=iqfils(ifils,iq)
     127        DO l=1,llm
     128         DO ij=1,ip1jmp1
     129           q(ij,l,iq2)=zq(ij,l,iq2)
     130         ENDDO
     131         DO ij=1,ip1jm+1,iip1
     132            q(ij+iim,l,iq2)=q(ij,l,iq2)
     133         ENDDO
     134        ENDDO
     135      enddo !do ifils=1,nqdesc(iq)   
     136      endif ! if (nqdesc(iq).gt.0) then   
    112137
    113138      RETURN
    114139      END
    115       SUBROUTINE vlx(q,pente_max,masse,u_m)
     140      RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
     141      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    116142
    117143c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    133159c   Arguments:
    134160c   ----------
    135       REAL masse(ip1jmp1,llm),pente_max
     161      REAL masse(ip1jmp1,llm,nqtot),pente_max
    136162      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
    137       REAL q(ip1jmp1,llm)
     163      REAL q(ip1jmp1,llm,nqtot)
    138164      REAL w(ip1jmp1,llm)
     165      INTEGER iq ! CRisi
    139166c
    140167c      Local
    … …  
    149176      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    150177      REAL u_mq(ip1jmp1,llm)
     178
     179      ! CRisi
     180      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
     181      INTEGER ifils,iq2 ! CRisi
    151182
    152183      Logical extremum,first,testcpu
    … …  
    182213         DO l = 1, llm
    183214            DO ij=iip2,ip1jm-1
    184                dxqu(ij)=q(ij+1,l)-q(ij,l)
     215               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    185216c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    186 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     217c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    187218            ENDDO
    188219            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    237268         DO l = 1, llm
    238269            DO ij=iip2,ip1jm-1
    239                dxqu(ij)=q(ij+1,l)-q(ij,l)
     270               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    240271            ENDDO
    241272            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    279310      DO l=1,llm
    280311       DO ij=iip2,ip1jm-1
    281           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    282      ,                     1.+u_m(ij,l)/masse(ij+1,l),
     312          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     313     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
    283314     ,                     u_m(ij,l))
    284315          zdum(ij,l)=0.5*zdum(ij,l)
    285316          u_mq(ij,l)=cvmgp(
    286      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    287      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     317     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     318     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    288319     ,                u_m(ij,l))
    289320          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    297328      DO l=1,llm
    298329       DO ij=iip2,ip1jm-1
    299 c       print*,'masse(',ij,')=',masse(ij,l)
     330c       print*,'masse(',ij,')=',masse(ij,l,iq)
    300331          IF (u_m(ij,l).gt.0.) THEN
    301              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
    302              u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
     332             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
     333             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l))
    303334          ELSE
    304              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
    305              u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
     335             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
     336             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
     337     &           -0.5*zdum(ij,l)*dxq(ij+1,l))
    306338          ENDIF
    307339       ENDDO
    … …  
    373405                     i=ijq-(j-1)*iip1
    374406c   accumulation pour les mailles completements advectees
    375                      do while(zu_m.gt.masse(ijq,l))
    376                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    377                         zu_m=zu_m-masse(ijq,l)
     407                     do while(zu_m.gt.masse(ijq,l,iq))
     408                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
     409     &                          *masse(ijq,l,iq)
     410                        zu_m=zu_m-masse(ijq,l,iq)
    378411                        i=mod(i-2+iim,iim)+1
    379412                        ijq=(j-1)*iip1+i
    … …  
    381414c   ajout de la maille non completement advectee
    382415                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    383      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     416     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
     417     &                  *dxq(ijq,l))
    384418                  ELSE
    385419                     ijq=ij+1
    386420                     i=ijq-(j-1)*iip1
    387421c   accumulation pour les mailles completements advectees
    388                      do while(-zu_m.gt.masse(ijq,l))
    389                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    390                         zu_m=zu_m+masse(ijq,l)
     422                     do while(-zu_m.gt.masse(ijq,l,iq))
     423                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     424     &                          *masse(ijq,l,iq)
     425                        zu_m=zu_m+masse(ijq,l,iq)
    391426                        i=mod(i,iim)+1
    392427                        ijq=(j-1)*iip1+i
    393428                     ENDDO
    394429c   ajout de la maille non completement advectee
    395                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    396      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     430                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     431     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    397432                  ENDIF
    398433               ENDDO
    … …  
    411446      ENDDO
    412447
     448! CRisi: appel récursif de l'advection sur les fils.
     449! Il faut faire ça avant d'avoir mis à jour q et masse
     450      !write(*,*) 'vlsplt 326: iq,nqfils(iq)=',iq,nqfils(iq)
     451     
     452      if (nqdesc(iq).gt.0) then 
     453       do ifils=1,nqdesc(iq)
     454         iq2=iqfils(ifils,iq)
     455         DO l=1,llm
     456          DO ij=iip2,ip1jm
     457           ! On a besoin de q et masse seulement entre iip2 et ip1jm
     458           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     459           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     460          enddo   
     461         enddo
     462        enddo !do ifils=1,nqdesc(iq)
     463        do ifils=1,nqfils(iq)
     464         iq2=iqfils(ifils,iq)
     465         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
     466        enddo !do ifils=1,nqfils(iq)
     467      endif !if (nqfils(iq).gt.0) then
     468! end CRisi
     469
    413470
    414471c   calcul des tENDances
    … …  
    416473      DO l=1,llm
    417474         DO ij=iip2+1,ip1jm
    418             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    419             q(ij,l)=(q(ij,l)*masse(ij,l)+
     475            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     476            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    420477     &      u_mq(ij-1,l)-u_mq(ij,l))
    421478     &      /new_m
    422             masse(ij,l)=new_m
     479            masse(ij,l,iq)=new_m
    423480         ENDDO
    424481c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    425482         DO ij=iip1+iip1,ip1jm,iip1
    426             q(ij-iim,l)=q(ij,l)
    427             masse(ij-iim,l)=masse(ij,l)
    428          ENDDO
    429       ENDDO
     483            q(ij-iim,l,iq)=q(ij,l,iq)
     484            masse(ij-iim,l,iq)=masse(ij,l,iq)
     485         ENDDO
     486      ENDDO
     487
     488      ! retablir les fils en rapport de melange par rapport a l'air:
     489      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
     490      ! puis on boucle en longitude
     491      if (nqdesc(iq).gt.0) then 
     492       do ifils=1,nqdesc(iq)
     493         iq2=iqfils(ifils,iq) 
     494         DO l=1,llm
     495          DO ij=iip2+1,ip1jm
     496            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     497          enddo
     498          DO ij=iip1+iip1,ip1jm,iip1
     499             q(ij-iim,l,iq2)=q(ij,l,iq2)
     500          enddo ! DO ij=ijb+iip1-1,ije,iip1
     501         enddo !DO l=1,llm
     502        enddo !do ifils=1,nqdesc(iq)
     503      endif !if (nqfils(iq).gt.0) then
     504
    430505c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    431506c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
    … …  
    434509      RETURN
    435510      END
    436       SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
     511      RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
     512      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
     513      USE comconst_mod, ONLY: pi
    437514c
    438515c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    446523c
    447524c   --------------------------------------------------------------------
    448       USE comconst_mod, ONLY: pi
    449 
    450525      IMPLICIT NONE
    451526c
    … …  
    457532c   Arguments:
    458533c   ----------
    459       REAL masse(ip1jmp1,llm),pente_max
     534      REAL masse(ip1jmp1,llm,nqtot),pente_max
    460535      REAL masse_adv_v( ip1jm,llm)
    461       REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
     536      REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm)
     537      INTEGER iq ! CRisi
    462538c
    463539c      Local
    … …  
    484560      SAVE sinlon,coslon,sinlondlon,coslondlon
    485561      SAVE airej2,airejjm
     562
     563      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     564      INTEGER ifils,iq2 ! CRisi
     565
    486566c
    487567c
    … …  
    490570      DATA first,testcpu/.true.,.false./
    491571      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
     572
     573      !write(*,*) 'vly 578: entree, iq=',iq
    492574
    493575      IF(first) THEN
    … …  
    522604
    523605      DO i = 1, iim
    524       airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
    525       airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     606      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
     607      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    526608      ENDDO
    527609      qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    531613
    532614      DO ij=1,ip1jm
    533          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     615         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    534616         adyqv(ij)=abs(dyqv(ij))
    535617      ENDDO
    … …  
    546628
    547629      DO ij=1,iip1
    548          dyq(ij,l)=qpns-q(ij+iip1,l)
    549          dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     630         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
     631         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    550632      ENDDO
    551633
    … …  
    668750      ENDDO
    669751
     752      !write(*,*) 'vly 756'
    670753      DO l=1,llm
    671754       DO ij=1,ip1jm
    672755          IF(masse_adv_v(ij,l).gt.0) THEN
    673               qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
    674      ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
     756              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
     757     ,                   0.5*(1.-masse_adv_v(ij,l)
     758     ,                   /masse(ij+iip1,l,iq))
    675759          ELSE
    676               qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
    677      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
     760              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
     761     ,                   0.5*(1.+masse_adv_v(ij,l)
     762     ,                   /masse(ij,l,iq))
    678763          ENDIF
    679764          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    681766      ENDDO
    682767
     768! CRisi: appel récursif de l'advection sur les fils.
     769! Il faut faire ça avant d'avoir mis à jour q et masse
     770      !write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq)
     771   
     772      if (nqfils(iq).gt.0) then 
     773       do ifils=1,nqdesc(iq)
     774         iq2=iqfils(ifils,iq)
     775         DO l=1,llm
     776         DO ij=1,ip1jmp1
     777           ! attention, chaque fils doit avoir son masseq, sinon, le 1er
     778           ! fils ecrase le masseq de ses freres.
     779           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     780           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     781          enddo   
     782         enddo
     783        enddo !do ifils=1,nqdesc(iq)
     784
     785        do ifils=1,nqfils(iq)
     786         iq2=iqfils(ifils,iq)
     787         call vly(Ratio,pente_max,masseq,qbyv,iq2)
     788        enddo !do ifils=1,nqfils(iq)
     789      endif !if (nqfils(iq).gt.0) then
    683790
    684791      DO l=1,llm
    685792         DO ij=iip2,ip1jm
    686             newmasse=masse(ij,l)
     793            newmasse=masse(ij,l,iq)
    687794     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    688             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    689      &         /newmasse
    690             masse(ij,l)=newmasse
     795            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
     796     &         -qbyv(ij-iip1,l))/newmasse
     797            masse(ij,l,iq)=newmasse
    691798         ENDDO
    692799c.-. ancienne version
    … …  
    696803         convpn=SSUM(iim,qbyv(1,l),1)
    697804         convmpn=ssum(iim,masse_adv_v(1,l),1)
    698          massepn=ssum(iim,masse(1,l),1)
     805         massepn=ssum(iim,masse(1,l,iq),1)
    699806         qpn=0.
    700807         do ij=1,iim
    701             qpn=qpn+masse(ij,l)*q(ij,l)
     808            qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
    702809         enddo
    703810         qpn=(qpn+convpn)/(massepn+convmpn)
    704811         do ij=1,iip1
    705             q(ij,l)=qpn
     812            q(ij,l,iq)=qpn
    706813         enddo
    707814
    … …  
    711818         convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    712819         convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    713          masseps=ssum(iim, masse(ip1jm+1,l),1)
     820         masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
    714821         qps=0.
    715822         do ij = ip1jm+1,ip1jmp1-1
    716             qps=qps+masse(ij,l)*q(ij,l)
     823            qps=qps+masse(ij,l,iq)*q(ij,l,iq)
    717824         enddo
    718825         qps=(qps+convps)/(masseps+convmps)
    719826         do ij=ip1jm+1,ip1jmp1
    720             q(ij,l)=qps
     827            q(ij,l,iq)=qps
    721828         enddo
    722829
    … …  
    732839c        DO ij = 1,iip1
    733840c           q(ij,l)=newq
    734 c           masse(ij,l)=newmasse*aire(ij)
     841c           masse(ij,l,iq)=newmasse*aire(ij)
    735842c        ENDDO
    736843c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    … …  
    742849c        DO ij = ip1jm+1,ip1jmp1
    743850c           q(ij,l)=newq
    744 c           masse(ij,l)=newmasse*aire(ij)
     851c           masse(ij,l,iq)=newmasse*aire(ij)
    745852c        ENDDO
    746853c._. fin nouvelle version
    747854      ENDDO
     855 
     856! retablir les fils en rapport de melange par rapport a l'air:
     857      if (nqfils(iq).gt.0) then 
     858       do ifils=1,nqdesc(iq)
     859         iq2=iqfils(ifils,iq) 
     860         DO l=1,llm
     861          DO ij=1,ip1jmp1
     862            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     863          enddo
     864         enddo
     865        enddo !do ifils=1,nqdesc(iq)
     866      endif !if (nqfils(iq).gt.0) then
     867
     868      !write(*,*) 'vly 853: sortie'
    748869
    749870      RETURN
    750871      END
    751       SUBROUTINE vlz(q,pente_max,masse,w)
     872      RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)
     873      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    752874c
    753875c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    769891c   Arguments:
    770892c   ----------
    771       REAL masse(ip1jmp1,llm),pente_max
    772       REAL q(ip1jmp1,llm)
     893      REAL masse(ip1jmp1,llm,nqtot),pente_max
     894      REAL q(ip1jmp1,llm,nqtot)
    773895      REAL w(ip1jmp1,llm+1)
     896      INTEGER iq
    774897c
    775898c      Local
    … …  
    782905      REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
    783906      REAL sigw
     907
     908      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     909      INTEGER ifils,iq2 ! CRisi
    784910
    785911      LOGICAL testcpu
    … …  
    795921c    On oriente tout dans le sens de la pression c'est a dire dans le
    796922c    sens de W
     923
     924      !write(*,*) 'vlz 923: entree'
    797925
    798926#ifdef BIDON
    … …  
    803931      DO l=2,llm
    804932         DO ij=1,ip1jmp1
    805             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
     933            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
    806934            adzqw(ij,l)=abs(dzqw(ij,l))
    807935         ENDDO
    … …  
    825953      ENDDO
    826954
     955      !write(*,*) 'vlz 954'
    827956      DO ij=1,ip1jmp1
    828957         dzq(ij,1)=0.
    … …  
    841970c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
    842971
     972       !write(*,*) 'vlz 969'
    843973       DO l = 1,llm-1
    844974         do  ij = 1,ip1jmp1
    845975          IF(w(ij,l+1).gt.0.) THEN
    846              sigw=w(ij,l+1)/masse(ij,l+1)
    847              wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
     976             sigw=w(ij,l+1)/masse(ij,l+1,iq)
     977             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq)
     978     &           +0.5*(1.-sigw)*dzq(ij,l+1))
    848979          ELSE
    849              sigw=w(ij,l+1)/masse(ij,l)
    850              wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
     980             sigw=w(ij,l+1)/masse(ij,l,iq)
     981             wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l))
    851982          ENDIF
    852983         ENDDO
    … …  
    858989       ENDDO
    859990
     991! CRisi: appel récursif de l'advection sur les fils.
     992! Il faut faire ça avant d'avoir mis à jour q et masse
     993      !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
     994      if (nqfils(iq).gt.0) then 
     995       do ifils=1,nqdesc(iq)
     996         iq2=iqfils(ifils,iq)
     997         DO l=1,llm
     998          DO ij=1,ip1jmp1
     999           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     1000           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)       
     1001          enddo   
     1002         enddo
     1003        enddo !do ifils=1,nqdesc(iq)
     1004       
     1005        do ifils=1,nqfils(iq)
     1006         iq2=iqfils(ifils,iq)         
     1007         call vlz(Ratio,pente_max,masseq,wq,iq2)
     1008        enddo !do ifils=1,nqfils(iq)
     1009      endif !if (nqfils(iq).gt.0) then
     1010! end CRisi 
     1011
    8601012      DO l=1,llm
    8611013         DO ij=1,ip1jmp1
    862             newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
    863             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
     1014            newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l)
     1015            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l))
    8641016     &         /newmasse
    865             masse(ij,l)=newmasse
    866          ENDDO
    867       ENDDO
    868 
     1017            masse(ij,l,iq)=newmasse
     1018         ENDDO
     1019      ENDDO
     1020
     1021! retablir les fils en rapport de melange par rapport a l'air:
     1022      if (nqfils(iq).gt.0) then 
     1023       do ifils=1,nqdesc(iq)
     1024         iq2=iqfils(ifils,iq) 
     1025         DO l=1,llm
     1026          DO ij=1,ip1jmp1
     1027            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     1028          enddo
     1029         enddo
     1030        enddo !do ifils=1,nqdesc(iq)
     1031      endif !if (nqfils(iq).gt.0) then
     1032      !write(*,*) 'vlsplt 1032'
    8691033
    8701034      RETURN

  • trunk/LMDZ.COMMON/libf/dyn3d/vlspltqs.F

    r1422 r1508  
    1 !
    2 ! $Header$
    3 !
     1c
     2c $Id: vlspltqs.F 2286 2015-05-20 13:27:07Z emillour $
     3c
    44       SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
    5      ,                                  p,pk,teta                 )
     5     ,                                  p,pk,teta,iq             )
     6       USE infotrac, ONLY: nqtot,nqdesc,iqfils
    67c
    78c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
    … …  
    3334      REAL masse(ip1jmp1,llm),pente_max
    3435      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    35       REAL q(ip1jmp1,llm)
     36      REAL q(ip1jmp1,llm,nqtot)
    3637      REAL w(ip1jmp1,llm),pdt
    3738      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
     39      INTEGER iq ! CRisi
    3840c
    3941c      Local
    … …  
    4143c
    4244      INTEGER i,ij,l,j,ii
     45      INTEGER ifils,iq2 ! CRisi
    4346c
    4447      REAL qsat(ip1jmp1,llm)
    45       REAL zm(ip1jmp1,llm)
     48      REAL zm(ip1jmp1,llm,nqtot)
    4649      REAL mu(ip1jmp1,llm)
    4750      REAL mv(ip1jm,llm)
    4851      REAL mw(ip1jmp1,llm+1)
    49       REAL zq(ip1jmp1,llm)
     52      REAL zq(ip1jmp1,llm,nqtot)
    5053      REAL temps1,temps2,temps3
    5154      REAL zzpbar, zzw
    … …  
    6366      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
    6467      REAL ptarg,pdelarg,foeew,zdelta
    65 ! ADAPTATION GCM POUR CP(T)
    6668      REAL tempe(ip1jmp1,llm)
    6769
    … …  
    115117      ENDDO
    116118
    117       CALL SCOPY(ijp1llm,q,1,zq,1)
    118       CALL SCOPY(ijp1llm,masse,1,zm,1)
     119      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
     120      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
     121      if (nqdesc(iq).gt.0) then 
     122       do ifils=1,nqdesc(iq)
     123        iq2=iqfils(ifils,iq)
     124        CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
     125       enddo 
     126      endif !if (nqfils(iq).gt.0) then
    119127
    120128c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
    121       call vlxqs(zq,pente_max,zm,mu,qsat)
    122 
     129      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
    123130
    124131c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
    125132
    126       call vlyqs(zq,pente_max,zm,mv,qsat)
    127 
     133      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
    128134
    129135c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
    130136
    131       call vlz(zq,pente_max,zm,mw)
    132 
     137      call vlz(zq,pente_max,zm,mw,iq)
    133138
    134139c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
    135140c     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')
    136141
    137       call vlyqs(zq,pente_max,zm,mv,qsat)
    138 
     142      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
    139143
    140144c     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
    141145c     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')
    142146
    143       call vlxqs(zq,pente_max,zm,mu,qsat)
     147      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
    144148
    145149c     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')
    … …  
    149153      DO l=1,llm
    150154         DO ij=1,ip1jmp1
    151            q(ij,l)=zq(ij,l)
     155           q(ij,l,iq)=zq(ij,l,iq)
    152156         ENDDO
    153157         DO ij=1,ip1jm+1,iip1
    154             q(ij+iim,l)=q(ij,l)
    155          ENDDO
    156       ENDDO
     158            q(ij+iim,l,iq)=q(ij,l,iq)
     159         ENDDO
     160      ENDDO
     161      ! CRisi: aussi pour les fils
     162      if (nqdesc(iq).gt.0) then
     163      do ifils=1,nqdesc(iq)
     164        iq2=iqfils(ifils,iq)
     165        DO l=1,llm
     166         DO ij=1,ip1jmp1
     167           q(ij,l,iq2)=zq(ij,l,iq2)
     168         ENDDO
     169         DO ij=1,ip1jm+1,iip1
     170            q(ij+iim,l,iq2)=q(ij,l,iq2)
     171         ENDDO
     172        ENDDO
     173      enddo !do ifils=1,nqdesc(iq) 
     174      endif ! if (nqfils(iq).gt.0) then
     175      !write(*,*) 'vlspltqs 183: fin de la routine'
    157176
    158177      RETURN
    159178      END
    160       SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
     179      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
     180      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
     181
    161182c
    162183c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    175196c   Arguments:
    176197c   ----------
    177       REAL masse(ip1jmp1,llm),pente_max
     198      REAL masse(ip1jmp1,llm,nqtot),pente_max
    178199      REAL u_m( ip1jmp1,llm )
    179       REAL q(ip1jmp1,llm)
     200      REAL q(ip1jmp1,llm,nqtot)
    180201      REAL qsat(ip1jmp1,llm)
     202      INTEGER iq ! CRisi
    181203c
    182204c      Local
    … …  
    191213      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    192214      REAL u_mq(ip1jmp1,llm)
     215
     216      ! CRisi
     217      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
     218      INTEGER ifils,iq2 ! CRisi
    193219
    194220      Logical first,testcpu
    … …  
    223249         DO l = 1, llm
    224250            DO ij=iip2,ip1jm-1
    225                dxqu(ij)=q(ij+1,l)-q(ij,l)
     251               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    226252c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    227 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     253c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    228254            ENDDO
    229255            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    277303         DO l = 1, llm
    278304            DO ij=iip2,ip1jm-1
    279                dxqu(ij)=q(ij+1,l)-q(ij,l)
     305               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    280306            ENDDO
    281307            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    319345      DO l=1,llm
    320346       DO ij=iip2,ip1jm-1
    321           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    322      ,                     1.+u_m(ij,l)/masse(ij+1,l),
     347          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     348     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
    323349     ,                     u_m(ij,l))
    324350          zdum(ij,l)=0.5*zdum(ij,l)
    325351          u_mq(ij,l)=cvmgp(
    326      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    327      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     352     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     353     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    328354     ,                u_m(ij,l))
    329355          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    337363       DO ij=iip2,ip1jm-1
    338364          IF (u_m(ij,l).gt.0.) THEN
    339              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
     365             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
    340366             u_mq(ij,l)=u_m(ij,l)*
    341      $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
     367     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
    342368          ELSE
    343              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
     369             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
    344370             u_mq(ij,l)=u_m(ij,l)*
    345      $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
     371     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
    346372          ENDIF
    347373       ENDDO
    … …  
    412438                     i=ijq-(j-1)*iip1
    413439c   accumulation pour les mailles completements advectees
    414                      do while(zu_m.gt.masse(ijq,l))
    415                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    416                         zu_m=zu_m-masse(ijq,l)
     440                     do while(zu_m.gt.masse(ijq,l,iq))
     441                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
     442     &                          *masse(ijq,l,iq)
     443                        zu_m=zu_m-masse(ijq,l,iq)
    417444                        i=mod(i-2+iim,iim)+1
    418445                        ijq=(j-1)*iip1+i
    … …  
    420447c   ajout de la maille non completement advectee
    421448                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    422      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     449     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
     450     &                  *dxq(ijq,l))
    423451                  ELSE
    424452                     ijq=ij+1
    425453                     i=ijq-(j-1)*iip1
    426454c   accumulation pour les mailles completements advectees
    427                      do while(-zu_m.gt.masse(ijq,l))
    428                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    429                         zu_m=zu_m+masse(ijq,l)
     455                     do while(-zu_m.gt.masse(ijq,l,iq))
     456                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     457     &                          *masse(ijq,l,iq)
     458                        zu_m=zu_m+masse(ijq,l,iq)
    430459                        i=mod(i,iim)+1
    431460                        ijq=(j-1)*iip1+i
    432461                     ENDDO
    433462c   ajout de la maille non completement advectee
    434                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    435      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     463                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     464     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    436465                  ENDIF
    437466               ENDDO
    … …  
    450479      ENDDO
    451480
     481! CRisi: appel récursif de l'advection sur les fils.
     482! Il faut faire ça avant d'avoir mis à jour q et masse
     483      !write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq)
     484     
     485      if (nqfils(iq).gt.0) then 
     486       do ifils=1,nqdesc(iq)
     487         iq2=iqfils(ifils,iq)
     488         DO l=1,llm
     489          DO ij=iip2,ip1jm
     490           ! On a besoin de q et masse seulement entre iip2 et ip1jm
     491           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     492           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     493          enddo   
     494         enddo
     495        enddo !do ifils=1,nqdesc(iq)
     496        do ifils=1,nqfils(iq)
     497         iq2=iqfils(ifils,iq)
     498         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
     499        enddo !do ifils=1,nqfils(iq)
     500      endif !if (nqfils(iq).gt.0) then
     501! end CRisi
    452502
    453503c   calcul des tendances
    … …  
    455505      DO l=1,llm
    456506         DO ij=iip2+1,ip1jm
    457             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    458             q(ij,l)=(q(ij,l)*masse(ij,l)+
     507            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     508            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    459509     &      u_mq(ij-1,l)-u_mq(ij,l))
    460510     &      /new_m
    461             masse(ij,l)=new_m
     511            masse(ij,l,iq)=new_m
    462512         ENDDO
    463513c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    464514         DO ij=iip1+iip1,ip1jm,iip1
    465             q(ij-iim,l)=q(ij,l)
    466             masse(ij-iim,l)=masse(ij,l)
    467          ENDDO
    468       ENDDO
     515            q(ij-iim,l,iq)=q(ij,l,iq)
     516            masse(ij-iim,l,iq)=masse(ij,l,iq)
     517         ENDDO
     518      ENDDO
     519
     520      ! retablir les fils en rapport de melange par rapport a l'air:
     521      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
     522      ! puis on boucle en longitude
     523      if (nqdesc(iq).gt.0) then 
     524       do ifils=1,nqdesc(iq)
     525         iq2=iqfils(ifils,iq) 
     526         DO l=1,llm
     527          DO ij=iip2+1,ip1jm
     528            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     529          enddo
     530          DO ij=iip1+iip1,ip1jm,iip1
     531             q(ij-iim,l,iq2)=q(ij,l,iq2)
     532          enddo ! DO ij=ijb+iip1-1,ije,iip1
     533         enddo !DO l=1,llm
     534        enddo !do ifils=1,nqdesc(iq)
     535      endif !if (nqfils(iq).gt.0) then
    469536
    470537c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    … …  
    474541      RETURN
    475542      END
    476       SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
     543      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
     544      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
     545      USE comconst_mod, ONLY: pi
    477546c
    478547c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    486555c
    487556c   --------------------------------------------------------------------
    488       USE comconst_mod, ONLY: pi
    489 
    490557      IMPLICIT NONE
    491558c
    … …  
    497564c   Arguments:
    498565c   ----------
    499       REAL masse(ip1jmp1,llm),pente_max
     566      REAL masse(ip1jmp1,llm,nqtot),pente_max
    500567      REAL masse_adv_v( ip1jm,llm)
    501       REAL q(ip1jmp1,llm)
     568      REAL q(ip1jmp1,llm,nqtot)
    502569      REAL qsat(ip1jmp1,llm)
     570      INTEGER iq ! CRisi
    503571c
    504572c      Local
    … …  
    524592      SAVE sinlon,coslon,sinlondlon,coslondlon
    525593      SAVE airej2,airejjm
     594
     595      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     596      INTEGER ifils,iq2 ! CRisi
    526597c
    527598c
    … …  
    562633
    563634      DO i = 1, iim
    564       airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
    565       airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     635      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
     636      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    566637      ENDDO
    567638      qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    571642
    572643      DO ij=1,ip1jm
    573          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     644         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    574645         adyqv(ij)=abs(dyqv(ij))
    575646      ENDDO
    … …  
    586657
    587658      DO ij=1,iip1
    588          dyq(ij,l)=qpns-q(ij+iip1,l)
    589          dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     659         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
     660         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    590661      ENDDO
    591662
    … …  
    705776       DO ij=1,ip1jm
    706777         IF( masse_adv_v(ij,l).GT.0. ) THEN
    707            qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
    708      ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
     778           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
     779     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
     780     ,      /masse(ij+iip1,l,iq)))
    709781         ELSE
    710               qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
    711      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
     782              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
     783     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
    712784         ENDIF
    713785          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    716788
    717789
     790! CRisi: appel récursif de l'advection sur les fils.
     791! Il faut faire ça avant d'avoir mis à jour q et masse
     792      !write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
     793   
     794      if (nqfils(iq).gt.0) then 
     795       do ifils=1,nqdesc(iq)
     796         iq2=iqfils(ifils,iq)
     797         DO l=1,llm
     798         DO ij=1,ip1jmp1
     799           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     800           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     801          enddo   
     802         enddo
     803        enddo !do ifils=1,nqdesc(iq)
     804
     805        do ifils=1,nqfils(iq)
     806         iq2=iqfils(ifils,iq)
     807         !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
     808         call vly(Ratio,pente_max,masseq,qbyv,iq2)
     809        enddo !do ifils=1,nqfils(iq)
     810      endif !if (nqfils(iq).gt.0) then
     811
    718812      DO l=1,llm
    719813         DO ij=iip2,ip1jm
    720             newmasse=masse(ij,l)
     814            newmasse=masse(ij,l,iq)
    721815     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    722             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    723      &         /newmasse
    724             masse(ij,l)=newmasse
     816            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
     817     &         -qbyv(ij-iip1,l))/newmasse
     818            masse(ij,l,iq)=newmasse
    725819         ENDDO
    726820c.-. ancienne version
    … …  
    728822         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
    729823         DO ij = 1,iip1
    730             newmasse=masse(ij,l)+convmpn*aire(ij)
    731             q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
     824            newmasse=masse(ij,l,iq)+convmpn*aire(ij)
     825            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
    732826     &               newmasse
    733             masse(ij,l)=newmasse
     827            masse(ij,l,iq)=newmasse
    734828         ENDDO
    735829         convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
    736830         convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
    737831         DO ij = ip1jm+1,ip1jmp1
    738             newmasse=masse(ij,l)+convmps*aire(ij)
    739             q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
     832            newmasse=masse(ij,l,iq)+convmps*aire(ij)
     833            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
    740834     &               newmasse
    741             masse(ij,l)=newmasse
     835            masse(ij,l,iq)=newmasse
    742836         ENDDO
    743837c.-. fin ancienne version
    … …  
    752846c        DO ij = 1,iip1
    753847c           q(ij,l)=newq
    754 c           masse(ij,l)=newmasse*aire(ij)
     848c           masse(ij,l,iq)=newmasse*aire(ij)
    755849c        ENDDO
    756850c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    … …  
    762856c        DO ij = ip1jm+1,ip1jmp1
    763857c           q(ij,l)=newq
    764 c           masse(ij,l)=newmasse*aire(ij)
     858c           masse(ij,l,iq)=newmasse*aire(ij)
    765859c        ENDDO
    766860c._. fin nouvelle version
    767861      ENDDO
    768862
     863      !write(*,*) 'vly 866'
     864
     865! retablir les fils en rapport de melange par rapport a l'air:
     866      if (nqdesc(iq).gt.0) then 
     867       do ifils=1,nqdesc(iq)
     868         iq2=iqfils(ifils,iq) 
     869         DO l=1,llm
     870          DO ij=1,ip1jmp1
     871            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     872          enddo
     873         enddo
     874        enddo !do ifils=1,nqdesc(iq)
     875      endif !if (nqfils(iq).gt.0) then
     876      !write(*,*) 'vly 879'
     877
    769878      RETURN
    770879      END
Note: See TracChangeset for help on using the changeset viewer.