Ignore:
Timestamp:
May 7, 2015, 5:45:04 PM (9 years ago)
Author:
crisi
Message:

Adding isotopes in the dynamics and more generally tracers of tracers.
CRisi

Location:
LMDZ5/trunk/libf/dyn3d
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d/advtrac.F90

    r2239 r2270  
    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
    … …  
    223223     !     Appel des sous programmes d'advection
    224224     !-----------------------------------------------------------
    225      do iq=1,nqtot
     225
     226     if (ok_iso_verif) then
     227           call check_isotopes_seq(q,1,ip1jmp1,'advtrac 162')
     228     endif !if (ok_iso_verif) then
     229
     230     do iq=1,nqperes
    226231        !        call clock(t_initial)
    227232        if(iadv(iq) == 0) cycle
    … …  
    230235        !   ----------------------------------------------------------------
    231236        if(iadv(iq).eq.10) THEN
    232            call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
    233 
     237           ! CRisi: on fait passer tout q pour avoir acces aux fils
     238           
     239           write(*,*) 'advtrac 239: iq,q(1721,19,:)=',iq,q(1721,19,:)     
     240           call vlsplt(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq)
     241           
    234242           !   ----------------------------------------------------------------
    235243           !   Schema "pseudo amont" + test sur humidite specifique
    … …  
    238246        else if(iadv(iq).eq.14) then
    239247           !
    240            CALL vlspltqs( q(1,1,1), 2., massem, wg , &
    241                 pbarug,pbarvg,dtvr,p,pk,teta )
     248           write(*,*) 'advtrac 248: iq,q(1721,19,:)=',iq,q(1721,19,:)
     249           CALL vlspltqs( q, 2., massem, wg , &
     250                pbarug,pbarvg,dtvr,p,pk,teta,iq)
     251           
    242252           !   ----------------------------------------------------------------
    243253           !   Schema de Frederic Hourdin
    … …  
    388398     end DO
    389399
     400     if (ok_iso_verif) then
     401           call check_isotopes_seq(q,ip1jmp1,'advtrac 397')
     402     endif !if (ok_iso_verif) then
    390403
    391404     !------------------------------------------------------------------

  • LMDZ5/trunk/libf/dyn3d/caladvtrac.F

    r1907 r2270  
    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
    … …  
    8184           ENDDO
    8285          ENDDO
    83          
    84           CALL qminimum( q, 2, finmasse )
     86
     87          CALL qminimum( q, nqtot, finmasse )
    8588
    8689          CALL SCOPY   ( ip1jmp1*llm, masse, 1, finmasse,       1 )
    … …  
    9295          dtvrtrac = iapp_tracvl * dtvr
    9396c
    94            DO iq = 1 , 2
     97           DO iq = 1 , nqtot
    9598            DO l = 1 , llm
    9699             DO ij = 1,ip1jmp1
    … …  
    105108        if (planet_type.eq."earth") then
    106109! Earth-specific treatment for the first 2 tracers (water)
    107           dq(:,:,1:2)=0.
     110          dq(:,:,1:nqtot)=0.
    108111        endif ! of if (planet_type.eq."earth")
    109112      ENDIF ! of IF( iapptrac.EQ.iapp_tracvl )

  • LMDZ5/trunk/libf/dyn3d/dynetat0.F

    r1930 r2270  
    297297           write(lunout,*)"          Il est donc initialise a zero"
    298298           q(:,:,:,iq)=0.
     299
     300           ! CRisi: pour les isotopes, on peut faire init théorique
     301           ! distill de Rayleigh très simplifiée
     302           if (ok_isotopes) then
     303              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.0)) then
     304                q(:,:,:,iq)=q(:,:,:,iqpere(iq))                         
     305     &                   *tnat(iso_num(iq))                             
     306     &                   *(q(:,:,:,iqpere(iq))/30.e-3)                 
     307     &                   **(alpha_ideal(iso_num(iq))-1)
     308              endif
     309              if ((iso_num(iq).gt.0).and.(zone_num(iq).eq.1)) then
     310                  q(:,:,:,iq)=q(:,:,:,iqiso(iso_indnum(iq),
     311     &                   phase_num(iq)))
     312              endif 
     313           endif !if (ok_isotopes) then
    299314        ELSE
    300315           ierr = NF90_GET_VAR(nid, nvarid, q(:,:,:,iq))

  • LMDZ5/trunk/libf/dyn3d/iniacademic.F90

    r2087 r2270  
    55
    66  USE filtreg_mod, ONLY: inifilr
    7   USE infotrac, ONLY : nqtot
     7  USE infotrac
    88  USE control_mod, ONLY: day_step,planet_type
    99#ifdef CPP_IOIPSL
    … …  
    262262              if (i == 2) q(:,:,i)=1.e-15
    263263              if (i.gt.2) q(:,:,i)=0.
     264
     265              ! CRisi: init des isotopes
     266              ! distill de Rayleigh très simplifiée
     267              if (ok_isotopes) then
     268                if ((iso_num(i).gt.0).and.(zone_num(i).eq.0)) then         
     269                   q(:,:,i)=q(:,:,iqpere(i))             &
     270      &                  *tnat(iso_num(i))               &
     271      &                  *(q(:,:,iqpere(i))/30.e-3)      &
     272      &                  **(alpha_ideal(iso_num(i))-1)
     273                endif               
     274                if ((iso_num(i).gt.0).and.(zone_num(i).eq.1)) then
     275                  q(:,:,i)=q(:,:,iqiso(iso_indnum(i),phase_num(i)))
     276                endif
     277              endif !if (ok_isotopes) then
     278
    264279           enddo
    265280        else
    266281           q(:,:,:)=0
    267282        endif ! of if (planet_type=="earth")
     283
     284        if (ok_iso_verif) then
     285           call check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
     286        endif !if (ok_iso_verif) then
    268287
    269288        ! add random perturbation to temperature

  • LMDZ5/trunk/libf/dyn3d/leapfrog.F

    r2239 r2270  
    1111      use IOIPSL
    1212#endif
    13       USE infotrac, ONLY: nqtot
     13      USE infotrac, ONLY: nqtot,ok_iso_verif
    1414      USE guide_mod, ONLY : guide_main
    1515      USE write_field, ONLY: writefield
    … …  
    235235      jH_cur = jH_cur - int(jH_cur)
    236236
     237        if (ok_iso_verif) then
     238           call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
     239        endif !if (ok_iso_verif) then
    237240
    238241#ifdef CPP_IOIPSL
    … …  
    265268!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
    266269!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
     270
     271        if (ok_iso_verif) then
     272           call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
     273        endif !if (ok_iso_verif) then
    267274
    268275   2  CONTINUE ! Matsuno backward or leapfrog step begins here
    … …  
    305312      endif
    306313
     314
     315        if (ok_iso_verif) then
     316           call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
     317        endif !if (ok_iso_verif) then
     318
    307319c-----------------------------------------------------------------------
    308320c   calcul des tendances dynamiques:
    … …  
    321333c   calcul des tendances advection des traceurs (dont l'humidite)
    322334c   -------------------------------------------------------------
     335
     336        if (ok_iso_verif) then
     337           call check_isotopes_seq(q,ip1jmp1,
     338     &           'leapfrog 686: avant caladvtrac')
     339        endif !if (ok_iso_verif) then
    323340
    324341      IF( forward. OR . leapf )  THEN
    … …  
    327344     *        p, masse, dq,  teta,
    328345     .        flxw, pk)
     346
    329347         
    330348         IF (offline) THEN
    … …  
    346364c   ----------------------------------
    347365
    348 
    349        CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
     366        if (ok_iso_verif) then
     367           write(*,*) 'leapfrog 720'
     368           call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
     369        endif !if (ok_iso_verif) then
     370       
     371       CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
    350372     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
    351373!     $              finvmaold                                    )
    352374
     375       if (ok_iso_verif) then
     376          write(*,*) 'leapfrog 724'
     377           call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
     378        endif !if (ok_iso_verif) then
    353379
    354380c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
    … …  
    516542        CALL massdair(p,masse)
    517543
     544        if (ok_iso_verif) then
     545           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
     546        endif !if (ok_iso_verif) then
    518547
    519548c-----------------------------------------------------------------------
    … …  
    600629c   preparation du pas d'integration suivant  ......
    601630
     631        if (ok_iso_verif) then
     632           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
     633        endif !if (ok_iso_verif) then
     634
    602635      IF ( .NOT.purmats ) THEN
    603636c       ........................................................
    … …  
    657690            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
    658691
     692        if (ok_iso_verif) then
     693           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
     694        endif !if (ok_iso_verif) then
     695
    659696c-----------------------------------------------------------------------
    660697c   ecriture de la bande histoire:
    … …  
    735772      ELSE ! of IF (.not.purmats)
    736773
     774        if (ok_iso_verif) then
     775           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
     776        endif !if (ok_iso_verif) then
     777
    737778c       ........................................................
    738779c       ..............       schema  matsuno        ...............
    … …  
    757798
    758799            ELSE ! of IF(forward) i.e. backward step
     800 
     801        if (ok_iso_verif) then
     802           call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
     803        endif !if (ok_iso_verif) then 
    759804
    760805              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN

  • LMDZ5/trunk/libf/dyn3d/qminimum.F

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

  • LMDZ5/trunk/libf/dyn3d/vlsplt.F

    r1907 r2270  
    33c
    44
    5       SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt)
     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
    … …  
    3233c      REAL masse(iip1,jjp1,llm),pente_max
    3334      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    34       REAL q(ip1jmp1,llm)
     35      REAL q(ip1jmp1,llm,nqtot)
    3536c      REAL q(iip1,jjp1,llm)
    3637      REAL w(ip1jmp1,llm),pdt
     38      INTEGER iq ! CRisi
    3739c
    3840c      Local
    … …  
    4244      INTEGER ijlqmin,iqmin,jqmin,lqmin
    4345c
    44       REAL zm(ip1jmp1,llm),newmasse
     46      REAL zm(ip1jmp1,llm,nqtot),newmasse
    4547      REAL mu(ip1jmp1,llm)
    4648      REAL mv(ip1jm,llm)
    4749      REAL mw(ip1jmp1,llm+1)
    48       REAL zq(ip1jmp1,llm),zz
     50      REAL zq(ip1jmp1,llm,nqtot),zz
    4951      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
    5052      REAL second,temps0,temps1,temps2,temps3
    … …  
    5557      SAVE temps1,temps2,temps3
    5658      INTEGER iminn,imaxx
     59      INTEGER ifils,iq2 ! CRisi
    5760
    5861      REAL qmin,qmax
    … …  
    7982         mw(ij,llm+1)=0.
    8083      ENDDO
    81      
    82       CALL SCOPY(ijp1llm,q,1,zq,1)
    83       CALL SCOPY(ijp1llm,masse,1,zm,1)
     84           
     85      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
     86      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
     87       
     88      if (nqdesc(iq).gt.0) then 
     89        do ifils=1,nqdesc(iq)
     90          iq2=iqfils(ifils,iq)
     91          CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
     92        enddo 
     93      endif !if (nqfils(iq).gt.0) then
    8494
    8595cprint*,'Entree vlx1'
    8696c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
    87       call vlx(zq,pente_max,zm,mu)
     97      call vlx(zq,pente_max,zm,mu,iq)
    8898cprint*,'Sortie vlx1'
    8999c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
    90100
    91101c print*,'Entree vly1'
    92       call vly(zq,pente_max,zm,mv)
     102
     103      call vly(zq,pente_max,zm,mv,iq)
    93104c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
    94105cprint*,'Sortie vly1'
    95       call vlz(zq,pente_max,zm,mw)
     106      call vlz(zq,pente_max,zm,mw,iq)
    96107c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
    97108
    98109
    99       call vly(zq,pente_max,zm,mv)
     110      call vly(zq,pente_max,zm,mv,iq)
    100111c       call minmaxq(zq,qmin,qmax,'apres vly     ')
    101112
    102113
    103       call vlx(zq,pente_max,zm,mu)
     114      call vlx(zq,pente_max,zm,mu,iq)
    104115c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
    105116       
    … …  
    107118      DO l=1,llm
    108119         DO ij=1,ip1jmp1
    109            q(ij,l)=zq(ij,l)
     120           q(ij,l,iq)=zq(ij,l,iq)
    110121         ENDDO
    111122         DO ij=1,ip1jm+1,iip1
    112             q(ij+iim,l)=q(ij,l)
    113          ENDDO
    114       ENDDO
     123            q(ij+iim,l,iq)=q(ij,l,iq)
     124         ENDDO
     125      ENDDO
     126      ! CRisi: aussi pour les fils
     127      if (nqdesc(iq).gt.0) then
     128      do ifils=1,nqdesc(iq)
     129        iq2=iqfils(ifils,iq)
     130        DO l=1,llm
     131         DO ij=1,ip1jmp1
     132           q(ij,l,iq2)=zq(ij,l,iq2)
     133         ENDDO
     134         DO ij=1,ip1jm+1,iip1
     135            q(ij+iim,l,iq2)=q(ij,l,iq2)
     136         ENDDO
     137        ENDDO
     138      enddo !do ifils=1,nqdesc(iq)   
     139      endif ! if (nqdesc(iq).gt.0) then   
    115140
    116141      RETURN
    117142      END
    118       SUBROUTINE vlx(q,pente_max,masse,u_m)
     143      SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
     144      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    119145
    120146c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    139165c   Arguments:
    140166c   ----------
    141       REAL masse(ip1jmp1,llm),pente_max
     167      REAL masse(ip1jmp1,llm,nqtot),pente_max
    142168      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
    143       REAL q(ip1jmp1,llm)
     169      REAL q(ip1jmp1,llm,nqtot)
    144170      REAL w(ip1jmp1,llm)
     171      INTEGER iq ! CRisi
    145172c
    146173c      Local
    … …  
    155182      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    156183      REAL u_mq(ip1jmp1,llm)
     184
     185      ! CRisi
     186      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
     187      INTEGER ifils,iq2 ! CRisi
    157188
    158189      Logical extremum,first,testcpu
    … …  
    188219         DO l = 1, llm
    189220            DO ij=iip2,ip1jm-1
    190                dxqu(ij)=q(ij+1,l)-q(ij,l)
     221               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    191222c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    192 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     223c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    193224            ENDDO
    194225            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    243274         DO l = 1, llm
    244275            DO ij=iip2,ip1jm-1
    245                dxqu(ij)=q(ij+1,l)-q(ij,l)
     276               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    246277            ENDDO
    247278            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    285316      DO l=1,llm
    286317       DO ij=iip2,ip1jm-1
    287           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    288      ,                     1.+u_m(ij,l)/masse(ij+1,l),
     318          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     319     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
    289320     ,                     u_m(ij,l))
    290321          zdum(ij,l)=0.5*zdum(ij,l)
    291322          u_mq(ij,l)=cvmgp(
    292      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    293      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     323     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     324     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    294325     ,                u_m(ij,l))
    295326          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    303334      DO l=1,llm
    304335       DO ij=iip2,ip1jm-1
    305 c       print*,'masse(',ij,')=',masse(ij,l)
     336c       print*,'masse(',ij,')=',masse(ij,l,iq)
    306337          IF (u_m(ij,l).gt.0.) THEN
    307              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
    308              u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
     338             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
     339             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l))
    309340          ELSE
    310              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
    311              u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
     341             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
     342             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
     343     &           -0.5*zdum(ij,l)*dxq(ij+1,l))
    312344          ENDIF
    313345       ENDDO
    … …  
    379411                     i=ijq-(j-1)*iip1
    380412c   accumulation pour les mailles completements advectees
    381                      do while(zu_m.gt.masse(ijq,l))
    382                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    383                         zu_m=zu_m-masse(ijq,l)
     413                     do while(zu_m.gt.masse(ijq,l,iq))
     414                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
     415     &                          *masse(ijq,l,iq)
     416                        zu_m=zu_m-masse(ijq,l,iq)
    384417                        i=mod(i-2+iim,iim)+1
    385418                        ijq=(j-1)*iip1+i
    … …  
    387420c   ajout de la maille non completement advectee
    388421                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    389      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     422     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
     423     &                  *dxq(ijq,l))
    390424                  ELSE
    391425                     ijq=ij+1
    392426                     i=ijq-(j-1)*iip1
    393427c   accumulation pour les mailles completements advectees
    394                      do while(-zu_m.gt.masse(ijq,l))
    395                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    396                         zu_m=zu_m+masse(ijq,l)
     428                     do while(-zu_m.gt.masse(ijq,l,iq))
     429                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     430     &                          *masse(ijq,l,iq)
     431                        zu_m=zu_m+masse(ijq,l,iq)
    397432                        i=mod(i,iim)+1
    398433                        ijq=(j-1)*iip1+i
    399434                     ENDDO
    400435c   ajout de la maille non completement advectee
    401                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    402      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     436                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     437     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    403438                  ENDIF
    404439               ENDDO
    … …  
    417452      ENDDO
    418453
     454! CRisi: appel récursif de l'advection sur les fils.
     455! Il faut faire ça avant d'avoir mis à jour q et masse
     456      write(*,*) 'vlsplt 326: iq,nqfils(iq)=',iq,nqfils(iq)
     457     
     458      if (nqdesc(iq).gt.0) then 
     459       do ifils=1,nqdesc(iq)
     460         iq2=iqfils(ifils,iq)
     461         DO l=1,llm
     462          DO ij=iip2,ip1jm
     463           ! On a besoin de q et masse seulement entre iip2 et ip1jm
     464           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     465           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     466          enddo   
     467         enddo
     468        enddo !do ifils=1,nqdesc(iq)
     469        do ifils=1,nqfils(iq)
     470         iq2=iqfils(ifils,iq)
     471         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
     472        enddo !do ifils=1,nqfils(iq)
     473      endif !if (nqfils(iq).gt.0) then
     474! end CRisi
     475
    419476
    420477c   calcul des tENDances
    … …  
    422479      DO l=1,llm
    423480         DO ij=iip2+1,ip1jm
    424             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    425             q(ij,l)=(q(ij,l)*masse(ij,l)+
     481            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     482            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    426483     &      u_mq(ij-1,l)-u_mq(ij,l))
    427484     &      /new_m
    428             masse(ij,l)=new_m
     485            masse(ij,l,iq)=new_m
    429486         ENDDO
    430487c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    431488         DO ij=iip1+iip1,ip1jm,iip1
    432             q(ij-iim,l)=q(ij,l)
    433             masse(ij-iim,l)=masse(ij,l)
    434          ENDDO
    435       ENDDO
     489            q(ij-iim,l,iq)=q(ij,l,iq)
     490            masse(ij-iim,l,iq)=masse(ij,l,iq)
     491         ENDDO
     492      ENDDO
     493
     494      ! retablir les fils en rapport de melange par rapport a l'air:
     495      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
     496      ! puis on boucle en longitude
     497      if (nqdesc(iq).gt.0) then 
     498       do ifils=1,nqdesc(iq)
     499         iq2=iqfils(ifils,iq) 
     500         DO l=1,llm
     501          DO ij=iip2+1,ip1jm
     502            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     503          enddo
     504          DO ij=iip1+iip1,ip1jm,iip1
     505             q(ij-iim,l,iq2)=q(ij,l,iq2)
     506          enddo ! DO ij=ijb+iip1-1,ije,iip1
     507         enddo !DO l=1,llm
     508        enddo !do ifils=1,nqdesc(iq)
     509      endif !if (nqfils(iq).gt.0) then
     510
    436511c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    437512c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
    … …  
    440515      RETURN
    441516      END
    442       SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
     517      SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
     518      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    443519c
    444520c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    464540c   Arguments:
    465541c   ----------
    466       REAL masse(ip1jmp1,llm),pente_max
     542      REAL masse(ip1jmp1,llm,nqtot),pente_max
    467543      REAL masse_adv_v( ip1jm,llm)
    468       REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
     544      REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm)
     545      INTEGER iq ! CRisi
    469546c
    470547c      Local
    … …  
    491568      SAVE sinlon,coslon,sinlondlon,coslondlon
    492569      SAVE airej2,airejjm
     570
     571      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     572      INTEGER ifils,iq2 ! CRisi
     573
    493574c
    494575c
    … …  
    497578      DATA first,testcpu/.true.,.false./
    498579      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
     580
     581      write(*,*) 'vly 578: entree, iq=',iq
    499582
    500583      IF(first) THEN
    … …  
    529612
    530613      DO i = 1, iim
    531       airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
    532       airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     614      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
     615      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    533616      ENDDO
    534617      qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    538621
    539622      DO ij=1,ip1jm
    540          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     623         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    541624         adyqv(ij)=abs(dyqv(ij))
    542625      ENDDO
    … …  
    553636
    554637      DO ij=1,iip1
    555          dyq(ij,l)=qpns-q(ij+iip1,l)
    556          dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     638         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
     639         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    557640      ENDDO
    558641
    … …  
    675758      ENDDO
    676759
     760      write(*,*) 'vly 756'
    677761      DO l=1,llm
    678762       DO ij=1,ip1jm
    679763          IF(masse_adv_v(ij,l).gt.0) THEN
    680               qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
    681      ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
     764              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
     765     ,                   0.5*(1.-masse_adv_v(ij,l)
     766     ,                   /masse(ij+iip1,l,iq))
    682767          ELSE
    683               qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
    684      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
     768              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
     769     ,                   0.5*(1.+masse_adv_v(ij,l)
     770     ,                   /masse(ij,l,iq))
    685771          ENDIF
    686772          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    688774      ENDDO
    689775
     776! CRisi: appel récursif de l'advection sur les fils.
     777! Il faut faire ça avant d'avoir mis à jour q et masse
     778      write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq)
     779   
     780      if (nqfils(iq).gt.0) then 
     781       do ifils=1,nqdesc(iq)
     782         iq2=iqfils(ifils,iq)
     783         DO l=1,llm
     784         DO ij=1,ip1jmp1
     785           ! attention, chaque fils doit avoir son masseq, sinon, le 1er
     786           ! fils ecrase le masseq de ses freres.
     787           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     788           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     789          enddo   
     790         enddo
     791        enddo !do ifils=1,nqdesc(iq)
     792
     793        do ifils=1,nqfils(iq)
     794         iq2=iqfils(ifils,iq)
     795         call vly(Ratio,pente_max,masseq,qbyv,iq2)
     796        enddo !do ifils=1,nqfils(iq)
     797      endif !if (nqfils(iq).gt.0) then
    690798
    691799      DO l=1,llm
    692800         DO ij=iip2,ip1jm
    693             newmasse=masse(ij,l)
     801            newmasse=masse(ij,l,iq)
    694802     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    695             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    696      &         /newmasse
    697             masse(ij,l)=newmasse
     803            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
     804     &         -qbyv(ij-iip1,l))/newmasse
     805            masse(ij,l,iq)=newmasse
    698806         ENDDO
    699807c.-. ancienne version
    … …  
    703811         convpn=SSUM(iim,qbyv(1,l),1)
    704812         convmpn=ssum(iim,masse_adv_v(1,l),1)
    705          massepn=ssum(iim,masse(1,l),1)
     813         massepn=ssum(iim,masse(1,l,iq),1)
    706814         qpn=0.
    707815         do ij=1,iim
    708             qpn=qpn+masse(ij,l)*q(ij,l)
     816            qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
    709817         enddo
    710818         qpn=(qpn+convpn)/(massepn+convmpn)
    711819         do ij=1,iip1
    712             q(ij,l)=qpn
     820            q(ij,l,iq)=qpn
    713821         enddo
    714822
    … …  
    718826         convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    719827         convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    720          masseps=ssum(iim, masse(ip1jm+1,l),1)
     828         masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
    721829         qps=0.
    722830         do ij = ip1jm+1,ip1jmp1-1
    723             qps=qps+masse(ij,l)*q(ij,l)
     831            qps=qps+masse(ij,l,iq)*q(ij,l,iq)
    724832         enddo
    725833         qps=(qps+convps)/(masseps+convmps)
    726834         do ij=ip1jm+1,ip1jmp1
    727             q(ij,l)=qps
     835            q(ij,l,iq)=qps
    728836         enddo
    729837
    … …  
    739847c        DO ij = 1,iip1
    740848c           q(ij,l)=newq
    741 c           masse(ij,l)=newmasse*aire(ij)
     849c           masse(ij,l,iq)=newmasse*aire(ij)
    742850c        ENDDO
    743851c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    … …  
    749857c        DO ij = ip1jm+1,ip1jmp1
    750858c           q(ij,l)=newq
    751 c           masse(ij,l)=newmasse*aire(ij)
     859c           masse(ij,l,iq)=newmasse*aire(ij)
    752860c        ENDDO
    753861c._. fin nouvelle version
    754862      ENDDO
     863 
     864! retablir les fils en rapport de melange par rapport a l'air:
     865      if (nqfils(iq).gt.0) then 
     866       do ifils=1,nqdesc(iq)
     867         iq2=iqfils(ifils,iq) 
     868         DO l=1,llm
     869          DO ij=1,ip1jmp1
     870            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     871          enddo
     872         enddo
     873        enddo !do ifils=1,nqdesc(iq)
     874      endif !if (nqfils(iq).gt.0) then
     875
     876      write(*,*) 'vly 853: sortie'
    755877
    756878      RETURN
    757879      END
    758       SUBROUTINE vlz(q,pente_max,masse,w)
     880      SUBROUTINE vlz(q,pente_max,masse,w,iq)
     881      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    759882c
    760883c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    779902c   Arguments:
    780903c   ----------
    781       REAL masse(ip1jmp1,llm),pente_max
    782       REAL q(ip1jmp1,llm)
     904      REAL masse(ip1jmp1,llm,nqtot),pente_max
     905      REAL q(ip1jmp1,llm,nqtot)
    783906      REAL w(ip1jmp1,llm+1)
     907      INTEGER iq
    784908c
    785909c      Local
    … …  
    792916      REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
    793917      REAL sigw
     918
     919      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     920      INTEGER ifils,iq2 ! CRisi
    794921
    795922      LOGICAL testcpu
    … …  
    805932c    On oriente tout dans le sens de la pression c'est a dire dans le
    806933c    sens de W
     934
     935      write(*,*) 'vlz 923: entree'
    807936
    808937#ifdef BIDON
    … …  
    813942      DO l=2,llm
    814943         DO ij=1,ip1jmp1
    815             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
     944            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
    816945            adzqw(ij,l)=abs(dzqw(ij,l))
    817946         ENDDO
    … …  
    835964      ENDDO
    836965
     966      write(*,*) 'vlz 954'
    837967      DO ij=1,ip1jmp1
    838968         dzq(ij,1)=0.
    … …  
    851981c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
    852982
     983       write(*,*) 'vlz 969'
    853984       DO l = 1,llm-1
    854985         do  ij = 1,ip1jmp1
    855986          IF(w(ij,l+1).gt.0.) THEN
    856              sigw=w(ij,l+1)/masse(ij,l+1)
    857              wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
     987             sigw=w(ij,l+1)/masse(ij,l+1,iq)
     988             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq)
     989     &           +0.5*(1.-sigw)*dzq(ij,l+1))
    858990          ELSE
    859              sigw=w(ij,l+1)/masse(ij,l)
    860              wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
     991             sigw=w(ij,l+1)/masse(ij,l,iq)
     992             wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l))
    861993          ENDIF
    862994         ENDDO
    … …  
    8681000       ENDDO
    8691001
     1002! CRisi: appel récursif de l'advection sur les fils.
     1003! Il faut faire ça avant d'avoir mis à jour q et masse
     1004      write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
     1005      if (nqfils(iq).gt.0) then 
     1006       do ifils=1,nqdesc(iq)
     1007         iq2=iqfils(ifils,iq)
     1008         DO l=1,llm
     1009          DO ij=1,ip1jmp1
     1010           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     1011           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)       
     1012          enddo   
     1013         enddo
     1014        enddo !do ifils=1,nqdesc(iq)
     1015       
     1016        do ifils=1,nqfils(iq)
     1017         iq2=iqfils(ifils,iq)         
     1018         call vlz(Ratio,pente_max,masseq,wq,iq2)
     1019        enddo !do ifils=1,nqfils(iq)
     1020      endif !if (nqfils(iq).gt.0) then
     1021! end CRisi 
     1022
    8701023      DO l=1,llm
    8711024         DO ij=1,ip1jmp1
    872             newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
    873             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
     1025            newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l)
     1026            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l))
    8741027     &         /newmasse
    875             masse(ij,l)=newmasse
    876          ENDDO
    877       ENDDO
    878 
     1028            masse(ij,l,iq)=newmasse
     1029         ENDDO
     1030      ENDDO
     1031
     1032! retablir les fils en rapport de melange par rapport a l'air:
     1033      if (nqfils(iq).gt.0) then 
     1034       do ifils=1,nqdesc(iq)
     1035         iq2=iqfils(ifils,iq) 
     1036         DO l=1,llm
     1037          DO ij=1,ip1jmp1
     1038            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     1039          enddo
     1040         enddo
     1041        enddo !do ifils=1,nqdesc(iq)
     1042      endif !if (nqfils(iq).gt.0) then
     1043      write(*,*) 'vlsplt 1032'
    8791044
    8801045      RETURN

  • LMDZ5/trunk/libf/dyn3d/vlspltqs.F

    r1907 r2270  
    33c
    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
    … …  
    3536      REAL masse(ip1jmp1,llm),pente_max
    3637      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    37       REAL q(ip1jmp1,llm)
     38      REAL q(ip1jmp1,llm,nqtot)
    3839      REAL w(ip1jmp1,llm),pdt
    3940      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
     41      INTEGER iq ! CRisi
    4042c
    4143c      Local
    … …  
    4345c
    4446      INTEGER i,ij,l,j,ii
     47      INTEGER ifils,iq2 ! CRisi
    4548c
    4649      REAL qsat(ip1jmp1,llm)
    47       REAL zm(ip1jmp1,llm)
     50      REAL zm(ip1jmp1,llm,nqtot)
    4851      REAL mu(ip1jmp1,llm)
    4952      REAL mv(ip1jm,llm)
    5053      REAL mw(ip1jmp1,llm+1)
    51       REAL zq(ip1jmp1,llm)
     54      REAL zq(ip1jmp1,llm,nqtot)
    5255      REAL temps1,temps2,temps3
    5356      REAL zzpbar, zzw
    … …  
    116119      ENDDO
    117120
    118       CALL SCOPY(ijp1llm,q,1,zq,1)
    119       CALL SCOPY(ijp1llm,masse,1,zm,1)
     121      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
     122      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
     123      if (nqdesc(iq).gt.0) then 
     124       do ifils=1,nqdesc(iq)
     125        iq2=iqfils(ifils,iq)
     126        CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
     127       enddo 
     128      endif !if (nqfils(iq).gt.0) then
    120129
    121130c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
    122       call vlxqs(zq,pente_max,zm,mu,qsat)
    123 
     131      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
    124132
    125133c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
    126134
    127       call vlyqs(zq,pente_max,zm,mv,qsat)
    128 
     135      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
    129136
    130137c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
    131138
    132       call vlz(zq,pente_max,zm,mw)
    133 
     139      call vlz(zq,pente_max,zm,mw,iq)
    134140
    135141c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
    136142c     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')
    137143
    138       call vlyqs(zq,pente_max,zm,mv,qsat)
    139 
     144      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
    140145
    141146c     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
    142147c     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')
    143148
    144       call vlxqs(zq,pente_max,zm,mu,qsat)
     149      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
    145150
    146151c     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')
    … …  
    150155      DO l=1,llm
    151156         DO ij=1,ip1jmp1
    152            q(ij,l)=zq(ij,l)
     157           q(ij,l,iq)=zq(ij,l,iq)
    153158         ENDDO
    154159         DO ij=1,ip1jm+1,iip1
    155             q(ij+iim,l)=q(ij,l)
    156          ENDDO
    157       ENDDO
     160            q(ij+iim,l,iq)=q(ij,l,iq)
     161         ENDDO
     162      ENDDO
     163      ! CRisi: aussi pour les fils
     164      if (nqdesc(iq).gt.0) then
     165      do ifils=1,nqdesc(iq)
     166        iq2=iqfils(ifils,iq)
     167        DO l=1,llm
     168         DO ij=1,ip1jmp1
     169           q(ij,l,iq2)=zq(ij,l,iq2)
     170         ENDDO
     171         DO ij=1,ip1jm+1,iip1
     172            q(ij+iim,l,iq2)=q(ij,l,iq2)
     173         ENDDO
     174        ENDDO
     175      enddo !do ifils=1,nqdesc(iq) 
     176      endif ! if (nqfils(iq).gt.0) then
     177      write(*,*) 'vlspltqs 183: fin de la routine'
    158178
    159179      RETURN
    160180      END
    161       SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
     181      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
     182      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
     183
    162184c
    163185c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    179201c   Arguments:
    180202c   ----------
    181       REAL masse(ip1jmp1,llm),pente_max
     203      REAL masse(ip1jmp1,llm,nqtot),pente_max
    182204      REAL u_m( ip1jmp1,llm )
    183       REAL q(ip1jmp1,llm)
     205      REAL q(ip1jmp1,llm,nqtot)
    184206      REAL qsat(ip1jmp1,llm)
     207      INTEGER iq ! CRisi
    185208c
    186209c      Local
    … …  
    195218      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    196219      REAL u_mq(ip1jmp1,llm)
     220
     221      ! CRisi
     222      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
     223      INTEGER ifils,iq2 ! CRisi
    197224
    198225      Logical first,testcpu
    … …  
    227254         DO l = 1, llm
    228255            DO ij=iip2,ip1jm-1
    229                dxqu(ij)=q(ij+1,l)-q(ij,l)
     256               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    230257c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    231 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     258c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    232259            ENDDO
    233260            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    281308         DO l = 1, llm
    282309            DO ij=iip2,ip1jm-1
    283                dxqu(ij)=q(ij+1,l)-q(ij,l)
     310               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    284311            ENDDO
    285312            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    323350      DO l=1,llm
    324351       DO ij=iip2,ip1jm-1
    325           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    326      ,                     1.+u_m(ij,l)/masse(ij+1,l),
     352          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     353     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
    327354     ,                     u_m(ij,l))
    328355          zdum(ij,l)=0.5*zdum(ij,l)
    329356          u_mq(ij,l)=cvmgp(
    330      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    331      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     357     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     358     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    332359     ,                u_m(ij,l))
    333360          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    341368       DO ij=iip2,ip1jm-1
    342369          IF (u_m(ij,l).gt.0.) THEN
    343              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
     370             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
    344371             u_mq(ij,l)=u_m(ij,l)*
    345      $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
     372     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
    346373          ELSE
    347              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
     374             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
    348375             u_mq(ij,l)=u_m(ij,l)*
    349      $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
     376     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
    350377          ENDIF
    351378       ENDDO
    … …  
    416443                     i=ijq-(j-1)*iip1
    417444c   accumulation pour les mailles completements advectees
    418                      do while(zu_m.gt.masse(ijq,l))
    419                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    420                         zu_m=zu_m-masse(ijq,l)
     445                     do while(zu_m.gt.masse(ijq,l,iq))
     446                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
     447     &                          *masse(ijq,l,iq)
     448                        zu_m=zu_m-masse(ijq,l,iq)
    421449                        i=mod(i-2+iim,iim)+1
    422450                        ijq=(j-1)*iip1+i
    … …  
    424452c   ajout de la maille non completement advectee
    425453                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    426      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     454     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
     455     &                  *dxq(ijq,l))
    427456                  ELSE
    428457                     ijq=ij+1
    429458                     i=ijq-(j-1)*iip1
    430459c   accumulation pour les mailles completements advectees
    431                      do while(-zu_m.gt.masse(ijq,l))
    432                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    433                         zu_m=zu_m+masse(ijq,l)
     460                     do while(-zu_m.gt.masse(ijq,l,iq))
     461                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     462     &                          *masse(ijq,l,iq)
     463                        zu_m=zu_m+masse(ijq,l,iq)
    434464                        i=mod(i,iim)+1
    435465                        ijq=(j-1)*iip1+i
    436466                     ENDDO
    437467c   ajout de la maille non completement advectee
    438                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    439      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     468                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     469     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    440470                  ENDIF
    441471               ENDDO
    … …  
    454484      ENDDO
    455485
     486! CRisi: appel récursif de l'advection sur les fils.
     487! Il faut faire ça avant d'avoir mis à jour q et masse
     488      write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq)
     489     
     490      if (nqfils(iq).gt.0) then 
     491       do ifils=1,nqdesc(iq)
     492         iq2=iqfils(ifils,iq)
     493         DO l=1,llm
     494          DO ij=iip2,ip1jm
     495           ! On a besoin de q et masse seulement entre iip2 et ip1jm
     496           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     497           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     498          enddo   
     499         enddo
     500        enddo !do ifils=1,nqdesc(iq)
     501        do ifils=1,nqfils(iq)
     502         iq2=iqfils(ifils,iq)
     503         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
     504        enddo !do ifils=1,nqfils(iq)
     505      endif !if (nqfils(iq).gt.0) then
     506! end CRisi
    456507
    457508c   calcul des tendances
    … …  
    459510      DO l=1,llm
    460511         DO ij=iip2+1,ip1jm
    461             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    462             q(ij,l)=(q(ij,l)*masse(ij,l)+
     512            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     513            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    463514     &      u_mq(ij-1,l)-u_mq(ij,l))
    464515     &      /new_m
    465             masse(ij,l)=new_m
     516            masse(ij,l,iq)=new_m
    466517         ENDDO
    467518c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    468519         DO ij=iip1+iip1,ip1jm,iip1
    469             q(ij-iim,l)=q(ij,l)
    470             masse(ij-iim,l)=masse(ij,l)
    471          ENDDO
    472       ENDDO
     520            q(ij-iim,l,iq)=q(ij,l,iq)
     521            masse(ij-iim,l,iq)=masse(ij,l,iq)
     522         ENDDO
     523      ENDDO
     524
     525      ! retablir les fils en rapport de melange par rapport a l'air:
     526      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
     527      ! puis on boucle en longitude
     528      if (nqdesc(iq).gt.0) then 
     529       do ifils=1,nqdesc(iq)
     530         iq2=iqfils(ifils,iq) 
     531         DO l=1,llm
     532          DO ij=iip2+1,ip1jm
     533            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     534          enddo
     535          DO ij=iip1+iip1,ip1jm,iip1
     536             q(ij-iim,l,iq2)=q(ij,l,iq2)
     537          enddo ! DO ij=ijb+iip1-1,ije,iip1
     538         enddo !DO l=1,llm
     539        enddo !do ifils=1,nqdesc(iq)
     540      endif !if (nqfils(iq).gt.0) then
    473541
    474542c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    … …  
    478546      RETURN
    479547      END
    480       SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
     548      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
     549      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    481550c
    482551c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    502571c   Arguments:
    503572c   ----------
    504       REAL masse(ip1jmp1,llm),pente_max
     573      REAL masse(ip1jmp1,llm,nqtot),pente_max
    505574      REAL masse_adv_v( ip1jm,llm)
    506       REAL q(ip1jmp1,llm)
     575      REAL q(ip1jmp1,llm,nqtot)
    507576      REAL qsat(ip1jmp1,llm)
     577      INTEGER iq ! CRisi
    508578c
    509579c      Local
    … …  
    529599      SAVE sinlon,coslon,sinlondlon,coslondlon
    530600      SAVE airej2,airejjm
     601
     602      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     603      INTEGER ifils,iq2 ! CRisi
    531604c
    532605c
    … …  
    567640
    568641      DO i = 1, iim
    569       airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
    570       airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     642      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
     643      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    571644      ENDDO
    572645      qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    576649
    577650      DO ij=1,ip1jm
    578          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     651         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    579652         adyqv(ij)=abs(dyqv(ij))
    580653      ENDDO
    … …  
    591664
    592665      DO ij=1,iip1
    593          dyq(ij,l)=qpns-q(ij+iip1,l)
    594          dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     666         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
     667         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    595668      ENDDO
    596669
    … …  
    710783       DO ij=1,ip1jm
    711784         IF( masse_adv_v(ij,l).GT.0. ) THEN
    712            qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
    713      ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
     785           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
     786     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
     787     ,      /masse(ij+iip1,l,iq)))
    714788         ELSE
    715               qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
    716      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
     789              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
     790     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
    717791         ENDIF
    718792          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    721795
    722796
     797! CRisi: appel récursif de l'advection sur les fils.
     798! Il faut faire ça avant d'avoir mis à jour q et masse
     799      write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
     800   
     801      if (nqfils(iq).gt.0) then 
     802       do ifils=1,nqdesc(iq)
     803         iq2=iqfils(ifils,iq)
     804         DO l=1,llm
     805         DO ij=1,ip1jmp1
     806           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     807           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     808          enddo   
     809         enddo
     810        enddo !do ifils=1,nqdesc(iq)
     811
     812        do ifils=1,nqfils(iq)
     813         iq2=iqfils(ifils,iq)
     814         write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
     815         call vly(Ratio,pente_max,masseq,qbyv,iq2)
     816        enddo !do ifils=1,nqfils(iq)
     817      endif !if (nqfils(iq).gt.0) then
     818
    723819      DO l=1,llm
    724820         DO ij=iip2,ip1jm
    725             newmasse=masse(ij,l)
     821            newmasse=masse(ij,l,iq)
    726822     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    727             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    728      &         /newmasse
    729             masse(ij,l)=newmasse
     823            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
     824     &         -qbyv(ij-iip1,l))/newmasse
     825            masse(ij,l,iq)=newmasse
    730826         ENDDO
    731827c.-. ancienne version
    … …  
    733829         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
    734830         DO ij = 1,iip1
    735             newmasse=masse(ij,l)+convmpn*aire(ij)
    736             q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
     831            newmasse=masse(ij,l,iq)+convmpn*aire(ij)
     832            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
    737833     &               newmasse
    738             masse(ij,l)=newmasse
     834            masse(ij,l,iq)=newmasse
    739835         ENDDO
    740836         convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
    741837         convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
    742838         DO ij = ip1jm+1,ip1jmp1
    743             newmasse=masse(ij,l)+convmps*aire(ij)
    744             q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
     839            newmasse=masse(ij,l,iq)+convmps*aire(ij)
     840            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
    745841     &               newmasse
    746             masse(ij,l)=newmasse
     842            masse(ij,l,iq)=newmasse
    747843         ENDDO
    748844c.-. fin ancienne version
    … …  
    757853c        DO ij = 1,iip1
    758854c           q(ij,l)=newq
    759 c           masse(ij,l)=newmasse*aire(ij)
     855c           masse(ij,l,iq)=newmasse*aire(ij)
    760856c        ENDDO
    761857c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    … …  
    767863c        DO ij = ip1jm+1,ip1jmp1
    768864c           q(ij,l)=newq
    769 c           masse(ij,l)=newmasse*aire(ij)
     865c           masse(ij,l,iq)=newmasse*aire(ij)
    770866c        ENDDO
    771867c._. fin nouvelle version
    772868      ENDDO
    773869
     870      write(*,*) 'vly 866'
     871
     872! retablir les fils en rapport de melange par rapport a l'air:
     873      if (nqdesc(iq).gt.0) then 
     874       do ifils=1,nqdesc(iq)
     875         iq2=iqfils(ifils,iq) 
     876         DO l=1,llm
     877          DO ij=1,ip1jmp1
     878            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     879          enddo
     880         enddo
     881        enddo !do ifils=1,nqdesc(iq)
     882      endif !if (nqfils(iq).gt.0) then
     883      write(*,*) 'vly 879'
     884
    774885      RETURN
    775886      END
Note: See TracChangeset for help on using the changeset viewer.