Changeset 2296 for trunk


Ignore:
Timestamp:
Apr 24, 2020, 6:55:59 PM (5 years ago)
Author:
mvals
Message:

Mars GCM:
dyn3dpar: Implementation in the parallel version of the dynamics of the dynamical transport of the isotopes using the same scheme as the one
implemented by Camille Risi in the Earth GCM (see LMDZ6/libf/dyn3dmem, and Risi et al. 2009: genealogy of the tracers+transport of the isotopic
ratio). This is added to prepare the future implementation of the HDO cycle in the GCM. These changes had been already added in the sequential part.
dyn3d: corrections to prevent from dividing by zero in the case of the use of the isotopes scheme.
traceur.def.isotopes: example of how to write the traceur.def if you want to use the new dynamics with the genealogy of the isotopes: air is the
father of all, H2O is the father of HDO.
MV

Location:
trunk
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

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

    r1508 r2296  
    456456          DO ij=iip2,ip1jm
    457457           ! 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)
     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           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     461           masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
     462           if (q(ij,l,iq).gt.1e-16) then
     463             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     464           else
     465             Ratio(ij,l,iq2)=0.
     466           endif
    460467          enddo   
    461468         enddo
    … …  
    473480      DO l=1,llm
    474481         DO ij=iip2+1,ip1jm
    475             new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     482            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     483            new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),1e-16)
    476484            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    477485     &      u_mq(ij-1,l)-u_mq(ij,l))
    … …  
    777785           ! attention, chaque fils doit avoir son masseq, sinon, le 1er
    778786           ! 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)     
     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           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     790           masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
     791           if (q(ij,l,iq).gt.1e-16) then
     792             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     793           else
     794             Ratio(ij,l,iq2)=0.
     795           endif
    781796          enddo   
    782797         enddo
    … …  
    9971012         DO l=1,llm
    9981013          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)       
     1014           !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     1015           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     1016           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     1017           masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
     1018           if (q(ij,l,iq).gt.1e-16) then
     1019             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     1020           else
     1021             Ratio(ij,l,iq2)=0.
     1022           endif     
    10011023          enddo   
    10021024         enddo

  • trunk/LMDZ.COMMON/libf/dyn3d_common/infotrac.F90

    r1971 r2296  
    109109 
    110110    INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
    111     INTEGER :: iq, new_iq, iiq, jq, ierr, ierr2, ierr3
     111    INTEGER :: iq, new_iq, iiq, jq, ierr, ierr1, ierr2, ierr3, ierr4 !MVals: add ierr1 et ierr4 for the isotope case
    112112    INTEGER :: ifils,ipere,generation ! CRisi
    113113    LOGICAL :: continu,nouveau_traceurdef
    … …  
    533533            ! try to be smart when reading traceur.def
    534534            read(90,'(80a)') line ! store the line from traceur.def
     535            ! if format is hadv, vadv, tnom_0, tnom_transp
     536            read(line,*,iostat=ierr1) hadv(iq),vadv(iq),tnom_0(iq),tnom_transp(iq)
     537            if (ierr1.ne.0) then
    535538            ! assume format is hadv,vadv,tnom_0
    536             read(line,*,iostat=ierr2) hadv(iq),vadv(iq),tnom_0(iq)
    537             if (ierr2.ne.0) then
     539             read(line,*,iostat=ierr2) hadv(iq),vadv(iq),tnom_0(iq)
     540             if (ierr2.ne.0) then
    538541              ! maybe format is tnom0,hadv,vadv
    539542              read(line,*,iostat=ierr3) tnom_0(iq),hadv(iq),vadv(iq)
    540               if (ierr3.ne.0) then
    541                 ! assume only tnom0 is provided (havd and vad default to 10)
    542                 read(line,*) tnom_0(iq)
     543               if (ierr3.ne.0) then
     544                !assuming values of hadv et vadv
    543545                hadv(iq)=10
    544546                vadv(iq)=10
     547                read(line,*, iostat=ierr4) tnom_0(iq), tnom_transp(iq)
     548                if (ierr4.ne.0) then
     549                ! assume only tnom0 is provided (havd and vad default to 10)
     550                  read(line,*) tnom_0(iq)
     551                  tnom_transp(iq)='air'
     552                endif
    545553              endif
     554             endif
     555            !tnom_transp(iq)='air' ! no isotopes... for now... !LRossi: we now add the possibility to use isotopes (HDO cycle)
    546556            endif ! of if(ierr2.ne.0)
    547             tnom_transp(iq)='air' ! no isotopes... for now...
    548557          END DO ! of DO iq=1,nqtrue
    549558          CLOSE(90) 

  • trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90

    r2126 r2296  
    1717  USE Vampir
    1818  USE times
    19   USE infotrac, ONLY: nqtot, iadv
     19  USE infotrac, ONLY: nqtot, iadv, ok_iso_verif
    2020  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type, &
    2121                         force_conserv_tracer
    … …  
    271271
    272272     !ym          call massbar_p(massem,massebx,masseby) 
    273 
    274273     call vlspltgen_p( q,iadv, 2., massem, wg , &
    275274          pbarug,pbarvg,dtvr,p_tmp,pk_tmp,teta_tmp )
    276275
     276          !write(*,*) 'advtrac 162: apres appel vlspltgen_loc'
     277      if (ok_iso_verif) then
     278           call check_isotopes(q,1,ip1jmp1,'advtrac 162')
     279      endif !if (ok_iso_verif) then
    277280
    278281!!! ATTENTION !!!! TOUT CE QUI EST ENTRE ICI ET 1234 EST OBSOLETE !!!!!!!
    … …  
    291294        if(iadv(iq).eq.10) THEN
    292295
    293            call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
     296!           call vlsplt_p(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
     297           call vlsplt_p(q,2.,massem,wg,pbarug,pbarvg,dtvr,iq) !MVals
    294298
    295299           !   ----------------------------------------------------------------

  • trunk/LMDZ.COMMON/libf/dyn3dpar/vlsplt_p.F

    r1422 r2296  
    1       SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt)
     1      SUBROUTINE vlsplt_p(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
    22c
    33c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    1717      USE mod_hallo
    1818      USE Vampir
     19      USE infotrac, ONLY : nqtot,nqdesc,iqfils ! CRisi
    1920      IMPLICIT NONE
    2021c
    … …  
    2829c      REAL masse(iip1,jjp1,llm),pente_max
    2930      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    30       REAL q(ip1jmp1,llm)
    31 c      REAL q(iip1,jjp1,llm)
     31      REAL q(ip1jmp1,llm,nqtot) !CRisi: add dimension nqtot
    3232      REAL w(ip1jmp1,llm),pdt
     33      INTEGER iq !CRisi
    3334c
    3435c      Local
    … …  
    3839      INTEGER ijlqmin,iqmin,jqmin,lqmin
    3940c
    40       REAL zm(ip1jmp1,llm),newmasse
     41      REAL zm(ip1jmp1,llm,nqtot),newmasse
    4142      REAL mu(ip1jmp1,llm)
    4243      REAL mv(ip1jm,llm)
    43       REAL mw(ip1jmp1,llm+1)
    44       REAL zq(ip1jmp1,llm),zz
     44      REAL mw(ip1jmp1,llm+1,nqtot)
     45      REAL zq(ip1jmp1,llm,nqtot),zz
    4546      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
    4647      REAL second,temps0,temps1,temps2,temps3
    … …  
    5152      SAVE temps1,temps2,temps3
    5253      INTEGER iminn,imaxx
     54      INTEGER ifils,iq2 ! CRisi
    5355
    5456      REAL qmin,qmax
    … …  
    9395      DO l=1,llm
    9496        DO ij=ijb,ije
    95            mw(ij,l)=w(ij,l) * zzw
     97           mw(ij,l,iq)=w(ij,l) * zzw
    9698        ENDDO
    9799      ENDDO
    98100
    99101      DO ij=ijb,ije
    100          mw(ij,llm+1)=0.
     102         mw(ij,llm+1,iq)=0.
    101103      ENDDO
    102104     
    … …  
    106108       ijb=ij_begin
    107109       ije=ij_end
    108        zq(ijb:ije,:)=q(ijb:ije,:)
    109        zm(ijb:ije,:)=masse(ijb:ije,:)
     110       zq(ijb:ije,:,iq)=q(ijb:ije,:,iq)
     111       zm(ijb:ije,:,iq)=masse(ijb:ije,:)
    110112     
    111113     
    112114c       print*,'Entree vlx1'
    113115c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
    114       call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1)
    115       call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end)
     116      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_begin+2*iip1-1,iq)
     117      call vlx_p(zq,pente_max,zm,mu,ij_end-2*iip1+1,ij_end,iq)
    116118      call VTb(VTHallo)
    117119      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest1)
    … …  
    119121      call SendRequest(MyRequest1)
    120122      call VTe(VTHallo)
    121       call vlx_p(zq,pente_max,zm,mu,ij_begin+2*iip1,ij_end-2*iip1)
     123      call vlx_p(zq,pente_max,zm,mu,ij_begin+2*iip1,ij_end-2*iip1,iq)
    122124c      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
    123125      call VTb(VTHallo)
    … …  
    133135c      call exchange_hallo(zm,ip1jmp1,llm,1,1)
    134136     
    135       call vly_p(zq,pente_max,zm,mv)
     137      call vly_p(zq,pente_max,zm,mv,iq)
    136138c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
    137139c       print*,'Sortie vly1'
    138       call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1)
    139       call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end)
     140      call vlz_p(zq,pente_max,zm,mw,ij_begin,ij_begin+2*iip1-1,iq)
     141      call vlz_p(zq,pente_max,zm,mw,ij_end-2*iip1+1,ij_end,iq)
    140142      call VTb(VTHallo)
    141143      call Register_Hallo(zq,ip1jmp1,llm,2,2,2,2,MyRequest2)
    … …  
    143145      call SendRequest(MyRequest2)
    144146      call VTe(VTHallo)
    145       call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1)
     147      call vlz_p(zq,pente_max,zm,mw,ij_begin+2*iip1,ij_end-2*iip1,iq)
    146148      call VTb(VTHallo)
    147149      call WaitRecvRequest(MyRequest2)
    … …  
    154156     
    155157     
    156       call vly_p(zq,pente_max,zm,mv)
     158      call vly_p(zq,pente_max,zm,mv,iq)
    157159c       call minmaxq(zq,qmin,qmax,'apres vly     ')
    158160
    159161
    160       call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end)
     162      call vlx_p(zq,pente_max,zm,mu,ij_begin,ij_end,iq)
    161163c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
    162164
    … …  
    167169      DO l=1,llm
    168170         DO ij=ijb,ije
    169            q(ij,l)=zq(ij,l)
     171           q(ij,l,iq)=zq(ij,l,iq)
    170172         ENDDO
    171173      ENDDO
    … …  
    174176      DO l=1,llm
    175177         DO ij=ijb,ije-iip1+1,iip1
    176             q(ij+iim,l)=q(ij,l)
    177          ENDDO
    178       ENDDO
     178            q(ij+iim,l,iq)=q(ij,l,iq)
     179         ENDDO
     180      ENDDO
     181
     182      ! CRisi: aussi pour les fils
     183      if (nqdesc(iq).gt.0) then
     184      do ifils=1,nqdesc(iq)
     185        iq2=iqfils(ifils,iq)
     186        DO l=1,llm
     187         DO ij=ijb,ije
     188           q(ij,l,iq2)=zq(ij,l,iq2)
     189         ENDDO
     190         DO ij=ijb,ije-iip1+1,iip1
     191            q(ij+iim,l,iq2)=q(ij,l,iq2)
     192         ENDDO
     193        ENDDO
     194      enddo !do ifils=1,nqdesc(iq)   
     195      endif ! if (nqdesc(iq).gt.0) then
    179196
    180197      call WaitSendRequest(MyRequest1)
    … …  
    185202     
    186203     
    187       SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x)
     204      RECURSIVE SUBROUTINE vlx_p(q,pente_max,masse,u_m,ijb_x,ije_x,iq)
    188205
    189206c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    197214c   --------------------------------------------------------------------
    198215      USE Parallel_lmdz
     216      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    199217      IMPLICIT NONE
    200218c
    … …  
    205223c   Arguments:
    206224c   ----------
    207       REAL masse(ip1jmp1,llm),pente_max
     225      REAL masse(ip1jmp1,llm,nqtot),pente_max
    208226      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
    209       REAL q(ip1jmp1,llm)
    210       REAL w(ip1jmp1,llm)
     227      REAL q(ip1jmp1,llm,nqtot) ! CRisi: ajout dimension nqtot
     228      !REAL w(ip1jmp1,llm,nqtot) !MVals seems not useful in this case
     229      INTEGER iq ! CRisi
    211230c
    212231c      Local
    … …  
    222241      REAL u_mq(ip1jmp1,llm)
    223242
     243      REAL Ratio(ip1jmp1,llm,nqtot) ! CRisi,MVals: Ratio zq(fils)/zq(pere)
     244      INTEGER ifils,iq2 ! CRisi,MVals: indices pour les traceurs fils
     245
    224246      Logical extremum
    225247
    … …  
    238260      if (pole_nord.and.ijb==1) ijb=ijb+iip1
    239261      if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
    240          
     262       
    241263      IF (pente_max.gt.-1.e-5) THEN
    242264c       IF (pente_max.gt.10) THEN
    … …  
    250272           
    251273            DO ij=ijb,ije-1
    252                dxqu(ij)=q(ij+1,l)-q(ij,l)
     274               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    253275c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    254 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     276c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    255277            ENDDO
    256278            DO ij=ijb+iip1-1,ije,iip1
    … …  
    306328         DO l = 1, llm
    307329            DO ij=ijb,ije-1
    308                dxqu(ij)=q(ij+1,l)-q(ij,l)
     330               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    309331            ENDDO
    310332            DO ij=ijb+iip1-1,ije,iip1
    … …  
    326348c$OMP END DO NOWAIT
    327349      ENDIF ! (pente_max.lt.-1.e-5)
    328 
    329350c   bouclage de la pente en iip1:
    330351c   -----------------------------
    … …  
    348369      DO l=1,llm
    349370       DO ij=ijb,ije-1
    350           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    351      ,                     1.+u_m(ij,l)/masse(ij+1,l),
    352      ,                     u_m(ij,l))
     371          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     372     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
     373     ,                     u_m(ij,l,iq))
    353374          zdum(ij,l)=0.5*zdum(ij,l)
    354375          u_mq(ij,l)=cvmgp(
    355      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    356      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     376     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     377     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    357378     ,                u_m(ij,l))
    358379          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    367388      DO l=1,llm
    368389       DO ij=ijb,ije-1
    369 c       print*,'masse(',ij,')=',masse(ij,l)
     390c       print*,'masse(',ij,')=',masse(ij,l,iq)
    370391          IF (u_m(ij,l).gt.0.) THEN
    371              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
    372              u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
     392             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
     393             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l))
    373394          ELSE
    374              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
    375              u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
     395             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
     396             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)-0.5*zdum(ij,l)
     397     &                                              *dxq(ij+1,l))
    376398          ENDIF
    377399       ENDDO
    … …  
    426448c      PRINT*,'Nombre de points pour lesquels on advect plus que le'
    427449c     &       ,'contenu de la maille : ',n0
     450
    428451c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    429452         DO l=1,llm
    … …  
    450473                     i=ijq-(j-1)*iip1
    451474c   accumulation pour les mailles completements advectees
    452                      do while(zu_m.gt.masse(ijq,l))
    453                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    454                         zu_m=zu_m-masse(ijq,l)
     475                     do while(zu_m.gt.masse(ijq,l,iq))
     476                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)*
     477     &                                    masse(ijq,l,iq)
     478                        zu_m=zu_m-masse(ijq,l,iq)
    455479                        i=mod(i-2+iim,iim)+1
    456480                        ijq=(j-1)*iip1+i
    … …  
    458482c   ajout de la maille non completement advectee
    459483                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    460      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     484     &               (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
     485     &                                             *dxq(ijq,l))
    461486                  ELSE
    462487                     ijq=ij+1
    463488                     i=ijq-(j-1)*iip1
    464489c   accumulation pour les mailles completements advectees
    465                      do while(-zu_m.gt.masse(ijq,l))
    466                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    467                         zu_m=zu_m+masse(ijq,l)
     490                     do while(-zu_m.gt.masse(ijq,l,iq))
     491                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     492     &                                         *masse(ijq,l,iq)
     493                        zu_m=zu_m+masse(ijq,l,iq)
    468494                        i=mod(i,iim)+1
    469495                        ijq=(j-1)*iip1+i
    470496                     ENDDO
    471497c   ajout de la maille non completement advectee
    472                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    473      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     498                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     499     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    474500                  ENDIF
    475501               ENDDO
    … …  
    491517c$OMP END DO NOWAIT
    492518
     519! CRisi: appel récursif de l'advection sur les fils.
     520! Il faut faire ça avant d'avoir mis à jour q et masse
     521      if (nqfils(iq).gt.0) then 
     522       do ifils=1,nqdesc(iq)
     523         iq2=iqfils(ifils,iq)
     524c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     525         DO l=1,llm
     526          DO ij=ijb,ije
     527           ! On a besoin de q et masse seulement entre ijb et ije. On ne
     528           ! les calcule donc que de ijb à ije
     529           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     530           masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
     531           if (q(ij,l,iq).gt.1e-16) then
     532             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     533           else
     534             Ratio(ij,l,iq2)=0.
     535           endif
     536          enddo   
     537         enddo
     538c$OMP END DO NOWAIT
     539        enddo !do ifils=1,nqdesc(iq)
     540        do ifils=1,nqfils(iq)
     541         iq2=iqfils(ifils,iq)
     542         call vlx_p(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
     543        enddo !do ifils=1,nqfils(iq)
     544      endif !if (nqfils(iq).gt.0) then
     545! end CRisi
    493546c   calcul des tENDances
    494547c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    495548      DO l=1,llm
    496549         DO ij=ijb+1,ije
    497             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    498             q(ij,l)=(q(ij,l)*masse(ij,l)+
     550            !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     551            new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),1e-16)
     552            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    499553     &      u_mq(ij-1,l)-u_mq(ij,l))
    500554     &      /new_m
    501             masse(ij,l)=new_m
     555            masse(ij,l,iq)=new_m
    502556         ENDDO
    503557c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    504558         DO ij=ijb+iip1-1,ije,iip1
    505             q(ij-iim,l)=q(ij,l)
    506             masse(ij-iim,l)=masse(ij,l)
     559            q(ij-iim,l,iq)=q(ij,l,iq)
     560            masse(ij-iim,l,iq)=masse(ij,l,iq)
    507561         ENDDO
    508562      ENDDO
    … …  
    511565c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
    512566
     567! retablir les fils en rapport de melange par rapport a l'air:
     568      ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio
     569      ! puis on boucle en longitude
     570      if (nqfils(iq).gt.0) then 
     571       do ifils=1,nqdesc(iq)
     572         iq2=iqfils(ifils,iq) 
     573c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     574         DO l=1,llm
     575          DO ij=ijb+1,ije
     576            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     577          enddo
     578          DO ij=ijb+iip1-1,ije,iip1
     579             q(ij-iim,l,iq2)=q(ij,l,iq2)
     580          enddo ! DO ij=ijb+iip1-1,ije,iip1
     581         enddo !DO l=1,llm
     582c$OMP END DO NOWAIT
     583        enddo !do ifils=1,nqdesc(iq)
     584      endif !if (nqfils(iq).gt.0) then
    513585
    514586      RETURN
    … …  
    516588
    517589
    518       SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v)
     590      RECURSIVE SUBROUTINE vly_p(q,pente_max,masse,masse_adv_v,iq)
    519591c
    520592c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    529601c   --------------------------------------------------------------------
    530602      USE parallel_lmdz
     603      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    531604      USE comconst_mod, ONLY: pi
    532605
    … …  
    540613c   Arguments:
    541614c   ----------
    542       REAL masse(ip1jmp1,llm),pente_max
     615      REAL masse(ip1jmp1,llm,nqtot),pente_max
    543616      REAL masse_adv_v( ip1jm,llm)
    544       REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
     617      REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm) ! CRisi: ajout dimension nqtot
     618      INTEGER iq ! CRisi
    545619c
    546620c      Local
    … …  
    571645      SAVE airej2,airejjm
    572646c$OMP THREADPRIVATE(airej2,airejjm)
     647
     648      REAL Ratio(ip1jmp1,llm,nqtot) ! CRisi,MVals: Ratio zq(fils)/zq(pere)
     649      INTEGER ifils,iq2 ! CRisi,MVals: indices pour les traceurs fils
    573650c
    574651c
    … …  
    613690      if (pole_nord) then
    614691        DO i = 1, iim
    615           airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
     692          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
    616693        ENDDO
    617694        qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    620697      if (pole_sud) then
    621698        DO i = 1, iim
    622           airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     699          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    623700        ENDDO
    624701        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
    … …  
    635712     
    636713      DO ij=ijb,ije
    637          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     714         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    638715         adyqv(ij)=abs(dyqv(ij))
    639716      ENDDO
    … …  
    654731      IF (pole_nord) THEN
    655732        DO ij=1,iip1
    656            dyq(ij,l)=qpns-q(ij+iip1,l)
     733           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
    657734        ENDDO
    658735       
    … …  
    676753
    677754        DO ij=1,iip1
    678            dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     755           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    679756        ENDDO
    680757
    … …  
    812889       DO ij=ijb,ije
    813890          IF(masse_adv_v(ij,l).gt.0) THEN
    814               qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
    815      ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
     891              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
     892     ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l,iq))
    816893          ELSE
    817               qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
    818      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
     894              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
     895     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq))
    819896          ENDIF
    820897          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    822899      ENDDO
    823900c$OMP END DO NOWAIT
     901
     902! CRisi: appel récursif de l'advection sur les fils.
     903! Il faut faire ça avant d'avoir mis à jour q et masse
     904      !write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq)
     905
     906      ijb=ij_begin-2*iip1
     907      ije=ij_end+2*iip1
     908      if (pole_nord) ijb=ij_begin
     909      if (pole_sud)  ije=ij_end
     910   
     911      if (nqfils(iq).gt.0) then 
     912       do ifils=1,nqdesc(iq)
     913         iq2=iqfils(ifils,iq)
     914c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     915         DO l=1,llm
     916         DO ij=ijb,ije
     917           !masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     918           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     919           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     920           masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
     921           if (q(ij,l,iq).gt.1e-16) then
     922             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     923           else
     924             Ratio(ij,l,iq2)=0.
     925           endif   
     926          enddo   
     927         enddo
     928c$OMP END DO NOWAIT
     929        enddo !do ifils=1,nqdesc(iq)
     930
     931        do ifils=1,nqfils(iq)
     932         iq2=iqfils(ifils,iq)
     933         call vly_p(Ratio,pente_max,masse,qbyv,iq2)
     934        enddo !do ifils=1,nqfils(iq)
     935      endif !if (nqfils(iq).gt.0) then
     936! end CRisi
    824937     
    825938      ijb=ij_begin
    … …  
    831944      DO l=1,llm
    832945         DO ij=ijb,ije
    833             newmasse=masse(ij,l)
     946            newmasse=masse(ij,l,iq)
    834947     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    835948     
    836             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    837      &         /newmasse
    838             masse(ij,l)=newmasse
     949            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)
     950     &         +qbyv(ij,l)-qbyv(ij-iip1,l))/newmasse
     951            masse(ij,l,iq)=newmasse
    839952         ENDDO
    840953c.-. ancienne version
    … …  
    844957           convpn=SSUM(iim,qbyv(1,l),1)
    845958           convmpn=ssum(iim,masse_adv_v(1,l),1)
    846            massepn=ssum(iim,masse(1,l),1)
     959           massepn=ssum(iim,masse(1,l,iq),1)
    847960           qpn=0.
    848961           do ij=1,iim
    849               qpn=qpn+masse(ij,l)*q(ij,l)
     962              qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
    850963           enddo
    851964           qpn=(qpn+convpn)/(massepn+convmpn)
    852965           do ij=1,iip1
    853               q(ij,l)=qpn
     966              q(ij,l,iq)=qpn
    854967           enddo
    855968         endif
    … …  
    862975           convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    863976           convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    864            masseps=ssum(iim, masse(ip1jm+1,l),1)
     977           masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
    865978           qps=0.
    866979           do ij = ip1jm+1,ip1jmp1-1
    867               qps=qps+masse(ij,l)*q(ij,l)
     980              qps=qps+masse(ij,l,iq)*q(ij,l,iq)
    868981           enddo
    869982           qps=(qps+convps)/(masseps+convmps)
    870983           do ij=ip1jm+1,ip1jmp1
    871               q(ij,l)=qps
     984              q(ij,l,iq)=qps
    872985           enddo
    873986         endif
    … …  
    8991012c$OMP END DO NOWAIT
    9001013
     1014! CRisi: retablir les fils en rapport de melange par rapport a l'air:
     1015      ijb=ij_begin
     1016      ije=ij_end
     1017!      if (pole_nord) ijb=ij_begin
     1018!      if (pole_sud)  ije=ij_end
     1019
     1020      if (nqfils(iq).gt.0) then 
     1021       do ifils=1,nqdesc(iq)
     1022         iq2=iqfils(ifils,iq) 
     1023c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     1024         DO l=1,llm
     1025          DO ij=ijb,ije
     1026            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     1027          enddo
     1028         enddo
     1029c$OMP END DO NOWAIT
     1030        enddo !do ifils=1,nqdesc(iq)
     1031      endif !if (nqfils(iq).gt.0) then
     1032
    9011033      RETURN
    9021034      END
    … …  
    9041036     
    9051037     
    906       SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x)
     1038      RECURSIVE SUBROUTINE vlz_p(q,pente_max,masse,w,ijb_x,ije_x,iq)
    9071039c
    9081040c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    9171049c   --------------------------------------------------------------------
    9181050      USE Parallel_lmdz
     1051      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    9191052      IMPLICIT NONE
    9201053c
    … …  
    9251058c   Arguments:
    9261059c   ----------
    927       REAL masse(ip1jmp1,llm),pente_max
    928       REAL q(ip1jmp1,llm)
    929       REAL w(ip1jmp1,llm+1)
     1060      REAL masse(ip1jmp1,llm,nqtot),pente_max
     1061      REAL q(ip1jmp1,llm,nqtot)
     1062      REAL w(ip1jmp1,llm+1,nqtot)
     1063      INTEGER iq
    9301064c
    9311065c      Local
    … …  
    9341068      INTEGER i,ij,l,j,ii
    9351069c
    936       REAL,SAVE :: wq(ip1jmp1,llm+1)
     1070!      REAL,SAVE :: wq(ip1jmp1,llm+1)
     1071      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: wq !MVals
     1072     
    9371073      REAL newmasse
    9381074
    … …  
    9561092c    On oriente tout dans le sens de la pression c'est a dire dans le
    9571093c    sens de W
     1094      !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
     1095      ! Ces varibles doivent être déclarées en pointer et en save dans
     1096      ! vlz_loc si on veut qu'elles soient vues par tous les threads.
     1097      !REAL Ratio(ip1jmp1,llm,nqtot) ! CRisi
     1098      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: Ratio ! CRisi
     1099      INTEGER ifils,iq2 ! CRisi
     1100      LOGICAL, SAVE :: first=.TRUE.!MVals
     1101!$OMP THREADPRIVATE(first)
     1102
     1103         IF (first) THEN !MVals
     1104            first=.FALSE.
     1105!$OMP MASTER
     1106            ALLOCATE(wq(ip1jmp1,llm+1,nqtot)) !MVals
     1107            ALLOCATE(Ratio(ip1jmp1,llm,nqtot)) !MVals
     1108!$OMP END MASTER
     1109!$OMP BARRIER
     1110         END IF !MVals
    9581111
    9591112#ifdef BIDON
    … …  
    9691122      DO l=2,llm
    9701123         DO ij=ijb,ije
    971             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
     1124            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
    9721125            adzqw(ij,l)=abs(dzqw(ij,l))
    9731126         ENDDO
    … …  
    9991152         dzq(ij,llm)=0.
    10001153      ENDDO
     1154!      ALLOCATE(wq(ip1jmp1,llm+1,nqtot)) !MVals
    10011155c$OMP END MASTER
    10021156c$OMP BARRIER
    … …  
    10151169       DO l = 1,llm-1
    10161170         do  ij = ijb,ije
    1017           IF(w(ij,l+1).gt.0.) THEN
    1018              sigw=w(ij,l+1)/masse(ij,l+1)
    1019              wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
     1171          IF(w(ij,l+1,iq).gt.0.) THEN
     1172             sigw=w(ij,l+1,iq)/masse(ij,l+1,iq)
     1173             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l+1,iq)+0.5*(1.-sigw)
     1174     &                                                    *dzq(ij,l+1))
    10201175          ELSE
    1021              sigw=w(ij,l+1)/masse(ij,l)
    1022              wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
     1176             sigw=w(ij,l+1,iq)/masse(ij,l,iq)
     1177             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l,iq)-0.5*(1.+sigw)
     1178     &                                                  *dzq(ij,l))
    10231179          ENDIF
    10241180         ENDDO
    … …  
    10281184c$OMP MASTER
    10291185       DO ij=ijb,ije
    1030           wq(ij,llm+1)=0.
    1031           wq(ij,1)=0.
     1186          wq(ij,llm+1,iq)=0.
     1187          wq(ij,1,iq)=0.
    10321188       ENDDO
    10331189c$OMP END MASTER
    10341190c$OMP BARRIER
    10351191
     1192! CRisi: appel récursif de l'advection sur les fils.
     1193! Il faut faire ça avant d'avoir mis à jour q et masse
     1194      !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
     1195      if (nqfils(iq).gt.0) then 
     1196       do ifils=1,nqdesc(iq)
     1197         iq2=iqfils(ifils,iq)
     1198         !print*,"margaux:vlsplt,PERE",iq
     1199         !print*,"margaux:vlsplt,FILS",iq2
     1200c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     1201         DO l=1,llm
     1202          DO ij=ijb,ije
     1203           !masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     1204           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     1205           !MVals: veiller a ce qu'on n'ait pas de denominateur nul
     1206           masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),1e-16)
     1207           if (q(ij,l,iq).gt.1e-16) then
     1208             Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     1209           else
     1210             Ratio(ij,l,iq2)=0.
     1211           endif
     1212           !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
     1213           w(ij,l,iq2)=wq(ij,l,iq)
     1214          enddo   
     1215         enddo
     1216c$OMP END DO NOWAIT
     1217        enddo !do ifils=1,nqdesc(iq)
     1218c$OMP BARRIER
     1219
     1220        do ifils=1,nqfils(iq)
     1221         iq2=iqfils(ifils,iq)
     1222         call vlz_p(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2)
     1223        enddo !do ifils=1,nqfils(iq)
     1224      endif !if (nqfils(iq).gt.0) then
     1225! end CRisi 
     1226
     1227! CRisi: On rajoute ici une barrière car on veut être sur que tous les
     1228! wq soient synchronisés
     1229
     1230c$OMP BARRIER
    10361231c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    10371232      DO l=1,llm
    10381233         DO ij=ijb,ije
    1039             newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
    1040             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
    1041      &         /newmasse
    1042             masse(ij,l)=newmasse
    1043          ENDDO
    1044       ENDDO
    1045 c$OMP END DO NOWAIT
    1046 
     1234            newmasse=masse(ij,l,iq)+w(ij,l+1,iq)-w(ij,l,iq)
     1235            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)
     1236     &         +wq(ij,l+1,iq)-wq(ij,l,iq))/newmasse
     1237            masse(ij,l,iq)=newmasse
     1238         ENDDO
     1239      ENDDO
     1240c$OMP END DO NOWAIT
     1241
     1242! retablir les fils en rapport de melange par rapport a l'air:
     1243      if (nqfils(iq).gt.0) then
     1244       do ifils=1,nqdesc(iq)
     1245         iq2=iqfils(ifils,iq) 
     1246c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     1247         DO l=1,llm
     1248          DO ij=ijb,ije
     1249            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     1250          enddo
     1251         enddo
     1252c$OMP END DO NOWAIT
     1253        enddo !do ifils=1,nqdesc(iq)
     1254      endif !if (nqfils(iq).gt.0) then
    10471255
    10481256      RETURN

  • trunk/LMDZ.COMMON/libf/dyn3dpar/vlspltgen_p.F

    r1422 r2296  
    2626      USE Write_Field_p
    2727      USE VAMPIR
    28       USE infotrac, ONLY : nqtot
     28      USE infotrac, ONLY : nqtot,nqperes,nqdesc,nqfils,iqfils,
     29     &    ok_iso_verif
    2930      USE comconst_mod, ONLY: cpp
    3031      IMPLICIT NONE
    … …  
    5354      REAL,SAVE :: mu(ip1jmp1,llm)
    5455      REAL,SAVE :: mv(ip1jm,llm)
    55       REAL,SAVE :: mw(ip1jmp1,llm+1)
     56      !REAL,SAVE :: mw(ip1jmp1,llm+1)
     57      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: mw !MVals
    5658      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: zq
    5759      REAL zzpbar, zzw
    … …  
    6567      REAL ptarg,pdelarg,foeew,zdelta
    6668      REAL tempe(ip1jmp1)
    67       INTEGER ijb,ije,iq
     69      INTEGER ijb,ije,iq,iq2,ifils
    6870      LOGICAL, SAVE :: firstcall=.TRUE.
    6971!$OMP THREADPRIVATE(firstcall)
    … …  
    9193!$OMP MASTER
    9294            ALLOCATE(zm(ip1jmp1,llm,nqtot))
     95            ALLOCATE(mw(ip1jmp1,llm+1,nqtot)) !MVals
    9396            ALLOCATE(zq(ip1jmp1,llm,nqtot))
    9497!$OMP END MASTER
    … …  
    155158      ije=ij_end
    156159
     160      DO iq=1,nqtot
    157161c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    158162      DO l=1,llm
    159163         DO ij=ijb,ije
    160             mw(ij,l)=w(ij,l) * zzw
     164            mw(ij,l,iq)=w(ij,l) * zzw
    161165         ENDDO
    162166      ENDDO
    163167c$OMP END DO NOWAIT
    164 
     168      ENDDO
     169
     170      DO iq=1,nqtot
    165171c$OMP MASTER
    166172      DO ij=ijb,ije
    167          mw(ij,llm+1)=0.
    168       ENDDO
    169 c$OMP END MASTER
     173         mw(ij,llm+1,iq)=0.
     174      ENDDO
     175c$OMP END MASTER
     176      ENDDO
    170177
    171178c      CALL SCOPY(ijp1llm,q,1,zq,1)
    … …  
    184191      ENDDO
    185192
     193      ! verif temporaire
     194      ijb=ij_begin
     195      ije=ij_end
     196      if (ok_iso_verif) then
     197        call check_isotopes(zq,ijb,ije,'vlspltgen_loc 197')
     198      endif !if (ok_iso_verif) then
    186199
    187200c$OMP BARRIER           
    188       DO iq=1,nqtot
    189 
     201!      DO iq=1,nqtot
     202      DO iq=1,nqperes ! CRisi: on ne boucle que sur les pères= ceux qui sont transportés directement par l'air
    190203        if(iadv(iq) == 0) then
    191204       
    … …  
    193206       
    194207        else if (iadv(iq)==10) then
    195 
    196 #ifdef _ADV_HALO       
    197           call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
    198      &               ij_begin,ij_begin+2*iip1-1)
    199           call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
    200      &               ij_end-2*iip1+1,ij_end)
     208#ifdef _ADV_HALO     
     209          call vlx_p(zq,pente_max,zm,mu,
     210     &               ij_begin,ij_begin+2*iip1-1,iq)
     211          call vlx_p(zq,pente_max,zm,mu,
     212     &               ij_end-2*iip1+1,ij_end,iq)
    201213#else
    202           call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
    203      &               ij_begin,ij_end)
     214          call vlx_p(zq,pente_max,zm,mu,
     215     &               ij_begin,ij_end,iq)
    204216#endif
    205217
    … …  
    209221          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
    210222          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
    211 
     223! CRisi
     224          do ifils=1,nqdesc(iq)
     225            iq2=iqfils(ifils,iq)
     226         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest1)
     227         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest1)
     228          enddo
    212229c$OMP MASTER
    213230          call VTe(VTHallo)
    … …  
    232249          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest1)
    233250          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest1)
    234 
     251!CRisi
     252          do ifils=1,nqdesc(iq)
     253            iq2=iqfils(ifils,iq)
     254         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest1)
     255         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest1)
     256          enddo
    235257c$OMP MASTER
    236258          call VTe(VTHallo)
    … …  
    242264        endif
    243265     
    244       enddo
     266      enddo !DO iq=1,nqperes
    245267     
    246268     
    … …  
    256278c$OMP END MASTER       
    257279c$OMP BARRIER
    258       do iq=1,nqtot
     280
     281!      ! verif temporaire
     282!      ijb=ij_begin-2*iip1
     283!      ije=ij_end+2*iip1 
     284!      if (pole_nord) ijb=ij_begin
     285!      if (pole_sud)  ije=ij_end 
     286      if (ok_iso_verif) then
     287           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 293')
     288      endif !if (ok_iso_verif) then
     289
     290!      do iq=1,nqtot
     291      do iq=1,nqperes !CRisi
    259292
    260293        if(iadv(iq) == 0) then
    … …  
    265298
    266299#ifdef _ADV_HALLO
    267           call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
    268      &               ij_begin+2*iip1,ij_end-2*iip1)
     300          call vlx_p(zq,pente_max,zm,mu,
     301     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
    269302#endif       
    270303        else if (iadv(iq)==14) then
    … …  
    295328c$OMP END MASTER
    296329c$OMP BARRIER
    297  
    298       do iq=1,nqtot
     330
     331      if (ok_iso_verif) then
     332           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 326')
     333      endif !if (ok_iso_verif) then       
     334      if (ok_iso_verif) then
     335           ijb=ij_begin-2*iip1
     336           ije=ij_end+2*iip1
     337           if (pole_nord) ijb=ij_begin
     338           if (pole_sud)  ije=ij_end
     339           call check_isotopes(zq,ijb,ije,'vlspltgen_loc 336')
     340      endif !if (ok_iso_verif) then
     341
     342!      do iq=1,nqtot
     343      do iq=1,nqperes !CRisi
    299344
    300345        if(iadv(iq) == 0) then
    … …  
    304349        else if (iadv(iq)==10) then
    305350       
    306           call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
     351          call vly_p(zq,pente_max,zm,mv,iq)
    307352 
    308353        else if (iadv(iq)==14) then
    … …  
    318363       enddo
    319364
    320 
    321       do iq=1,nqtot
     365      if (ok_iso_verif) then
     366           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 357')
     367      endif !if (ok_iso_verif) then
     368
     369!      do iq=1,nqtot
     370      do iq=1,nqperes !CRisi
    322371
    323372        if(iadv(iq) == 0) then
    … …  
    329378c$OMP BARRIER       
    330379#ifdef _ADV_HALLO
    331           call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
    332      &               ij_begin,ij_begin+2*iip1-1)
    333           call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
    334      &               ij_end-2*iip1+1,ij_end)
     380          call vlz_p(zq,pente_max,zm,mw,
     381     &               ij_begin,ij_begin+2*iip1-1,iq)
     382          call vlz_p(zq,pente_max,zm,mw,
     383     &               ij_end-2*iip1+1,ij_end,iq)
    335384#else
    336           call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
    337      &               ij_begin,ij_end)
     385          call vlz_p(zq,pente_max,zm,mw,
     386     &               ij_begin,ij_end,iq)
    338387#endif
    339388c$OMP BARRIER
    … …  
    345394          call Register_Hallo(zq(1,1,iq),ip1jmp1,llm,2,2,2,2,MyRequest2)
    346395          call Register_Hallo(zm(1,1,iq),ip1jmp1,llm,1,1,1,1,MyRequest2)
    347 
     396          ! CRisi
     397          do ifils=1,nqdesc(iq)
     398            iq2=iqfils(ifils,iq)
     399         call Register_Hallo(zq(1,1,iq2),ip1jmp1,llm,2,2,2,2,MyRequest2)
     400         call Register_Hallo(zm(1,1,iq2),ip1jmp1,llm,1,1,1,1,MyRequest2)
     401          enddo
    348402c$OMP MASTER
    349403          call VTe(VTHallo)
    … …  
    369423c$OMP END MASTER       
    370424
    371 c$OMP BARRIER
    372       do iq=1,nqtot
     425      if (ok_iso_verif) then
     426           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 429')
     427      endif !if (ok_iso_verif) then
     428
     429c$OMP BARRIER
     430!      do iq=1,nqtot
     431      do iq=1,nqperes !CRisi
    373432
    374433        if(iadv(iq) == 0) then
    … …  
    380439
    381440#ifdef _ADV_HALLO
    382           call vlz_p(zq(1,1,iq),pente_max,zm(1,1,iq),mw,
    383      &               ij_begin+2*iip1,ij_end-2*iip1)
     441          call vlz_p(zq,pente_max,zm,mw,
     442     &               ij_begin+2*iip1,ij_end-2*iip1,iq)
    384443#endif
    385444
    … …  
    408467c$OMP BARRIER
    409468
    410 
    411       do iq=1,nqtot
     469      !write(*,*) 'vlspltgen_loc 494'
     470      if (ok_iso_verif) then
     471           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 461')
     472      endif !if (ok_iso_verif) then
     473
     474!      do iq=1,nqtot
     475      do iq=1,nqperes !CRisi
    412476
    413477        if(iadv(iq) == 0) then
    … …  
    417481        else if (iadv(iq)==10) then
    418482       
    419           call vly_p(zq(1,1,iq),pente_max,zm(1,1,iq),mv)
     483          call vly_p(zq,pente_max,zm,mv,iq)
    420484 
    421485        else if (iadv(iq)==14) then
    … …  
    429493        endif
    430494       
    431        enddo
    432 
    433       do iq=1,nqtot
     495       enddo !do iq=1,nqperes
     496
     497      if (ok_iso_verif) then
     498           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 493')
     499      endif !if (ok_iso_verif) then
     500
     501!      do iq=1,nqtot
     502      do iq=1,nqperes !CRisi
    434503
    435504        if(iadv(iq) == 0) then
    … …  
    439508        else if (iadv(iq)==10) then
    440509       
    441           call vlx_p(zq(1,1,iq),pente_max,zm(1,1,iq),mu,
    442      &               ij_begin,ij_end)
     510          call vlx_p(zq,pente_max,zm,mu,
     511     &               ij_begin,ij_end,iq)
    443512 
    444513        else if (iadv(iq)==14) then
    … …  
    453522        endif
    454523       
    455        enddo
    456 
     524       enddo !do iq=1,nqperes
     525
     526      !write(*,*) 'vlspltgen 550: apres derniere serie de call vlx'
     527      if (ok_iso_verif) then
     528           call check_isotopes(zq,ij_begin,ij_end,'vlspltgen_loc 521')
     529      endif !if (ok_iso_verif) then
    457530     
    458531      ijb=ij_begin
    … …  
    483556      ENDDO
    484557       
    485      
     558      if (ok_iso_verif) then
     559           call check_isotopes(q,ij_begin,ij_end,'vlspltgen_loc 557')
     560      endif !if (ok_iso_verif) then
     561
    486562c$OMP BARRIER
    487563

  • trunk/LMDZ.MARS/README

    r2294 r2296  
    29642964Now default emissivity/albedo is 0.95/0.1
    29652965Flags to change these fields are now more explicit: surfemis_without_startfi instead of surfemis (idem for surfalbedo)
     2966
     2967== 24/04/2020 == MV
     2968dyn3dpar: Implementation in the parallel version of the dynamics of the dynamical transport of the isotopes using the same scheme as the one implemented by Camille Risi in the Earth GCM (see LMDZ6/libf/dyn3dmem, and Risi et al. 2009: genealogy of the tracers+transport of the isotopic ratio). This is added to prepare the future implementation of the HDO cycle in the GCM. These changes had been already added in the sequential part.
     2969dyn3d: corrections to prevent from dividing by zero in the case of the use of the isotopes scheme.
     2970traceur.def.isotopes: example of how to write the traceur.def if you want to use the new dynamics with the genealogy of the isotopes: air is the father of all, H2O is the father of HDO.
Note: See TracChangeset for help on using the changeset viewer.